Repository: CliMT/climt (language: python)

Issue:
Merge branch with dycore removed



Discussion / hints:

Fix patch:
diff --git a/.github/workflows/python-package.yml b/.github/workflows/python-package.yml
index 46cbe37c..97f712bb 100644
--- a/.github/workflows/python-package.yml
+++ b/.github/workflows/python-package.yml
@@ -5,7 +5,7 @@ name: Python package
 
 on:
   push:
-    branches: [ master ]
+    branches: [ master, without-dycore ]
   pull_request:
     branches: [ master, develop ]
 
@@ -15,24 +15,29 @@ jobs:
     strategy:
       matrix:
         os: [ubuntu-18.04, macos-10.15]
-        python-version: [py36 , py37, py38]
+        python-version: [py37, py38, py39, py310]
     steps:
     - uses: actions/checkout@v2
-    - name: Set up Python 36
-      if: ${{ matrix.python-version == 'py36' }} 
-      uses: actions/setup-python@v2
-      with:
-        python-version: 3.6
     - name: Set up Python 37
       if: ${{ matrix.python-version == 'py37' }}
       uses: actions/setup-python@v2
       with:
-        python-version: 3.7
+        python-version: "3.7"
     - name: Set up Python 38
       if: ${{ matrix.python-version == 'py38' }}
       uses: actions/setup-python@v2
       with:
-        python-version: 3.8
+        python-version: "3.8"
+    - name: Set up Python 39
+      if: ${{ matrix.python-version == 'py39' }}
+      uses: actions/setup-python@v2
+      with:
+        python-version: "3.9"
+    - name: Set up Python 310
+      if: ${{ matrix.python-version == 'py310' }}
+      uses: actions/setup-python@v2
+      with:
+        python-version: "3.10"
 
     - name: Install dependencies
       run: |
@@ -42,14 +47,16 @@ jobs:
     - name: Lint with flake8
       run: |
         # stop the build if there are Python syntax errors or undefined names
-        flake8 climt/ --count --show-source --statistics
+        flake8 --ignore=E501,E226,W504,W605 climt/ --count --show-source --statistics
     - name: Build climt
       run: |
             if [ "$RUNNER_OS" == "Linux" ]; then
               gfortran --version
-              CC=gcc-7 FC=gfortran-7 CLIMT_ARCH=Linux tox -e ${{ matrix.python-version }}
+              CLIMT_ARCH=Linux tox -e ${{ matrix.python-version }}
             elif [ "$RUNNER_OS" == "macOS" ]; then
-              CC=gcc-8 FC=gfortran-8 CLIMT_ARCH=Darwin tox -e ${{ matrix.python-version }}
+              # which gcc #/usr/bin/gcc
+              # which gfortran
+              CLIMT_ARCH=Darwin tox -e ${{ matrix.python-version }}
             else 
               echo "$RUNNER_OS not supported"
               exit 1
diff --git a/.github/workflows/windows-package.yml b/.github/workflows/windows-package.yml
index 9e9bb38f..37edcb05 100644
--- a/.github/workflows/windows-package.yml
+++ b/.github/workflows/windows-package.yml
@@ -4,7 +4,7 @@ name: windows package
 
 on:
   push:
-    branches: [ master ]
+    branches: [ master, without-dycore ]
   pull_request:
     branches: [ develop ]
 
@@ -16,7 +16,7 @@ jobs:
         shell: bash -l {0}
     strategy:
       matrix:
-        python-version: ["3.6","3.7", "3.8"]
+        python-version: ["3.6"]
     steps:
     - uses: actions/checkout@v2
     - uses: conda-incubator/setup-miniconda@v2
@@ -29,7 +29,7 @@ jobs:
           export CLIMT_OPT_FLAGS=-fno-asynchronous-unwind-tables
           conda activate test
           which python
-          conda install -c msys2 m2w64-gcc m2w64-gcc-fortran m2w64-winpthreads-git m2w64-openblas
+          # conda install -c msys2 m2w64-gcc m2w64-gcc-fortran m2w64-winpthreads-git m2w64-openblas
           conda install libpython
           which gcc
           which gfortran
@@ -38,8 +38,8 @@ jobs:
           ls $HOME; cat $HOME/pydistutils.cfg
           python -m pip install -r requirements_dev.txt
           touch libmsvcr140.dll
-          cp /c/Miniconda/envs/test/Library/bin/vcruntime140.dll /c/Miniconda/envs/test/Library/bin/msvcr140.dll
-          ls -l
+          # cp /c/Miniconda/envs/test/Library/bin/vcruntime140.dll /c/Miniconda/envs/test/Library/bin/msvcr140.dll
+          # ls -l
           python setup.py develop
           py.test -v
 
diff --git a/.gitignore b/.gitignore
index ca4d37e3..6a24258e 100644
--- a/.gitignore
+++ b/.gitignore
@@ -72,3 +72,5 @@ target/
 
 # pyenv python configuration file
 .python-version
+
+.vscode/
\ No newline at end of file
diff --git a/HISTORY.rst b/HISTORY.rst
index e2aa6a14..955f969f 100644
--- a/HISTORY.rst
+++ b/HISTORY.rst
@@ -2,6 +2,11 @@
 History
 =======
 
+Latest
+-------
+
+* Removed dycore to move it to independent package
+
 v.0.16.15
 ---------
 
diff --git a/climt/__init__.py b/climt/__init__.py
index f5ae7011..f11fda47 100644
--- a/climt/__init__.py
+++ b/climt/__init__.py
@@ -11,7 +11,7 @@
     Frierson06LongwaveOpticalDepth, GrayLongwaveRadiation, HeldSuarez,
     GridScaleCondensation, BergerSolarInsolation, SimplePhysics, RRTMGLongwave,
     RRTMGShortwave,
-    EmanuelConvection, SlabSurface, GFSDynamicalCore,
+    EmanuelConvection, SlabSurface,
     DcmipInitialConditions, IceSheet, Instellation, DryConvectiveAdjustment, BucketHydrology)
 
 
@@ -25,7 +25,7 @@
     Frierson06LongwaveOpticalDepth, GrayLongwaveRadiation, HeldSuarez,
     GridScaleCondensation, BergerSolarInsolation, SimplePhysics, RRTMGLongwave,
     RRTMGShortwave,
-    EmanuelConvection, SlabSurface, GFSDynamicalCore, DcmipInitialConditions,
+    EmanuelConvection, SlabSurface, DcmipInitialConditions,
     IceSheet, Instellation, DryConvectiveAdjustment, BucketHydrology)
 
 __version__ = '0.16.25'
diff --git a/climt/_components/__init__.py b/climt/_components/__init__.py
index d69265d8..7e6fe9ff 100644
--- a/climt/_components/__init__.py
+++ b/climt/_components/__init__.py
@@ -7,7 +7,6 @@
 from .emanuel import EmanuelConvection
 from .slab_surface import SlabSurface
 from .surface_ice import IceSheet
-from .gfs import GFSDynamicalCore
 from .dcmip import DcmipInitialConditions
 from .instellation import Instellation
 from .dry_convection import DryConvectiveAdjustment
@@ -17,5 +16,5 @@
     Frierson06LongwaveOpticalDepth, GrayLongwaveRadiation,
     HeldSuarez, GridScaleCondensation, BergerSolarInsolation, SimplePhysics,
     RRTMGLongwave, RRTMGShortwave, EmanuelConvection, SlabSurface,
-    GFSDynamicalCore, DcmipInitialConditions, IceSheet,
+    DcmipInitialConditions, IceSheet,
     Instellation, DryConvectiveAdjustment, BucketHydrology)
diff --git a/climt/_components/bucket_hydrology/component.py b/climt/_components/bucket_hydrology/component.py
index 3efd43a5..bebfa584 100644
--- a/climt/_components/bucket_hydrology/component.py
+++ b/climt/_components/bucket_hydrology/component.py
@@ -1,14 +1,13 @@
-from math import sqrt
 from sympl import Stepper
+from sympl import initialize_numpy_arrays_with_properties
+import numpy as np
 
 
 class BucketHydrology(Stepper):
     """
     Manages surface energy and moisture balance
-
     This component assumes that the surface is a slab with some heat capacity and moisture holding capacity.
     Calculates the sensible and latent heat flux, takes precipitation values as input.
-
     """
 
     input_properties = {
@@ -76,6 +75,10 @@ class BucketHydrology(Stepper):
             'dims': ['mid_levels', '*'],
             'units': 'm s^-1',
         },
+        'area_type': {
+            'dims': ['*'],
+            'units': 'dimensionless',
+        },
     }
 
     diagnostic_properties = {
@@ -98,21 +101,25 @@ class BucketHydrology(Stepper):
     }
 
     output_properties = {
-        'surface_temperature': {'units': 'degK'},
-        'lwe_thickness_of_soil_moisture_content': {'units': 'm'},
+        'surface_temperature': {
+            'dims': ['*'],
+            'units': 'degK',
+        },
+        'lwe_thickness_of_soil_moisture_content': {
+            'dims': ['*'],
+            'units': 'm'},
     }
 
     def __init__(self, soil_moisture_max=0.15, beta_parameter=0.75,
                  specific_latent_heat_of_water=2260000, bulk_coefficient=0.0011, **kwargs):
         """
         Args:
-
         soil_moisture_max:
             The maximum moisture that can be held by the surface_temperature
         beta_parameter:
             A constant value that is used in the beta_factor calculation.
         bulk_coefficient:
-            The bulk transfer coeffiecient that is used to calculate
+            The bulk transfer coefficient that is used to calculate
             maximum evaporation rate and sensible heat flux
         """
         self._smax = soil_moisture_max
@@ -129,33 +136,43 @@ def array_call(self, state, timestep):
 
         beta_factor = 0
 
-        north_wind_speed = state['northward_wind'][0][0]
-        east_wind_speed = state['eastward_wind'][0][0]
+        new_state = initialize_numpy_arrays_with_properties(
+            self.output_properties, state, self.input_properties
+        )
+
+        diagnostics = initialize_numpy_arrays_with_properties(
+            self.diagnostic_properties, state, self.input_properties
+        )
+
+        north_wind_speed = state['northward_wind'][0]
+        east_wind_speed = state['eastward_wind'][0]
+
+        wind_speed = np.sqrt(np.power(north_wind_speed, 2) +
+                             np.power(east_wind_speed, 2))
 
-        wind_speed = sqrt(pow(north_wind_speed, 2) +
-                          pow(east_wind_speed, 2))
         potential_evaporation = self._c * wind_speed * (state['surface_specific_humidity'] - state['specific_humidity'][0])
 
         precipitation_rate = state['convective_precipitation_rate'] + state['stratiform_precipitation_rate']
+        diagnostics['precipitation_rate'] = precipitation_rate
 
-        soil_moisture = state['lwe_thickness_of_soil_moisture_content'][0]
+        soil_moisture = state['lwe_thickness_of_soil_moisture_content']
 
-        soil_moisture_tendency = 0
+        beta_factor = np.ones(soil_moisture.shape)
 
-        if soil_moisture >= self._g * self._smax:
-            beta_factor = 1
-        else:
-            beta_factor = soil_moisture/(self._g*self._smax)
+        mask = soil_moisture <= self._g * self._smax
+        beta_factor[mask] = soil_moisture[mask]/(self._g*self._smax)
 
         evaporation_rate = beta_factor * potential_evaporation
+        diagnostics['evaporation_rate'] = evaporation_rate
 
-        if soil_moisture < self._smax or precipitation_rate <= evaporation_rate:
-            soil_moisture_tendency = precipitation_rate - evaporation_rate
-        else:
-            soil_moisture_tendency = 0
+        soil_moisture_tendency = np.zeros(soil_moisture.shape)
+        mask = np.logical_or(soil_moisture < self._smax, precipitation_rate <= evaporation_rate)
+        soil_moisture_tendency[mask] = precipitation_rate[mask] - evaporation_rate[mask]
 
         surface_upward_latent_heat_flux = self._l * evaporation_rate
         surface_upward_sensible_heat_flux = self._c * wind_speed * (state['surface_temperature'] - state['air_temperature'][0])
+        diagnostics['surface_upward_sensible_heat_flux'] = surface_upward_sensible_heat_flux
+        diagnostics['surface_upward_latent_heat_flux'] = surface_upward_latent_heat_flux
 
         net_heat_flux = (
             state['downwelling_shortwave_flux_in_air'][:, 0] +
@@ -163,30 +180,15 @@ def array_call(self, state, timestep):
             state['upwelling_shortwave_flux_in_air'][:, 0] -
             state['upwelling_longwave_flux_in_air'][:, 0] -
             surface_upward_sensible_heat_flux -
-            surface_upward_latent_heat_flux
-        )
+            surface_upward_latent_heat_flux)
 
         mass_surface_slab = state['surface_material_density'] * \
             state['soil_layer_thickness']
         heat_capacity_surface = mass_surface_slab * state['heat_capacity_of_soil']
 
-        new_surface_temperature = state['surface_temperature'] + \
+        new_state['surface_temperature'] = state['surface_temperature'] + \
             (net_heat_flux/heat_capacity_surface * timestep.total_seconds())
-        new_soil_moisture = state['lwe_thickness_of_soil_moisture_content'] + \
-            (soil_moisture_tendency * timestep.total_seconds())
-
-        new_soil_moisture[new_soil_moisture > 0.15] = 0.15
-
-        new_state = {
-            'surface_temperature': new_surface_temperature,
-            'lwe_thickness_of_soil_moisture_content': new_soil_moisture,
-        }
-
-        diagnostics = {
-            'precipitation_rate': precipitation_rate,
-            'surface_upward_sensible_heat_flux': surface_upward_sensible_heat_flux,
-            'surface_upward_latent_heat_flux': surface_upward_latent_heat_flux,
-            'evaporation_rate': evaporation_rate,
-        }
+        new_state['lwe_thickness_of_soil_moisture_content'] = np.minimum(state['lwe_thickness_of_soil_moisture_content'] +
+                                                                         (soil_moisture_tendency * timestep.total_seconds()), 0.15)
 
         return diagnostics, new_state
diff --git a/climt/_components/gfs/__init__.py b/climt/_components/gfs/__init__.py
deleted file mode 100644
index 2f0fdef8..00000000
--- a/climt/_components/gfs/__init__.py
+++ /dev/null
@@ -1,3 +0,0 @@
-from .component import GFSDynamicalCore
-
-__all__ = (GFSDynamicalCore)
diff --git a/climt/_components/gfs/_gfs_dynamics.pyx b/climt/_components/gfs/_gfs_dynamics.pyx
deleted file mode 100644
index e668b457..00000000
--- a/climt/_components/gfs/_gfs_dynamics.pyx
+++ /dev/null
@@ -1,587 +0,0 @@
-cimport numpy as cnp
-import numpy as np
-import cython
-
-# Typedef for function pointer returning void and taking no arguments (for now)
-
-ctypedef void (*pyPhysicsCallback)(double *, double *, double *, double *,
-                                   double *, double *, double *,
-                              double *, double *, double *, double *,
-                                   double *, double *, double *)
-
-# Variables that control the behaviour of the model
-# bint is Cython's boolean type
-cdef extern:
-    bint dry, adiabatic, heldsuarez
-
-cdef extern:
-    void gfs_lnps_to_spectral()
-    void gfs_vrt_div_to_spectral()
-    void gfs_virtemp_to_spectral()
-    void gfs_tracer_to_spectral()
-
-#Function to set global mass of dry air
-cdef extern:
-    void gfs_set_dry_pressure(double *pdry_init)
-
-#Function to set top of model pressure
-cdef extern:
-    void gfs_set_model_top_pressure(double *p_toa)
-
-#Function to set time step
-cdef extern:
-    void gfs_set_time_step(double *dt)
-
-#Function to get longitude/latitude map from shtns
-cdef extern:
-    void gfs_get_lon_lat(double *longitudes, double *latitudes)
-
-#Function to set physical constants in the dycore
-cdef extern:
-    void gfs_set_constants(double *planetary_radius,
-                           double *omega,
-                           double *gas_constant,
-                           double *gas_constant_dry_air,
-                           double *gas_constant_condensible,
-                           double *gravity,
-                           double *specific_heat_dry_air,
-                           double *specific_heat_condensible)
-
-#Function to set grid dimensions in dycore
-cdef extern:
-    void gfs_set_model_dimensions(int *num_lats,
-                                  int *num_lons,
-                                  int *num_levs,
-                                  int *truncation,
-                                  int *spectral_dim,
-                                  int *num_tracers)
-
-#Function to set model time (might be eliminated if I have time!
-#since it is not really used)
-cdef extern:
-    void gfs_set_model_time(double *time)
-
-# Function to init dynamics
-cdef extern:
-    void gfs_init_dynamics(int *num_damp_levs, double *tau_damping)
-
-# Function to init physics
-cdef extern:
-    void gfs_init_physics()
-
-#Function to calculate pressure fields after other fields have been updated
-cdef extern:
-    void gfs_calculate_pressure()
-
-# Function to step fields by one dt
-cdef extern:
-    void gfs_take_one_step()
-
-#Function to register an external callback
-#cdef extern:
-#    void gfsRegisterPhysicsCallback(pyPhysicsCallback callback)
-
-#Function to add u,v tendency to vrt,div tendency
-cdef extern:
-    void gfs_uv_to_vrtdiv(
-        double *pyUg,
-        double *pyVg,
-        double *pyVrtg,
-        double *pyDivg)
-
-#Function to set tendencies within the fortran code, which will be called by the
-# dynamical core after the semi-implicit step
-cdef extern:
-    void gfs_set_tendencies(
-        double *pyVrtTend,
-        double *pyDivTend,
-        double *pyVirtTempTend,
-        double *pyLnpsTend,
-        double *pyTracerTend)
-
-#Function to set model topography
-cdef extern:
-    void gfs_set_topography(double *surface_geopotential)
-
-#Function to set initial atmospheric mass
-cdef extern:
-    void gfs_set_dry_pressure(double *dry_pressure)
-
-#Function to deallocate arrays in physics, etc.,
-cdef extern:
-    void gfs_finalise()
-
-#Function to convert the input grid arrays to spectral arrays
-# The data must be available in py{Ug,Vg,VirtTempg,Tracerg,Lnpsg}
-cdef extern:
-    void gfs_convert_to_spectral()
-
-#Function to convert the spectral arrays to grid arrays
-# The result will be available in py{Ug,Vg,VirtTempg,Tracerg,Lnpsg}
-cdef extern:
-    void gfs_convert_to_grid()
-
-#Function to initialise the arrays for computation
-# These are allocated in Python and passed in
-cdef extern:
-    void gfs_initialise_spectral_arrays(
-        complex *pyVrtSpec,
-        complex *pyDivSpec,
-        complex *pyVirtTempSpec,
-        complex *pyTracerSpec,
-        complex *pyTopoSpec,
-        complex *pyLnPsSpec,
-        complex *pyDissSpec,
-        complex *pyDmpProf,
-        complex *pyDiffProf)
-
-cdef extern:
-    void gfs_assign_grid_arrays(
-        double *pyUg,
-        double *pyVg,
-        double *pyVirtTempg,
-        double *pyLnPsg,
-        double *pyTracerg,
-        double *pyVrtg,
-        double *pyDivg)
-
-cdef extern:
-    void gfs_initialise_grid_arrays(
-        double *pyDlnpdtg,
-        double *pyEtaDotg,
-        double *pyPhisg,
-        double *pyDPhisdxg,
-        double *pyDPhisdyg,
-        double *pyDlnpsdtg)
-
-cdef extern:
-    void gfs_assign_pressure_arrays(
-        double *pySurfPressure,
-        double *pyPressGrid,
-        double *pyInterfacePressure)
-
-cdef extern:
-    void gfs_init_coordinate_arrays(
-        double *ak, double *bk, double *ck,
-        double *dbk, double *si, double *sl,
-        double *bkl, double *alfa, double *rlnp, double *dpk)
-
-cdef extern:
-    void gfs_initialise_spectral_physics_arrays(
-        complex *pyVrtSpecTend,
-        complex *pyDivSpecTend,
-        complex *pyVirtTempSpecTend,
-        complex *pyTracerSpecTend,
-        complex *pyLnPsSpec)
-
-cdef extern:
-    void gfs_initialise_grid_physics_arrays(
-        double *pyPwat)
-
-cdef extern:
-    void gfs_initialise_timestepper_arrays(
-        double *pyDmhyb, double *pyAmhyb, double *pyBmhyb,
-        double *pyTref, double *pyPkref, double *pyDpkref,
-        double *pyAlfaref, double *pySvhyb, double *pyTorhyb)
-
-cdef extern:
-    void gfs_assign_sht_arrays(
-        double *pyLap, double *pyInvLap, double *pyGaussWeights,
-        double *pyLons, double *pyLats, double *pyAreaWeights,
-        int *pyDegree, int *pyOrder, int *pyNlons, int *pyNlats,
-        int *pyNlm)
-
-
-# Grid space arrays, cython for initialisation
-cdef cnp.double_t[:,:,:,:] pyTracerg
-cdef cnp.double_t[:,:,:] pyUg, pyVg, pyVrtg, pyDivg,\
-        pyVirtTempg, pyDlnpdtg, pyEtaDotg
-cdef cnp.double_t[:,:] pyLnPsg, pyPhis, pyDPhisdx,\
-        pyDPhisdy, pyDlnpsdt, pyPwat
-
-# Pressure arrays, in grid space
-cdef cnp.double_t[:,:] pySurfPressure
-cdef cnp.double_t[:,:,:] pyPressGrid
-cdef cnp.double_t[:,:,:] pyInterfacePressure
-
-# Vertical Coordinate arrays
-cdef cnp.double_t[:] pyAk, pyBk, pyCk, pyDbk, pySi, pySl, pyBkl
-cdef cnp.double_t[:,:,:] pyAlfa, pyRlnp, pyDpk
-
-# Timestepper arrays
-cdef cnp.double_t[:, :] pyAmhyb, pyBmhyb
-cdef cnp.double_t[:, :, :, :] pyDmhyb
-cdef cnp.double_t[:] pyTref, pyPkref, pyDpkref, pyAlfaref
-cdef cnp.double_t[:] pySvhyb, pyTorhyb
-
-# Spectral arrays, using cython for declaration
-cdef cnp.complex128_t[:,:,:] pyTracerSpec
-
-cdef cnp.complex128_t[:,:] pyVrtSpec, pyDivSpec, pyVirtTempSpec
-
-cdef cnp.complex128_t[:] pyTopoSpec, pyLnPsSpec,\
-            pyDissSpec, pyDmpProf, pyDiffProf
-
-cdef cnp.complex128_t[:,:,:] pyTracerSpecTend
-
-cdef cnp.complex128_t[:,:] pyVrtSpecTend, pyDivSpecTend, pyVirtTempSpecTend
-
-cdef cnp.complex128_t[:] pyLnPsSpecTend
-
-# Temporary arrays for setting tendency terms
-cdef cnp.double_t[:,:,:] tempVrtTend, tempDivTend, tempVirtTempTend, tempUTend,tempVTend
-cdef cnp.double_t[:,:,:,:] tempTracerTend
-cdef cnp.double_t[:,:] tempLnpsTend
-
-# Shtns arrays
-cdef cnp.double_t[:, :] pyLats, pyLons, pyAreaWeights
-cdef cnp.double_t[:] pyLap, pyInvLap, pyGaussWeights
-cdef cnp.int_t[:] pyDegree, pyOrder
-
-
-# Grid size
-#cdef public int numLats, numLons, numTrunc, numLevs, spectralDim, numTracers
-cdef int nlm
-
-# Model state
-cdef int modelIsInitialised, physicsEnabled
-cdef double zero_model_time=0
-
-# Physics subroutine pointer
-cdef object physicsCallback
-
-#TODO get constant from physcons.f90
-cdef double fv
-
-def lnps_to_spectral():
-    gfs_lnps_to_spectral()
-
-def vrt_div_to_spectral():
-    gfs_vrt_div_to_spectral()
-
-def virtemp_to_spectral():
-    gfs_virtemp_to_spectral()
-
-def tracer_to_spectral():
-    gfs_tracer_to_spectral()
-
-def set_constants(double radius, double omega,
-                  double R, double Rd, double Rv,
-                  double g, double Cp, double C_cond):
-
-    gfs_set_constants(&radius, &omega, &R, &Rd, &Rv, &g, &Cp, &C_cond)
-
-#@cython.boundscheck(False)
-def set_model_grid(
-        int nlats, int nlons,
-        int nlevs, int ntrunc, int ndimspec,
-        int ntracers,
-        cnp.double_t[:] ak_in,
-        cnp.double_t[:] bk_in,
-    ):
-
-    global pyDmhyb, pyAmhyb, pyBmhyb,\
-            pyTref, pyPkref, pyDpkref,\
-            pyAlfaref, pySvhyb, pyTorhyb,\
-            pyAk, pyBk, pyCk, pySi, pySl, pyDbk,\
-            pyBkl, pyAlfa, pyRlnp, pyDpk,\
-            pyLap, pyInvLap, pyGaussWeights, pyLats, pyLons, pyAreaWeights,\
-            pyDegree, pyOrder, nlm
-
-
-    nlm = int((ntrunc+1)*(ntrunc+2)/2)
-
-    pyLap = np.zeros(nlm, dtype=np.double)
-    pyInvLap = np.zeros(nlm, dtype=np.double)
-    pyGaussWeights = np.zeros(nlats, dtype=np.double)
-
-    pyDegree = np.zeros(nlm, dtype=np.int)
-    pyOrder = np.zeros(nlm, dtype=np.int)
-
-    pyLats = np.zeros((nlats, nlons), dtype=np.double)
-    pyLons = np.zeros((nlats, nlons), dtype=np.double)
-    pyAreaWeights = np.zeros((nlats, nlons), dtype=np.double)
-
-    pyDmhyb = np.zeros((3, ntrunc+1, nlevs, nlevs),
-                       dtype=np.double)
-
-    pyAmhyb = np.zeros((nlevs, nlevs), dtype=np.double)
-    pyBmhyb = np.zeros((nlevs, nlevs), dtype=np.double)
-
-    pyTref = np.zeros((nlevs), dtype=np.double)
-    pyPkref = np.zeros((nlevs+1), dtype=np.double)
-    pyDpkref = np.zeros((nlevs), dtype=np.double)
-    pyAlfaref = np.zeros((nlevs), dtype=np.double)
-    pySvhyb = np.zeros((nlevs), dtype=np.double)
-    pyTorhyb = np.zeros((nlevs), dtype=np.double)
-
-    pyAk = np.zeros((nlevs+1), dtype=np.double)
-    pyBk = np.zeros((nlevs+1), dtype=np.double)
-    cdef int i
-    for i in range(nlevs+1):
-        pyAk[i] = ak_in[i]
-        pyBk[i] = bk_in[i]
-    pySi = np.zeros((nlevs+1), dtype=np.double)
-    pyCk = np.zeros((nlevs), dtype=np.double)
-    pyDbk = np.zeros((nlevs), dtype=np.double)
-    pySl = np.zeros((nlevs), dtype=np.double)
-    pyBkl = np.zeros((nlevs), dtype=np.double)
-
-    pyAlfa = np.zeros((nlevs, nlats, nlons), dtype=np.double)
-    pyRlnp = np.zeros((nlevs, nlats, nlons), dtype=np.double)
-    pyDpk = np.zeros((nlevs, nlats, nlons), dtype=np.double)
-
-    gfs_set_model_dimensions(&nlats, &nlons, &nlevs, &ntrunc,
-                             &ndimspec, &ntracers)
-
-    gfs_initialise_timestepper_arrays(
-        <double *>&pyDmhyb[0,0,0,0],
-        <double *>&pyAmhyb[0,0],
-        <double *>&pyBmhyb[0,0],
-        <double *>&pyTref[0],
-        <double *>&pyPkref[0],
-        <double *>&pyDpkref[0],
-        <double *>&pyAlfaref[0],
-        <double *>&pySvhyb[0],
-        <double *>&pyTorhyb[0])
-
-    gfs_init_coordinate_arrays(
-        <double *>&pyAk[0],
-        <double *>&pyBk[0],
-        <double *>&pyCk[0],
-        <double *>&pyDbk[0],
-        <double *>&pySi[0],
-        <double *>&pySl[0],
-        <double *>&pyBkl[0],
-        <double *>&pyAlfa[0,0,0],
-        <double *>&pyRlnp[0,0,0],
-        <double *>&pyDpk[0,0,0])
-
-    gfs_assign_sht_arrays(
-        <double *>&pyLap[0],
-        <double *>&pyInvLap[0],
-        <double *>&pyGaussWeights[0],
-        <double *>&pyLons[0, 0],
-        <double *>&pyLats[0, 0],
-        <double *>&pyAreaWeights[0, 0],
-        <int *>&pyDegree[0],
-        <int *>&pyOrder[0],
-        &nlons, &nlats, &nlm)
-
-    init_grid_arrays(nlons, nlats, nlevs, ntracers)
-    init_spectral_arrays(ndimspec, nlevs, ntracers)
-
-
-# Initialise dynamics and physics
-def init_model(
-    double dry_pressure, int num_damp_levels, double tau_damping,
-    double model_top_pressure):
-
-    # Now that the arrays are initialised, call dynamics and physics init
-
-
-    gfs_set_dry_pressure(&dry_pressure)
-    gfs_set_model_top_pressure(&model_top_pressure)
-    gfs_set_model_time(&zero_model_time)
-    gfs_init_dynamics(&num_damp_levels, &tau_damping)
-    #gfs_get_lon_lat(<double *>&__longitudes[0,0], <double *>&__latitudes[0,0])
-
-    longitudes = np.ascontiguousarray(pyLons).copy()
-    latitudes = np.ascontiguousarray(pyLats).copy()
-    area_weights = np.ascontiguousarray(pyAreaWeights).copy()
-    gaussian_weights = np.ascontiguousarray(pyGaussWeights).copy()
-    # sigma_levels = np.ascontiguousarray(pySl).copy()
-    # sigma_inteface_levels = np.ascontiguousarray(pySi).copy()
-
-    return gaussian_weights, area_weights, latitudes, longitudes#, sigma_levels, sigma_interface_levels
-
-
-# Create the spectral arrays (defined in spectral_data.f90)
-@cython.boundscheck(False)
-def init_spectral_arrays(spectral_dim, num_levs, num_tracers):
-
-    global pyTracerSpec, pyTracerSpecTend, pyTopoSpec, \
-            pyLnPsSpec, pyLnPsSpecTend,\
-            pyDissSpec, pyDmpProf, pyDiffProf,\
-            pyVrtSpec, pyVrtSpecTend, pyDivSpec, pyDivSpecTend,\
-            pyVirtTempSpec, pyVirtTempSpecTend
-
-    pyTracerSpec = np.zeros(
-        (num_tracers, num_levs, spectral_dim), dtype=complex)
-
-    pyVrtSpec = np.zeros((num_levs, spectral_dim), dtype=complex)
-    pyVirtTempSpec = np.zeros(
-        (spectral_dim, num_levs), dtype=complex)
-
-    pyDivSpec = np.zeros((num_levs, spectral_dim), dtype=complex)
-
-
-    pyTopoSpec = np.zeros(spectral_dim, dtype=complex)
-    pyLnPsSpec = np.zeros(spectral_dim, dtype=complex)
-    pyDissSpec = np.zeros(spectral_dim, dtype=complex)
-
-    pyDmpProf = np.zeros(num_levs, dtype=complex)
-    pyDiffProf = np.zeros(num_levs, dtype=complex)
-
-    pyTracerSpecTend = np.zeros(
-        (num_tracers, num_levs, spectral_dim), dtype=complex)
-
-    pyVrtSpecTend = np.zeros(
-        (num_levs, spectral_dim), dtype=complex)
-
-    pyVirtTempSpecTend = np.zeros(
-        (num_levs, spectral_dim), dtype=complex)
-
-    pyDivSpecTend = np.zeros(
-        (num_levs, spectral_dim), dtype=complex)
-
-    pyLnPsSpecTend = np.zeros(spectral_dim, dtype=complex)
-
-    gfs_initialise_spectral_arrays(
-        <double complex *>&pyVrtSpec[0,0],
-        <double complex *>&pyDivSpec[0,0],
-        <double complex *>&pyVirtTempSpec[0,0],
-        <double complex *>&pyTracerSpec[0,0,0],
-        <double complex *>&pyTopoSpec[0],
-        <double complex *>&pyLnPsSpec[0],
-        <double complex *>&pyDissSpec[0],
-        <double complex *>&pyDmpProf[0],
-        <double complex *>&pyDiffProf[0])
-
-    gfs_initialise_spectral_physics_arrays(
-        <double complex *>&pyVrtSpecTend[0,0],
-        <double complex *>&pyDivSpecTend[0,0],
-        <double complex *>&pyVirtTempSpecTend[0,0],
-        <double complex *>&pyTracerSpecTend[0,0,0],
-        <double complex *>&pyLnPsSpecTend[0])
-
-
-
-# Create the grid arrays (defined in grid_data.f90)
-def init_grid_arrays(num_lons, num_lats, num_levs, num_tracers):
-
-    global pyDlnpdtg, pyEtaDotg, pyPhis, pyDPhisdx, pyDPhisdy, \
-            pyDlnpsdt, pyPwat, tempVrtTend, tempDivTend
-
-    pyDlnpdtg = np.zeros((num_levs, num_lats, num_lons),
-                         dtype=np.double)
-    pyEtaDotg = np.zeros((num_levs+1, num_lats, num_lons),
-                         dtype=np.double)
-
-
-    pyPhis = np.zeros((num_lats, num_lons), dtype=np.double)
-    pyDPhisdx = np.zeros((num_lats, num_lons), dtype=np.double)
-    pyDPhisdy = np.zeros((num_lats, num_lons), dtype=np.double)
-    pyDlnpsdt = np.zeros((num_lats, num_lons), dtype=np.double)
-    pyPwat = np.zeros((num_lats, num_lons), dtype=np.double)
-
-
-
-    tempVrtTend = np.zeros((num_levs, num_lats, num_lons),
-                           dtype=np.double)
-    tempDivTend = np.zeros((num_levs, num_lats, num_lons),
-                           dtype=np.double)
-
-
-    gfs_initialise_grid_arrays(
-        <double *>&pyDlnpdtg[0,0,0],
-        <double *>&pyEtaDotg[0,0,0],
-        <double *>&pyPhis[0,0],
-        <double *>&pyDPhisdx[0,0],
-        <double *>&pyDPhisdy[0,0],
-        <double *>&pyDlnpsdt[0,0])
-
-    gfs_initialise_grid_physics_arrays(
-        <double *>&pyPwat[0,0])
-
-#Intialise pressure arrays
-def assign_pressure_arrays(
-    cnp.double_t[:, :] surface_pressure,
-    cnp.double_t[:, :, :] midlevel_pressure,
-    cnp.double_t[:, :, :] interface_pressure):
-
-
-    gfs_assign_pressure_arrays(
-        <double *>&surface_pressure[0,0],
-        <double *>&midlevel_pressure[0,0,0],
-        <double *>&interface_pressure[0,0,0])
-
-# Set tendencies for dynamical core to use in physics
-@cython.boundscheck(False)
-def assign_tendencies(
-    cnp.double_t[:, :, :] u_tend,
-    cnp.double_t[:, :, :] v_tend,
-    cnp.double_t[:, :, :] virtemp_tend,
-    cnp.double_t[:, :] lnps_tend,
-    cnp.double_t[:, :, :, :] tracer_tend):
-
-    global tempVrtTend, tempDivTend
-
-    gfs_uv_to_vrtdiv(
-        <double *>&u_tend[0,0,0],
-        <double *>&v_tend[0,0,0],
-        <double *>&tempVrtTend[0,0,0],
-        <double *>&tempDivTend[0,0,0])
-
-    gfs_set_tendencies(
-        <double *>&tempVrtTend[0,0,0],
-        <double *>&tempDivTend[0,0,0],
-        <double *>&virtemp_tend[0,0,0],
-        <double *>&lnps_tend[0,0],
-        <double *>&tracer_tend[0,0,0,0])
-
-# Register a callback which calculates the physics (to be used in stand-alone
-# mode only)
-
-#cdef setPhysicsCallback(physicsFnPtr):
-
-#    physicsCallback = physicsFnPtr
-#    physicsEnabled = True
-    #        gfsRegisterPhysicsCallback(testFunc)
-
-
-@cython.boundscheck(False)
-def assign_grid_arrays(
-    cnp.double_t[:, :, :] u,
-    cnp.double_t[:, :, :] v,
-    cnp.double_t[:, :, :] tvirt,
-    cnp.double_t[:, :] lnpsg,
-    cnp.double_t[:, :, :, :] tracers,
-    cnp.double_t[:, :, :] vorticity,
-    cnp.double_t[:, :, :] divergence):
-
-    gfs_assign_grid_arrays(
-        <double *>&u[0,0,0],
-        <double *>&v[0,0,0],
-        <double *>&tvirt[0,0,0],
-        <double *>&lnpsg[0,0],
-        <double *>&tracers[0,0,0,0],
-        <double *>&vorticity[0,0,0],
-        <double *>&divergence[0,0,0])
-
-
-def update_spectral_arrays():
-    gfs_convert_to_spectral()
-
-
-def take_one_step():
-    gfs_take_one_step()
-
-def convert_to_grid():
-    gfs_convert_to_grid()
-
-def calculate_pressure():
-    gfs_calculate_pressure()
-
-def set_topography(cnp.double_t[:, :] topography):
-
-    gfs_set_topography(<double *>&topography[0,0])
-
-def set_time_step(double time_step):
-
-    gfs_set_time_step(&time_step)
-
-def shut_down_model():
-    #Remember to set time to zero!
-    gfs_finalise()
-    gfs_set_model_time(&zero_model_time)
diff --git a/climt/_components/gfs/component.py b/climt/_components/gfs/component.py
deleted file mode 100644
index 1ba3fcdf..00000000
--- a/climt/_components/gfs/component.py
+++ /dev/null
@@ -1,599 +0,0 @@
-from __future__ import division
-from ..._core import ensure_contiguous_state
-from sympl import (
-    get_constant, TendencyStepper, initialize_numpy_arrays_with_properties,
-    get_numpy_arrays_with_properties, ImplicitTendencyComponentComposite,
-    get_tracer_names, restore_data_arrays_with_properties,
-)
-import numpy as np
-import sys
-import logging
-try:
-    from . import _gfs_dynamics
-except ImportError as error:
-    logging.warning(
-        'Import failed. GFS dynamical core is likely not compiled and will not '
-        'be available.'
-    )
-    print(error)
-
-
-class GFSError(Exception):
-    pass
-
-
-def get_valid_properties(gfs_properties, prognostic_properties, property_type):
-    return_dict = {}
-    for name, properties in prognostic_properties.items():
-        if name not in gfs_properties:
-            return_dict[name] = properties
-        elif 'dims' in prognostic_properties.keys():
-            extra_dims = set(prognostic_properties['dims']).difference(['*'] + list(gfs_properties[name]['dims']))
-            if len(extra_dims) != 0:
-                raise GFSError(
-                    'Cannot handle TendencyComponent with {} {} '
-                    'that has extra dimensions {} not used by GFS'.format(
-                        property_type, name, extra_dims)
-                )
-    return return_dict
-
-
-class GFSDynamicalCore(TendencyStepper):
-    """
-    Climt interface to the GFS dynamical core. The GFS
-    code is available on `github`_.
-
-    .. _github:
-       https://github.com/jswhit/gfs-dycore
-    """
-
-    _gfs_input_properties = {
-        'latitude': {
-            'units': 'degrees_N',
-            'dims': ['lat', 'lon'],
-        },
-        'longitude': {
-            'units': 'degrees_E',
-            'dims': ['lat', 'lon'],
-        },
-        'air_temperature': {
-            'units': 'degK',
-            'dims': ['mid_levels', 'lat', 'lon'],
-        },
-        'atmosphere_hybrid_sigma_pressure_a_coordinate_on_interface_levels': {
-            'units': 'dimensionless',
-            'dims': ['interface_levels'],
-            'alias': 'a_coord',
-        },
-        'atmosphere_hybrid_sigma_pressure_b_coordinate_on_interface_levels': {
-            'units': 'dimensionless',
-            'dims': ['interface_levels'],
-            'alias': 'b_coord',
-        },
-        'air_pressure': {
-            'dims': ['mid_levels', 'lat', 'lon'],
-            'units': 'Pa'
-        },
-        'air_pressure_on_interface_levels': {
-            'dims': ['interface_levels', 'lat', 'lon'],
-            'units': 'Pa'
-        },
-        'surface_air_pressure': {
-            'units': 'Pa',
-            'dims': ['lat', 'lon'],
-        },
-        'eastward_wind': {
-            'units': 'm s^-1',
-            'dims': ['mid_levels', 'lat', 'lon'],
-        },
-        'northward_wind': {
-            'units': 'm s^-1',
-            'dims': ['mid_levels', 'lat', 'lon'],
-        },
-        'divergence_of_wind': {
-            'units': 's^-1',
-            'dims': ['mid_levels', 'lat', 'lon'],
-        },
-        'atmosphere_relative_vorticity': {
-            'units': 's^-1',
-            'dims': ['mid_levels', 'lat', 'lon'],
-        },
-        'surface_geopotential': {
-            'units': 'm^2 s^-2',
-            'dims': ['lat', 'lon'],
-        },
-    }
-
-    _gfs_output_properties = {
-        'air_temperature': {'units': 'degK'},
-        'air_pressure': {
-            'dims': ['mid_levels', 'lat', 'lon'],
-            'units': 'Pa'
-        },
-        'air_pressure_on_interface_levels': {
-            'dims': ['interface_levels', 'lat', 'lon'],
-            'units': 'Pa'
-        },
-        'surface_air_pressure': {'units': 'Pa'},
-        'eastward_wind': {'units': 'm s^-1'},
-        'northward_wind': {'units': 'm s^-1'},
-        'divergence_of_wind': {'units': 's^-1'},
-        'atmosphere_relative_vorticity': {'units': 's^-1'},
-    }
-
-    _gfs_diagnostic_properties = {}
-
-    input_properties = None
-    output_properties = None
-    diagnostic_properties = None
-
-    uses_tracers = True
-    tracer_dims = ('tracer', 'mid_levels', 'lat', 'lon')
-    prepend_tracers = (('specific_humidity', 'kg/kg'),)
-
-    @property
-    def spectral_names(self):
-        return (
-            'eastward_wind', 'northward_wind', 'air_temperature',
-            'surface_air_pressure') + get_tracer_names()
-
-    def __init__(
-            self,
-            tendency_component_list=None,
-            number_of_damped_levels=0,
-            damping_timescale=2.*86400,
-            moist=True,
-            zero_negative_moisture=True,
-    ):
-        """
-        Initialise the GFS dynamical core.
-
-        Args:
-            tendency_component_list (list of TendencyComponent, optional):
-                The TendencyComponent objects to use for spectral time stepping.
-
-            number_of_damped_levels (int, optional):
-                The number of levels from the model top which are Rayleigh damped.
-
-            damping_timescale (float, optional):
-                The damping timescale in :math:`s` to use for top-of-model Rayleigh damping.
-
-            zero_negative_moisture (bool, optional):
-                If True, all negative values of moisture will be set to zero
-                before tracers are returned.
-        """
-        tendency_component_list = tendency_component_list or []
-        self._tendency_component = ImplicitTendencyComponentComposite(*tendency_component_list)
-        bad_diagnostics = set(
-            self._tendency_component.diagnostic_properties.keys()).intersection(
-            self.spectral_names)
-        if len(bad_diagnostics) > 0:
-            raise GFSError(
-                'Cannot use TendencyComponent components that produce {} as diagnostic '
-                'output as these must only be stepped spectrally.'.format(
-                    bad_diagnostics,
-                )
-            )
-        bad_diagnostics = set(
-            self._tendency_component.diagnostic_properties.keys()).intersection(
-            self._gfs_diagnostic_properties.keys())
-        if len(bad_diagnostics) > 0:
-            raise GFSError(
-                'Cannot use TendencyComponent components that produce {} as diagnostic'
-                ' output as these are already diagnosed by GFS.'.format(
-                    bad_diagnostics,
-                )
-            )
-
-        self._update_constants()
-
-        self._damping_levels = number_of_damped_levels
-        self._tau_damping = damping_timescale
-        self._zero_negative_moisture = zero_negative_moisture
-
-        self.initialized = False
-
-        self.input_properties = self._gfs_input_properties.copy()
-        self.output_properties = self._gfs_output_properties.copy()
-        self.diagnostic_properties = self._gfs_diagnostic_properties.copy()
-        super(GFSDynamicalCore, self).__init__()
-        self.input_properties.update(
-            get_valid_properties(
-                self._gfs_input_properties, self._tendency_component.input_properties, 'input')
-        )
-        self.output_properties.update(
-            get_valid_properties(
-                self._gfs_output_properties, self._tendency_component.tendency_properties, 'output')
-        )
-        self.diagnostic_properties.update(
-            get_valid_properties(
-                self._gfs_diagnostic_properties, self._tendency_component.diagnostic_properties, 'diagnostic')
-        )
-
-    def _update_constants(self):
-        self._radius = get_constant('planetary_radius', 'm')
-        self._omega = get_constant('planetary_rotation_rate', 's^-1')
-        self._R = get_constant('universal_gas_constant', 'J/mole/K')
-        self._Rd = get_constant('gas_constant_of_dry_air', 'J/kg/K')
-        self._Rv = get_constant('gas_constant_of_vapor_phase', 'J/kg/K')
-        self._g = get_constant('gravitational_acceleration', 'm/s^2')
-        self._Cp = get_constant('heat_capacity_of_dry_air_at_constant_pressure', 'J/kg/K')
-        self._Cvap = get_constant('heat_capacity_of_vapor_phase', 'J/kg/K')
-        self._fvirt = (1 - self._Rd/self._Rv)/(self._Rd/self._Rv)
-        self._dry_pressure = get_constant('reference_air_pressure', 'Pa')
-        self._toa_pressure = get_constant('top_of_model_pressure', 'Pa')
-
-    def _initialize_model(self, state, timestep):
-        assert not self.initialized
-        self._time_step = timestep
-        self._num_tracers = state['tracers'].shape[0]
-        self._num_levs, self._num_lats, self._num_lons = state['air_temperature'].shape
-
-        self._truncation = int(self._num_lons/3 - 2)
-
-        self._spectral_dim = int(
-            (self._truncation + 1)*(self._truncation + 2)/2)
-
-        _gfs_dynamics.set_time_step(
-            self._time_step.total_seconds())
-
-        _gfs_dynamics.set_constants(
-            self._radius, self._omega, self._R, self._Rd, self._Rv, self._g,
-            self._Cp, self._Cvap)
-
-        _gfs_dynamics.set_model_grid(
-            self._num_lats,
-            self._num_lons,
-            self._num_levs,
-            self._truncation,
-            self._spectral_dim,
-            self._num_tracers,
-            state['a_coord'][::-1],
-            state['b_coord'][::-1],
-        )
-
-        logging.info('Initialising dynamical core, this could take some time...')
-
-        gaussian_weights, area_weights, latitudes, longitudes = \
-            _gfs_dynamics.init_model(
-                self._dry_pressure,
-                self._damping_levels,
-                self._tau_damping,
-                self._toa_pressure)
-
-        np.testing.assert_allclose(latitudes*180./np.pi, state['latitude'])
-        np.testing.assert_allclose(longitudes*180./np.pi, state['longitude'])
-
-        logging.info('Done!')
-
-        self._hash_u = 1000
-        self._hash_v = 1000
-        self._hash_temperature = 1000
-        self._hash_surface_pressure = 1000
-        self._hash_tracers = 1000
-
-        self.initialized = True
-
-    def __call__(self, state, timestep):
-        """
-        Gets diagnostics from the current model state and steps the state
-        forward in time according to the timestep.
-
-        Args
-        ----
-        state : dict
-            A model state dictionary satisfying the input_properties of this
-            object.
-        timestep : timedelta
-            The amount of time to step forward.
-
-        Returns
-        -------
-        diagnostics : dict
-            Diagnostics from the timestep of the input state.
-        new_state : dict
-            A dictionary whose keys are strings indicating
-            state quantities and values are the value of those quantities
-            at the timestep after input state.
-
-        Raises
-        ------
-        KeyError
-            If a required quantity is missing from the state.
-        InvalidStateError
-            If state is not a value input for the Stepper instance
-            for other reasons.
-        """
-        self._check_self_is_initialized()
-        self._input_checker.check_inputs(state)
-        raw_state = get_numpy_arrays_with_properties(state, self.input_properties)
-        if self.uses_tracers:
-            raw_state['tracers'] = self._tracer_packer.pack(state)
-        raw_state['time'] = state['time']
-
-        tendencies, diagnostics = self._tendency_component(state, timestep)
-        for name, value in tendencies.items():
-            if name in self.input_properties.keys():
-                tendencies[name] = value.to_units(self.input_properties[name]['units'] + ' s^-1')
-        raw_diagnostics, raw_new_state = self.array_call(
-            raw_state, timestep, prognostic_tendencies=tendencies)
-        if self.uses_tracers:
-            new_state = self._tracer_packer.unpack(raw_new_state.pop('tracers'), state)
-        else:
-            new_state = {}
-        if self.tendencies_in_diagnostics:
-            self._insert_tendencies_to_diagnostics(
-                raw_state, raw_new_state, timestep, raw_diagnostics)
-
-        diagnostics.update(restore_data_arrays_with_properties(
-            raw_diagnostics, self.diagnostic_properties,
-            state, self.input_properties, ignore_missing=True))
-        new_state.update(restore_data_arrays_with_properties(
-            raw_new_state, self.output_properties,
-            state, self.input_properties, ignore_missing=True))
-
-        gfs_output_quantities = list(self._gfs_output_properties.keys())
-        for tracer in self.prepend_tracers:
-            gfs_output_quantities.append(tracer[0])
-        remaining = set(tendencies.keys()).difference(gfs_output_quantities)
-
-        for name in remaining:
-            new_state[name] = state[name] + tendencies[name]*timestep.total_seconds()
-        for key in state.keys():
-            if key not in new_state:
-                new_state[key] = state[key]
-
-        check_new_state = {
-            name: quantity
-            for (name, quantity) in new_state.items()
-            if name in self.output_properties.keys()
-        }
-        self._diagnostic_checker.check_diagnostics(diagnostics)
-        self._output_checker.check_outputs(check_new_state)
-        return diagnostics, new_state
-
-    @ensure_contiguous_state
-    def array_call(self, state, timestep, prognostic_tendencies=None):
-        """ Step the dynamical core by one step
-
-        Args:
-            state (dict): The state dictionary of numpy arrays.
-
-        Returns:
-            new_state, diagnostics (dict):
-                The new state and associated diagnostics.
-        """
-        prognostic_tendencies = prognostic_tendencies or {}
-        self._update_constants()
-        nlev, nlat, nlon = state['air_temperature'].shape
-        if nlat < 16:
-            raise GFSError('GFS requires at least 16 latitudes.')
-        if nlon < 12:
-            raise GFSError('GFS requires at least 12 longitudes')
-        if not self.initialized:
-            self._initialize_model(state, timestep)
-        if nlev != self._num_levs:
-            raise GFSError(
-                'Number of vertical levels may not change between successive '
-                'calls to GFS. Last time was {}, this time is {}'.format(
-                    self._num_levs, nlev)
-            )
-        if nlat != self._num_lats:
-            raise GFSError(
-                'Number of latitudes may not change between successive '
-                'calls to GFS. Last time was {}, this time is {}'.format(
-                    self._num_lats, nlat)
-            )
-        if nlon != self._num_lons:
-            raise GFSError(
-                'Number of longitudes may not change between successive '
-                'calls to GFS. Last time was {}, this time is {}'.format(
-                    self._num_lons, nlon)
-            )
-        if timestep.total_seconds() != self._time_step.total_seconds():
-            raise GFSError(
-                'GFSDynamicalCore can only be run with a constant timestep.'
-            )
-
-        outputs = initialize_numpy_arrays_with_properties(
-            self.output_properties, state, self.input_properties,
-            prepend_tracers=self.prepend_tracers,
-            tracer_dims=self.tracer_dims,
-        )
-
-        lnsp = np.log(state['surface_air_pressure'])
-        t_virt = state['air_temperature']*(
-            1 + self._fvirt*state['tracers'][0, :, :, :])
-
-        outputs['air_pressure_on_interface_levels'][:] = (
-            state['air_pressure_on_interface_levels'][::-1, :, :])
-        for name in (
-                'air_pressure', 'tracers', 'eastward_wind', 'northward_wind',
-                'divergence_of_wind', 'atmosphere_relative_vorticity',
-                'surface_air_pressure',):
-            if np.product(outputs[name].shape) > 0:
-                outputs[name][:] = state[name]
-
-        _gfs_dynamics.assign_grid_arrays(
-            outputs['eastward_wind'],
-            outputs['northward_wind'],
-            t_virt,
-            lnsp,
-            outputs['tracers'],
-            outputs['atmosphere_relative_vorticity'],
-            outputs['divergence_of_wind'])
-
-        _gfs_dynamics.assign_pressure_arrays(
-            outputs['surface_air_pressure'],
-            outputs['air_pressure'],
-            outputs['air_pressure_on_interface_levels'])
-
-        _gfs_dynamics.set_topography(state['surface_geopotential'])
-
-        _gfs_dynamics.calculate_pressure()
-
-        np.testing.assert_allclose(outputs['air_pressure'], state['air_pressure'])
-        np.testing.assert_allclose(
-            outputs['air_pressure_on_interface_levels'][::-1, :, :],
-            state['air_pressure_on_interface_levels']
-        )
-
-        self._update_spectral_arrays(state)
-
-        tendency_arrays = \
-            self._get_tendency_arrays(
-                prognostic_tendencies, state['air_temperature'].shape)
-
-        # see Pg. 12 in gfsModelDoc.pdf
-        virtual_temp_tend = tendency_arrays['air_temperature']*(
-            1 + self._fvirt*state['tracers'][0, :, :, :]) + \
-            self._fvirt*t_virt*tendency_arrays['tracers'][0, :, :, :]
-
-        # dlnps/dt = (1/ps)*dps/dt
-        lnps_tend = ((1. / state['surface_air_pressure']) *
-                     tendency_arrays['surface_air_pressure'])
-
-        _gfs_dynamics.assign_tendencies(
-            tendency_arrays['eastward_wind'],
-            tendency_arrays['northward_wind'],
-            virtual_temp_tend,
-            lnps_tend,
-            tendency_arrays['tracers'])
-
-        _gfs_dynamics.take_one_step()
-        _gfs_dynamics.convert_to_grid()
-        _gfs_dynamics.calculate_pressure()
-
-        if self._zero_negative_moisture:
-            set_negatives_to_zero(outputs['tracers'][0, :, :, :])
-        outputs['air_temperature'][:] = t_virt/(
-            1 + self._fvirt*outputs['tracers'][0, :, :, :])
-
-        outputs['air_pressure_on_interface_levels'][:] = \
-            outputs['air_pressure_on_interface_levels'][::-1, :, :]
-
-        outputs['time'] = state['time']
-
-        return {}, outputs
-
-    def _get_tendency_arrays(self, tendencies, T_shape):
-        out_arrays = {}
-        out_arrays.update(self._get_3d_tendencies(tendencies, T_shape))
-        out_arrays.update(self._get_2d_tendencies(tendencies, T_shape))
-        out_arrays['tracers'] = self._get_tracer_tendencies(tendencies, T_shape)
-        return out_arrays
-
-    def _get_3d_tendencies(self, tendencies, T_shape):
-        return self._get_tendencies(
-            tendencies,
-            ['air_temperature', 'eastward_wind', 'northward_wind'],
-            ['mid_levels', 'lat', 'lon'],
-            T_shape
-        )
-
-    def _get_2d_tendencies(self, tendencies, T_shape):
-        return self._get_tendencies(
-            tendencies,
-            ['surface_air_pressure'],
-            ['lat', 'lon'],
-            T_shape[1:]
-        )
-
-    def _get_tendencies(self, tendencies, quantity_names, dims, shape):
-        out_arrays = {}
-        for name in quantity_names:
-            if name in tendencies:
-                property_dict = {
-                    name: {
-                        'dims': dims,
-                        'units': tendencies[name].attrs['units'],
-                    }
-                }
-                out_arrays.update(
-                    get_numpy_arrays_with_properties(tendencies, property_dict))
-            else:
-                out_arrays[name] = np.zeros(shape)
-        # must broadcast any singleton dimensions
-        for name, current_array in out_arrays.items():
-            if current_array.shape != shape:
-                new_array = np.empty(shape)
-                new_array[:] = current_array
-                out_arrays[name] = new_array
-        return out_arrays
-
-    def _get_tracer_tendencies(self, tendencies, T_shape):
-        return_array = np.empty([self._num_tracers] + list(T_shape))
-        for i, name in enumerate(self._tracer_packer.tracer_names):
-            if name in tendencies:
-                property_dict = {
-                    name: {
-                        'dims': ['mid_levels', 'lat', 'lon'],
-                        'units': tendencies[name].attrs['units'],
-                    }
-                }
-                tend = get_numpy_arrays_with_properties(
-                    tendencies, property_dict)[name]
-                return_array[i, :, :, :] = tend
-            else:
-                return_array[i, :, :, :] = 0.
-        return return_array
-
-    def _update_spectral_arrays(self, state):
-        """
-        Checks the state to see if any arrays with spectral counterparts have
-        been modified since they were last returned, and if they have will
-        update the specctral counterpart to reflect the new state array.
-        """
-        u_hash = get_hash(state['eastward_wind'])
-        v_hash = get_hash(state['northward_wind'])
-        if (u_hash != self._hash_u or v_hash != self._hash_v):
-            _gfs_dynamics.vrt_div_to_spectral()
-            self._hash_u = u_hash
-            self._hash_v = v_hash
-        T_hash = get_hash(state['air_temperature'])
-        if T_hash != self._hash_temperature:
-            _gfs_dynamics.virtemp_to_spectral()
-            self._hash_temperature = T_hash
-        tracer_hash = get_hash(state['tracers'])
-        if tracer_hash != self._hash_tracers:
-            _gfs_dynamics.tracer_to_spectral()
-            self._hash_tracers = tracer_hash
-        p_surf_hash = get_hash(state['surface_air_pressure'])
-        if p_surf_hash != self._hash_surface_pressure:
-            _gfs_dynamics.lnps_to_spectral()
-            self._hash_surface_pressure = p_surf_hash
-
-    def __del__(self):
-        """ call shutdown in fortran code """
-        logging.info("Cleaning up dynamical core...")
-        _gfs_dynamics.shut_down_model()
-        logging.info("Done!")
-
-
-def get_hash(array):
-    if sys.version_info > (3, 0):
-        return hash(array.data.tobytes())
-    else:
-        array.flags.writeable = False
-        return hash(array.data)
-
-
-def set_negatives_to_zero(array):
-    array[array < 0] = 0
-    return array
-
-
-'''
-def return_tendency_arrays_or_zeros(quantity_list, state, tendencies):
-
-    tendency_list = []
-    for quantity in quantity_list:
-        if quantity in tendencies.keys():
-            tendency_list.append(np.asfortranarray(tendencies[quantity].values))
-        elif quantity in state.keys():
-            tendency_list.append(np.zeros(state[quantity].shape, order='F'))
-        else:
-            raise IndexError("{} not found in input state or tendencies".format(quantity))
-
-    return tendency_list
-'''
diff --git a/climt/_components/gfs/gfsModelDoc.pdf b/climt/_components/gfs/gfsModelDoc.pdf
deleted file mode 100644
index b03313da..00000000
Binary files a/climt/_components/gfs/gfsModelDoc.pdf and /dev/null differ
diff --git a/climt/_lib/GFS/Makefile b/climt/_lib/GFS/Makefile
deleted file mode 100644
index 3f26f61c..00000000
--- a/climt/_lib/GFS/Makefile
+++ /dev/null
@@ -1,89 +0,0 @@
-OMPFLAG		 = -fopenmp
-DEBUGFLAGS   =-fimplicit-none
-GFS_FFLAGS	 = $(FFLAGS) -march=native -ffast-math -pipe -fomit-frame-pointer -fbacktrace
-SHAREFLAGS 	 = #-static
-
-.SUFFIXES: .o .f90 .F90
-
-OBJS=	 params.o \
-	 kinds.o\
-	 dyn_init.o\
-	 dyn_run.o\
-	 physcons.o\
-	 phy_data.o\
-	 spectral_data.o\
-	 grid_data.o\
-	 pressure_data.o\
-	 finalize.o\
-	 dyn_finalize.o\
-	 semimp_data.o\
-	 run.o\
-         phy_run.o\
-         phy_init.o\
-	 phy_finalize.o\
-	 shtns.o
-	 #init.o\
-
-all: $(python_lib)
-
-python_lib : $(OBJS)
-	ar rcs libgfs_dycore.a $(OBJS) 
-
-force_look:
-	true
-
-clean: force_look
-	rm -f $(OBJS) *.mod *genmod* libgfs_dycore.a
-
-
-.f90.o:
-	$(FC) -Wall $(GFS_FFLAGS) $(DEBUGFLAGS) $(SHAREFLAGS) $(OMPFLAG) -c $<
-
-.F90.o:
-	$(FC) -Wall $(GFS_FFLAGS) $(DEBUGFLAGS)  $(SHAREFLAGS) $(OMPFLAG) -c $<
-
-.f.o:
-	$(FC) -Wall $(GFS_FFLAGS) $(DEBUGFLAGS) $(SHAREFLAGS)  $(OMPFLAG) -free -c $<
-
-%.o : %.mod
-
-kinds.o:	kinds.f90
-
-physcons.o:	physcons.f90 kinds.o
-
-constants.o:	constants.f90 kinds.o
-
-shtns.o:	shtns.f90 kinds.o
-
-params.o:	params.f90 kinds.o 
-
-init.o:  	init.f90 params.o dyn_init.o phy_init.o phy_data.o physcons.o kinds.o 
-
-finalize.o:  	finalize.f90 params.o dyn_finalize.o phy_finalize.o 
-
-spectral_data.o:  	spectral_data.f90 kinds.o
-
-phy_data.o:  	phy_data.f90 kinds.o params.o physcons.o
-
-semimp_data.o: semimp_data.f90 kinds.o shtns.o pressure_data.o physcons.o
-
-grid_data.o:  	grid_data.f90 kinds.o params.o
-
-pressure_data.o:  	pressure_data.f90 kinds.o params.o physcons.o grid_data.o
-
-dyn_init.o:  	dyn_init.f90 params.o shtns.o spectral_data.o pressure_data.o physcons.o grid_data.o semimp_data.o 
-		$(FC) -Wall $(GFS_FFLAGS) $(OMPFLAG) -c dyn_init.f90
-
-phy_init.o:     phy_init.f90 phy_data.o
-
-phy_finalize.o:     phy_finalize.f90 
-
-phy_run.o:      phy_run.f90 pressure_data.o shtns.o spectral_data.o physcons.o kinds.o params.o grid_data.o 
-		$(FC) -Wall $(GFS_FFLAGS) $(OMPFLAG) -c phy_run.f90
-	        
-dyn_finalize.o:  dyn_finalize.f90 shtns.o spectral_data.o pressure_data.o grid_data.o 
-
-dyn_run.o:	dyn_run.f90 pressure_data.o shtns.o spectral_data.o physcons.o kinds.o params.o grid_data.o 
-
-run.o:   	run.f90 dyn_run.o phy_run.o kinds.o params.o  
-
diff --git a/climt/_lib/GFS/Makefile.conf b/climt/_lib/GFS/Makefile.conf
deleted file mode 120000
index 0e7359e9..00000000
--- a/climt/_lib/GFS/Makefile.conf
+++ /dev/null
@@ -1,1 +0,0 @@
-Makefile.conf.gfortran
\ No newline at end of file
diff --git a/climt/_lib/GFS/Makefile.conf.gfortran b/climt/_lib/GFS/Makefile.conf.gfortran
deleted file mode 100644
index 83cb5bb1..00000000
--- a/climt/_lib/GFS/Makefile.conf.gfortran
+++ /dev/null
@@ -1,17 +0,0 @@
-#gnu fortran 
-#FC = gfortran
-OMPFLAG = -fopenmp
-###FFLAGS=-O0 -pedantic -fimplicit-none -fbounds-check -fbacktrace -Wall -fcheck-array-temporaries -g
-DEBUGFLAGS=-fimplicit-none
-FFLAGS=-O3 -march=native -ffast-math -pipe -fomit-frame-pointer -fbacktrace -fno-range-check -fPIC
-#FFLAGS=-O3 -march=core2 -ffast-math -pipe -fomit-frame-pointer -fbacktrace -fno-range-check
-FREEFLAG=-ffree-form
-AUTOFLAG=-frecursive
-EXEC          = ../bin/gfs
-LDFLAGS       = $(OMPFLAG)
-SHAREFLAGS = #-static
-#LIBS          = -L/home/$(USER)/lib -lshtns_omp -L/usr/lib/ -llapack -L/usr/local/lib -lfftw3_omp -lfftw3 -lrt -lm  libgfs/libgfs.a w3lib/libw3.a 
-# JOY removing libw3 and libgfs (physics)
-LIBS          =  -lshtns_omp -L/usr/lib/ -llapack -lfftw3_omp -lfftw3 -lrt -lm 
-CFLAGS  = -O2# -DLINUX
-#CC      = gcc
diff --git a/climt/_lib/GFS/Makefile.conf.ifort b/climt/_lib/GFS/Makefile.conf.ifort
deleted file mode 100644
index 65a92fb9..00000000
--- a/climt/_lib/GFS/Makefile.conf.ifort
+++ /dev/null
@@ -1,19 +0,0 @@
-#intel fortran 
-FC            = ifort
-OMPFLAG       = -openmp
-FREEFLAG      = -FR
-AUTOFLAG      = -auto
-#AUTOFLAG      =  -openmp
-FFLAGS        = -O3 -xHost -convert big_endian  
-DEBUGFLAGS    = -warn all -implicitnone -traceback 
-#DEBUGFLAGS    = 
-EXEC          = ../bin/gfs
-LDFLAGS       = -openmp -mkl=sequential 
-# jet
-LIBS          = -L/pan2/projects/gfsenkf/whitaker/lib -lshtns -lfftw3 libgfs/libgfs.a w3lib/libw3.a
-# zeus
-LIBS          = -L/scratch1/portfolios/BMC/fim/whitaker/lib -lshtns-icc -L/scratch1/portfolios/BMC/fim/whitaker/intel/lib -lfftw3 libgfs/libgfs.a w3lib/libw3.a
-# mac
-#LIBS          = -L/Users/$(USER)/lib -lshtns -L/opt/local/lib -lfftw3 libgfs/libgfs.a w3lib/libw3.a
-CFLAGS  = -O2 -DLINUX
-CC      = gcc
diff --git a/climt/_lib/GFS/driver.f90 b/climt/_lib/GFS/driver.f90
deleted file mode 100644
index a236617b..00000000
--- a/climt/_lib/GFS/driver.f90
+++ /dev/null
@@ -1,19 +0,0 @@
-program driver
-
-  ! driver for gfs dynamical core
-
-  use init_mod, only: init
-  use run_mod, only: run
-  use finalize_mod, only: finalize 
-
-  implicit none
-
-  ! read in namelist, initial conditions, 
-  ! initialize physics and dynamics, allocate arrays.
-  call init() 
-  ! run model, write out at specified intervals.
-  call run()
-  ! deallocate arrays.
-  call finalize()
-
-end program driver
diff --git a/climt/_lib/GFS/dyn_finalize.f90 b/climt/_lib/GFS/dyn_finalize.f90
deleted file mode 100644
index b6d23ac0..00000000
--- a/climt/_lib/GFS/dyn_finalize.f90
+++ /dev/null
@@ -1,35 +0,0 @@
-module dyn_finalize
-! finalize dynamics (provides finalize_dyn subroutine)
-
- use shtns, only: shtns_destroy
- use spectral_data, only: destroy_specdata
- use pressure_data, only: destroy_pressdata
- use grid_data, only: destroy_griddata
-! JOY removing stochastic routines 
-! use stoch_data, only: destroy_stochdata
- use semimp_data, only: destroy_semimpdata
- use params, only: explicit, iau
-! JOY remove Incremental Analysis Unit (IAU) forcing                          
-! use iau_module, only: destroy_iau
-
- implicit none
- private
- public :: finalize_dyn
-
- contains
-
- subroutine finalize_dyn()
-    ! call routines to deallocate arrays
-    call destroy_specdata()
-    call destroy_griddata()
-    call destroy_pressdata()
-    if (.not. explicit) call destroy_semimpdata()
-! JOY removing stochastic routines 
-!    call destroy_stochdata()
-    call shtns_destroy()
-! JOY remove Incremental Analysis Unit (IAU) forcing                          
-!    if (iau) call destroy_iau()
- end subroutine finalize_dyn
-
-
-end module dyn_finalize
diff --git a/climt/_lib/GFS/dyn_init.f90 b/climt/_lib/GFS/dyn_init.f90
deleted file mode 100644
index 140a3346..00000000
--- a/climt/_lib/GFS/dyn_init.f90
+++ /dev/null
@@ -1,331 +0,0 @@
-module dyn_init
-! initialize dynamics (including reading in initial conditions).
-! also does IO of spectral data.
-! public subroutines:
-! init_dyn:  allocate spectral arrays (init_specdata),
-!  read in initial conditions, initalize pressure variables (init_pressdata,
-!  get ak,bk from IC file header), compute gradient of orography,
-!  set up linear damping operators (disspec,diff_prof,damp_prof), initialize arrays
-!  for semi-implicit time stepping.
-! wrtout_sig: write out spectral data.
-! readin_sig: read in spectral data.
- use kinds, only: r_kind, r_single, r_double
-! JOY removing sigio dependency
-! use sigio_module, only: sigio_sclose,sigio_swohdc,&
-!  sigio_srohdc,sigio_aldata,sigio_data,sigio_head,sigio_sropen,sigio_srdata,sigio_axdata
- use params, only: &
- nlons,nlats,nlevs,ndimspec,ntrunc,initfile,dry,ndiss,efold,dcmip,polar_opt,&
- heldsuarez,explicit,tstart,idate_start,hdif_fac,hdif_fac2,fshk,ntrac,taustratdamp,&
- ntoz,ntclw,pdryini,massfix,toa_pressure
- !sighead
- !JOY removed sighead to remove sigio dependency
- use shtns, only: shtns_init, spectogrd, grdtospec, getgrad, getvrtdivspec, lap, lats, lons
- use spectral_data, only: vrtspec,divspec,virtempspec,tracerspec,topospec,lnpsspec,&
-                          disspec,diff_prof,dmp_prof,init_specdata
- use pressure_data, only: ak,bk,ck,dbk,bkl,sl,si,psg,prs,pk,dpk,calc_pressdata
- use grid_data, only: init_griddata, dphisdx, dphisdy, phis, ug, vg, virtempg, &
- tracerg,lnpsg,dlnpdtg
-! JOY removing stochastic routines
-! use stoch_data, only: init_stochdata
- use physcons, only: rerth => con_rerth, rd => con_rd, cp => con_cp, &
-                     omega => con_omega, grav => con_g, pi => con_pi, &
-                     fv => con_fvirt
- use semimp_data, only: init_semimpdata
- use iso_c_binding
-
-
- implicit none
-!namelist/dynamics_list/ak,bk
- private
- !public :: init_dyn, wrtout_sig, readin_sig, wrtout_gfsgrb
- !! removing all output calls
- ! JOY remove routines which read initial conditions from file
- !public :: init_dyn, wrtout_sig, readin_sig
- public :: init_dyn
-
- real(r_kind) :: slrd0=0
- contains
- !JOY subroutine to set topography from arbitrary array
- subroutine set_topography(surf_geop) bind(c, name='gfs_set_topography')
-
-    real(r_kind), intent(in), dimension(nlons,nlats) :: surf_geop
-
-    call grdtospec(surf_geop/grav, topospec)
-    call spectogrd(grav*topospec, phis)
-    call getgrad(grav*topospec, dphisdx, dphisdy, rerth)
-
- end subroutine set_topography
-
- subroutine init_dyn(py_damping_level, py_taustratdamp) bind(c, name='gfs_init_dynamics')
-    integer k
-    ! allocate arrays
-    ! JOY data array inits will be done from python
-
-    !call init_specdata()
-    !call init_griddata()
-    integer(c_int), intent(in):: py_damping_level
-    real(c_double), intent(in):: py_taustratdamp
-
-    ! initialize spherical harmonic lib
-    !print *,'Initialising shtns'
-    call shtns_init(nlons,nlats,ntrunc,nthreads=1,polar_opt=polar_opt)
-        ! convert to ln(ps) in Pa.
-    !print *,'calling spectogrd'
-    !print *, 'Grid size', nlats, nlons, ntrunc
-    !print *, shape(lnpsspec)
-    !call spectogrd(lnpsspec, psg)
-    !psg = 1000.*exp(psg) ! convert to Pa
-    !call grdtospec(log(psg), lnpsspec) ! back to spectral.
-
-
-    do k=1,nlevs
-       dbk(k) = bk(k+1)-bk(k)
-       bkl(k) = 0.5*(bk(k+1)+bk(k))
-       ck(k)  = ak(k+1)*bk(k)-ak(k)*bk(k+1)
-    enddo
-         ! compute sigma coordinate quantities (bottom to top).
-    do k=1,nlevs+1
-        si(nlevs+2-k)= (ak(k) - toa_pressure)/(pdryini - toa_pressure)+bk(k) ! si are now sigmas
-    enddo
-    do k=1,nlevs
-        sl(k) = 0.5*(si(k)+si(k+1))
-         !sl(k) = ((si(k)**(rk+1.) - si(k+1)**(rk+1.))/&
-         !        ((rk+1.)*(si(k)-si(k+1))))**(1./rk)
-    enddo
-
-    slrd0 = si(nlevs+1 - py_damping_level)
-    taustratdamp = py_taustratdamp
-
-    !call spectogrd(grav*topospec, phis)
-    !print *,'min/max surface geopotential',minval(phis),maxval(phis)
-    !call spectogrd(lnpsspec, psg)
-    !psg = exp(psg)
-    !print *,'min/max sfc pressure (hPa)',minval(psg/100.),maxval(psg/100.)
-    ! initialize model interface and level pressures, related variables.
-    !call calc_pressdata(lnpsg)
-    ! compute gradient of surface orography
-    !call getgrad(grav*topospec, dphisdx, dphisdy, rerth)
-    ! hyper-diffusion operator (plus rayleigh damping in upper stratosphere)
-
-    call setdampspec(ndiss,efold,hdif_fac,hdif_fac2,fshk,disspec,diff_prof,dmp_prof)
-
-    ! initialize arrays for semi-implicit adjustments.
-    if (.not. explicit) call init_semimpdata()
-    ! initialize stochastic data.
-    ! JOY removing stochastic routines
-    !call init_stochdata()
- end subroutine init_dyn
-
- subroutine copyspecin(rspecdata,cspecdata)
-    real(r_single), intent(in) :: rspecdata(2*ndimspec)
-    complex(r_kind), intent(out) :: cspecdata(ndimspec)
-    integer n,nn
-    nn = 1
-    ! factor of sqrt(2.*pi) accounts for difference in normalization
-    ! between ncep libs and shtns (which uses orthonormalized norm)
-    do n=1,ndimspec
-       cspecdata(n) = sqrt(2.*pi)*cmplx(rspecdata(nn),rspecdata(nn+1))
-       nn = nn + 2
-    enddo
- end subroutine copyspecin
-
- subroutine copyspecout(cspecdata, rspecdata)
-    real(r_single), intent(out) :: rspecdata(2*ndimspec)
-    complex(r_kind), intent(in) :: cspecdata(ndimspec)
-    integer n,nn
-    nn = 1
-    ! factor of sqrt(2.*pi) accounts for difference in normalization
-    ! between ncep libs and shtns (which uses orthonormalized norm)
-    do n=1,ndimspec
-       rspecdata(nn) = real(cspecdata(n))/sqrt(2.*pi)
-       rspecdata(nn+1) = imag(cspecdata(n))/sqrt(2.*pi)
-       nn = nn + 2
-    enddo
- end subroutine copyspecout
-
-
- subroutine twodtooned(data2,data1)
-   real(r_kind), intent(in) :: data2(nlons,nlats)
-   real(r_kind), intent(out) :: data1(nlons*nlats)
-   integer i,j,n
-   do n=1,nlons*nlats
-      j = 1+(n-1)/nlons
-      i = n-(j-1)*nlons
-      data1(n) = data2(i,j)
-   enddo
- end subroutine twodtooned
-
-! subroutine dcmip_ics(testcase,dry)
-!   ! DCMIP 2012 test case initial conditions (test cases 4 and 5)
-!   ! (idealized baroclinic instability or trop cyclone test case initial conditions)
-!   use dcmip_initial_conditions_test_4, only: test4_baroclinic_wave
-!   use dcmip_initial_conditions_test_5, only: test5_tropical_cyclone
-!   integer, intent(in) :: testcase
-!   logical, intent(in) :: dry
-!   real(8) lon,lat,z,X,p,u,v,w,t,phi_s,ps,rho,q,q1,q2
-!   integer moist,zcoords,i,j,k,nt
-!
-!   if (testcase/10 .eq. 5) then
-!      if (dry) then
-!         print *,'cannot run dcmip trop cyclone test case with dry=.true.'
-!         stop
-!      endif
-!      print *,'replacing initial conds with dcmip tropical cyclone case..'
-!   else if (testcase/10 .eq. 4) then
-!      print *,'replacing initial conds with dcmip baroclinic wave test case..'
-!   else
-!      print *,'dcmip must be 40 something or 50 something'
-!      stop
-!   endif
-!   zcoords = 0 ! p, not z, is specified.
-!   if (dry) then
-!     moist = 0
-!   else
-!     moist = 1
-!   endif
-!   if (testcase/10 .eq. 4) psg = 1.e5
-!   if (testcase/10 .eq. 5) then
-!      psg = 1.015e5
-!      lnpsg = log(psg)
-!      call grdtospec(lnpsg,lnpsspec)
-!      call calc_pressdata(lnpsg)
-!      ! run once, just to calculate ps
-!      ! ps then used to compute pressure on model levels for subsequent call
-!      do k=1,nlevs
-!      do j=1,nlats
-!      do i=1,nlons
-!         lon = lons(i,j); lat = lats(i,j); p = prs(i,j,k)
-!         call test5_tropical_cyclone (lon,lat,p,z,zcoords,u,v,w,t,phi_s,ps,rho,q)
-!         psg(i,j) = ps
-!      enddo
-!      enddo
-!      enddo
-!   endif
-!   ! calculate pressure on model levels given ps.
-!   lnpsg = log(psg)
-!   call grdtospec(lnpsg,lnpsspec)
-!   call calc_pressdata(lnpsg)
-!   ! calculate u,v,t,q,phis and other tracers.
-!   do k=1,nlevs
-!      do j=1,nlats
-!      do i=1,nlons
-!         lon = lons(i,j); lat = lats(i,j); p = prs(i,j,k)
-!         if (testcase/10 .eq. 4) then
-!            call test4_baroclinic_wave (moist,X,lon,lat,p,z,zcoords,u,v,w,t,phi_s,ps,rho,q,q1,q2)
-!         else if (testcase/10 .eq. 5) then
-!            call test5_tropical_cyclone (lon,lat,p,z,zcoords,u,v,w,t,phi_s,ps,rho,q)
-!         endif
-!         ug(i,j,k) = u; vg(i,j,k) = v; phis(i,j) = phi_s
-!         if (dry) then
-!            virtempg(i,j,k) = t
-!            tracerg(i,j,k,1) = 0.
-!            tracerg(i,j,k,2) = q1
-!            tracerg(i,j,k,3) = q2
-!         else
-!            virtempg(i,j,k) = t*(1. + fv*q)
-!            tracerg(i,j,k,1) = q
-!            if (testcase/10 .eq. 4) then
-!               tracerg(i,j,k,2) = q1
-!               tracerg(i,j,k,3) = q2
-!            else if (testcase .eq. 52) then
-!               tracerg(i,j,k,2:ntrac) = 0.
-!            endif
-!         endif
-!       enddo
-!       enddo
-!       call getvrtdivspec(ug(:,:,k),vg(:,:,k),vrtspec(:,k),divspec(:,k),rerth)
-!       call grdtospec(virtempg(:,:,k),virtempspec(:,k))
-!       do nt=1,ntrac
-!          call grdtospec(tracerg(:,:,k,nt),tracerspec(:,k,nt))
-!       enddo
-!   enddo
-!   call grdtospec(phis/grav, topospec)
-!
-! end subroutine dcmip_ics
-
-
- subroutine heldsuarez_ics(toa_pressure)
-   ! jablonowski and williamson (2006, QJR, p. 2943, doi: 10.1256/qj.06.12)
-   ! initial perturbation on an isothermal state.
-   real(r_kind), dimension(nlons,nlats) :: rnh,xnh,rsh,xsh, rand_vel
-   real(r_kind) :: lonc,latc,up,pertrad
-   real(c_double), intent(in):: toa_pressure
-   integer k
-   print *,'replacing initial conds with held and suarez test case..'
-   lonc = pi/9.
-   up = 1.
-   pertrad = rerth/10.
-   !JOY changing ICs
-   latc = 0.
-   !latc = 2.*pi/9.
-   !xnh = sin(latc)*sin(lats) + cos(latc)*cos(lats)*cos(lons-lonc)
-   xnh = sin(latc)*sin(lats) + cos(latc)*cos(lats)
-   latc = -2.*pi/9.
-   !xsh = sin(latc)*sin(lats) + cos(latc)*cos(lats)*cos(lons-lonc)
-   xsh = sin(latc)*sin(lats) + cos(latc)*cos(lats)
-   rnh = rerth*acos(xnh)
-   rsh = 0.*rerth*acos(xsh)
-   virtempg = 300.  ! isothermal state.
-   vg = 0.
-   psg = 1.e5
-   lnpsg = log(psg)
-   call grdtospec(lnpsg,lnpsspec)
-   call calc_pressdata(lnpsg)
-   do k=1,nlevs
-      ! add a barotropic zonal wind perturbation (opp sign in each hemisphere)
-      !JOY changed to single tropical jet
-      call random_number(rand_vel)
-      ug(:,:,k) = 0.8*rand_vel
-      !ug(:,:,k) = up*(exp(-(rnh/pertrad)**2)+exp(-(rsh/pertrad)**2))
-      call getvrtdivspec(ug(:,:,k),vg(:,:,k),vrtspec(:,k),divspec(:,k),rerth)
-      call grdtospec(virtempg(:,:,k),virtempspec(:,k))
-   enddo
-   topospec = 0.
- end subroutine heldsuarez_ics
-
- subroutine setdampspec(ndiss,efold,hdif_fac,hdif_fac2,fshk,disspec,diff_prof,dmp_prof)
-   ! set hyper-diffusion and upper level rayleigh damping parameters/structure.
-   ! if efold and/or ndiss are zero, default GFS parameters are used.
-   integer, intent(inout) :: ndiss
-   real(r_kind), intent(inout) :: efold,fshk
-   real(r_kind), intent(in) :: hdif_fac,hdif_fac2
-   real(r_kind), intent(out),dimension(ndimspec) :: disspec
-   real(r_kind), dimension(nlevs),intent(out) :: diff_prof,dmp_prof
-   integer k
-   real(r_kind) dmp_prof1
-   ! if ndiss=0, use default value of 8
-   if (ndiss == 0) ndiss = 8
-   ! if efold <= 0, use GFS defaults.
-   if (efold .le. 0) then
-      if (ntrunc > 170) then
-         efold = 3600./(hdif_fac2*(ntrunc/170.)**4*1.1)
-      else if (ntrunc == 126 .or. ntrunc == 170) then
-         efold = 1./(hdif_fac2*12.E15/(RERTH**4)*(80.*81.)**2)
-      else
-         efold = 1./(hdif_fac2*3.E15/(RERTH**4)*(80.*81.)**2)
-      end if
-      efold = 2.*efold ! mysterious factor 2 in deldifs.f
-      !print *,'ndiss,efold =',ndiss,efold
-   end if
-   if (fshk .le. 0) then
-      ! factor to multiply diffusion per scale height.
-      fshk = 1.0*hdif_fac
-      if (ntrunc > 170) fshk = 2.2*hdif_fac
-      if (ntrunc == 126) fshk = 1.5*hdif_fac
-   end if
-   !slrd0=0.00        ! SIGMA LEVEL AT WHICH TO BEGIN RAYLEIGH MOMTUM DAMPING
-   dmp_prof1=1./taustratdamp ! RECIPROCAL OF TIME SCALE PER SCALE HEIGHT
-                      ! ABOVE BEGINNING SIGMA LEVEL FOR RAYLEIGH DAMPING
-   dmp_prof = 0.
-   diff_prof = sl**(log(1./fshk))
-   !print *,'profiles for diffusion and linear damping:'
-   !print *,'(level, sigma, diffusion enhancment, linear drag coeff)'
-   do k=1,nlevs
-      if (sl(k) .lt. slrd0) dmp_prof(k)=dmp_prof1*log(slrd0/sl(k))
-   !   print *,k,sl(k),diff_prof(k),dmp_prof(k)
-   enddo
-   disspec(:) = -(1./efold)*(lap(:)/minval(lap))**(ndiss/2)
-   return
- end subroutine setdampspec
-
-end module dyn_init
diff --git a/climt/_lib/GFS/dyn_run.f90 b/climt/_lib/GFS/dyn_run.f90
deleted file mode 100644
index 143f5885..00000000
--- a/climt/_lib/GFS/dyn_run.f90
+++ /dev/null
@@ -1,748 +0,0 @@
- module dyn_run
-! calculation of model tendencies.
-! Public subroutines:
-! getdyntend:  dynamics tendencies (not including linear damping and diffusion,
-! called in run.f90)
-! used by getdyntend:
-! getomega: calculate dlnpdt (omega/p)
-! getpresgrad: calculate pressure gradient terms in momentum eqns.
-! getvadv: calculate vertical advection terms.
-
- use params, only: nlevs,ntrunc,nlons,nlats,ndimspec,dt,ntrac,pdryini,&
-   dcmip,dry,vcamp,svc
- use kinds, only: r_kind,r_double
- use shtns, only: degree,order,areawts,&
- lap,invlap,lats,grdtospec,spectogrd,getuv,getvrtdivspec,getgrad,areawts
- use spectral_data, only:  lnpsspec, vrtspec, divspec, virtempspec,&
- tracerspec, topospec, disspec, dmp_prof, diff_prof
- use grid_data, only: ug,vg,vrtg,virtempg,divg,tracerg,dlnpdtg,etadot,lnpsg,&
- phis,dphisdx,dphisdy,dlnpsdt
- use pressure_data, only:  ak, bk, ck, dbk, dpk, rlnp, pk, alfa, dpk, psg,&
- calc_pressdata
-! JOY removing stochastic routines
-! use stoch_data, only:  grd_svc, vfact_svc
- use physcons, only: rerth => con_rerth, rd => con_rd, cp => con_cp, &
-               eps => con_eps, omega => con_omega, cvap => con_cvap, &
-               grav => con_g, fv => con_fvirt, kappa => con_rocp
- use iso_c_binding
- implicit none
- private
-
- public :: getdyntend, getomega, getpresgrad, getvadv, dry_mass_fixer,&
-            gfsConvertToGrid, gfsConvertToSpec, gfs_uv_to_vrtdiv
-
- contains
-
-!JOY adding subroutine to convert u,v grids to vrt,div grids
-subroutine gfs_uv_to_vrtdiv(ug1,vg1,vrtg1,divg1) bind(c,name='gfs_uv_to_vrtdiv')
-
-    real(c_double), dimension(nlons,nlats,nlevs), intent(in) :: ug1, vg1
-    real(c_double), dimension(nlons,nlats,nlevs), intent(out) :: vrtg1, divg1
-    complex(c_double), dimension(:,:), allocatable :: vrtspec1, divspec1
-    integer k
-
-    allocate(vrtspec1(ndimspec,nlevs),divspec1(ndimspec,nlevs))
-!$omp parallel do private(k)
-    do k=1,nlevs
-
-        call getvrtdivspec(ug1(:,:,k),vg1(:,:,k),vrtspec1(:,k),divspec1(:,k),rerth);
-    enddo
-
-!$omp end parallel do
-!$omp parallel do private(k)
-    do k=1,nlevs
-
-        call spectogrd(vrtspec1(:,k),vrtg1(:,:,k))
-        call spectogrd(divspec1(:,k),divg1(:,:,k))
-    enddo
-
-!$omp end parallel do
-
-
-
-    deallocate(vrtspec1,divspec1)
-
-end subroutine gfs_uv_to_vrtdiv
-
-!JOY adding wrapper functions to convert to and from spectral/grid arrays
- subroutine gfsConvertToGrid() bind(c,name='gfs_convert_to_grid')
-
-! JOY adding dummy variables for wrapper functions. Not efficient to allocate
-! these variables every time, but it is the cleanest way.
-
-    complex(r_kind), dimension(:,:), allocatable :: dummy3d
-    complex(r_kind), dimension(:,:,:), allocatable :: dummy4d
-    complex(r_kind), dimension(:), allocatable :: dummy2d
-
-    allocate(dummy3d(ndimspec,nlevs))
-    allocate(dummy4d(ndimspec,nlevs,ntrac))
-    allocate(dummy2d(ndimspec))
-
-    call getdyntend(dummy3d,dummy3d,dummy3d,dummy4d,dummy2d,1,.True.)
-
-    deallocate(dummy2d,dummy3d,dummy4d)
-
- end subroutine gfsConvertToGrid
-
- subroutine gfsConvertToSpec() bind(c,name='gfs_convert_to_spectral')
-
-
-    integer i,k
-
-!$omp parallel do private(k,i) schedule(dynamic)
-    do k=1,nlevs
-
-        call getvrtdivspec(ug(:,:,k),vg(:,:,k),vrtspec(:,k),divspec(:,k),rerth)
-        call grdtospec(virtempg(:,:,k),virtempspec(:,k))
-
-        do i=1,ntrac
-            call grdtospec(tracerg(:,:,k,i),tracerspec(:,k,i))
-        enddo
-    enddo
-!$omp end parallel do
- call grdtospec(lnpsg,lnpsspec)
- end subroutine
-
- subroutine lnpsToSpec() bind(c, name='gfs_lnps_to_spectral')
-    call grdtospec(lnpsg,lnpsspec)
- end subroutine
-
- subroutine vrtdivToSpec() bind(c, name='gfs_vrt_div_to_spectral')
-     integer k
-!$omp parallel do private(k) schedule(dynamic)
-     do k=1, nlevs
-        call getvrtdivspec(ug(:,:,k),vg(:,:,k),vrtspec(:,k),divspec(:,k),rerth)
-     enddo
-!$omp end parallel do
- end subroutine
-
- subroutine virtempToSpec() bind(c, name='gfs_virtemp_to_spectral')
-     integer k
-!$omp parallel do private(k) schedule(dynamic)
-     do k=1, nlevs
-        call grdtospec(virtempg(:,:,k),virtempspec(:,k))
-     enddo
-!$omp end parallel do
- end subroutine
-
- subroutine tracerToSpec() bind(c, name='gfs_tracer_to_spectral')
-     integer k, i
-!$omp parallel do private(k,i) schedule(dynamic)
-     do k=1, nlevs
-        call grdtospec(virtempg(:,:,k),virtempspec(:,k))
-        do i=1,ntrac
-            call grdtospec(tracerg(:,:,k,i),tracerspec(:,k,i))
-        enddo
-     enddo
-!$omp end parallel do
- end subroutine
-
- subroutine getdyntend(dvrtspecdt,ddivspecdt,dvirtempspecdt,&
-                       dtracerspecdt,dlnpsspecdt,rkstage,just_do_inverse_transform)
-   ! compute dynamics tendencies (not including hyper-diffusion.and linear drag,
-   ! which are treated implicitly)
-   ! based on hybrid sigma-pressure dynamical core described in
-   ! http://www.emc.ncep.noaa.gov/officenotes/newernotes/on462.pdf
-   ! The only difference is that here we use a forward in time
-   ! runge-kutta scheme instead of robert-assellin time filtered leap-frog.
-   logical, optional, intent(in) :: just_do_inverse_transform
-   logical :: early_return = .false. ! if true. spectral-> grid only
-   complex(r_kind), intent(out), dimension(ndimspec,nlevs) :: &
-   dvrtspecdt,ddivspecdt,dvirtempspecdt
-   complex(r_kind), intent(out), dimension(ndimspec,nlevs,ntrac) :: &
-   dtracerspecdt
-   complex(r_kind), intent(out), dimension(ndimspec) :: dlnpsspecdt
-! local variables.
-   complex(r_kind), dimension(:,:),allocatable :: workspec
-   real(r_kind), dimension(ndimspec) :: smoothfact
-   real(r_kind), dimension(:,:), allocatable :: dlnpsdx,dlnpsdy
-   real(r_kind), dimension(:,:,:), allocatable :: &
-   prsgx,prsgy,vadvu,vadvv,vadvt,vadvq,dvirtempdx,dvirtempdy
-   real(r_kind), dimension(:,:,:), allocatable :: dvrtdx,dvrtdy
-   integer, intent(in) :: rkstage ! runge-kutta stage (1,2,3)
-   integer k,nt,ntrac_use
-   logical :: profile = .false. ! print out timing stats
-   real(r_kind) pmoist,pdry,epstiny,rnn1,rnn0
-   real(8) t1,t2,t0
-   integer(8) count, count_rate, count_max
-
-   epstiny = tiny(epstiny)
-
-   ! should tendencies be computed, or just spectral -> grid?
-   if (present(just_do_inverse_transform)) then
-     early_return = .true.
-   else
-     early_return = .false.
-   endif
-
-   ! use alloctable arrays instead of automatic arrays
-   ! for these to avoid segfaults in openmp.
-   allocate(dlnpsdx(nlons,nlats))
-   allocate(dlnpsdy(nlons,nlats))
-   allocate(prsgx(nlons,nlats,nlevs))
-   allocate(prsgy(nlons,nlats,nlevs))
-   allocate(vadvu(nlons,nlats,nlevs))
-   allocate(vadvv(nlons,nlats,nlevs))
-   allocate(vadvt(nlons,nlats,nlevs))
-   allocate(vadvq(nlons,nlats,nlevs))
-   allocate(dvirtempdx(nlons,nlats,nlevs))
-   allocate(dvirtempdy(nlons,nlats,nlevs))
-   allocate(workspec(ndimspec,nlevs))
-   if (vcamp > epstiny .or. svc > epstiny) then
-      ! smoothing factors used to compute vorticity confinement
-      ! anti-diffusion of vorticity.
-      rnn0 = 0.5*ntrunc*(0.5*ntrunc+1)
-      do k=1,ndimspec
-         rnn1 = degree(k)*(degree(k)+1)
-         smoothfact(k)=exp(-(rnn1/rnn0))
-      enddo
-      allocate(dvrtdx(nlons,nlats,nlevs))
-      allocate(dvrtdy(nlons,nlats,nlevs))
-   endif
-
-   ! compute u,v,virt temp, vorticity, divergence, ln(ps)
-   ! and specific humidity on grid from spectral coefficients.
-   call system_clock(count, count_rate, count_max)
-   t1 = count*1.d0/count_rate
-   t0 = t1
-!$omp parallel do private(k,nt) schedule(dynamic)
-   do k=1,nlevs
-      call getuv(vrtspec(:,k),divspec(:,k),ug(:,:,k),vg(:,:,k),rerth)
-      call spectogrd(vrtspec(:,k),vrtg(:,:,k))
-      call spectogrd(divspec(:,k),divg(:,:,k))
-      call spectogrd(virtempspec(:,k),virtempg(:,:,k))
-      ! gradient of virtual temperature on grid.
-      call getgrad(virtempspec(:,k),dvirtempdx(:,:,k),dvirtempdy(:,:,k),rerth)
-      ! specific humidity, other tracers on grid.
-      do nt=1,ntrac
-         call spectogrd(tracerspec(:,k,nt),tracerg(:,:,k,nt))
-      enddo
-      ! vorticity gradient only needed for vorticity confinement.
-      if (vcamp > epstiny .or. svc > epstiny) then
-         call getgrad(vrtspec(:,k),dvrtdx(:,:,k),dvrtdy(:,:,k),rerth)
-      end if
-   enddo
-!$omp end parallel do
-   !print *,'min/max ug',minval(ug),maxval(ug)
-   !print *,'min/max virtempg',minval(virtempg),maxval(virtempg)
-   call system_clock(count, count_rate, count_max)
-   t2 = count*1.d0/count_rate
-   if (profile) print *,'1 time=',t2-t1
-
-   ! lnps on grid.
-   call spectogrd(lnpsspec, lnpsg)
-
-   ! compute pressure on interfaces and model layers using lnpsg
-   ! i.e. calculate alpha,delp,rlnp,dpk,psg,pk from ak,bk,lnpsg
-   ! results stored in module pressure_data
-   call calc_pressdata(lnpsg)
-   !print *,'min/max lnpsg',minval(lnpsg),maxval(lnpsg)
-   !print *,'min/max psg',minval(psg),maxval(psg)
-
-! print out global mean dry ps (only first step in RK scheme).
-   ! mcgibbon: This should be unnecessary and makes a bad assumption that all
-   ! tracers are moisture, so let's just get rid of this.
-   !if (rkstage .eq. 1) then
-   !   pmoist = 0.
-   !   if (.not. dry) then
-   !      ntrac_use = ntrac
-   !      !if (dcmip > 0) ntrac_use = 1
-   !      do k=1,nlevs
-   !         do nt=1,ntrac_use
-   !            pmoist = pmoist + sum(areawts*dpk(:,:,nlevs-k+1)*tracerg(:,:,k,nt))
-   !         enddo
-   !      enddo
-   !   endif
-   !   pdry = sum(areawts*psg) - pmoist
-   !   if (pdryini .eq. 0) then
-   !      pdryini = pdry
-   !      print *,'pdryini = ',pdryini
-   !   else
-   !      !print *, "SUM OF PSG: ", sum(psg);
-   !      !print *,'pdry = ',pdry
-   !   endif
-   !endif
-
-   ! gradient of surface pressure.
-   ! gradient of surface geopotential precomputed in dyn_init,
-   ! saved in module grid_data (dphisdx,dphisdy).
-   call getgrad(lnpsspec, dlnpsdx, dlnpsdy, rerth)
-
-   !print *,'min/max ps',minval(psg/100.),maxval(psg/100.)
-
-   ! get pressure vertical velocity divided by pressure (dlnpdtg),
-   ! tendency of lnps, etadot.
-   ! etadot goes top to bottom, dlnpdtg goes bottom to top
-   call getomega(ug,vg,divg,dlnpsdx,dlnpsdy,dlnpsdt,dlnpdtg,etadot,&
-                 vadvu,vadvv,vadvt,vadvq,prsgx) ! work storage
-
-   ! return before computing tendencies (after updating
-   ! grid data).
-   if (early_return) then
-      deallocate(vadvq,workspec,dvirtempdx,dvirtempdy)
-      deallocate(prsgx,prsgy,vadvu,vadvv,vadvt)
-      if (vcamp > epstiny) deallocate(dvrtdx,dvrtdy)
-      deallocate(dlnpsdx,dlnpsdy)
-      return
-   endif
-
-   ! get tendency of lnps on grid.
-   call grdtospec(dlnpsdt,dlnpsspecdt)
-   ! get pressure gradient terms (prsgx,prsgy)
-   call getpresgrad(virtempg,dvirtempdx,dvirtempdy,dphisdx,dphisdy,dlnpsdx,dlnpsdy,&
-                    prsgx,prsgy,&
-                    vadvu,vadvv,vadvt,vadvq) ! work storage
-
-   ! get vertical advection terms  (input etadot is top to bottom)
-   call getvadv(ug,etadot,vadvu)
-   call getvadv(vg,etadot,vadvv)
-   call getvadv(virtempg,etadot,vadvt)
-
-   ! temporary storage of planetary vorticity
-   dlnpsdx = 2.*omega*sin(lats)
-
-   ! compute energy conversion term, store in vadvq.
-   if (ntrac > 0) then
-      !$omp parallel do private(k)
-      do k=1,nlevs
-         ! section 1.5 in ON 461 (eqn 1.0.3).
-         !vadvq(:,:,k) = &
-         !kappa*dlnpdtg(:,:,k)*virtempg(:,:,k)*&
-         !(1.+fv*tracerg(:,:,k,1))/(1.+(cvap/cp-1.)*tracerg(:,:,k,1))
-         ! GFS appears to missing term in numerator (1. + fv*q)?
-         ! GFS has (in gfidi_hyb.f).  This form is consistent with
-         ! ECMWF IFS documentation.
-         vadvq(:,:,k) = &
-         kappa*dlnpdtg(:,:,k)*virtempg(:,:,k)/(1.+(cvap/cp-1.)*tracerg(:,:,k,1))
-      enddo
-      !$omp end parallel do
-   else
-      !$omp parallel do private(k)
-      do k=1,nlevs
-         vadvq(:,:,k) = kappa*dlnpdtg(:,:,k)*virtempg(:,:,k)
-      enddo
-      !$omp end parallel do
-   endif
-   call system_clock(count, count_rate, count_max)
-   t1 = count*1.d0/count_rate
-   if (profile) print *,'2 time=',t1-t2
-
-   ! compute tendencies of virt temp, vort and div in spectral space
-!$omp parallel do private(k) schedule(dynamic)
-   do k=1,nlevs
-      ! add pressure gradient force to vertical advection terms
-      ! (so prsgy and prsgx can be re-used)
-      vadvv(:,:,k) = vadvv(:,:,k) - prsgy(:,:,k)
-      vadvu(:,:,k) = vadvu(:,:,k) - prsgx(:,:,k)
-      ! virtual temp tendency
-      prsgx(:,:,k) = -ug(:,:,k)*dvirtempdx(:,:,k) - vg(:,:,k)*dvirtempdy(:,:,k) - &
-                      vadvt(:,:,k) + vadvq(:,:,k)
-      call grdtospec(prsgx(:,:,k), dvirtempspecdt(:,k))
-      ! flux terms for vort, div eqns
-      prsgx(:,:,k) = ug(:,:,k)*(vrtg(:,:,k) + dlnpsdx(:,:)) + vadvv(:,:,k)
-      prsgy(:,:,k) = vg(:,:,k)*(vrtg(:,:,k) + dlnpsdx(:,:)) - vadvu(:,:,k)
-      if (vcamp > epstiny .or. svc > epstiny) then
-         ! abs(grad(vrt)) - stored in vadvu
-         vadvu(:,:,k) = sqrt(dvrtdx(:,:,k)**2 + dvrtdy(:,:,k)**2)
-         ! upgradient advection velocity is vcamp*grad(vrt)/abs(grad(vrt))
-         ! (unit normal vector pointing up vorticity gradient scaled by
-         ! vcamp, which has units of velocity)
-         ! vcamp stored in vadvv
-         vadvv(:,:,k) = vcamp
-         ! vcamp modulated by stochastic pattern.
-         ! JOY removing stochastic routines
-         !if (svc > epstiny) vadvv(:,:,k) = (1.+vfact_svc(k)*grd_svc)*vadvv(:,:,k)
-         where (vadvu(:,:,k) > epstiny)
-            dvrtdx(:,:,k) = dvrtdx(:,:,k)/vadvu(:,:,k)
-            dvrtdy(:,:,k) = dvrtdy(:,:,k)/vadvu(:,:,k)
-            prsgx(:,:,k) = prsgx(:,:,k) + vadvv(:,:,k)*dvrtdx(:,:,k)*abs(vrtg(:,:,k))
-            prsgy(:,:,k) = prsgy(:,:,k) + vadvv(:,:,k)*dvrtdy(:,:,k)*abs(vrtg(:,:,k))
-         end where
-      end if
-      call getvrtdivspec(prsgx(:,:,k),prsgy(:,:,k),ddivspecdt(:,k),dvrtspecdt(:,k),rerth)
-      ! flip sign of vort tend.
-      dvrtspecdt(:,k) = -dvrtspecdt(:,k)
-      !if (vcamp > epstiny .or. svc > epstiny) then
-      !   ! add simplified vorticity confinement (anti-diffusion of vorticity)
-      !   ! operates on smoothed vorticity field.
-      !   ! multiply anti-diffusion by stochastic random pattern if svc > 0
-      !   if (svc > epstiny) then
-      !      workspec(:,k) = (lap/rerth**2)*smoothfact*vrtspec(:,k)
-      !      call spectogrd(workspec(:,k), prsgx(:,:,k))
-      !      prsgx(:,:,k) = (vcamp+vfact_svc(k)*grd_svc)*prsgx(:,:,k)
-      !      call grdtospec(prsgx(:,:,k),workspec(:,k))
-      !      dvrtspecdt(:,k) = dvrtspecdt(:,k) - workspec(:,k)
-      !   else ! if not stochastic, no transforms needed (VC is cost-free)
-      !      dvrtspecdt(:,k) = dvrtspecdt(:,k) - vcamp*(lap/rerth**2)*smoothfact*vrtspec(:,k)
-      !   endif
-      !endif
-      ! add laplacian(KE) term to div tendency
-      prsgx(:,:,k) = 0.5*(ug(:,:,k)**2+vg(:,:,k)**2)
-      call grdtospec(prsgx(:,:,k),workspec(:,k))
-      ddivspecdt(:,k) = ddivspecdt(:,k) - &
-      (lap(:)/rerth**2)*workspec(:,k)
-   enddo
-!$omp end parallel do
-   call system_clock(count, count_rate, count_max)
-   t2 = count*1.d0/count_rate
-   if (profile) print *,'3 time=',t2-t1
-
-   ! compute tendency of tracers (including specific humidity) in spectral space.
-   do nt=1,ntrac
-   ! use positive-definite vertical advection.
-   call getvadv_tracers(tracerg(:,:,:,nt),etadot,vadvq)
-!$omp parallel do private(k) schedule(dynamic)
-   do k=1,nlevs
-      ! gradient of specific humidity on grid.
-      call getgrad(tracerspec(:,k,nt),dvirtempdx(:,:,k),dvirtempdy(:,:,k),rerth)
-      ! specific humidity tendency
-      prsgx(:,:,k) = &
-      -ug(:,:,k)*dvirtempdx(:,:,k)-vg(:,:,k)*dvirtempdy(:,:,k)-vadvq(:,:,k)
-      call grdtospec(prsgx(:,:,k), dtracerspecdt(:,k,nt))
-   enddo
-!$omp end parallel do
-   enddo
-   call system_clock(count, count_rate, count_max)
-   t1 = count*1.d0/count_rate
-   if (profile) print *,'4 time=',t1-t2
-   !print *,'getdyntend time=',t1-t0
-
-   ! deallocate storage.
-   deallocate(vadvq,workspec,dvirtempdx,dvirtempdy)
-   deallocate(prsgx,prsgy,vadvu,vadvv,vadvt)
-   deallocate(dlnpsdx,dlnpsdy)
-   if (vcamp > epstiny) deallocate(dvrtdx,dvrtdy)
-
-   return
- end subroutine getdyntend
-
- subroutine getomega(ug,vg,divg,dlnpsdx,dlnpsdy,dlnpsdt,dlnpdtg,etadot,&
-! pass in work storage so it can be re-used, saving memory.
-  workb,workc,cg,cb,db)
-    ! compute omega, etadot, tendency of lnps
-    ! all input and output arrays oriented bottom to top (k=1 is near ground)
-    real(r_kind), intent(in), dimension(nlons,nlats,nlevs) :: ug,vg,divg
-    real(r_kind), intent(in), dimension(nlons,nlats) :: dlnpsdx,dlnpsdy
-    ! omega (pressure vertical velocity divided by pressure) on model layers.
-    real(r_kind), intent(out), dimension(nlons,nlats,nlevs) :: dlnpdtg
-    ! etadot (vertical velocity in hybrid coords) on layer interfaces.
-    real(r_kind), intent(out), dimension(nlons,nlats,nlevs+1) :: etadot
-    real(r_kind), intent(inout), dimension(nlons,nlats) :: dlnpsdt
-! work space:
-    real(r_kind), intent(inout), dimension(nlons,nlats,nlevs) :: &
-    workb,workc,cg,cb,db
-! local scalars
-    integer k
-
-!$omp parallel do private(k)
-    do k=1,nlevs
-     cg(:,:,k)=ug(:,:,nlevs+1-k)*dlnpsdx(:,:)+vg(:,:,nlevs+1-k)*dlnpsdy(:,:)
-    enddo
-!$omp end parallel do
-
-    db(:,:,1)=divg(:,:,nlevs)*dpk(:,:,1)
-    cb(:,:,1)=cg(:,:,1)*dbk(1)
-
-    do k=1,nlevs-1
-     db(:,:,k+1)=db(:,:,k)+divg(:,:,nlevs-k)*dpk(:,:,k+1)
-     cb(:,:,k+1)=cb(:,:,k)+cg(:,:,k+1)*dbk(k+1)
-    enddo
-
-    dlnpsdt(:,:) = -db(:,:,nlevs)/psg(:,:)-cb(:,:,nlevs)
-    etadot(:,:,1) = 0.; etadot(:,:,nlevs+1) = 0.
-!$omp parallel do private(k)
-    do k=1,nlevs-1
-       etadot(:,:,k+1)=-psg(:,:)*(bk(k+1)*dlnpsdt(:,:)+cb(:,:,k)) - db(:,:,k)
-    enddo
-!$omp end parallel do
-
-    workb(:,:,1)=alfa(:,:,1)* &
-                ( divg(:,:,nlevs)*dpk(:,:,1)+psg(:,:)*cb(:,:,1)*dbk(1) )
-
-!$omp parallel do private(k)
-    do k=2,nlevs
-      workb(:,:,k)=rlnp(:,:,k)*( db(:,:,k-1)+psg(:,:)*cb(:,:,k-1) )+&
-      alfa(:,:,k)*( divg(:,:,nlevs+1-k)*dpk(:,:,k)+psg(:,:)*cg(:,:,k)*dbk(k) )
-    enddo
-!$omp end parallel do
-
-    workc(:,:,1)=psg(:,:)*cg(:,:,1)*dbk(1)
-
-!$omp parallel
-!$omp do private(k)
-    do k=2,nlevs
-      workc(:,:,k)=psg(:,:)*cg(:,:,k)*( dbk(k)+ck(k)*rlnp(:,:,k)/dpk(:,:,k) )
-    enddo
-!$omp end do
-!$omp do private(k)
-    do k=1,nlevs
-      dlnpdtg(:,:,nlevs+1-k)=(workc(:,:,k)-workb(:,:,k))/dpk(:,:,k)
-    enddo
-!$omp end do
-!$omp end parallel
-
-    return
- end subroutine getomega
-
- subroutine getpresgrad(virtempg,dvirtempdx,dvirtempdy,dphisdx,dphisdy,dlnpsdx,dlnpsdy,&
-                        prsgx,prsgy,&
-! pass in work storage so it can be re-used, saving memory.
-                        cofa,cofb,px3u,px3v)
-    ! compute pressure gradient terms
-    ! all input and output arrays oriented bottom to top (k=1 is near ground)
-    real(r_kind), intent(in), dimension(nlons,nlats,nlevs) :: &
-    virtempg,dvirtempdx,dvirtempdy
-    real(r_kind), intent(in), dimension(nlons,nlats) :: &
-    dlnpsdx,dlnpsdy,dphisdx,dphisdy
-    ! pressure gradient terms
-    real(r_kind), intent(out), dimension(nlons,nlats,nlevs) :: &
-    prsgx,prsgy
-! work storage
-    real(r_kind), dimension(nlons,nlats,nlevs), intent(inout) :: &
-    cofa,cofb,px3u,px3v
-    integer k
-
-    cofb(:,:,1)=-(1./dpk(:,:,1))*(                 alfa(:,:,1)*dbk(1))
-
-!$omp parallel
-!$omp do private(k)
-    do k=2,nlevs
-       cofb(:,:,k)=-(1./dpk(:,:,k))*(bk(k)*rlnp(:,:,k)+alfa(:,:,k)*dbk(k))
-    enddo
-!$omp end do
-!$omp do private(k)
-    do k=1,nlevs
-        prsgx(:,:,nlevs-k+1)=cofb(:,:,k)*rd*virtempg(:,:,nlevs+1-k)*psg(:,:)*dlnpsdx(:,:)
-        prsgy(:,:,nlevs-k+1)=cofb(:,:,k)*rd*virtempg(:,:,nlevs+1-k)*psg(:,:)*dlnpsdy(:,:)
-    enddo
-!$omp end do
-!$omp do private(k)
-    do k=1,nlevs
-       cofa(:,:,k)=-(1./dpk(:,:,k))*( &
-        bk(k+1)*pk(:,:,k)/pk(:,:,k+1) - bk(k) &
-       +rlnp(:,:,k)*( bk(k)-pk(:,:,k)*dbk(k)/dpk(:,:,k) )  )
-    enddo
-!$omp end do
-!$omp end parallel
-
-    cofb(:,:,nlevs)=0.
-    cofb(:,:,nlevs-1)= &
-      -rd*( bk(nlevs+1)/pk(:,:,nlevs+1)-bk(nlevs)/pk(:,:,nlevs) )*virtempg(:,:,1)
-
-
-    do k=2,nlevs-1
-       cofb(:,:,nlevs-k)=cofb(:,:,nlevs+1-k) &
-       -rd*(bk(nlevs+2-k)/pk(:,:,nlevs+2-k)-bk(nlevs+1-k)/pk(:,:,nlevs+1-k))*&
-                                                        virtempg(:,:,k)
-    enddo
-
-    px3u(:,:,nlevs)=0.
-    px3v(:,:,nlevs)=0.
-    px3u(:,:,nlevs-1)=-rd*rlnp(:,:,nlevs)*dvirtempdx(:,:,1)
-    px3v(:,:,nlevs-1)=-rd*rlnp(:,:,nlevs)*dvirtempdy(:,:,1)
-    do k=2,nlevs-1
-       px3u(:,:,nlevs-k)=px3u(:,:,nlevs+1-k)-rd*rlnp(:,:,nlevs+1-k)*dvirtempdx(:,:,k)
-       px3v(:,:,nlevs-k)=px3v(:,:,nlevs+1-k)-rd*rlnp(:,:,nlevs+1-k)*dvirtempdy(:,:,k)
-    enddo
-
-!$omp parallel do private(k)
-    do k=1,nlevs
-       prsgx(:,:,nlevs-k+1)=prsgx(:,:,nlevs-k+1)-dphisdx(:,:)& !px1u
-       +cofb(:,:,k)*psg(:,:)*dlnpsdx(:,:)& !px2u
-       +px3u(:,:,k)&
-       -rd*alfa(:,:,k)*dvirtempdx(:,:,nlevs+1-k)& !px4u
-       -cofa(:,:,k)*rd*virtempg(:,:,nlevs+1-k)*psg(:,:)*dlnpsdx(:,:) !px5u
-       prsgy(:,:,nlevs-k+1)=prsgy(:,:,nlevs-k+1)-dphisdy(:,:)& !px1v
-       +cofb(:,:,k)*psg(:,:)*dlnpsdy(:,:)& !px2v
-       +px3v(:,:,k)&
-       -rd*alfa(:,:,k)*dvirtempdy(:,:,nlevs+1-k)& !px4v
-       -cofa(:,:,k)*rd*virtempg(:,:,nlevs+1-k)*psg(:,:)*dlnpsdy(:,:) !px5v
-    enddo
-!$omp end parallel do
-
- end subroutine getpresgrad
-
- subroutine getvadv(datag,etadot,vadv)
-   ! compute vertical advection of datag, using etadot - result in vadv
-   ! datag, vadv bottom to top, etadot top to bottom.
-   real(r_kind), intent(in), dimension(nlons,nlats,nlevs) :: datag
-   real(r_kind), intent(in), dimension(nlons,nlats,nlevs+1) :: etadot
-   real(r_kind), intent(out), dimension(nlons,nlats,nlevs) :: vadv
-   integer k
-
-   vadv(:,:,nlevs)= &
-   (0.5/dpk(:,:,1))*etadot(:,:,2)*( datag(:,:,nlevs-1)-datag(:,:,nlevs))
-
-   vadv(:,:,1)= &
-   (0.5/dpk(:,:,nlevs))*etadot(:,:,nlevs)*( datag(:,:,1)-datag(:,:,2) )
-
-!$omp parallel do private(k)
-   do k=2,nlevs-1
-      vadv(:,:,nlevs+1-k)= &
-   (0.5/dpk(:,:,k))*( etadot(:,:,k+1)*( datag(:,:,nlevs  -k)-datag(:,:,nlevs+1-k) )+&
-                      etadot(:,:,k  )*( datag(:,:,nlevs+1-k)-datag(:,:,nlevs+2-k) ) )
-   enddo
-!$omp end parallel do
-
-   return
- end subroutine getvadv
-
- subroutine getvadv_tracers(datag,etadot,vadv)
-   ! compute vertical advection of datag, using etadot - result in vadv
-   ! datag, vadv bottom to top, etadot top to bottom.
-   ! uses positive definite scheme of Thuburn (1993, DOI:
-   ! 10.1002/qj.49711951107)
-   real(r_kind), intent(in), dimension(nlons,nlats,nlevs) :: datag
-   real(r_kind), intent(in), dimension(nlons,nlats,nlevs+1) :: etadot
-   real(r_kind), intent(out), dimension(nlons,nlats,nlevs) :: vadv
-   ! local variables
-   real(r_kind), dimension(:,:,:), allocatable :: datag_half, datag_d
-   integer i,j,k
-   real(r_kind) epstiny,phi
-   !real(r_kind) phi(nlons,nlats)
-
-   allocate(datag_half(nlons,nlats,0:nlevs))
-   allocate(datag_d(nlons,nlats,0:nlevs))
-
-!$omp parallel do private(k)
-   do k=1,nlevs-1
-      datag_half(:,:,k) = 0.5*(datag(:,:,nlevs-k)+datag(:,:,nlevs+1-k))
-      datag_d(:,:,k) = datag(:,:,nlevs-k) - datag(:,:,nlevs+1-k)
-   enddo
-!$omp end parallel do
-
-   datag_half(:,:,0) = datag(:,:,nlevs)
-   datag_half(:,:,nlevs) = datag(:,:,1)
-
-!$omp parallel do private(i)
-   do i=1,nlons
-      where (datag(i,:,nlevs) >= 0)
-         datag_d(i,:,0) = datag(i,:,nlevs) - &
-         max(0.,2.*datag(i,:,nlevs)-datag(i,:,nlevs-1))
-      elsewhere
-         datag_d(i,:,0) = datag(i,:,nlevs) - &
-         min(0.,2.*datag(i,:,nlevs)-datag(i,:,nlevs-1))
-      end where
-      where (datag(i,:,1) >= 0)
-         datag_d(i,:,nlevs) = max(0.,2.*datag(i,:,1)-datag(i,:,2)) - datag(i,:,1)
-      elsewhere
-         datag_d(i,:,nlevs) = min(0.,2.*datag(i,:,1)-datag(i,:,2)) - datag(i,:,1)
-      end where
-   enddo
-!$omp end parallel do
-
-! same as above, but without where blocks inside parallel region.
-!!$omp parallel do private(i,j)
-!   do i=1,nlons
-!      do j=1,nlats
-!         if (datag(i,j,nlevs) >= 0) then
-!            datag_d(i,j,0) = datag(i,j,nlevs) - &
-!            max(0.,2.*datag(i,j,nlevs)-datag(i,j,nlevs-1))
-!         else
-!            datag_d(i,j,0) = datag(i,j,nlevs) - &
-!            min(0.,2.*datag(i,j,nlevs)-datag(i,j,nlevs-1))
-!         endif
-!         if (datag(i,j,1) >= 0) then
-!            datag_d(i,j,nlevs) = max(0.,2.*datag(i,j,1)-datag(i,j,2)) - datag(i,j,1)
-!         else
-!            datag_d(i,j,nlevs) = min(0.,2.*datag(i,j,1)-datag(i,j,2)) - datag(i,j,1)
-!         endif
-!      enddo
-!   enddo
-!!$omp end parallel do
-
-   ! to prevent NaNs in computation of van leer limiter
-   epstiny = tiny(epstiny)
-!$omp parallel do private(i,j,k)
-   do k=1,nlevs-1
-      do j=1,nlats
-      do i=1,nlons
-         if (abs(datag_d(i,j,k)) .lt. epstiny) then
-            datag_d(i,j,k) = sign(epstiny, datag_d(i,j,k))
-         endif
-      enddo
-      enddo
-   enddo
-!$omp end parallel do
-
-!   ! apply flux limiter.
-!!note: this segfaults with intel 12.1 at T574
-!!$omp parallel do private(k,phi)
-!   do k=1,nlevs-1
-!      where(etadot(:,:,k+1) > 0.)
-!         phi = datag_d(:,:,k-1)/datag_d(:,:,k)
-!         datag_half(:,:,k) = datag(:,:,nlevs+1-k) + &
-!         (phi+abs(phi))/(1.+abs(phi))*(datag_half(:,:,k)-datag(:,:,nlevs+1-k))
-!      elsewhere
-!         phi = datag_d(:,:,k+1)/datag_d(:,:,k)
-!         datag_half(:,:,k) = datag(:,:,nlevs-k) + &
-!         (phi+abs(phi))/(1.+abs(phi))*(datag_half(:,:,k)-datag(:,:,nlevs-k))
-!      end where
-!   enddo
-!!$omp end parallel do
-
-!    ! apply flux limiter.
-! to avoid t574 segfault, get rid of where statements inside parallel region.
-!$omp parallel do private(i,j,k,phi)
-    do k=1,nlevs-1
-       do j=1,nlats
-       do i=1,nlons
-          if (etadot(i,j,k+1) > 0.) then
-             phi = datag_d(i,j,k-1)/datag_d(i,j,k)
-             datag_half(i,j,k) = datag(i,j,nlevs+1-k) + &
-             (phi+abs(phi))/(1.+abs(phi))*(datag_half(i,j,k)-datag(i,j,nlevs+1-k))
-          else
-             phi = datag_d(i,j,k+1)/datag_d(i,j,k)
-             datag_half(i,j,k) = datag(i,j,nlevs-k) + &
-             (phi+abs(phi))/(1.+abs(phi))*(datag_half(i,j,k)-datag(i,j,nlevs-k))
-          endif
-       enddo
-       enddo
-     enddo
-!$omp end parallel do
-
-!$omp parallel do private(k)
-   do k=1,nlevs
-      vadv(:,:,nlevs+1-k) = (1./dpk(:,:,k))*(&
-      (datag_half(:,:,k)*etadot(:,:,k+1) - datag_half(:,:,k-1)*etadot(:,:,k))+&
-      (datag(:,:,nlevs+1-k)*(etadot(:,:,k)-etadot(:,:,k+1))))
-   enddo
-!$omp end parallel do
-
-   deallocate(datag_half, datag_d)
-
-   return
- end subroutine getvadv_tracers
-
- subroutine dry_mass_fixer(psg,pwat,dlnpsspecdt,dtx)
-! adjust psg to global integral of dry mass is equal
-! to pdryini.
-! input: psg,pwatg (gridded surface pressure, precipitable water)
-! input/output: dlnpsspecdt (spectral lnps tendency) - modified
-!  to include correction term for dry mass adjustment.
-! input: dtx (time step)
-      real(r_kind), intent(in), dimension(nlons,nlats) ::&
-          psg(nlons,nlats),pwat(nlons,nlats)
-      complex(r_kind), intent(inout) :: dlnpsspecdt(ndimspec)
-      real(r_double), intent(in) :: dtx
-      real(r_kind) pmean,pdry,pcorr,pwatg
-      real(r_kind) dpsdt(nlons,nlats)
-      complex(r_kind) workspec(ndimspec)
-! compute global mean dry ps.
-      pwatg = sum(areawts*pwat)
-      !print *,'global mean pwat = ',pwatg
-      pmean = sum(areawts*psg)
-      pdry = pmean - grav*pwatg
-      !print *,'pdry after physics update',pdry
-! implied ps correction needed to return dry mass to initial value
-      pcorr = (pdryini + grav*pwatg)/pmean
-! apply correction as a multiplication to provisional value of ps so as
-! not to change gradients of ln(ps).
-      dpsdt = psg*pcorr
-! compute implied lnps tendency in spectral space.
-      dpsdt = log(dpsdt)
-      call grdtospec(dpsdt,workspec)
-      dlnpsspecdt = (workspec - lnpsspec)/dtx
-  end subroutine dry_mass_fixer
-
-end module dyn_run
diff --git a/climt/_lib/GFS/finalize.f90 b/climt/_lib/GFS/finalize.f90
deleted file mode 100644
index 2202c97a..00000000
--- a/climt/_lib/GFS/finalize.f90
+++ /dev/null
@@ -1,22 +0,0 @@
-module finalize_mod
-! finalization module (provides finalize driver subroutine)
-private
-public :: finalize
-contains
-
-subroutine finalize() bind(c,name='gfs_finalise');
-  use dyn_finalize ,only: finalize_dyn
-  use phy_finalize ,only: finalize_phy
-  !JOY this was not being called at all
-  use shtns, only: shtns_destroy
-
-  implicit none
-
-  call finalize_dyn ()
-  call finalize_phy ()
-  !JOY this was not being called at all
-  call shtns_destroy()
-  
-  return
-end subroutine finalize
-end module finalize_mod
diff --git a/climt/_lib/GFS/gfs_namelist b/climt/_lib/GFS/gfs_namelist
deleted file mode 100644
index e69de29b..00000000
diff --git a/climt/_lib/GFS/grid_data.f90 b/climt/_lib/GFS/grid_data.f90
deleted file mode 100644
index a801cb64..00000000
--- a/climt/_lib/GFS/grid_data.f90
+++ /dev/null
@@ -1,86 +0,0 @@
-module grid_data
-! gaussian grid model state variable data.
-! public subroutines:
-! init_griddata: allocate arrays.
-! destroy_griddata: deallocate arrays.
- use kinds, only: r_kind
- use params, only: nlons,nlats,nlevs,ntrac
- use iso_c_binding, only: c_double
- implicit none
- private
- public :: init_griddata, destroy_griddata
- real(r_kind), pointer, public, dimension(:,:,:) :: ug,vg,vrtg,divg,&
- virtempg,etadot,dlnpdtg
- real(r_kind), pointer, public, dimension(:,:,:,:) :: tracerg
-! (nlons,nlons,nlevs) arrays (bottom to top unless otherwise noted)
-! they are transformed to the grid from spectral space in subroutine
-! getdyntend in module dyn_run.
-! ug: zonal wind
-! vg: meridional wind
-! vrtg: vorticity
-! divg: divergence
-! tracerg: tracers (first one is specific humidity)
-! virtempg: virtual temperature
-! the following last two are computed in subroutine omega from module dyn_run
-! (called by dyntend):
-! etadot: vertical motion in hybrid sigma-pressure (top to bottom)
-! dlnpdtg: d(lnp)/dt = omega/p (compute in subroutine omega from module dyn_run)
-! dlnpsdt: local tendency of ln(ps)
- real(r_kind), pointer, public, dimension(:,:) :: &
- dlnpsdt,lnpsg,dlnpdst,dphisdx,dphisdy,phis
-
- contains
-
-     subroutine assign_grid_arrays(py_ug, py_vg, py_virt_tempg, py_lnpsg, &
-             py_tracerg, py_vortg, py_divg) bind(c, name='gfs_assign_grid_arrays')
-
-         real(c_double), target, intent(inout):: py_ug(nlons, nlats, nlevs)
-         real(c_double), target, intent(inout):: py_vg(nlons, nlats, nlevs)
-         real(c_double), target, intent(inout):: py_virt_tempg(nlons, nlats, nlevs)
-         real(c_double), target, intent(inout):: py_vortg(nlons, nlats, nlevs)
-         real(c_double), target, intent(inout):: py_divg(nlons, nlats, nlevs)
-
-         real(c_double), target, intent(inout):: py_lnpsg(nlons, nlats)
-         real(c_double), target, intent(inout):: py_tracerg(nlons, nlats, nlevs, ntrac)
-
-         ug => py_ug
-         vg => py_vg
-         vrtg => py_vortg
-         divg => py_divg
-         virtempg => py_virt_tempg
-         tracerg => py_tracerg
-         lnpsg => py_lnpsg
-
-     end subroutine
-
-     subroutine init_griddata(pyDlnpdtg, pyEtaDotg, pyPhis, &
-             pyDPhisdx, pyDPhisdy, pyDlnpsdt) bind(c,name='gfs_initialise_grid_arrays')
-
-         real(c_double), target, intent(inout) :: pyDlnpdtg(nlons, nlats, nlevs)
-
-         real(c_double), target, intent(inout) :: pyEtaDotg(nlons,nlats,nlevs+1)
-
-         real(c_double), target, intent(inout) :: pyPhis(nlons,nlats), pyDPhisdx(nlons,nlats), &
-             pyDPhisdy(nlons,nlats), pyDlnpsdt(nlons,nlats)
-
-         !print *, 'Initialising grid data'
-
-         !print *,  shape(pyDlnpdtg);
-         dlnpdtg => pyDlnpdtg
-         etadot => pyEtaDotg
-         phis => pyPhis
-         dphisdy => pyDPhisdy
-         dphisdx => pyDPhisdx
-         dlnpsdt => pyDlnpsdt
-         !print *, shape(dlnpdtg);
-         !print *, ' Done Initialising grid data'
-     end subroutine init_griddata
-     
-     subroutine destroy_griddata()
-         !print *, 'Deallocating grid data'
-         nullify(ug,vg,vrtg,divg,virtempg,dlnpdtg,etadot)
-         nullify(tracerg)
-         nullify(lnpsg,phis,dphisdx,dphisdy,dlnpsdt)
-     end subroutine destroy_griddata
-
- end module grid_data
diff --git a/climt/_lib/GFS/init.f90 b/climt/_lib/GFS/init.f90
deleted file mode 100644
index 7142370b..00000000
--- a/climt/_lib/GFS/init.f90
+++ /dev/null
@@ -1,20 +0,0 @@
-module init_mod
-! initialization module (provides init driver subroutine)
-private
-!public :: init
-contains
-
-!subroutine init() bind(c,name='dycoreInit')
-!!  use params, only: read_namelist
-!  use dyn_init ,only: init_dyn
-!  use phy_init ,only: init_phy
-!
-!  implicit none
-!  print *, 'Im here!!!'
-!!  call read_namelist ()
-!  call init_dyn ()
-!  call init_phy ()
-!  
-!  return
-!end subroutine init
-end module init_mod
diff --git a/climt/_lib/GFS/kinds.f90 b/climt/_lib/GFS/kinds.f90
deleted file mode 100644
index 8403991e..00000000
--- a/climt/_lib/GFS/kinds.f90
+++ /dev/null
@@ -1,107 +0,0 @@
-! this module was extracted from the GSI version operational
-! at NCEP in Dec. 2007.
-module kinds
-!$$$  module documentation block
-!                .      .    .                                       .
-! module:   kinds
-!   prgmmr: treadon          org: np23                date: 2004-08-15
-!
-! abstract:  Module to hold specification kinds for variable declaration.
-!            This module is based on (copied from) Paul vanDelst's 
-!            type_kinds module found in the community radiative transfer
-!            model
-!
-! module history log:
-!   2004-08-15  treadon
-!
-! Subroutines Included:
-!
-! Functions Included:
-!
-! remarks:
-!   The numerical data types defined in this module are:
-!      i_byte    - specification kind for byte (1-byte) integer variable
-!      i_short   - specification kind for short (2-byte) integer variable
-!      i_long    - specification kind for long (4-byte) integer variable
-!      i_llong   - specification kind for double long (8-byte) integer variable
-!      r_single  - specification kind for single precision (4-byte) real variable
-!      r_double  - specification kind for double precision (8-byte) real variable
-!      r_quad    - specification kind for quad precision (16-byte) real variable
-!
-!      i_kind    - generic specification kind for default integer
-!      r_kind    - generic specification kind for default floating point
-!
-!
-! attributes:
-!   language: f90
-!   machine:  ibm RS/6000 SP
-!
-!$$$ end documentation block
-  implicit none
-  private
-
-! Integer type definitions below
-
-! Integer types
-  integer, parameter, public  :: i_byte  = selected_int_kind(1)      ! byte  integer
-  integer, parameter, public  :: i_short = selected_int_kind(4)      ! short integer
-  integer, parameter, public  :: i_long  = selected_int_kind(8)      ! long  integer
-  integer, parameter, private :: llong_t = selected_int_kind(16)     ! llong integer
-  integer, parameter, public  :: i_llong = max( llong_t, i_long )
-
-! Expected 8-bit byte sizes of the integer kinds
-  integer, parameter, public :: num_bytes_for_i_byte  = 1
-  integer, parameter, public :: num_bytes_for_i_short = 2
-  integer, parameter, public :: num_bytes_for_i_long  = 4
-  integer, parameter, public :: num_bytes_for_i_llong = 8
-
-! Define arrays for default definition
-  integer, parameter, private :: num_i_kinds = 4
-  integer, parameter, dimension( num_i_kinds ), private :: integer_types = (/ &
-       i_byte, i_short, i_long,  i_llong  /) 
-  integer, parameter, dimension( num_i_kinds ), private :: integer_byte_sizes = (/ &
-       num_bytes_for_i_byte, num_bytes_for_i_short, &
-       num_bytes_for_i_long, num_bytes_for_i_llong  /)
-
-! Default values
-! **** CHANGE THE FOLLOWING TO CHANGE THE DEFAULT INTEGER TYPE KIND ***
-  integer, parameter, public :: default_integer = 3  ! 1=byte, 
-                                                      ! 2=short, 
-                                                      ! 3=long, 
-                                                      ! 4=llong
-  integer, parameter, public  :: i_kind = integer_types( default_integer )
-  integer, parameter, public  :: num_bytes_for_i_kind = &
-       integer_byte_sizes( default_integer )
-
-
-! Real definitions below
-
-! Real types
-  integer, parameter, public  :: r_single = selected_real_kind(6)  ! single precision
-  integer, parameter, public  :: r_double = selected_real_kind(15) ! double precision
-  integer, parameter, private :: quad_t   = selected_real_kind(20) ! quad precision
-  integer, parameter, public  :: r_quad   = max( quad_t, r_double )
-
-! Expected 8-bit byte sizes of the real kinds
-  integer, parameter, public :: num_bytes_for_r_single = 4
-  integer, parameter, public :: num_bytes_for_r_double = 8
-  integer, parameter, public :: num_bytes_for_r_quad   = 16
-
-! Define arrays for default definition
-  integer, parameter, private :: num_r_kinds = 3
-  integer, parameter, dimension( num_r_kinds ), private :: real_kinds = (/ &
-       r_single, r_double, r_quad    /) 
-  integer, parameter, dimension( num_r_kinds ), private :: real_byte_sizes = (/ &
-       num_bytes_for_r_single, num_bytes_for_r_double, &
-       num_bytes_for_r_quad    /)
-
-! Default values
-! **** CHANGE THE FOLLOWING TO CHANGE THE DEFAULT REAL TYPE KIND ***
-  integer, parameter, public :: default_real = 2  ! 1=single, 
-                                                  ! 2=double, 
-                                                  ! 3=quad
-  integer, parameter, public  :: r_kind = real_kinds( default_real )
-  integer, parameter, public  :: num_bytes_for_r_kind = &
-       real_byte_sizes( default_real )
-
-end module kinds
diff --git a/climt/_lib/GFS/params.f90 b/climt/_lib/GFS/params.f90
deleted file mode 100644
index a6a68140..00000000
--- a/climt/_lib/GFS/params.f90
+++ /dev/null
@@ -1,251 +0,0 @@
-module params
-! holds model parameters
-! Public subroutines:
-! read_namelist: read namelist
- use kinds, only: r_kind,r_double
-! JOY removing sigio dependency
-! use sigio_module, only: sigio_head, sigio_srhead, sigio_sropen, sigio_sclose
-
-!JOY adding iso_c_binding required for relevant variables
-use iso_c_binding, only: c_double,c_int,c_bool
- implicit none
- private
-
- public :: initfile,sfcinitfile,dt,ntmax,ndimspec,nlons,nlats,&
- tstart,ndiss,efold,nlevs,ntrunc,dry,explicit,heldsuarez,dcmip,&
- ntout,fhdfi,fhzer,idate_start,adiabatic,hdif_fac,hdif_fac2,fshk,ntrac,ntoz,ntclw,&
- pdryini, toa_pressure, massfix,timestepsperhr,ncw,taustratdamp,polar_opt,ntdfi,gfsio_out,sigio_out,&
-! gfs phys parameters.
- nmtvr,fhlwr,fhswr,ictm,isol,ico2,iaer,ialb,iems,isubc_sw,isubc_lw,&
- iovr_sw,iovr_lw,newsas,ras,sashal,num_p3d,num_p2d,crick_proof,ccnorm,&
- norad_precip,crtrh,cdmbgwd,ccwf,dlqf,ctei_rm,psautco,prautco,evpco,wminco,flgmin,&
- old_monin,cnvgwd,mom4ice,shal_cnv,cal_pre,trans_trac,nst_fcst,moist_adj,mstrat,&
- ngptc,pre_rad,bkgd_vdif_m,bkgd_vdif_s,bkgd_vdif_h,gloopb_filter,&
-! iau parameters
- iau,iaufiles_fg,iaufiles_anl,iaufhrs,iau_delthrs,&
-! vorticity confinement parameters
- vcamp,svc,svc_tau,svc_lscale,iseed_svc,&
-! stochastic physics tendency parameters
- sppt,sppt_tau,sppt_lscale,iseed_sppt, &
-! additive stochastic humidity perturbations
- shum,shum_tau,shum_lscale,iseed_shum,clipsupersat,&
-! additive noise perts
- addnoise,addnoise_tau,addnoise_lscale,addnoise_vfilt,iseed_addnoise,&
- addnoise_kenorm, addnoise_dissfact,addnoise_vrtonly
-
-
- !! JOY adding bind(c) arguments wherever needed
-
- character(len=500) :: initfile ! init cond filename
- character(len=500) :: sfcinitfile ! surface init cond filename
-! integer(c_int)     :: fhmax ! hours to run
-! integer(c_int)     :: fhout ! interval for IO
- integer            :: fhzer ! interval to zero accumulated arrays
-! half window length (hrs) for digital filter launch (=0 mean no dfi)
- integer            :: fhdfi=0
- real(c_double)     :: deltim=0    ! namelist input time step (secs)
-
- real(c_double)     :: dt    ! time step (secs) (=deltim or 3600/timestepsperhr)
-! bind(c) :: fhmax,fhout
- bind(c) :: deltim,dt
-
- real(c_double) :: pdryini ! initial dry ps
- bind(c) :: pdryini
-
- real(c_double) :: toa_pressure !top of model pressure
- bind(c) :: toa_pressure
-
- logical(c_bool)  :: massfix=.true. ! apply dry mass 'fixer'
- bind(c) :: massfix
-
- logical    :: gfsio_out=.false. ! write out 'gfsio' grib 1 files.
- logical    :: sigio_out=.true. ! write out 'sigma' spectral binary files.
- integer    :: ntmax ! time steps to run
- integer    :: ntdfi ! number of time steps in dfi window is 2*ntdfi+1
-
- integer(c_int)    :: nlons=1 ! number of longitudes on grid
- integer(c_int)    :: nlats=1 ! number of latitudes on grid
- integer(c_int)    :: nlevs=1 ! number of levels on grid
- integer(c_int)    :: ntrunc=1 ! spectral truncation
- integer(c_int)    :: ndimspec=1 ! spectral array dimension
- 
-
-! JOY removing sigio dependency
-! type(sigio_head),save  :: sighead ! header struct from initfile
- logical(c_bool)    :: dry = .false. ! no moisture, cloud condensate or ozone.
- logical(c_bool)    :: adiabatic = .false. ! don't call physics
- ! held-suarez forcing
- logical(c_bool)    :: heldsuarez = .false.
- bind(c) :: dry,adiabatic,heldsuarez
-
- logical    :: iau = .false. ! iau forcing included
- integer    :: iau_delthrs = 6 ! iau time interval (to scale increments)
- character(len=120), dimension(7) ::  iaufiles_fg,iaufiles_anl
- real(r_kind), dimension(7) :: iaufhrs
- ! dcmip test cases
- integer    :: dcmip = -1 ! 4x for baroclinic wave, 5x for tropical cyclone.
- ! 3-stage, 2nd order RK time differencing is used.
- ! explicit=.true. means use explicit time differencing
- ! if .false., semi-implicit or additive RK (Ascher, Spiteri and Ruth 1997
- ! DOI:10.1016/S0168-9274(97)00056-1) is used.
- logical    :: explicit = .false. ! use semi-implicit RK
- ! starting forecast time in seconds (read in from initfile)
- real(r_kind) :: tstart
- integer    :: idate_start(4) ! starting date (hr,month,day,year)
- integer    :: ntout ! time step interval for IO
- integer    :: ndiss=0 ! hyperdiffusion order (0 for GFS defaults)
- real(r_double) :: polar_opt=1.e-10 ! polar optimization threshold for transforms
- ! efolding scale for smallest resolvable wave (0 for GFS defaults)
- real(r_kind) :: efold=0. 
- real(r_kind) :: hdif_fac=1.0 ! multiplier for height-dep part of hyper-diff
- real(r_kind) :: hdif_fac2=1.0 ! multiplier to increase hyper-diffusion
- ! amplitude of vertically varying part of hyper-diff (1 means no variation,
- ! zero gives GFS resolution dependent defaults)
- real(r_kind) :: fshk=0 
- integer(c_int) :: ntrac=3 ! number of tracers (including specific humidity)
-
-   
-!JOY bind the grid size variables   
- bind(c) :: nlons, nlats, nlevs, ntrunc, ndimspec, ntrac
-
- integer :: ntoz=2 ! ozone tracer number
- integer :: ntclw=3 ! cloud condensate tracer number
- integer :: nmtvr=14 ! number of fields in mtnvar file.
- integer :: ngptc=24 ! pass this many profiles to gbphys in one call.
- real(r_kind) :: taustratdamp=2.*86400. ! extra linear drag near top of model
-! parameters relevant for GFS physics
- ! interval in hours to call long-wave radiation (0 means every time step)
- real(r_kind) :: fhlwr=0 
- ! interval in hours to call short-wave radiation (0 means every time step)
- real(r_kind) :: fhswr=0 
- ! ictm controls source for controls source for time sensitive external data (e.g. CO2,
- ! solcon, aerosols, etc)
- integer :: ictm=0 ! use data at initial cond time, or latest available.
- ! isol controls solar constant data source.
- integer :: isol=0 ! use prescribed solar constant
- ! ico2 controls co2 data source for radiation
- integer :: ico2=0 ! prescribed global mean co2
- ! iaer controls aerosols scheme selection
- integer :: iaer=1 ! default aerosol
- ! ialb controls surface albedo for sw radiation
- integer :: ialb=0 ! use climo albedo based on sfc type
- ! iems controls surface emissivity and sfc air/ground temp for lw radiation
- integer :: iems=0 ! used fixed value of 1.0
- ! isubc_sw,lw controls sub-column cloud approximation in radiation
- integer :: isubc_sw=0 ! sw clouds without sub-column approximation
- integer :: isubc_lw=0 ! lw clouds without sub-column approximation
- ! iovr_sw,lw controls cloud overlapping method in radiation
- integer :: iovr_sw=1 ! max-random overlap clouds
- integer :: iovr_lw=1 ! max-random overlap clouds
- logical :: newsas=.true. ! 'new. SAS convection scheme
- logical :: ras=.false. ! RAS convection scheme.
- logical :: sashal=.true. ! 'new' mass-sflux based shallow convection
- logical :: mstrat=.false. ! flag for moorthi approach for stratus
- logical :: pre_rad=.false. ! flag for debugging in gbphys
- integer :: num_p3d=4 ! # of 3d microphysics fields (value for Zhao microphysics)
- integer :: num_p2d=3 ! # of 2d microphysics fields (value for Zhao)
- logical :: crick_proof = .false.
- logical :: ccnorm = .false.
- logical :: norad_precip = .false. ! only used for Ferrier microphysics
- ! range of droplet number concentrations for Ferrier scheme.
- integer :: ncw(2)=(/50,150/)
- ! critical relative humidity at the surface, PBL and atmosphere top   
- real(r_kind) :: crtrh(3)=(/0.85,0.85,0.85/)
- ! multiplication factors for cdmb and gwd  (Mtn Blking and GWD tuning factors)
- real(r_kind) :: cdmbgwd(2)=(/1.0,1.0/)
- ! RAS convection mult factor for critical cloud work function.
- real(r_kind) :: ccwf(2)=(/1.0,1.0/)
- ! factor for cloud condensate detrainment from cloud edges in RAS.
- real(r_kind) :: dlqf(2)=(/0.0,0.0/)
- ! critical cloud top entrainment instability criteria (for mstrat=.true.)
- real(r_kind) :: ctei_rm(2)=(/10.0,10.0/)
- ! auto conversion coeff from ice to snow for Zhao microphysics
- real(r_kind) :: psautco(2)=(/4.0e-4,4.0e-4/)
- ! auto conversion coeff from cloud to rain for Zhao microphysics
- real(r_kind) :: prautco(2)=(/1.0e-4,1.0e-4/)
- real(r_kind) :: evpco=2.e-5 ! Zhao scheme evap coefficient for lg-scale rain
- real(r_kind) :: wminco(2)=(/1.e-5,1.e-5/) !  water and ice minimum threshold for Zhao 
- real(r_kind) :: flgmin(2)=(/0.2,0.2/) ! (Ferrier only) range of minimum large ice fraction
- real(r_kind) :: vcamp=0. ! vorticity confinement amplitude
- real(r_kind) :: svc=0.   ! stochastic vorticity confinement amplitude
- real(r_kind) :: svc_tau=0.      ! stochastic vorticity confinement time scale
- real(r_kind) :: svc_lscale=0.   ! stochastic vorticity confinement length scale
- integer :: iseed_svc=0 ! random seed for stochastic vc (zero means use clock)
- real(r_kind) :: sppt=0.  ! stochastic physics tendency amplitude
- real(r_kind) :: sppt_tau=0.  ! stochastic physics tendency time scale
- real(r_kind) :: sppt_lscale=0.  ! stochastic dynamics tendency length scale
- integer :: iseed_sppt=0 ! random seed for sppt (0 means use system clock)
- real(r_kind) :: shum=0.  ! stochastic humidity pert amplitude
- real(r_kind) :: shum_tau=0.  ! stochastic humidity pert time scale
- real(r_kind) :: shum_lscale=0.  ! stochastic humidity pert length scale
- integer :: iseed_shum=0 ! random seed for stochastic humid pert (0 means use system clock)
- logical :: addnoise_kenorm=.false. ! add noise in ke norm (instead of psi norm)
- logical :: addnoise_dissfact=.false. ! modulate noise amplitude with ke diss
- logical :: addnoise_vrtonly=.true. ! include balanced t noise (if false)
- real(r_kind) :: addnoise=0.  ! add noise pert amplitude
- real(r_kind) :: addnoise_tau=0.  ! add noise pert time scale
- real(r_kind) :: addnoise_lscale=0.  ! add noise pert length scale
- integer :: iseed_addnoise=0 ! random seed for add noise pert (0 means use system clock)
- integer :: addnoise_vfilt=0 ! add noise 1-2-1 vertical filter parameter
- logical :: old_monin = .false. ! flag for old Monin-Obhukov surface layer
- logical :: cnvgwd = .false. ! flag for convective gravity wave drag
- logical :: mom4ice = .false. ! flag for MOM4 sea-ice scheme
- logical :: shal_cnv = .true. ! use shallow convection?
- logical :: cal_pre = .false. ! true for huiya's precip type algorithm
- logical :: trans_trac = .true. ! convective transport of tracers? (RAS only)
- integer :: nst_fcst=0 ! 0 - AM only, 1 - uncoupled, 2 - coupled
- logical :: moist_adj = .false. 
- logical :: gloopb_filter = .true. ! apply spectral filter to physics tendencies
- ! make sure stochastic perts don't create neg or supersat humidities.
- logical :: clipsupersat=.false. 
-
- real(r_kind) :: bkgd_vdif_m = 3.0 ! background vertical diffusion for momentum
- real(r_kind) :: bkgd_vdif_h = 1.0 ! background vertical diffusion for heat, q
- real(r_kind) :: bkgd_vdif_s = 0.2 ! sigma threshold for background mom. diffusn 
- ! if dt not given, but timestepsperhr is, dt=3600/timestepsperhr
- real(r_double) :: timestepsperhr = -1
- !JOY added to remove sigio dependency
- integer(c_int) :: ntracin = 0
- bind(c) :: ntracin
-
-contains
-
-    subroutine set_dry_pressure(py_pdryini)bind(c, name='gfs_set_dry_pressure')
-
-        real(c_double), intent(in):: py_pdryini
-
-        pdryini = py_pdryini
-
-    end subroutine
-
-    subroutine set_model_top_pressure(py_ptoa)bind(c, name='gfs_set_model_top_pressure')
-
-        real(c_double), intent(in):: py_ptoa
-
-        toa_pressure = py_ptoa
-
-    end subroutine
-
-    subroutine gfs_set_time_step(py_dt)bind(c, name='gfs_set_time_step')
-
-        real(c_double), intent(in):: py_dt
-
-        dt = py_dt
-        deltim = py_dt
-
-    end subroutine
-
-    subroutine gfs_set_model_dimensions(py_nlats, py_nlons, py_nlevs, py_ntrunc,&
-                    py_ndimspec, py_ntrac)bind(c,name='gfs_set_model_dimensions')
-
-        integer(c_int), intent(in):: py_nlats, py_nlons, py_nlevs, py_ntrunc, py_ndimspec, py_ntrac
-
-        nlats = py_nlats
-        nlons = py_nlons
-        nlevs = py_nlevs
-        ntrunc = py_ntrunc
-        ndimspec = py_ndimspec
-        ntrac = py_ntrac
-
-    end subroutine
-end module params
diff --git a/climt/_lib/GFS/phy_data.f90 b/climt/_lib/GFS/phy_data.f90
deleted file mode 100644
index 30314030..00000000
--- a/climt/_lib/GFS/phy_data.f90
+++ /dev/null
@@ -1,72 +0,0 @@
-module phy_data
- use kinds, only: r_kind,r_double
- use params, only: nlats,nlons,nlevs,ntrac,ndimspec
- use iso_c_binding
- ! stub
- ! public subroutines:
- ! init_phydata: allocate and populate arrays.
- ! destroy_phydata: deallocate arrays.
-
- implicit none
- private
- 
- !Temporary variable to store tendencies. This will eventually
- ! be replaced by a tendency function outside the code.
-
- complex(r_kind), dimension(:,:), pointer :: spectral_vort_tend, spectral_div_tend, spectral_virtemp_tend
-
- complex(r_kind), dimension(:,:,:), pointer :: spectral_tracer_tend
- complex(r_kind), dimension(:), pointer :: spectral_lnps_tend
-
- public :: init_phydata, destroy_phydata, wrtout_flx, wrtout_sfc, pwat
- public :: spectral_vort_tend, spectral_div_tend, spectral_virtemp_tend, spectral_tracer_tend, spectral_lnps_tend
-
- real(r_kind), pointer, dimension(:,:) :: pwat
-
- contains
-
-     subroutine init_phydata(py_spec_vort_tend, py_spec_div_tend, py_spec_theta_tend,&
-             py_spec_tracer_tend, py_spec_lnps_tend)bind(c, name='gfs_initialise_spectral_physics_arrays')
-
-         complex(c_double), dimension(:,:), target, intent(inout):: py_spec_vort_tend(ndimspec, nlevs),&
-            py_spec_div_tend(ndimspec, nlevs), py_spec_theta_tend(ndimspec, nlevs)
-         complex(c_double), target, intent(inout):: py_spec_tracer_tend(ndimspec, nlevs, ntrac)
-         complex(c_double), target, intent(inout):: py_spec_lnps_tend(ndimspec)
-
-         !print *, 'Initialising physics'
-         spectral_vort_tend => py_spec_vort_tend
-         spectral_div_tend => py_spec_div_tend
-         spectral_virtemp_tend => py_spec_theta_tend
-         spectral_lnps_tend => py_spec_lnps_tend
-         spectral_tracer_tend => py_spec_tracer_tend
-
-     end subroutine init_phydata
-
-     subroutine init_phydata_pwat(py_pwat)bind(c, name='gfs_initialise_grid_physics_arrays')
-         real(c_double), target, intent(inout):: py_pwat(nlons, nlats)
-
-         pwat => py_pwat
-         pwat = 0.
-     end subroutine
-
-     subroutine destroy_phydata()
-         nullify(pwat)
-         nullify(spectral_lnps_tend)
-         nullify(spectral_vort_tend,spectral_div_tend,spectral_virtemp_tend)
-         nullify(spectral_tracer_tend)
-     end subroutine destroy_phydata
-
-     subroutine wrtout_sfc(fhour,filename)
-         implicit none
-         real(r_kind), intent(in) :: fhour
-         character(len=120) filename
-     end subroutine wrtout_sfc
-
-     subroutine wrtout_flx(fhour,ta,filename)
-         implicit none
-         real(r_kind),intent(in) :: fhour
-         real(r_double),intent(in) :: ta
-         character(len=120),intent(in) :: filename
-     end subroutine wrtout_flx
-
- end module phy_data
diff --git a/climt/_lib/GFS/phy_finalize.f90 b/climt/_lib/GFS/phy_finalize.f90
deleted file mode 100644
index bd67013e..00000000
--- a/climt/_lib/GFS/phy_finalize.f90
+++ /dev/null
@@ -1,17 +0,0 @@
-module phy_finalize
-! stub that provides do nothing finalize_phy routine.
-
- use phy_data, only: destroy_phydata
- implicit none
- private
- public :: finalize_phy
-
- contains
-
- subroutine finalize_phy()
-    call destroy_phydata()
-    return
- end subroutine finalize_phy
-
-
-end module phy_finalize
diff --git a/climt/_lib/GFS/phy_init.f90 b/climt/_lib/GFS/phy_init.f90
deleted file mode 100644
index f8aad6b5..00000000
--- a/climt/_lib/GFS/phy_init.f90
+++ /dev/null
@@ -1,14 +0,0 @@
-module phy_init
-! stub that provides do nother init_phy routine.
-
- use phy_data, only: init_phydata
- implicit none
- private
- public :: init_phy
-
- contains
-
- subroutine init_phy() bind(c,name='gfs_init_physics')
- end subroutine init_phy
-
-end module phy_init
diff --git a/climt/_lib/GFS/phy_run.f90 b/climt/_lib/GFS/phy_run.f90
deleted file mode 100644
index 08c8137b..00000000
--- a/climt/_lib/GFS/phy_run.f90
+++ /dev/null
@@ -1,137 +0,0 @@
- module phy_run
-! compute physics tendencies for simplifed held-suarez forcing
-! getphytend: compute physics tendencies.
-! Public subroutines:
-! getphytend: compute newtonian damping of temperature, drag terms
-! in vort and div eqns.
-
- use params, only: nlevs,nlons,nlats,ndimspec,ntrac,pdryini
- use kinds, only: r_kind, r_double
- use shtns, only: grdtospec, lats, areawts, print_bool
- !use grid_data, only: vrtg,divg,virtempg
- !use pressure_data, only:  prs,psg
- !use spectral_data, only: lnpsspec
-
- use phy_data, only: spectral_vort_tend, spectral_div_tend, spectral_virtemp_tend, spectral_tracer_tend, spectral_lnps_tend
- use,intrinsic  :: iso_c_binding
-
-
- implicit none
- private
-
-
- ! Define interface of call-back routine.
-   !    ABSTRACT INTERFACE
-   !      SUBROUTINE callback (x)
- !         USE, INTRINSIC :: ISO_C_BINDING
- !          REAL(KIND=C_DOUBLE), INTENT(IN), VALUE :: x
-!         END SUBROUTINE callback
-!       END INTERFACE
- ! Define interface of call-back routine.
-       ABSTRACT INTERFACE
-         SUBROUTINE callback ()
-           USE, INTRINSIC :: ISO_C_BINDING
-         END SUBROUTINE callback
-       END INTERFACE
-
- public :: getphytend, set_tendencies
- !public :: getphytend, gfsRegisterPhysicsCallback, set_tendencies
-
- !procedure(callback), pointer :: pyPhysicsCallback => null()
-
- integer(c_int) :: phy_trac
- contains
-
- !subroutine gfsRegisterPhysicsCallback(callback) bind(c,name='gfs_register_physics_callback')
-
- !   type(c_funptr), intent(in), value :: callback
-
- !   print *, 'Registering callback'
- !   print *, callback
- !   call c_f_procpointer(callback, pyPhysicsCallback)
-
- !   call pyPhysicsCallback()
-
- !end subroutine gfsRegisterPhysicsCallback
-
- subroutine getphytend(dvrtspecdt,ddivspecdt,dvirtempspecdt,dtracerspecdt,dlnpsspecdt,t,dt)
-     complex(r_kind), intent(inout), dimension(ndimspec,nlevs) :: &
-         dvrtspecdt,ddivspecdt,dvirtempspecdt
-     complex(r_kind), intent(out), dimension(ndimspec,nlevs,ntrac) :: &
-         dtracerspecdt
-     real(r_double), intent(in) :: t,dt
-     complex(r_kind), intent(inout), dimension(ndimspec) :: dlnpsspecdt
-     integer k, trac
-
-
-     phy_trac = ntrac
-     !$omp parallel do private(k)
-     do k=1,nlevs
-     dvrtspecdt(:,k) = spectral_vort_tend(:,k)
-     ddivspecdt(:,k) = spectral_div_tend(:,k)
-     dvirtempspecdt(:,k) = spectral_virtemp_tend(:,k)
-     do trac=1,phy_trac
-         dtracerspecdt(:,k,trac) = spectral_tracer_tend(:,k,trac)
-     enddo
-     enddo
-     !$omp end parallel do
-     dlnpsspecdt(:) = spectral_lnps_tend(:)
-
-
- end subroutine getphytend
-
-
-subroutine set_tendencies(dvrtgdt,ddivgdt,dvirtempgdt,dlnpsgdt,dtracergdt) bind(c, name='gfs_set_tendencies')
-
-    real(c_double), intent(in), dimension(nlons,nlats,nlevs) :: &
-        dvrtgdt,ddivgdt,dvirtempgdt
-    real(c_double), intent(in), dimension(nlons,nlats,nlevs,ntrac) :: &
-        dtracergdt
-    real(c_double), intent(in), dimension(nlons,nlats) :: dlnpsgdt
-    !complex(r_kind), dimension(:,:),allocatable :: forcingspec
-    real(c_double), dimension(nlons, nlats) :: temp_grid_tracer
-    integer k,trac
-
-
-    !print *, dtracergdt(:, :, 1, 1)
-    phy_trac = ntrac
-    !print *, maxval(abs(dvrtgdt(:,:,9))), maxval(abs(ddivgdt)), maxval(abs(dvirtempgdt(:,:,9)))
-    !print *, maxval(abs(dtracergdt(:,:,:,1)))
-    !$omp parallel do private(k)
-    do k=1,nlevs
-    call grdtospec(dvrtgdt(:,:,k), spectral_vort_tend(:,k))
-    call grdtospec(ddivgdt(:,:,k), spectral_div_tend(:,k))
-    call grdtospec(dvirtempgdt(:,:,k), spectral_virtemp_tend(:,k))
-    if (k == 1) then
-        print_bool = 1
-    end if
-    print_bool = 0
-    !spectral_vort_tend(:,k) = dvrtspecdt(:,k)
-    !spectral_div_tend(:,k) = ddivspecdt(:,k)
-    !spectral_virtemp_tend(:,k) = dvirtempspecdt(:,k)
-    enddo
-    !$omp end parallel do
-    do trac=1,phy_trac
-        do k=1,nlevs
-        !print *, ''
-        !print *, 'Level: ', k
-        temp_grid_tracer(:,:) = dtracergdt(:,:,k,trac)
-        !print *, temp_grid_tracer
-        !print *,
-            call grdtospec(temp_grid_tracer, spectral_tracer_tend(:,k,trac))
-        enddo
-    enddo
-
-    call grdtospec(dlnpsgdt(:,:), spectral_lnps_tend(:))
-    !spectral_lnps_tend(:) = dlnpsspecdt(:)
-    !print *, maxval(abs(realpart((spectral_virtemp_tend(:,:)))))
-    !print *, maxval(abs(realpart((spectral_tracer_tend(:,:,:)))))
-    !print *, maxval(abs(imagpart((spectral_tracer_tend(:,:,1)))))
-
-    !deallocate(spectral_lnps_tend)
-    !deallocate(spectral_vort_tend,spectral_div_tend,spectral_virtemp_tend)
-    !deallocate(forcingspec)
-end subroutine set_tendencies
-
-
-end module phy_run
diff --git a/climt/_lib/GFS/phy_run_stub.f90 b/climt/_lib/GFS/phy_run_stub.f90
deleted file mode 100644
index e3f111b7..00000000
--- a/climt/_lib/GFS/phy_run_stub.f90
+++ /dev/null
@@ -1,25 +0,0 @@
- module phy_run
-! do-nothing stub. Provides a dummy getphytend for
-! adiabatic model test.
-
- use params, only: ndimspec, nlevs, ntrac
- use kinds, only: r_kind
- implicit none
- private
-
- public :: getphytend
-
- contains
-
- subroutine getphytend(dvrtspecdt,ddivspecdt,dvirtempspecdt,dtracerspecdt,dlnpsspecdt,t,dt)
-   ! compute physics tendencies for held-suarez test case.
-   complex(r_kind), intent(inout), dimension(ndimspec,nlevs) :: &
-   dvrtspecdt,ddivspecdt,dvirtempspecdt
-   complex(r_kind), intent(out), dimension(ndimspec,nlevs,ntrac) :: &
-   dtracerspecdt
-   real(r_kind), intent(in) :: t,dt
-   complex(r_kind), intent(inout), dimension(ndimspec) :: dlnpsspecdt
-   return
- end subroutine getphytend
-
-end module phy_run
diff --git a/climt/_lib/GFS/physcons.f90 b/climt/_lib/GFS/physcons.f90
deleted file mode 100644
index 68f44841..00000000
--- a/climt/_lib/GFS/physcons.f90
+++ /dev/null
@@ -1,99 +0,0 @@
-!   this module contains some the most frequently used math and        !
-!   physics constatns for gcm models.                                  !
-!                                                                      !
-!   references:                                                        !
-!     as set in NMC handbook from Smithsonian tables.                  !
-!                                                                      !
-          module physcons                
-!
-  use kinds, only : r_kind
-  use iso_c_binding
-!
-  implicit none
-!
-  public
-
-!  --- ...  Math constants
-
-  real(r_kind),parameter:: con_pi     =3.1415926535897931 ! pi
-  real(r_kind),parameter:: con_sqrt2  =1.414214e+0        ! square root of 2
-  real(r_kind),parameter:: con_sqrt3  =1.732051e+0        ! square root of 3
-
-!  --- ...  Geophysics/Astronomy constants
-
-  real(r_kind):: con_rerth  =6.3712e+6      ! radius of earth   (m)
-  real(r_kind):: con_g      =9.80665e+0     ! gravity           (m/s2)
-  real(r_kind):: con_omega  =7.2921e-5      ! ang vel of earth  (1/s)
-  !real(r_kind),parameter:: con_p0     =1.01325e5      ! std atms pressure (pa)
-  !real(r_kind),parameter:: con_solr   =1.3660e+3      ! solar constant    (W/m2)-liu(2002)
-
-!  --- ...  Thermodynamics constants
-
-  real(r_kind):: con_rgas   =8.314472       ! molar gas constant  (J/mol/K)
-  real(r_kind):: con_rd     =2.8705e+2      ! gas constant air    (J/kg/K)
-  real(r_kind):: con_rv     =4.6150e+2      ! gas constant H2O    (J/kg/K)
-  real(r_kind):: con_cp     =1.0046e+3      ! spec heat air @p    (J/kg/K)
-!  real(r_kind),parameter:: con_cv     =7.1760e+2      ! spec heat air @v    (J/kg/K)
-  real(r_kind):: con_cvap   =1.8460e+3      ! spec heat H2O gas   (J/kg/K)
-!  real(r_kind),parameter:: con_cliq   =4.1855e+3      ! spec heat H2O liq   (J/kg/K)
-!  real(r_kind),parameter:: con_csol   =2.1060e+3      ! spec heat H2O ice   (J/kg/K)
-!  real(r_kind),parameter:: con_hvap   =2.5000e+6      ! lat heat H2O cond   (J/kg)
-!  real(r_kind),parameter:: con_hfus   =3.3358e+5      ! lat heat H2O fusion (J/kg)
-!  real(r_kind),parameter:: con_psat   =6.1078e+2      ! pres at H2O 3pt     (Pa)  
-!  real(r_kind),parameter:: con_t0c    =2.7315e+2      ! temp at 0C          (K)
-!  real(r_kind),parameter:: con_ttp    =2.7316e+2      ! temp at H2O 3pt     (K)
-!  real(r_kind),parameter:: con_tice   =2.7120e+2      ! temp freezing sea     (K)
-!  real(r_kind),parameter:: con_jcal   =4.1855E+0      ! joules per calorie  ()
-!  real(r_kind),parameter:: con_rhw0   =1022.0         ! sea water reference density (kg/m^3)
-!  real(r_kind),parameter:: con_epsq   =1.0E-12        ! min q for computing precip type
-
-!  Secondary constants
-
-  real(r_kind):: con_rocp
-  real(r_kind):: con_cpor
-  real(r_kind):: con_rog
-  real(r_kind):: con_fvirt
-  real(r_kind):: con_eps
-  !real(r_kind),parameter:: con_epsm1  =con_rd/con_rv-1.
-  !real(r_kind),parameter:: con_dldt   =con_cvap-con_cliq
-  !real(r_kind),parameter:: con_xpona  =-con_dldt/con_rv
-  !real(r_kind),parameter:: con_xponb  =-con_dldt/con_rv+con_hvap/(con_rv*con_ttp)
-
-!  --- ...  Other Physics/Chemistry constants (source: 2002 CODATA)
-
-!  real(r_kind),parameter:: con_c      =2.99792458e+8  ! speed of light      (m/s)
-!  real(r_kind),parameter:: con_plnk   =6.6260693e-34  ! planck constatn     (J/s)
-!  real(r_kind),parameter:: con_boltz  =1.3806505e-23  ! boltzmann constant  (J/K)
-!  real(r_kind),parameter:: con_sbc    =5.670400e-8    ! stefan-boltzmann    (W/m2/K4)
-!  real(r_kind),parameter:: con_avgd   =6.0221415e23   ! avogadro constant   (1/mol)
-!  real(r_kind),parameter:: con_gasv   =22413.996e-6   ! vol of ideal gas at 273.15k, 101.325kpa (m3/mol)
-!  real(r_kind),parameter:: con_amd    =28.9644        ! molecular wght of dry air (g/mol)
-!  real(r_kind),parameter:: con_amw    =18.0154        ! molecular wght of water vapor (g/mol)
-!  real(r_kind),parameter:: con_amo3   =47.9982        ! molecular wght of o3  (g/mol)
-!  real(r_kind),parameter:: con_amco2  =44.011         ! molecular wght of co2 (g/mol)
-!  real(r_kind),parameter:: con_amo2   =31.9999        ! molecular wght of o2  (g/mol)
-!  real(r_kind),parameter:: con_amch4  =16.043         ! molecular wght of ch4 (g/mol)
-!  real(r_kind),parameter:: con_amn2o  =44.013         ! molecular wght of n2o (g/mol)
-
-  contains
-    subroutine set_physical_constants(radius, omega, R, Rd, Rv, g, Cp, C_cond)bind(c, name='gfs_set_constants')
-        real(c_double), intent(in):: radius, omega, R, Rd, Rv, g, Cp, C_cond
-
-        con_rerth = radius
-        con_omega = omega
-        con_rgas = R
-        con_rd = Rd
-        con_rv = Rv
-        con_g = g
-        con_cp = Cp
-        con_cvap = C_cond
-
-        con_rocp = con_rd/con_cp
-        con_cpor = con_cp/con_rd
-        con_rog = con_rd/con_g
-        con_fvirt = con_rv/con_rd-1.
-        con_eps = con_rd/con_rv
-
-    end subroutine
-
-end module physcons
diff --git a/climt/_lib/GFS/pressure_data.f90 b/climt/_lib/GFS/pressure_data.f90
deleted file mode 100644
index 453b343e..00000000
--- a/climt/_lib/GFS/pressure_data.f90
+++ /dev/null
@@ -1,135 +0,0 @@
-module pressure_data
-! data for model level pressures and pressure-related
-! quantities.
-! init_pressdata: allocates arrays.
-! calc_pressdata: computes pressure related vars given ln(psg).
-! (called by getdyntend in module dyn_run).
-! destroy_pressdata: deallocate arrays.
- use params, only: nlons,nlats,nlevs,toa_pressure
- use kinds, only: r_kind
- use physcons, only: con_rd,con_cp,rk => con_rocp
- use iso_c_binding
- implicit none
- private
- public :: ak,bk,ck,dbk,bkl,psg,pk,alfa,rlnp,dpk,&
-  si,sl,prs, calc_pressdata, destroy_pressdata
-! ak,bk: definition of hybrid sigma-pressure coordinates
-! si,sl: equivalent values of sigma coordinate at interfaces, mid-levels.
-! dbk(k) = bk(k+1)-bk(k)
-! bkl(k) = 0.5*(bk(k+1)+bk(k))
-! ck(k)  = ak(k+1)*bk(k)-ak(k)*bk(k+1)
-! above 1-d arrays computed in init_dyn, defined from top to bottom
-! of model.
-! the following are definied by call to calc_pressdata:
-! psg: surface pressure (Pa)
-! pk: model interface pressures = ak(k) + bk(k)*psg(:,:) (k=1,nlevs+1)
-! dpk = pk(:,:,k+1)-pk(:,:,k)
-! prs: model layer pressures
-! rlnp = log(pk(:,:,k+1)/pk(:,:,k))
-! alfa = 1.-(pk(:,:,k)/dpk(:,:,k))*rlnp(:,:,k) (for k=1 alfa(:,:,1)=log(2))
-! all of the above arrays are defined top to bottom, except
-! prs which is bottom to top.
- real(c_double), dimension(:), pointer :: ak => null()
- real(c_double), dimension(:), pointer :: bk => null()
- real(c_double), dimension(:), pointer :: ck => null()
- real(c_double), dimension(:), pointer :: dbk => null()
- real(c_double), dimension(:), pointer :: si => null()
- real(c_double), dimension(:), pointer :: sl => null()
- real(c_double), dimension(:), pointer :: bkl => null()
-! real(c_double), dimension(:,:), allocatable :: psg
- real(c_double), dimension(:,:,:), pointer :: alfa => null()
- real(c_double), dimension(:,:,:), pointer :: rlnp => null()
- real(c_double), dimension(:,:,:), pointer :: dpk => null()
-!JOY converted prs, pk and psg to pointer, memory allocated by python
- real(c_double), dimension(:,:,:), pointer :: prs => null()
- real(c_double), dimension(:,:,:), pointer :: pk => null()
- real(c_double), dimension(:,:), pointer :: psg => null()
-
- contains
-
-     subroutine assign_coordinate_arrays(py_ak, py_bk, py_ck, py_dbk,&
-             py_si, py_sl, py_bkl, py_alfa, py_rlnp, py_dpk)bind(c, name='gfs_init_coordinate_arrays')
-
-         real(c_double), intent(inout), target:: py_ak(nlevs+1), py_bk(nlevs+1), py_ck(nlevs),&
-             py_dbk(nlevs), py_si(nlevs+1), py_sl(nlevs), py_bkl(nlevs)
-         real(c_double), intent(inout), target:: py_alfa(nlons, nlats, nlevs), &
-             py_rlnp(nlons, nlats, nlevs), py_dpk(nlons, nlats, nlevs)
-
-         ak => py_ak
-         bk => py_bk
-         ck => py_ck
-         dbk => py_dbk
-         si => py_si
-         sl => py_sl
-         bkl => py_bkl
-
-         alfa => py_alfa
-         rlnp => py_rlnp
-         dpk => py_dpk
-
-     end subroutine
-
-     subroutine assignPressureArrays(pySurfPressure, pyPressGrid, pyInterfacePressure) bind(c, name='gfs_assign_pressure_arrays')
-
-         real(c_double), intent(inout), target :: pyPressGrid(nlons,nlats,nlevs)
-         real(c_double), intent(inout), target :: pyInterfacePressure(nlons,nlats,nlevs)
-         real(c_double), intent(inout), target :: pySurfPressure(nlons,nlats)
-
-         prs => pyPressGrid
-         pk => pyInterfacePressure
-         psg => pySurfPressure
-
-     end subroutine assignPressureArrays
-
-
-     subroutine destroy_pressdata()
-         !print *, 'Deallocating pressure data'
-         nullify(ak,bk,ck,dbk,si,sl,bkl)
-         nullify(alfa,rlnp,dpk)
-         !JOY these are just pointers, no deallocation needed
-         nullify(prs,psg,pk)
-     end subroutine destroy_pressdata
-
-     subroutine updatePressure() bind(c, name='gfs_calculate_pressure')
-         use grid_data, only: lnpsg
-
-         call calc_pressdata(lnpsg)
-
-     end subroutine updatePressure
-
-     subroutine calc_pressdata(lnpsg)
-         real(r_kind),  intent(in) :: lnpsg(nlons,nlats)
-         ! update pressure related variables using latest estimate of lnps
-         integer k
-         ! surface press from log(ps)
-         psg = exp(lnpsg)
-         ! interface pressure (psg is ps)
-         ! ak,bk go top to bottom, so does pk
-         !$omp parallel
-         !$omp do private(k)
-         do k=1,nlevs+1
-         pk(:,:,k)=ak(k) + bk(k)*(psg(:,:) - toa_pressure)
-         enddo
-         !$omp end do
-         !$omp do private(k)
-         do k=1,nlevs
-         ! layer pressure thickness
-         dpk(:,:,k)=    pk(:,:,k+1) - pk(:,:,k)
-         ! sela's layer pressure from hyb2press.f
-         ! (goes from bottom to top, unlike pk)
-         prs(:,:,nlevs-k+1) = ((pk(:,:,k+1)**rk*pk(:,:,k+1) - pk(:,:,k)**rk*pk(:,:,k))/&
-             ((rk+1.)*dpk(:,:,k))) ** (1./rk)
-         enddo
-         !$omp end do
-         !$omp end parallel
-         alfa(:,:,1)=log(2.)
-         rlnp(:,:,1)=99999.99 !doesn't matter, should never be used.
-         !$omp parallel do private(k)
-         do k=2,nlevs
-         rlnp(:,:,k)= log( pk(:,:,k+1)/pk(:,:,k) )
-         alfa(:,:,k)= 1.-( pk(:,:,k)/dpk(:,:,k) )*rlnp(:,:,k)
-         enddo
-         !$omp end parallel do
-     end subroutine calc_pressdata
-
- end module pressure_data
diff --git a/climt/_lib/GFS/run.f90 b/climt/_lib/GFS/run.f90
deleted file mode 100644
index 6048465c..00000000
--- a/climt/_lib/GFS/run.f90
+++ /dev/null
@@ -1,375 +0,0 @@
-module run_mod
-! time step loop for model run.
-! Public subroutines:
-! run: main time step loop (advances model state, writes out
-! data at specified intervals).
-
-use kinds, only: r_kind,r_double
-use params, only: ndimspec, nlevs, ntmax, tstart, dt, nlons, nlats, nlevs,&
-  fhzer,ntrac,ntrunc,ntout, explicit, adiabatic, ntrac, &
-  massfix,gfsio_out, sigio_out, sfcinitfile, ntdfi, shum, svc, sppt,&
-  addnoise, addnoise_vfilt,addnoise_kenorm, addnoise_lscale, addnoise_dissfact,&
-  addnoise_vrtonly
-use shtns, only: grdtospec, spectogrd, lats, areawts, lap, invlap, degree, getuv
-use dyn_run, only: getdyntend, dry_mass_fixer
-use phy_run, only: getphytend
-use phy_data, only: wrtout_sfc, wrtout_flx, init_phydata, pwat
-!! JOY removing calls to write data
-!use dyn_init, only: wrtout_sig, wrtout_gfsgrb
-use spectral_data, only:  lnpsspec, vrtspec, divspec, virtempspec,&
-                          tracerspec, disspec, dmp_prof, diff_prof
-! JOY remove Incremental Analysis Unit (IAU) forcing                          
-!use iau_module, only: getiauforcing, init_iau, init_iauialized
-! these arrays used to print diagnostics after each time step.
-use pressure_data, only: psg
-use grid_data, only: ug,vg,dlnpsdt
-use physcons, only: rerth => con_rerth
-
-use iso_c_binding
-implicit none
-private
-
-real(c_double) :: t
-bind(c) :: t
-public :: take_one_step, t
-
-contains
-
-subroutine set_model_time(py_time) bind(c, name='gfs_set_model_time')
-    real(c_double), intent(in):: py_time
-    t = py_time
-    !print *, 'model time is now: ', t
-end subroutine
-
-
-subroutine take_one_step() bind(c,name='gfs_take_one_step')
-
-    !integer nt,ntstart
-    !real(r_kind) fh,pstendmean,spdmax
-    real(r_kind) fh
-    real(8) t1,t2
-    !real(r_kind), dimension(nlons,nlats) :: pstend
-    integer(8) count, count_rate, count_max
-
-    call system_clock(count, count_rate, count_max)
-    t1 = count*1.d0/count_rate
-    ! advance solution with RK
-    call advance(t)
-    fh = t/3600.
-    t = t + dt ! absolute forecast time.
-    call system_clock(count, count_rate, count_max)
-    t2 = count*1.d0/count_rate
-    !spdmax = maxval(sqrt(ug**2+vg**2)) ! max wind speed
-    !pstend = (36.*psg*dlnpsdt)**2 ! ps tend variance (mb/hr)**2
-    !pstendmean = sqrt(sum(pstend*areawts))
-    !if (modulo(fh, 10.) .eq. 0.) then
-    !    write(6,8998) fh,spdmax,minval(psg/100.),maxval(psg/100.),pstendmean,t2-t1
-    !    8998 format('t = ',f0.3,' hrs, spdmax = ',f7.3,', min/max ps = ',f7.2,'/',f7.2,', pstend = ',f0.3,', cpu time = ',f0.3)
-    !endif
-
-end subroutine take_one_step
-
-subroutine advance(t)
-    ! advance model state to next time step using 3-stage 2nd-order additive
-    ! (i.e. implicit-explicit (IMEX) or semi-implicit) runge-kutta scheme.
-!
-! Double Butcher tableau looks like this:
-!  *explicit part*                *implicit part* 
-!
-!   0    | 0                       0  | 0
-!  a21   | a21 0             aa21+aa22| aa21 aa22
-! a31+a32| a31 a32 0    aa31+aa32+aa33| aa31 aa32 aa33
-!   1    | b1  b2  b3 0            1  | bb1  bb2   bb3 bb4
-! ---------------------  ----------------------------------
-!          b1  b2  b3 0                 bb1  bb2   bb3 bb4
-!
-! the values of the coefficients are defined in semimp_data.f90
-!
-! progression in time: to time t0+a21*dt, then t0+(a31+a32)*dt, then dt.
-! for progression in time to be the same in implicit and
-! explicit parts, a21=aa21+a22 and a31+a32=aa31+aa32+aa33
-! consistency condition:
-! b1 + b2 + b3 = bb1 + bb2 + bb3 + bb4 = 1
-! accuracy condition (2nd order):
-! a21*b2+(a31+a32)*b3 = 1/2
-! (aa21+aa22)*bb2+(aa31+aa32+aa33)*bb3+bb4 = 1/2
-
-! JOY remove patterngenerator and stoch references    
-!  use patterngenerator, only: patterngenerator_advance
-!  use stoch_data, only: rpattern_svc,rpattern_sppt,&
-!     spec_svc,spec_sppt,grd_svc,grd_sppt,&
-!     spec_shum,grd_shum,rpattern_shum,getstochforcing,&
-!     rpattern_addnoise,specpsi_addnoise,vfact_addnoise
-  use semimp_data, only: amhyb,bmhyb,svhyb,d_hyb_m,tor_hyb,&
-                         a21,a31,a32,aa21,aa22,aa31,aa32,aa33,b1,b2,b3,bb1,bb2,bb3,bb4
-  real(r_double), intent(in) :: t
-  complex(r_kind),dimension(ndimspec) :: &
-      lnpsspec_orig, dlnpsdtlin1, dlnpsdtlin2,&
-      dlnpsspecdt_orig,dlnpsdtlin_orig,&
-      dlnpsspecdt1, dlnpsspecdt2
-  complex(r_kind), dimension(ndimspec,nlevs) :: &
-      vrtspec_orig,divspec_orig,virtempspec_orig,&
-      dvrtspecdt_orig,ddivspecdt_orig,dvirtempspecdt_orig,&
-      dvrtspecdt1,ddivspecdt1,dvirtempspecdt1,&
-      dvrtspecdt2,ddivspecdt2,dvirtempspecdt2,&
-      psiforcing,tvforcing,ddivdtlin_orig, dtvdtlin_orig,&
-      ddivdtlin1, ddivdtlin2, dtvdtlin1, dtvdtlin2
-  complex(r_kind), dimension(ndimspec,nlevs,ntrac) :: &
-      dtracerspecdt_orig,&
-      tracerspec_orig,dtracerspecdt1,dtracerspecdt2
-  complex(r_kind) :: rhs(nlevs)
-  integer nt, k, n
-  !logical :: profile = .false. ! print out timing stats
-  integer(8) count, count_rate, count_max
-  real(8) t1,t2,t0,t4
-
-
-  ! save original fields.
-!$omp workshare
-  vrtspec_orig = vrtspec
-  divspec_orig = divspec
-  virtempspec_orig = virtempspec
-  tracerspec_orig = tracerspec
-!$omp end workshare
-  lnpsspec_orig = lnpsspec
-  call system_clock(count, count_rate, count_max)
-  t0 = count*1.d0/count_rate
-
-
-  ! stage 1
-  ! compute dynamics tendencies.
-  call system_clock(count, count_rate, count_max)
-  t1 = count*1.d0/count_rate
-  call getdyntend(dvrtspecdt_orig,ddivspecdt_orig,dvirtempspecdt_orig,dtracerspecdt_orig,dlnpsspecdt_orig,1)
-  call system_clock(count, count_rate, count_max)
-  t2 = count*1.d0/count_rate
-
-!$omp workshare
-  vrtspec=vrtspec_orig+a21*dt*dvrtspecdt_orig
-  tracerspec=tracerspec_orig+a21*dt*dtracerspecdt_orig
-!$omp end workshare
-  if (explicit) then
-!$omp workshare
-     divspec=divspec_orig+a21*dt*ddivspecdt_orig
-     virtempspec=virtempspec_orig+a21*dt*dvirtempspecdt_orig
-     lnpsspec=lnpsspec_orig+a21*dt*dlnpsspecdt_orig
-!$omp end workshare
-  else
-! solve for updated divergence.
-! back subsitution to get updated virt temp, lnps.
-!$omp parallel do private(n,rhs)
-  do n=1,ndimspec
-! first remove linear terms from computed tendencies.
-      ddivdtlin_orig(n,:) = -lap(n)*& 
-      (matmul(amhyb,virtempspec_orig(n,:)) + tor_hyb(:)*lnpsspec_orig(n))
-      ddivspecdt_orig(n,:) = ddivspecdt_orig(n,:) - ddivdtlin_orig(n,:)
-      dtvdtlin_orig(n,:) = -matmul(bmhyb,divspec_orig(n,:))
-      dvirtempspecdt_orig(n,:) = dvirtempspecdt_orig(n,:) - dtvdtlin_orig(n,:)
-      dlnpsdtlin_orig(n) = -sum(svhyb(:)*divspec_orig(n,:))
-      dlnpsspecdt_orig(n) = dlnpsspecdt_orig(n) - dlnpsdtlin_orig(n)
-! compute explicit part of update
-      divspec(n,:) = divspec_orig(n,:) + dt*(a21*ddivspecdt_orig(n,:)  +&
-                                             aa21*ddivdtlin_orig(n,:))
-      virtempspec(n,:) = virtempspec_orig(n,:) + dt*(a21*dvirtempspecdt_orig(n,:)  +&
-                                                     aa21*dtvdtlin_orig(n,:))
-      lnpsspec(n) = lnpsspec_orig(n) + dt*(a21*dlnpsspecdt_orig(n) +&
-                                           aa21*dlnpsdtlin_orig(n))
-! solve implicit part for updated divergence.
-      rhs(:)   = divspec(n,:) - aa22*lap(n)*dt*&
-                 (matmul(amhyb,virtempspec(n,:)) + tor_hyb(:)*lnpsspec(n))
-      divspec(n,:) = matmul(d_hyb_m(:,:,degree(n)+1,1),rhs(:))
-! back substitute to get implicit contribution to virtual temp and lnps update.
-      virtempspec(n,:) = virtempspec(n,:) - aa22*dt*matmul(bmhyb,divspec(n,:))
-      lnpsspec(n) = lnpsspec(n) - aa22*dt*sum(svhyb(:)*divspec(n,:))
-  enddo
-!$omp end parallel do 
-  endif
-
-  ! stage 2
-  ! compute dynamics tendencies.
-  call system_clock(count, count_rate, count_max)
-  t1 = count*1.d0/count_rate
-  call getdyntend(dvrtspecdt1,ddivspecdt1,dvirtempspecdt1,dtracerspecdt1,dlnpsspecdt1,2)
-  call system_clock(count, count_rate, count_max)
-  t2 = count*1.d0/count_rate
-  !if (profile) print *,'time in getdyntend (stage 2) =',t2-t1
-
-  ! update vorticity and tracers (always explicit)
-!$omp workshare
-  vrtspec=vrtspec_orig+dt*(a31*dvrtspecdt_orig+a32*dvrtspecdt1)
-  tracerspec=tracerspec_orig+dt*(a31*dtracerspecdt_orig+a32*dtracerspecdt1)
-!$omp end workshare
-  if (explicit) then
-!$omp workshare
-     divspec=divspec_orig+dt*(a31*ddivspecdt_orig+a32*ddivspecdt1)
-     virtempspec=virtempspec_orig+dt*(a31*dvirtempspecdt_orig+a32*dvirtempspecdt1)
-     lnpsspec=lnpsspec_orig+dt*(a31*dlnpsspecdt_orig+a32*dlnpsspecdt1)
-!$omp end workshare
-  else
-! solve for updated divergence.
-! back subsitution to get updated virt temp, lnps.
-!$omp parallel do private(n,rhs)
-  do n=1,ndimspec
-! remove linear terms from computed tendencies.
-      ddivdtlin1(n,:) = -lap(n)*& 
-      (matmul(amhyb,virtempspec(n,:)) + tor_hyb(:)*lnpsspec(n))
-      ddivspecdt1(n,:) = ddivspecdt1(n,:) - ddivdtlin1(n,:)
-      dtvdtlin1(n,:) = -matmul(bmhyb,divspec(n,:))
-      dvirtempspecdt1(n,:) = dvirtempspecdt1(n,:) - dtvdtlin1(n,:)
-      dlnpsdtlin1(n) = -sum(svhyb(:)*divspec(n,:))
-      dlnpsspecdt1(n) = dlnpsspecdt1(n) - dlnpsdtlin1(n)
-! compute explicit part of update
-      divspec(n,:) = divspec_orig(n,:) + dt*(a31*ddivspecdt_orig(n,:)  +&
-                                             aa31*ddivdtlin_orig(n,:)  +&
-                                             a32*ddivspecdt1(n,:)      +&
-                                             aa32*ddivdtlin1(n,:))   
-      virtempspec(n,:) = virtempspec_orig(n,:) + dt*(a31*dvirtempspecdt_orig(n,:)  +&
-                                                     aa31*dtvdtlin_orig(n,:)       +&
-                                                     a32*dvirtempspecdt1(n,:)      +&
-                                                     aa32*dtvdtlin1(n,:))
-      lnpsspec(n) = lnpsspec_orig(n) + dt*(a31*dlnpsspecdt_orig(n)     +&
-                                           aa31*dlnpsdtlin_orig(n)     +&
-                                           a32*dlnpsspecdt1(n)         +&
-                                           aa32*dlnpsdtlin1(n))     
-! solve implicit part for updated divergence.
-      rhs(:)   = divspec(n,:) - aa33*lap(n)*dt*&
-                 (matmul(amhyb,virtempspec(n,:)) + tor_hyb(:)*lnpsspec(n))
-      divspec(n,:) = matmul(d_hyb_m(:,:,degree(n)+1,2),rhs(:))
-! back substitute to get implicit contribution to virtual temp and lnps update.
-      virtempspec(n,:) = virtempspec(n,:) - aa33*dt*matmul(bmhyb,divspec(n,:))
-      lnpsspec(n) = lnpsspec(n) - aa33*dt*sum(svhyb(:)*divspec(n,:))
-  enddo
-!$omp end parallel do 
-  endif
-
-  ! stage 3
-  ! compute dynamics tendencies.
-  call system_clock(count, count_rate, count_max)
-  t1 = count*1.d0/count_rate
-  call getdyntend(dvrtspecdt2,ddivspecdt2,dvirtempspecdt2,dtracerspecdt2,dlnpsspecdt2,3)
-  call system_clock(count, count_rate, count_max)
-  t2 = count*1.d0/count_rate
-  !if (profile) print *,'time in getdyntend (stage 3) =',t2-t1
-
-  ! final update of vorticity and tracers (always explicit)
-!$omp workshare
-  vrtspec=vrtspec_orig+dt*(b1*dvrtspecdt_orig+b2*dvrtspecdt1+b3*dvrtspecdt2)
-  tracerspec=tracerspec_orig+dt*(b1*dtracerspecdt_orig+b2*dtracerspecdt1+b3*dtracerspecdt2)
-!$omp end workshare
-  if (explicit) then
-!$omp workshare
-     divspec=divspec_orig+dt*(b1*ddivspecdt_orig+b2*ddivspecdt1+b3*ddivspecdt2)
-     virtempspec=virtempspec_orig+dt*(b1*dvirtempspecdt_orig+b2*dvirtempspecdt1+b3*dvirtempspecdt2)
-     lnpsspec=lnpsspec_orig+dt*(b1*dlnpsspecdt_orig+b2*dlnpsspecdt1+b3*dlnpsspecdt2)
-!$omp end workshare
-  else
-!$omp parallel do private(n,rhs)
-  do n=1,ndimspec
-! remove linear terms from computed tendencies.
-      ddivdtlin2(n,:) = -lap(n)*& 
-      (matmul(amhyb,virtempspec(n,:)) + tor_hyb(:)*lnpsspec(n))
-      ddivspecdt2(n,:) = ddivspecdt2(n,:) - ddivdtlin2(n,:)
-      dtvdtlin2(n,:) = -matmul(bmhyb,divspec(n,:))
-      dvirtempspecdt2(n,:) = dvirtempspecdt2(n,:) - dtvdtlin2(n,:)
-      dlnpsdtlin2(n) = -sum(svhyb(:)*divspec(n,:))
-      dlnpsspecdt2(n) = dlnpsspecdt2(n) - dlnpsdtlin2(n)
-! compute final update
-      divspec(n,:) = divspec_orig(n,:) + dt*(b1*ddivspecdt_orig(n,:)  +&
-                                             bb1*ddivdtlin_orig(n,:)  +& 
-                                             b2*ddivspecdt1(n,:)      +&
-                                             bb2*ddivdtlin1(n,:)      +& 
-                                             b3*ddivspecdt2(n,:)      +&
-                                             bb3*ddivdtlin2(n,:))       
-      virtempspec(n,:) = virtempspec_orig(n,:) + dt*(b1*dvirtempspecdt_orig(n,:)  +&
-                                                     bb1*dtvdtlin_orig(n,:)       +& 
-                                                     b2*dvirtempspecdt1(n,:)      +&
-                                                     bb2*dtvdtlin1(n,:)           +& 
-                                                     b3*dvirtempspecdt2(n,:)      +&
-                                                     bb3*dtvdtlin2(n,:))       
-      lnpsspec(n) = lnpsspec_orig(n) + dt*(b1*dlnpsspecdt_orig(n) +&
-                                           bb1*dlnpsdtlin_orig(n) +&
-                                           b2*dlnpsspecdt1(n)     +&
-                                           bb2*dlnpsdtlin1(n)     +&
-                                           b3*dlnpsspecdt2(n)     +&
-                                           bb3*dlnpsdtlin2(n))
-! solve implicit part for updated divergence.
-      if (abs(bb4) .gt. 1.e-5) then
-         rhs(:)   = divspec(n,:) - bb4*lap(n)*dt*&
-                    (matmul(amhyb,virtempspec(n,:)) + tor_hyb(:)*lnpsspec(n))
-         divspec(n,:) = matmul(d_hyb_m(:,:,degree(n)+1,3),rhs(:))
-! back substitute to get implicit contribution to virtual temp and lnps update.
-         virtempspec(n,:) = virtempspec(n,:) - bb4*dt*matmul(bmhyb,divspec(n,:))
-         lnpsspec(n) = lnpsspec(n) - bb4*dt*sum(svhyb(:)*divspec(n,:))
-      endif
-  enddo
-!$omp end parallel do 
-  endif
-
-  ! forward implicit treatment of linear damping/diffusion
-  !$omp parallel do private(k,nt)
-  do k=1,nlevs
-     vrtspec(:,k) = vrtspec(:,k)/(1. - (disspec(:)*diff_prof(k) - dmp_prof(k))*dt)
-     divspec(:,k) = divspec(:,k)/(1. - (disspec(:)*diff_prof(k) - dmp_prof(k))*dt)
-     virtempspec(:,k) = virtempspec(:,k)/(1. - disspec(:)*diff_prof(k)*dt)
-     do nt=1,ntrac
-        tracerspec(:,k,nt) = tracerspec(:,k,nt)/(1. - disspec(:)*diff_prof(k)*dt)
-     enddo
-  enddo
-  !$omp end parallel do
-
-  call system_clock(count, count_rate, count_max)
-  t4 = count*1.d0/count_rate
-  !if (profile) print *,'time in dynamics update=',t4-t0
-  ! apply physics parameterizations as an adjustment to fields updated by dynamics.
-  ! (time-split approach)
-  if (.not. adiabatic) then
-     ! update variables on grid.
-     call getdyntend(dvrtspecdt1,ddivspecdt1,dvirtempspecdt1,&
-          dtracerspecdt1,dlnpsspecdt1,3,.true.) ! <- .true. for spectogrd only
-     ! compute stochastic backscatter forcing, if needed
-     ! if addnoise_dissfact is true,
-     ! streamfunctione noise forcing by smoothed estimate of ke dissipation
-     ! rate.
-     ! Joy removing stochastic forcing routines
-     !if (addnoise > tiny(addnoise)) &
-     !call getstochforcing(specpsi_addnoise, psiforcing, tvforcing)
-     ! compute physics tendencies, apply as an adjustment to updated state
-     call system_clock(count, count_rate, count_max)
-     t1 = count*1.d0/count_rate
-     call getphytend(dvrtspecdt1,ddivspecdt1,dvirtempspecdt1,dtracerspecdt1,dlnpsspecdt1,t,dt)
-     call system_clock(count, count_rate, count_max)
-     t2 = count*1.d0/count_rate
-     ! if (profile) print *,'time in getphytend=',t2-t1
-!$omp workshare
-     vrtspec=vrtspec+dt*dvrtspecdt1
-     divspec=divspec+dt*ddivspecdt1
-     virtempspec=virtempspec+dt*dvirtempspecdt1
-     tracerspec=tracerspec+dt*dtracerspecdt1
-!$omp end workshare
-! modify spectral lnps tendency to include contribution
-! from dry mass 'fixer' (adjusts dry surface pressure to
-! remain equal to pdryini).
-     if (massfix) then
-        call dry_mass_fixer(psg,pwat,dlnpsspecdt1,dt)
-     endif
-     lnpsspec=lnpsspec+dt*dlnpsspecdt1 ! only needed for dry mass fixer.
-     ! add stochastic backsscatter forcing in vorticity eqn.
-     ! psiforcing has units of streamfunction tendency
-     if (addnoise > tiny(addnoise)) then
-         !$omp parallel do private(k)
-         do k=1,nlevs
-            if (addnoise_kenorm .and. addnoise_vrtonly) then
-               !vrtspec(:,k) = vrtspec(:,k) + dt*(lap/rerth**2)*sqrt(-invlap*rerth**2)*psiforcing(:,k)
-               vrtspec(:,k) = vrtspec(:,k) - dt*(sqrt(-lap)/rerth)*psiforcing(:,k)
-            else
-               vrtspec(:,k) = vrtspec(:,k) + dt*(lap/rerth**2)*psiforcing(:,k)
-            endif
-            virtempspec(:,k) = virtempspec(:,k) + dt*tvforcing(:,k)
-         enddo
-         !$omp end parallel do 
-     end if
-  end if
-
-end subroutine advance
-
-end module run_mod
diff --git a/climt/_lib/GFS/semimp_data.f90 b/climt/_lib/GFS/semimp_data.f90
deleted file mode 100644
index ca444695..00000000
--- a/climt/_lib/GFS/semimp_data.f90
+++ /dev/null
@@ -1,284 +0,0 @@
-module semimp_data
-    ! static data arrays used in semi-implicit time integration scheme.
-    ! Public subroutines:
-    ! init_semimpdata: allocate and populate arrays.
-    ! destroy_semimpdata: deallocate arrays.
-    use kinds, only: r_kind, default_real, r_double
-    use params, only: dt,ntrunc,nlons,nlats,nlevs
-    use pressure_data, only: ak,bk
-    use physcons, only: rd => con_rd, cp => con_cp, rerth => con_rerth,&
-        kappa => con_rocp
-    use iso_c_binding
-    implicit none
-    private
-    public :: init_semimpdata, destroy_semimpdata
-    real(r_kind), pointer, public, dimension(:,:,:,:) :: d_hyb_m
-    real(r_kind), pointer, public, dimension(:,:) :: amhyb,bmhyb
-    real(r_kind), pointer, public, dimension(:) ::  &
-        tref,pkref,dpkref,alfaref,svhyb,tor_hyb
-
-    ! scheme derived from Ascher, Spiteri and Ruuth 1997.
-    ! Applied Numerical Mathematics DOI:10.1016/S0168-9274(97)00056-1
-    ! it's a combination of schemes in sxns 2.5 and 2.6
-    ! from that paper (tableaus from 2.5, value of delta from 2.6)
-    !  
-    ! Double tableau:  c = sqrt(2)/2 = 1/sqrt(2)
-    ! *explicit part*  *implicit part* 
-    !
-    !   0 | 0             0 | 0
-    ! 1-c | 1-c 0       1-c | 0 1-c
-    !   1 | -c 1+c 0      1 | 0 c 1-c
-    ! --------------   ---------------
-    !       0  c 1-c          0 c 1-c
-    !
-    ! when used to discretize 'fast'/'slow' double oscillator eqn
-    ! this scheme is:
-    ! stable for 'slow' courant numbers < 2*c = sqrt(2) ~ 1.414
-    ! stable for all 'fast' courant numbers
-
-
-    !real(r_kind),parameter :: halfroot2    = sqrt(2.)/2.
-    !real(r_kind),parameter, public :: a21  = 1.-halfroot2
-    !real(r_kind),parameter, public :: a31  = -halfroot2
-    !real(r_kind),parameter, public :: a32  = 1.+halfroot2
-    !real(r_kind),parameter, public :: aa21 = 0.
-    !real(r_kind),parameter, public :: aa22 = 1.-halfroot2 
-    !real(r_kind),parameter, public :: aa31 = 0.
-    !real(r_kind),parameter, public :: aa32 = halfroot2
-    !real(r_kind),parameter, public :: aa33 = 1.-halfroot2 
-    !real(r_kind),parameter, public :: b1   = 0.  
-    !real(r_kind),parameter, public :: b2   = halfroot2
-    !real(r_kind),parameter, public :: b3   = 1.-halfroot2
-    !real(r_kind),parameter, public :: bb1  = 0. 
-    !real(r_kind),parameter, public :: bb2  = halfroot2
-    !real(r_kind),parameter, public :: bb3  = 1.-halfroot2
-    !real(r_kind),parameter, public :: bb4  = 0.
-
-    ! Lorenz 3-cycle + trapezoidal
-    !
-    ! has very weak instability for purely imaginary double
-    ! oscillator.
-    ! Largest amp factor for 'slow' courant # < sqrt(3) ~ 1.73,
-    ! which is not a problem for models with weak damping.
-    !
-    ! Double tableau:
-    !
-    !  *explicit part*                *implicit part* 
-    !
-    !   0  | 0                       0  | 0
-    !  1/3 | 1/3  0                 1/3 | 1/6 1/6 0
-    !  2/3 | 1/6  1/2 0             2/3 | 1/3  0 1/3  0
-    !   1  | 1/2 -1/2 1  0           1  | 3/8  0 3/8 1/4 
-    ! ---------------------         ----------------------
-    !        1/2 -1/2 1  0                3/8  0 3/8 1/4
-
-    ! real(r_kind),parameter, public :: a21  = 1./3.
-    ! real(r_kind),parameter, public :: a31  = 1./6.
-    ! real(r_kind),parameter, public :: a32  = 1./2.
-    ! real(r_kind),parameter, public :: aa21 = 1./6.
-    ! real(r_kind),parameter, public :: aa22 = 1./6.
-    ! real(r_kind),parameter, public :: aa31 = 1./3.
-    ! real(r_kind),parameter, public :: aa32 = 0.
-    ! real(r_kind),parameter, public :: aa33 = 1./3.
-    ! real(r_kind),parameter, public :: b1   = 1./2.
-    ! real(r_kind),parameter, public :: b2   = -1./2.
-    ! real(r_kind),parameter, public :: b3   = 1.0 
-    ! real(r_kind),parameter, public :: bb1  = 3./8.
-    ! real(r_kind),parameter, public :: bb2  = 0.
-    ! real(r_kind),parameter, public :: bb3  = 3./8.  
-    ! real(r_kind),parameter, public :: bb4  = 1./4.
-
-    ! ssp/tvd rk3 (shu 1988, http://epubs.siam.org/doi/abs/10.1137/0909073) IMEX
-    ! minimizes spurious oscillations - good for advection of tracers
-    ! (no negative entries in Butcher tableau)
-    ! stability similar to Lorenz 3-cycle IMEX
-    !
-    ! Double tableau:
-    !
-    !  *explicit part*                *implicit part* 
-    !
-    !   0  | 0                       0  | 0
-    !   1  | 1   0                   1  | 0.635  0.365 
-    !  1/2 | 1/4 1/4 0              1/2 | 0.3175 0.0   0.1825  
-    !   1  | 1/6 1/6 2/3 0           1  | 0.35   0.0   0.3    0.35
-    ! ---------------------         -------------------------------
-    !        1/6 1/6 2/3 0                0.35   0.0   0.3    0.35
-    !
-    real(r_kind),parameter, public :: a21  = 1.0
-    real(r_kind),parameter, public :: a31  = 0.25  
-    real(r_kind),parameter, public :: a32  = 0.25  
-    real(r_kind),parameter, public :: aa21 = 0.635
-    real(r_kind),parameter, public :: aa22 = 0.365
-    real(r_kind),parameter, public :: aa31 = 0.3175 
-    real(r_kind),parameter, public :: aa32 = 0.
-    real(r_kind),parameter, public :: aa33 = 0.1825  
-    real(r_kind),parameter, public :: b1   = 1./6.
-    real(r_kind),parameter, public :: b2   = 1./6.
-    real(r_kind),parameter, public :: b3   = 2./3.
-    real(r_kind),parameter, public :: bb1  = 0.35  
-    real(r_kind),parameter, public :: bb2  = 0.0
-    real(r_kind),parameter, public :: bb3  = 0.3    
-    real(r_kind),parameter, public :: bb4  = 0.35  
-
-    real(r_kind),parameter, public :: ref_temp = 300.
-    real(r_kind),parameter, public :: ref_press = 800.e2
-
-contains
-
-    subroutine init_semimp_arrays(py_d_hyb_m, py_amhyb, py_bmhyb,&
-        py_tref, py_pkref, py_dpkref, py_alfaref, py_svhyb, py_tor_hyb)bind(c, name='gfs_initialise_timestepper_arrays')
-
-    real(c_double), target :: py_d_hyb_m(nlevs, nlevs, ntrunc+1, 3)
-    real(c_double), target:: py_amhyb(nlevs, nlevs), py_bmhyb(nlevs, nlevs)
-    real(c_double), target:: py_tref(nlevs), py_pkref(nlevs+1), py_dpkref(nlevs),&
-        py_alfaref(nlevs), py_svhyb(nlevs), py_tor_hyb(nlevs)
-
-    d_hyb_m => py_d_hyb_m
-    amhyb => py_amhyb
-    bmhyb => py_bmhyb
-    tref => py_tref
-    pkref => py_pkref
-    dpkref => py_dpkref
-    alfaref => py_alfaref
-    svhyb => py_svhyb
-    tor_hyb => py_tor_hyb
-
-    end subroutine
-
-    subroutine init_semimpdata()
-        real(r_double) rnn1
-        integer irow,icol,icolbeg,i,j,k,n,icolend,nn,iret
-        real(r_double), allocatable, dimension(:,:) :: yecm,tecm,ym,rim
-        real(r_double), allocatable, dimension(:) :: vecm
-        integer, allocatable, dimension(:) :: ipiv
-        real(r_kind) epstiny,const
-
-        ! make sure progression in time the same in explicit and implicit
-        ! parts of additive runge-kutta scheme.
-        epstiny = 1.e-5
-        if (abs(a21-(aa21+aa22)) > epstiny .or. abs(a31+a32-(aa31+aa32+aa33)) > epstiny) then
-            print *,a21,a31+a32,abs(a21-(aa21+aa22)) 
-            print *,aa22,aa32+aa33,abs(a31+a32-(aa31+aa32+aa33))
-            print *,'time progression different in explicit and implicit parts of RK'
-            print *,'fix in semimp_data.f90!'
-            stop
-        endif
-        ! check to make sure b1+b2+b3=bb1+bb2+bb3+bb4=1
-        if (abs(b1+b2+b3-1.) > epstiny .or. abs(bb1+bb2+bb3+bb4-1) > epstiny) then
-            print *,'b1+b2+b3 or bb1+bb2+bb3+bb4 != 1'
-            print *,'fix in semimp_data.f90!'
-            stop
-        endif
-        ! check to make sure scheme is 2nd order
-        if (abs(a21*b2+(a31+a32)*b3-0.5) > epstiny .or.&
-            abs((aa21+aa22)*bb2+(aa31+aa32+aa33)*bb3+bb4-0.5) > epstiny) then
-            print *,a21*b2 + (a31+a32)*b3,abs(a21*b2+(a31+a32)*b3+bb4-0.5)
-            print *,(aa21+aa22)*bb2+(aa31+aa32+aa33)*bb3 ,&
-                abs((aa21+aa22)*bb2+(aa31+aa32+aa33)*bb3+bb4-0.5)
-            print *,'scheme not 2nd order'
-            print *,'fix in semimp_data.f90!'
-            stop
-        endif
-
-        ! module vars
-        allocate(ipiv(nlevs))
-        !allocate(amhyb(nlevs,nlevs),bmhyb(nlevs,nlevs))
-        !allocate(tref(nlevs),pkref(nlevs+1),dpkref(nlevs),alfaref(nlevs))
-        !allocate(svhyb(nlevs),tor_hyb(nlevs))
-        !allocate(d_hyb_m(nlevs,nlevs,ntrunc+1,3))
-        ! temp storage.
-        allocate(yecm(nlevs,nlevs),tecm(nlevs,nlevs))
-        allocate(ym(nlevs,nlevs))
-        allocate(vecm(nlevs))
-        allocate(rim(nlevs,nlevs)) ! identity matrix
-
-        tref = ref_temp
-        do k=1,nlevs+1
-        pkref(k) = ak(k) + bk(k)*ref_press
-        enddo
-        do k=1,nlevs
-        dpkref(k) = pkref(k+1)-pkref(k)
-        enddo
-        alfaref(1) = log(2.)
-        do k=2,nlevs
-        alfaref(k) = 1.-(pkref(k)/dpkref(k))*log(pkref(k+1)/pkref(k))
-        enddo
-        yecm=0.
-        do irow=1,nlevs
-        yecm(irow,irow)=alfaref(irow)*rd
-        icolbeg=irow+1
-        if(icolbeg<=nlevs)then
-            do icol=icolbeg,nlevs
-            yecm(irow,icol)=rd*log( pkref(icol+1)/pkref(icol) )
-            enddo
-        endif
-        enddo
-        tecm=0.
-        do irow=1,nlevs
-        tecm(irow,irow)=kappa*tref(irow)*alfaref(irow)
-        icolend=irow-1
-        do icol=1,icolend
-        tecm(irow,icol)=(kappa*tref(irow)*dpkref(icol)/ &
-            dpkref(irow))*log(pkref(irow+1)/pkref(irow))
-        enddo
-        enddo
-        vecm=dpkref/ref_press
-        ! amhyb is operator for linearized pres grad term (geopot part)
-        ! (operatores on virt temp)
-        ! bmhyb is operator for linearized energy conv term
-        ! (operates on div)
-        do j=1,nlevs
-        svhyb(j)=vecm(nlevs+1-j) ! for linearized lnps tend term (div)
-        do k=1,nlevs
-        amhyb(k,j)=yecm(nlevs+1-k,nlevs+1-j)
-        bmhyb(k,j)=tecm(nlevs+1-k,nlevs+1-j)
-        enddo
-        enddo
-        amhyb=amhyb/rerth**2 
-        ! times lnps in linearized pgf term (lapacian in div eqn)
-        tor_hyb=rd*tref/rerth**2
-
-        rim = 0.
-        do k=1,nlevs
-        rim(k,k) = 1.
-        enddo
-        ! computations that do not depend on wavenumber
-        do i=1,nlevs
-        do j=1,nlevs
-        ym(i,j) = tor_hyb(i)*svhyb(j)
-        enddo
-        do k=1,nlevs
-        do j=1,nlevs
-        ym(i,j) = ym(i,j) + amhyb(i,k)*bmhyb(k,j)
-        enddo
-        enddo
-        enddo
-
-        ! computations that do depend on wavenumber
-        ! enabling openmp for this loop doesn't work with intel MKL
-        do k=1,3
-        if (k .eq. 1)  const = aa22
-        if (k .eq. 2)  const = aa33
-        if (k .eq. 3)  const = bb4
-        !!$omp parallel do private(nn,n,rnn1,yecm,ipiv,iret,vecm)
-        do nn=1,ntrunc+1
-        n = nn-1
-        rnn1 = n*(n+1)
-        yecm = rim + (const*dt)**2*rnn1*ym
-        ! invert matrix using LAPACK, save in d_hyb_m
-        call dgetrf(nlevs,nlevs,yecm,nlevs,ipiv,iret)
-        call dgetri(nlevs,yecm,nlevs,ipiv,vecm,nlevs,iret)
-        d_hyb_m(:,:,nn,k) = yecm
-        enddo
-        !!$omp end parallel do 
-        enddo
-        deallocate(rim,yecm,tecm,ym,vecm,ipiv)
-
-    end subroutine init_semimpdata
-
-    subroutine destroy_semimpdata()
-        nullify(amhyb,bmhyb,d_hyb_m)
-        nullify(tref,pkref,dpkref,alfaref,svhyb,tor_hyb)
-    end subroutine destroy_semimpdata
-
-end module semimp_data
diff --git a/climt/_lib/GFS/shtns.f90 b/climt/_lib/GFS/shtns.f90
deleted file mode 100644
index e9d7d629..00000000
--- a/climt/_lib/GFS/shtns.f90
+++ /dev/null
@@ -1,417 +0,0 @@
-      module shtns
-
-! fortran95 interface to the SHTNS library
-! (http://users.isterre.fr/nschaeff/SHTns/index.html).
-! this version hardwired for gaussian grids.
-
-! The available routines are:
-
-! SUBROUTINE GRDTOSPEC(datagrid,dataspec):
-! converts input gridded data array (datagrid) to complex spectral
-! coefficients (dataspec).
-
-! SUBROUTINE SPECTOGRD(dataspec,datagrid):
-! converts input spectral coefficient array (dataspec) to a grid (datagrid).
-
-! SUBROUTINE GETUV(vrtspec,divspec,ugrid,vgrid,rsphere):
-! given input spectral coefficients of vorticity and divergence
-! (vrtspec,divspec) calculates gridded winds (ugrid,vgrid).
-! rsphere is the radius of the sphere.
-
-! SUBROUTINE GETVRTDIVSPEC(ugrid,vgrid,vrtspec,divspec,rsphere):
-! given input gridded winds (ugrid,vgrid) calculates spectral coefficients
-! of vorticity and divergence (vrtspec,divspec). rsphere is the radius of
-! the sphere.
-
-! SUBROUTINE GETGRAD(dataspec,gradx,grady,rsphere):
-! calculates gridded vector gradient (gradx,grady) given input
-! spectral coefficients. rsphere is the radius of the sphere.
-
-! SUBROUTINE SHTNS_CLEANUP():
-! Garbage collection. Deallocates all memory.
-! Call this when you are done with the module to free up memory.
-
-! SUBROUTINE SHTNS_INIT(NLONS,NLATS,NTRUNC):
-! initialize transform machinery for a given grid and spectral
-! truncation.  Called automatically by other routines if
-! grid (as described by input arrays) changes.
-
-! module data:
-! lats,lons: 2-d arrays with lat/lon values (radians) of transform grid.
-! nlm: number of complex spectral coeffs.
-! current_nlon, current_nlat: current values of nlons,nlats
-! gauwts: 1-d array (length current_nlat) of gaussian quadrature weights.
-! current_ntrunc: current spectral truncation
-! degree: array with degree (n) for each spectral coefficient.
-! order: array with order (m) for each spectral coefficient.
-! lap: lapacian operator (-n*(n+1))
-! invlap: inverse of laplacian operator.
-
-! All of these routines use triangular truncation.
-! The spherical harmonic coefficients are normalized using an
-! 'orthonormalized' norm
-! (http://en.wikipedia.org/wiki/Spherical_harmonics#Conventions).
-
-      use kinds, only: r_kind, default_real, r_double
-      use iso_c_binding
-      implicit none
-      private
-      public :: shtns_init,grdtospec,spectogrd,getuv,getvrtdivspec,&
-                getgrad,shtns_destroy
-      public :: gauwts, lats, lons, nlm, degree, order, lap, invlap,&
-                current_nlon, current_nlat, current_ntrunc, areawts,&
-            print_bool
-      INTEGER, PARAMETER :: SHT_NATIVE_LAYOUT=0
-      INTEGER, PARAMETER :: SHT_THETA_CONTIGUOUS=256
-      INTEGER, PARAMETER :: SHT_PHI_CONTIGUOUS=512
-      INTEGER, PARAMETER :: SHT_NO_CS_PHASE=1024
-      INTEGER, PARAMETER :: SHT_REAL_NORM=2048
-      INTEGER, PARAMETER :: SHT_ORTHONORMAL=0
-      INTEGER, PARAMETER :: SHT_FOURPI=1
-      INTEGER, PARAMETER :: SHT_SCHMIDT=2
-      INTEGER, PARAMETER :: SHT_GAUSS=0
-      INTEGER, PARAMETER :: SHT_AUTO=1
-      INTEGER, PARAMETER :: SHT_REG_FAST=2
-      INTEGER, PARAMETER :: SHT_REG_DCT=3
-      INTEGER, PARAMETER :: SHT_QUICK_INIT=4
-      INTEGER, PARAMETER :: SHT_REG_POLES=5
-      INTEGER, PARAMETER :: SHT_GAUSS_FLY=6
-      REAL(r_double), PARAMETER :: SHT_DEFAULT_POLAR_OPT=1.d-10
-      INTEGER, SAVE      :: current_nlon = -1
-      INTEGER, SAVE      :: current_nlat = -1
-      INTEGER, SAVE      :: current_ntrunc = -1
-      INTEGER            :: nlm
-      INTEGER :: print_bool
-! arrays allocated when nlon or nlat or ntrunc change.
-      REAL(r_kind), DIMENSION(:), pointer :: lap, invlap, gauwts
-      REAL(r_kind), DIMENSION(:,:), pointer :: lats
-      REAL(r_kind), DIMENSION(:,:), pointer :: lons, areawts
-      INTEGER, DIMENSION(:), pointer :: degree, order
-      real(r_double) :: popt ! polar optimization thresh
-      integer :: nth ! number of threads to use
-
-      contains
-
-      subroutine sht_assign_arrays(py_lap, py_invlap, py_gauwts,&
-              py_lons, py_lats, py_areawts, py_degree,&
-              py_order, py_nlons, py_nlats, py_nlm)&
-              bind(c, name='gfs_assign_sht_arrays')
-
-          REAL(c_double), DIMENSION(py_nlm), target :: py_lap, py_invlap
-          REAL(c_double), DIMENSION(py_nlats), target :: py_gauwts
-          REAL(c_double), DIMENSION(py_nlons,py_nlats), target :: py_lats
-          REAL(c_double), DIMENSION(py_nlons,py_nlats), target :: py_lons
-          REAL(c_double), DIMENSION(py_nlons,py_nlats), target :: py_areawts
-          INTEGER(c_int), DIMENSION(py_nlm), target :: py_degree, py_order
-          integer(c_int), intent(in) :: py_nlons, py_nlats, py_nlm
-
-          lap => py_lap
-          invlap => py_invlap
-          gauwts => py_gauwts
-          lats => py_lats
-          lons => py_lons
-          areawts => py_areawts
-          degree => py_degree
-          order => py_order
-
-      end subroutine sht_assign_arrays
-
-      !subroutine get_lon_lat(longitudes, latitudes) &
-      !        bind(c,name='gfs_get_lon_lat')
-
-      !real(c_double), intent(out), dimension(current_nlon,current_nlat)&
-      !          :: longitudes, latitudes
-
-      !if(allocated(lats)) then
-      !    latitudes(:,:) = lats(:,:)
-      !else
-      !    print *, 'latitudes not allocated'
-      !endif
-
-      !if(allocated(lons)) then
-      !    longitudes(:,:) = lons(:,:)
-      !else
-      !    print *, 'longitudes not allocated'
-      !endif
-
-      !end subroutine get_lon_lat
-
-      subroutine shtns_init(nlon,nlat,ntrunc,nthreads,polar_opt)
-! initialize library, allocate arrays.
-      integer, intent(in) :: nlon,nlat,ntrunc
-      integer, intent(in), optional :: nthreads
-      real(r_double), intent(in), optional :: polar_opt
-      real(r_double), dimension(:), allocatable :: lats1,gauwts1
-      real(r_double) pi
-      integer m,n,i,j
-      print_bool = 0
-      if (present(nthreads)) then
-         nth = nthreads
-      else
-         nth = 1
-      endif
-      if (present(polar_opt)) then
-         popt = polar_opt
-      else
-         popt = SHT_DEFAULT_POLAR_OPT
-      endif
-      call shtns_destroy()
-      call shtns_use_threads(0) ! JOY auto compute number of openmp threads.
-      !print *, 'In set_size', ntrunc
-      call shtns_set_size(ntrunc,ntrunc,1,SHT_ORTHONORMAL+SHT_NO_CS_PHASE)
-      !print *, 'In precompute'
-      !print *, nlat, nlon
-      call shtns_precompute(SHT_GAUSS_FLY,SHT_PHI_CONTIGUOUS,popt,nlat,nlon)
-      call shtns_calc_nlm(nlm,ntrunc,ntrunc,1)
-      current_nlat = nlat
-      current_nlon = nlon
-      current_ntrunc = ntrunc
-      if (nlm .ne. (ntrunc+1)*(ntrunc+2)/2) then
-         print *,'error: nlm not what expected',nlm,ntrunc
-         stop
-      endif
-      !if (.not. allocated(lats)) allocate(lats(nlon,nlat))
-      !if (.not. allocated(gauwts)) allocate(gauwts(nlat))
-      !if (.not. allocated(lons)) allocate(lons(nlon,nlat))
-      allocate(lats1(nlat))
-      allocate(gauwts1(nlat/2))
-      call shtns_cos_array(lats1)
-      call shtns_gauss_wts(gauwts1)
-      ! lats is 2d-array (lon,lat), lats1 is 1d
-      do i=1,nlon
-         lats(i,:) = asin(lats1)
-      enddo
-      do j=1,nlat/2
-         gauwts(j) = gauwts1(j)
-         gauwts(nlat-j+1) = gauwts1(j)
-      enddo
-      !allocate(areawts(nlon,nlat))
-      do i=1,nlon
-         areawts(i,:) = gauwts(:)
-      enddo
-      areawts = areawts/sum(areawts)
-      pi = 4.*atan(1.0)
-      do j=1,nlat
-      do i=1,nlon
-         lons(i,j) = 2.*pi*real(i-1)/nlon
-      enddo
-      enddo
-      deallocate(lats1,gauwts1)
-      !if (.not. allocated(degree)) allocate(degree(nlm))
-      !if (.not. allocated(order)) allocate(order(nlm))
-      degree = (/((n,n=m,ntrunc),m=0,ntrunc)/)
-      order = (/((m,n=m,ntrunc),m=0,ntrunc)/)
-      !if (.not. allocated(lap)) allocate(lap(nlm))
-      !if (.not. allocated(invlap)) allocate(invlap(nlm))
-      lap = -degree*(degree+1.0)
-      invlap = 0
-      invlap(2:) = 1./lap(2:)
-      end subroutine shtns_init
-
-      subroutine shtns_destroy()
-! deallocate arrays.
-      call shtns_reset()
-      !if (allocated(lats)) deallocate(lats)
-      !if (allocated(gauwts)) deallocate(gauwts)
-      !if (allocated(areawts)) deallocate(areawts)
-      !if (allocated(lons)) deallocate(lons)
-      !if (allocated(degree)) deallocate(degree)
-      !if (allocated(order)) deallocate(order)
-      !if (allocated(lap)) deallocate(lap)
-      !if (allocated(invlap)) deallocate(invlap)
-      !nullify(lats,gauwts,areawts,lons,degree,order,lap,invlap)
-      end subroutine shtns_destroy
-
-      subroutine grdtospec(datagrid,dataspec)
-! converts gridded input array (datagrid) to complex spectral coefficients
-! (dataspec).
-      real(c_double), dimension(:,:), intent(in) :: datagrid
-      complex(r_kind), dimension(:), intent(out) :: dataspec
-      real(r_double), dimension(:,:), allocatable :: datagrid_tmp
-      complex(r_double), dimension(:), allocatable :: dataspec_tmp
-      integer nlon,nlat,ntrunc
-      ! compute array dimensions and infer truncation limit
-      ! from size of input arrays
-      nlon = size(datagrid,1)
-      nlat = size(datagrid,2)
-      ntrunc = nint((-1.+sqrt(1+8*float(size(dataspec))))/2.)-1
-
-      if (nlon .ne. current_nlon .or. &
-          nlat .ne. current_nlat .or. &
-          ntrunc .ne. current_ntrunc) then
-          call shtns_init(nlon,nlat,ntrunc,nth,popt)
-      end if
-      if (default_real .ne. 2) then
-         ! if inputs not double precision convert,
-         ! then convert back.
-         !if (print_bool == 1) then
-         !   print *, datagrid(4,1:3)
-         !   print *, 'grdtospec not default_real: ', loc(datagrid)
-         !   print *, ''
-         !end if
-
-         allocate(datagrid_tmp(nlon,nlat))
-         allocate(dataspec_tmp(nlm))
-         datagrid_tmp = datagrid
-         call shtns_spat_to_sh(datagrid_tmp, dataspec_tmp)
-         dataspec = dataspec_tmp
-         deallocate(datagrid_tmp,dataspec_tmp)
-      else
-         call shtns_spat_to_sh(datagrid, dataspec)
-      endif
-      end subroutine grdtospec
-
-      subroutine spectogrd(dataspec,datagrid)
-! converts complex spectral coefficients (dataspec) to
-! gridded data array (datagrid).
-      real(r_kind), dimension(:,:), intent(out) :: datagrid
-      complex(r_kind), dimension(:), intent(in) :: dataspec
-      real(r_double), dimension(:,:), allocatable :: datagrid_tmp
-      complex(r_double), dimension(:), allocatable :: dataspec_tmp
-      integer nlon,nlat,ntrunc
-      ! compute array dimensions and infer truncation limit
-      ! from size of input arrays
-      nlon = size(datagrid,1)
-      nlat = size(datagrid,2)
-      ntrunc = nint((-1.+sqrt(1+8*float(size(dataspec))))/2.)-1
-      if (nlon .ne. current_nlon .or. &
-          nlat .ne. current_nlat .or. &
-          ntrunc .ne. current_ntrunc) then
-          call shtns_init(nlon,nlat,ntrunc,nth,popt)
-      end if
-      if (default_real .ne. 2) then
-         allocate(datagrid_tmp(nlon,nlat))
-         allocate(dataspec_tmp(nlm))
-         dataspec_tmp = dataspec
-         call shtns_sh_to_spat(dataspec_tmp, datagrid_tmp)
-         datagrid = datagrid_tmp
-         deallocate(datagrid_tmp,dataspec_tmp)
-      else
-         call shtns_sh_to_spat(dataspec, datagrid)
-      endif
-      end subroutine spectogrd
-
-      subroutine getuv(vrtspec,divspec,ugrid,vgrid,rsphere)
-! given spectral coefficients of vorticity and divergence
-! (vrtspec,divspec) compute gridded winds (ugrid,vgrid).
-! rsphere is the radius of the sphere in meters.
-      real(r_kind), intent(in) :: rsphere
-      real(r_kind), dimension(:,:), intent(out) :: ugrid,vgrid
-      complex(r_kind), dimension(:), intent(in) :: vrtspec,divspec
-      real(r_double), dimension(:,:), allocatable :: ugrid_tmp,vgrid_tmp
-      complex(r_double), dimension(:), allocatable :: &
-      vrtspec_tmp,divspec_tmp
-      integer nlon,nlat,ntrunc
-      ! compute array dimensions and infer truncation limit
-      ! from size of input arrays
-      nlon = size(ugrid,1)
-      nlat = size(ugrid,2)
-      ntrunc = nint((-1.+sqrt(1+8*float(size(vrtspec))))/2.)-1
-      if (nlon .ne. current_nlon .or. &
-          nlat .ne. current_nlat .or. &
-          ntrunc .ne. current_ntrunc) then
-          call shtns_init(nlon,nlat,ntrunc,nth,popt)
-      end if
-      if (default_real .ne. 2) then
-         allocate(ugrid_tmp(nlon,nlat))
-         allocate(vrtspec_tmp(nlm))
-         allocate(vgrid_tmp(nlon,nlat))
-         allocate(divspec_tmp(nlm))
-         vrtspec_tmp = invlap*rsphere*vrtspec
-         divspec_tmp = invlap*rsphere*divspec
-         call shtns_sphtor_to_spat(vrtspec_tmp,divspec_tmp,ugrid_tmp,vgrid_tmp)
-         ugrid = ugrid_tmp
-         vgrid = vgrid_tmp
-         deallocate(ugrid_tmp,vgrid_tmp,vrtspec_tmp,divspec_tmp)
-      else
-         allocate(vrtspec_tmp(nlm))
-         allocate(divspec_tmp(nlm))
-         vrtspec_tmp = invlap*rsphere*vrtspec
-         divspec_tmp = invlap*rsphere*divspec
-         call shtns_sphtor_to_spat(vrtspec_tmp,divspec_tmp,ugrid,vgrid)
-         deallocate(vrtspec_tmp,divspec_tmp)
-      endif
-      end subroutine getuv
-
-      subroutine getvrtdivspec(ugrid,vgrid,vrtspec,divspec,rsphere)
-! calculate spectral coefficients of vorticity and divergence
-! (vrtspec,divspec) given input gridded winds (ugrid,vgrid).
-! rsphere is the radius of the sphere in meters.
-      real(r_kind), intent(in) :: rsphere
-      real(r_kind), dimension(:,:), intent(in) :: ugrid,vgrid
-      complex(r_kind), dimension(:), intent(out) :: vrtspec,divspec
-      real(r_double), dimension(:,:), allocatable :: ugrid_tmp,vgrid_tmp
-      complex(r_double), dimension(:), allocatable :: &
-      vrtspec_tmp,divspec_tmp
-      integer nlon,nlat,ntrunc
-      ! compute array dimensions and infer truncation limit
-      ! from size of input arrays.
-      nlon = size(ugrid,1)
-      nlat = size(ugrid,2)
-      ntrunc = nint((-1.+sqrt(1+8*float(size(vrtspec))))/2.)-1
-      if (nlon .ne. current_nlon .or. &
-          nlat .ne. current_nlat .or. &
-          ntrunc .ne. current_ntrunc) then
-          call shtns_init(nlon,nlat,ntrunc,nth,popt)
-      end if
-      if (default_real .ne. 2) then
-         allocate(ugrid_tmp(nlon,nlat))
-         allocate(vrtspec_tmp(nlm))
-         allocate(vgrid_tmp(nlon,nlat))
-         allocate(divspec_tmp(nlm))
-         ugrid_tmp = ugrid
-         vgrid_tmp = vgrid
-         call shtns_spat_to_sphtor(ugrid_tmp, vgrid_tmp, vrtspec_tmp, divspec_tmp)
-         vrtspec = vrtspec_tmp
-         divspec = divspec_tmp
-         deallocate(ugrid_tmp,vgrid_tmp,vrtspec_tmp,divspec_tmp)
-      else
-         call shtns_spat_to_sphtor(ugrid, vgrid, vrtspec, divspec)
-      endif
-      vrtspec = (lap/rsphere)*vrtspec; divspec = (lap/rsphere)*divspec
-      end subroutine getvrtdivspec
-
-      subroutine getgrad(divspec,ugrid,vgrid,rsphere)
-! given spectral coefficients compute gridded vector gradient.
-! rsphere is the radius of the sphere in meters.
-      real(r_kind), intent(in) :: rsphere
-      real(r_kind), dimension(:,:), intent(out) :: ugrid,vgrid
-      complex(r_kind), dimension(:), intent(in) :: divspec
-      complex(r_kind), dimension(:), allocatable :: vrtspec
-      real(r_double), dimension(:,:), allocatable :: ugrid_tmp,vgrid_tmp
-      complex(r_double), dimension(:), allocatable :: &
-      vrtspec_tmp,divspec_tmp
-      integer nlon,nlat,ntrunc
-      ! compute array dimensions and infer truncation limit
-      ! from size of input arrays
-      nlon = size(ugrid,1)
-      nlat = size(ugrid,2)
-      ntrunc = nint((-1.+sqrt(1+8*float(size(divspec))))/2.)-1
-      if (nlon .ne. current_nlon .or. &
-          nlat .ne. current_nlat .or. &
-          ntrunc .ne. current_ntrunc) then
-          call shtns_init(nlon,nlat,ntrunc,nth,popt)
-      end if
-      if (default_real .ne. 2) then
-         allocate(ugrid_tmp(nlon,nlat))
-         allocate(vrtspec_tmp(nlm))
-         allocate(vgrid_tmp(nlon,nlat))
-         allocate(divspec_tmp(nlm))
-         vrtspec_tmp = 0.
-         divspec_tmp = divspec
-         call shtns_sphtor_to_spat(vrtspec_tmp,divspec_tmp,ugrid_tmp,vgrid_tmp)
-         ugrid = ugrid_tmp
-         vgrid = vgrid_tmp
-         deallocate(ugrid_tmp,vgrid_tmp,vrtspec_tmp,divspec_tmp)
-      else
-         allocate(vrtspec(nlm))
-         vrtspec = 0.
-         call shtns_sphtor_to_spat(vrtspec,divspec,ugrid,vgrid)
-         deallocate(vrtspec)
-      endif
-      ugrid = ugrid/rsphere; vgrid = vgrid/rsphere
-      end subroutine getgrad
-
-! that's it!
-
-      end module shtns
diff --git a/climt/_lib/GFS/spectral_data.f90 b/climt/_lib/GFS/spectral_data.f90
deleted file mode 100644
index 4a5a4f06..00000000
--- a/climt/_lib/GFS/spectral_data.f90
+++ /dev/null
@@ -1,71 +0,0 @@
-module spectral_data
-! spectral model state variable data
-! Public subroutines:
-! init_specdata: allocate arrays.
-! destroy_specdata: deallocate arrays.
- use kinds, only: r_kind
- use params, only: ndimspec, nlevs, ntrac
- use iso_c_binding
- implicit none
- private
- public :: init_specdata,destroy_specdata
- complex(r_kind), public, dimension(:,:), pointer :: &
- vrtspec,divspec,virtempspec
- complex(r_kind), public, dimension(:,:,:), pointer :: tracerspec
-! ndimspec by nlevs complex arrays (ndimspec = (ntrunc+1)*(ntrunc)/2):
-! these arrays are update each timestep in subroutine advance in run.f90.
-! vrtspec: vorticity
-! divspec: divergence
-! virtempspec: virtural temp
-! tracerspec: tracers (ntrac=1 is  specific humidity, ntrac=ntoz ozone, ntrac=ntclw cloud
-! condensate)
- complex(r_kind), public, dimension(:), pointer :: lnpsspec,topospec
- real(r_kind), public, dimension(:), pointer :: disspec,diff_prof,dmp_prof
-! ndimspec 1-d arrays:
-! lnpsspec: log(ps) in Pa
-! topospec: orography in meters (static, created in init_dyn)
-! disspec: hyperdiffusion operator (static, created in init_dyn)
-! diff_prof: vertical profile of hyperdiffusion coefficient.
-! dmp_prof: vertical profile of linear momentum drag (active
-! at top of model).
-
- contains
-!
-! Modifying the specdata init routines to accept arrays from python   
- subroutine init_specdata(pyVrtSpec, pyDivSpec, pyVirtTempSpec, pyTracerSpec, &
-         pyTopoSpec, pyLnPsSpec, pyDissSpec, pyDmpProf, pyDiffProf) & 
-            bind(c, name='gfs_initialise_spectral_arrays')
-    
-    complex(c_double_complex), target, intent(inout) :: pyVrtSpec(ndimspec,nlevs), &
-    pyDivSpec(ndimspec, nlevs), pyVirtTempSpec(ndimspec, nlevs), &
-    pyTracerSpec(ndimspec, nlevs, ntrac), pyTopoSpec(ndimspec), &
-    pyLnPsSpec(ndimspec)
-
-    real(c_double), target, intent(inout) :: pyDmpProf(nlevs), &
-            pyDiffProf(nlevs), pyDissSpec(ndimspec)
-
-    !print *, "initialising Spectral arrays"
-    !print *, shape(pyVrtSpec), shape(pyTracerSpec)
-
-    vrtspec => pyVrtSpec
-    divspec => pyDivSpec
-    tracerspec => pyTracerSpec
-    virtempspec => pyVirtTempSpec
-
-    topospec => pyTopoSpec
-    lnpsspec => pyLnPsSpec
-
-    disspec => pyDissSpec
-    dmp_prof => pyDmpProf
-    diff_prof => pyDiffProf
-
- end subroutine init_specdata
- subroutine destroy_specdata()
-    !print *, 'Deallocating spectral data'
-    nullify(vrtspec,divspec,virtempspec)
-    nullify(tracerspec)
-    nullify(topospec,lnpsspec,disspec,dmp_prof,diff_prof)
- end subroutine destroy_specdata
-
-
-end module spectral_data
diff --git a/climt/_lib/Makefile b/climt/_lib/Makefile
index a2ccfd29..efa0f02f 100644
--- a/climt/_lib/Makefile
+++ b/climt/_lib/Makefile
@@ -1,7 +1,6 @@
 FFTW_DIR = fftw-3.3.6-pl2
 SHT_DIR = shtns-2.7-r557
 BLAS_DIR = OpenBLAS-0.2.19
-GFS_DIR = GFS
 SIM_PHYS_DIR = simple_physics
 EMAN_DIR = emanuel
 RRTMG_LW_DIR = rrtmg_lw
@@ -20,55 +19,15 @@ LIB_DIR= $(BASE_DIR)/$(CLIMT_ARCH)/
 CLIMT_CFLAGS =  "-I$(BASE_DIR) -L$(LIB_DIR) $(CFLAGS)"
 LDFLAGS = "-L$(LIB_DIR)"
 
-CLIMT_LIB_TARGETS = blas_lib fftw_lib
 
 ifeq ($(CLIMT_ARCH),Windows)
 	BASE_DIR = $(subst \,/,$(PWD))
 	LIB_DIR= $(BASE_DIR)/$(CLIMT_ARCH)/
 	CLIMT_FFTW_FLAGS = --with-our-malloc --build=x86_64-w64-mingw32 --host=x86_64-w64-mingw32 --enable-static --with-g77-wrappers --with-windows-f77-mangling --enable-sse2 --enable-avx
-	CLIMT_LIB_TARGETS =
 endif
 
 
-all: $(CLIMT_LIB_TARGETS) sht_lib gfs_lib simple_physics_lib emanuel_lib rrtmg_lw_lib rrtmg_sw_lib libdcmip
-
-# FFTW3 Configuration
-
-$(CLIMT_ARCH)/.configured_fftw:
-	if [ ! -d $(FFTW_DIR) ]; then tar -xvzf $(FFTW_GZ) > log 2>&1; fi; cd $(FFTW_DIR); ./configure CFLAGS=$(CLIMT_CFLAGS) LDFLAGS=$(LDFLAGS) --prefix="$(BASE_DIR)" $(CLIMT_FFTW_FLAGS) > log.config.fftw 2>&1; tail -n 100 log.config.fftw && touch ../$(CLIMT_ARCH)/.configured_fftw
-
-$(CLIMT_ARCH)/libfftw3_omp.a: $(CLIMT_ARCH)/.configured_fftw
-	cd $(FFTW_DIR); make install > log.install.fftw 2>&1; tail -n 100 log.install.fftw; cp ../lib/libfftw3_omp.a $(LIB_DIR)
-
-$(CLIMT_ARCH)/libfftw3.a: $(CLIMT_ARCH)/.configured_fftw
-	cp lib/libfftw3.a $(LIB_DIR)
-
-fftw_lib: $(CLIMT_ARCH)/libfftw3_omp.a $(CLIMT_ARCH)/libfftw3.a
-
-
-# SHTNS Configuration
-
-$(CLIMT_ARCH)/.configured_sht:
-	if [ ! -d $(SHT_DIR) ]; then tar -xvzf $(SHT_GZ) ; fi; cd $(SHT_DIR); ./configure CFLAGS=$(CLIMT_CFLAGS) LDFLAGS=$(LDFLAGS) --prefix=$(BASE_DIR) --disable-python --enable-openmp && touch ../$(CLIMT_ARCH)/.configured_sht
-
-$(CLIMT_ARCH)/libshtns_omp.a: $(CLIMT_ARCH)/.configured_sht
-	cd $(SHT_DIR); make install > log 2>&1; cat log; cp ../lib/libshtns_omp.a $(LIB_DIR); ls -ltrh $(LIB_DIR)
-
-sht_lib: $(CLIMT_ARCH)/libshtns_omp.a
-
-# BLAS Configuration
-
-blas_lib: $(CLIMT_ARCH)/libopenblas.a
-
-$(CLIMT_ARCH)/libopenblas.a:
-	if [ ! -d $(BLAS_DIR) ]; then tar -xvzf $(BLAS_GZ) > log 2>&1; fi; cd $(BLAS_DIR); CFLAGS=$(CLIMT_CFLAGS) make NO_SHARED=1 NO_LAPACK=0 > log 2>&1 ; make PREFIX=$(BASE_DIR) NO_SHARED=1 NO_LAPACK=0 install > log.install.blas 2>&1 ; tail -n 100 log.install.blas;  cp ../lib/libopenblas.a $(LIB_DIR)
-
-# GFS Configuration
-
-gfs_lib: $(CLIMT_ARCH)/libgfs_dycore.a
-
-$(CLIMT_ARCH)/libgfs_dycore.a: $(GFS_DIR)/*.f90
-	cd $(GFS_DIR); make python_lib; cp libgfs_dycore.a $(LIB_DIR); touch ../../_components/gfs/_gfs_dynamics.pyx
+all: simple_physics_lib emanuel_lib rrtmg_lw_lib rrtmg_sw_lib libdcmip
 
 # Simple Physics Configuration
 
@@ -111,7 +70,6 @@ clean:
 	if [ -d $(FFTW_DIR) ]; then rm -Rf $(FFTW_DIR); fi
 	if [ -d $(SHT_DIR) ]; then rm -Rf $(SHT_DIR); fi
 	if [ -d $(BLAS_DIR) ]; then rm -Rf $(BLAS_DIR); fi
-	cd $(GFS_DIR); make clean
 	cd $(SIM_PHYS_DIR); make clean
 	cd $(EMAN_DIR); make clean
 	cd $(RRTMG_LW_DIR); make clean
diff --git a/climt/_lib/OpenBLAS-0.2.19.tar.gz b/climt/_lib/OpenBLAS-0.2.19.tar.gz
deleted file mode 100644
index 6f7a6013..00000000
Binary files a/climt/_lib/OpenBLAS-0.2.19.tar.gz and /dev/null differ
diff --git a/climt/_lib/Windows/libfftw3.a b/climt/_lib/Windows/libfftw3.a
deleted file mode 100644
index 1a3ca569..00000000
Binary files a/climt/_lib/Windows/libfftw3.a and /dev/null differ
diff --git a/climt/_lib/Windows/libfftw3_omp.a b/climt/_lib/Windows/libfftw3_omp.a
deleted file mode 100644
index 77154c92..00000000
Binary files a/climt/_lib/Windows/libfftw3_omp.a and /dev/null differ
diff --git a/climt/_lib/fftw-3.3.6-pl2.tar.gz b/climt/_lib/fftw-3.3.6-pl2.tar.gz
deleted file mode 100644
index 1ec2d183..00000000
Binary files a/climt/_lib/fftw-3.3.6-pl2.tar.gz and /dev/null differ
diff --git a/climt/_lib/shtns-2.7-r557.tar.gz b/climt/_lib/shtns-2.7-r557.tar.gz
deleted file mode 100644
index 88ea07f3..00000000
Binary files a/climt/_lib/shtns-2.7-r557.tar.gz and /dev/null differ
diff --git a/setup.py b/setup.py
index 2f343df7..e8ffb034 100644
--- a/setup.py
+++ b/setup.py
@@ -84,7 +84,7 @@ def guess_compiler_name(env_name):
 
 if operating_system == 'Linux':
     libraries = ['m', 'gfortran', 'rt']
-    default_link_args = ['-lgfortran']
+    default_link_args = ['-lgfortran', '-lm']
 
 if operating_system == 'Windows':
     compiled_base_dir = 'climt\\_lib'
@@ -94,7 +94,6 @@ def guess_compiler_name(env_name):
 lib_path_list = [os.path.join(compiled_path, operating_system+'/')]
 inc_path = os.path.join(compiled_path, 'include')
 include_dirs.append(inc_path)
-openblas_path = lib_path_list[0]+'/libopenblas.a'
 
 
 # Compile libraries
@@ -102,8 +101,8 @@ def guess_compiler_name(env_name):
 if 'FC' not in os.environ:
     if operating_system == 'Darwin':
         # guess_compiler_name('FC')
-        os.environ['FC'] = 'gfortran-6'
-        os.environ['F77'] = 'gfortran-6'
+        os.environ['FC'] = 'gfortran-10'
+        os.environ['F77'] = 'gfortran-10'
     else:
         os.environ['FC'] = 'gfortran'
         os.environ['F77'] = 'gfortran'
@@ -111,19 +110,18 @@ def guess_compiler_name(env_name):
 if 'CC' not in os.environ:
     if operating_system == 'Darwin':
         # guess_compiler_name('CC')
-        os.environ['CC'] = 'gcc-6'
+        os.environ['CC'] = 'gcc-10'
     else:
         os.environ['CC'] = 'gcc'
 
 if 'CLIMT_OPT_FLAGS' not in os.environ:
     os.environ['CLIMT_OPT_FLAGS'] = '-O3'
 
-if operating_system == 'Windows' :
+if operating_system == 'Windows':
     os.environ['CC'] = 'gcc.exe'
     os.environ['FC'] = 'gfortran.exe'
     os.environ['AR'] = 'gcc-ar.exe'
     libraries = []
-    openblas_path = os.path.join(os.environ['COMPILER_PATH'], '../lib/libopenblas.a')
     default_link_args = ['-l:libgfortran.a', '-l:libquadmath.a', '-l:libm.a']
     default_compile_args = ['-DMS_WIN64']
 
@@ -132,12 +130,15 @@ def guess_compiler_name(env_name):
 
 if operating_system == 'Darwin':
     gcc_dir = find_homebrew_gcc()
+    # print('gcc_dir', gcc_dir)
     for root, dirs, files in os.walk(gcc_dir):
         for line in files:
             if re.match('libgfortran.a', line):
                 if not ('i386' in root):
                     lib_path_list.append(root)
 
+    # print(lib_path_list)
+
     os.environ['FFLAGS'] += ' -mmacosx-version-min=10.7'
     os.environ['CFLAGS'] += ' -mmacosx-version-min=10.7'
     default_link_args = []
@@ -214,19 +215,6 @@ def run(self):
             library_dirs=lib_path_list,
             extra_link_args=[lib_path_list[0]+'/librrtmg_lw.a', '-fopenmp'] + default_link_args),
 
-        Extension(
-            'climt._components.gfs._gfs_dynamics',
-            sources=['climt/_components/gfs/_gfs_dynamics.pyx'],
-            libraries=libraries,
-            include_dirs=include_dirs,
-            extra_compile_args=['-fopenmp'] + default_compile_args,
-            library_dirs=lib_path_list,
-            extra_link_args=['-fopenmp', lib_path_list[0]+'/libgfs_dycore.a',
-                             lib_path_list[0]+'/libshtns_omp.a', lib_path_list[0]+'/libfftw3_omp.a',
-                             lib_path_list[0]+'/libfftw3.a', openblas_path] + default_link_args),
-
-        # lib_path+'/libshtns_omp.a', openblas_path, os.environ['COMPILER_PATH']+'../lib/libfftw3.a'] + default_link_args),
-
         Extension(
             'climt._components.rrtmg.sw._rrtmg_sw',
             sources=['climt/_components/rrtmg/sw/_rrtmg_sw.pyx'],
diff --git a/tox.ini b/tox.ini
index 56c2cf79..9170d76e 100644
--- a/tox.ini
+++ b/tox.ini
@@ -13,8 +13,8 @@ commands=flake8 climt
 ignore = E501,E226,W504,W605
 
 [testenv]
-passenv = CC FC
-    PYTHONPATH = {toxinidir}:{toxinidir}/climt
+# passenv = CC FC
+PYTHONPATH = {toxinidir}:{toxinidir}/climt
 deps =
     -r{toxinidir}/requirements_dev.txt
 commands =


Test patch:
diff --git a/tests/test_bucket_hydrology.py b/tests/test_bucket_hydrology.py
index 507ac1a9..13e4d16a 100644
--- a/tests/test_bucket_hydrology.py
+++ b/tests/test_bucket_hydrology.py
@@ -2,13 +2,15 @@
 import numpy as np
 from climt import (BucketHydrology, get_default_state)
 
+
 def get_quantities(state):
-    heat = state['surface_material_density'].values*state['soil_layer_thickness'].values*\
-           state['surface_temperature'].values*state['heat_capacity_of_soil'].values
+    heat = state['surface_material_density'].values * state['soil_layer_thickness'].values * \
+        state['surface_temperature'].values*state['heat_capacity_of_soil'].values
     moisture = state['lwe_thickness_of_soil_moisture_content'].values
 
     return heat, moisture
 
+
 time_step = timedelta(seconds=1)
 state = get_default_state([BucketHydrology()])
 
@@ -69,13 +71,14 @@ def get_quantities(state):
 new_heat, new_moisture = get_quantities(state)
 
 
-#Test starts
+# Test starts
 
 def test_bucket_hydrology_heat_conservation():
     forcing_amount = surf_forcing * time_step.total_seconds()
     assert np.isclose(new_heat - old_heat, forcing_amount,
                       rtol=0, atol=1e-3)
 
+
 def test_bucket_hydrology_moisture_conservation():
     forcing_amount = mois_forcing * time_step.total_seconds()
     assert np.isclose(new_moisture - old_moisture, forcing_amount,
diff --git a/tests/test_components.py b/tests/test_components.py
index e733e9bc..bcb9cd30 100644
--- a/tests/test_components.py
+++ b/tests/test_components.py
@@ -11,7 +11,7 @@
     Frierson06LongwaveOpticalDepth, GridScaleCondensation,
     BergerSolarInsolation, SimplePhysics, RRTMGLongwave,
     RRTMGShortwave, SlabSurface, EmanuelConvection,
-    DcmipInitialConditions, GFSDynamicalCore, BucketHydrology,
+    DcmipInitialConditions, BucketHydrology,
     IceSheet, Instellation, DryConvectiveAdjustment,
     get_grid)
 import climt
@@ -46,7 +46,7 @@ def load_dictionary(filename):
 def state_3d_to_1d(state):
     return_state = {}
     for name, value in state.items():
-        if name is 'time':
+        if name == 'time':
             return_state[name] = value
         else:
             dim_list = []
@@ -62,7 +62,7 @@ def state_3d_to_1d(state):
 def transpose_state(state, dims=None):
     return_state = {}
     for name, value in state.items():
-        if name is 'time':
+        if name == 'time':
             return_state[name] = state[name]
         else:
             if dims is None:
@@ -517,104 +517,6 @@ class TestDryConvection(ComponentBaseColumn, ComponentBase3D):
     def get_component_instance(self):
         return DryConvectiveAdjustment()
 
-
-@pytest.mark.skip("fails on CI, no idea why")
-class TestFullMoistGFSDycoreWithPhysics(ComponentBase3D):
-
-    def get_component_instance(self):
-        # Create Radiation Prognostic
-        radiation = climt.RRTMGLongwave()
-        # Create Convection Prognostic
-        convection = climt.EmanuelConvection()
-        # Create a SimplePhysics Prognostic
-        boundary_layer = TimeDifferencingWrapper(
-            climt.SimplePhysics()
-        )
-        return GFSDynamicalCore(
-            [radiation, convection, boundary_layer]
-        )
-
-
-@pytest.mark.skip("fails on CI, no idea why")
-class TestGFSDycore(ComponentBase3D):
-
-    def get_component_instance(self):
-        return GFSDynamicalCore()
-
-
-@pytest.mark.skip("fails on CI, no idea why")
-class TestGFSDycoreWithDcmipInitialConditions(ComponentBase3D):
-
-    def get_component_instance(self):
-        return GFSDynamicalCore()
-
-    def get_3d_input_state(self):
-        state = climt.get_default_state(
-            [self.get_component_instance()], grid_state=get_grid(nx=32, ny=32, nz=28))
-        # state = super(TestGFSDycoreWithDcmipInitialConditions, self).get_3d_input_state()
-        state.update(climt.DcmipInitialConditions(add_perturbation=True)(state))
-        return state
-
-
-@pytest.mark.skipif(sys.platform == 'win32',
-                    reason="fails on appveyor, no idea why")
-class TestGFSDycoreWithImplicitTendency(ComponentBase3D):
-
-    def get_component_instance(self):
-        return GFSDynamicalCore([EmanuelConvection()])
-
-    def get_3d_input_state(self):
-        state = climt.get_default_state(
-            [self.get_component_instance()], grid_state=get_grid(nx=16, ny=16, nz=28))
-        return state
-
-
-@pytest.mark.skipif(sys.platform == 'win32',
-                    reason="fails on appveyor, no idea why")
-class TestGFSDycoreWithHeldSuarez(ComponentBase3D):
-    def test_inputs_are_dry(self):
-        component = self.get_component_instance()
-        assert 'specific_humidity' not in component.input_properties.keys()
-        assert 'specific_humidity_on_interface_levels' not in component.input_properties.keys()
-
-    def get_component_instance(self):
-        return GFSDynamicalCore([HeldSuarez()])
-
-    def get_3d_input_state(self):
-        state = climt.get_default_state(
-            [self.get_component_instance()],
-            grid_state=get_grid(nx=16, ny=16, nz=28))
-        return state
-
-
-@pytest.mark.skipif(sys.platform == 'win32',
-                    reason="fails on appveyor, no idea why")
-class TestGFSDycoreWithGrayLongwaveRadiation(ComponentBase3D):
-
-    def get_component_instance(self):
-        return GFSDynamicalCore([GrayLongwaveRadiation()])
-
-    def get_3d_input_state(self):
-        state = climt.get_default_state(
-            [self.get_component_instance()],
-            grid_state=get_grid(nx=16, ny=16, nz=28))
-        return state
-
-
-@pytest.mark.skipif(sys.platform == 'win32',
-                    reason="fails on appveyor, no idea why")
-class TestGFSDycoreWithRRTMGLongwave(ComponentBase3D):
-
-    def get_component_instance(self):
-        radiation = RRTMGLongwave()
-        return GFSDynamicalCore([radiation], moist=True)
-
-    def get_3d_input_state(self):
-        state = climt.get_default_state(
-            [self.get_component_instance()], grid_state=get_grid(nx=16, ny=16, nz=28))
-        return state
-
-
 def test_piecewise_constant_component():
 
     radiation = UpdateFrequencyWrapper(RRTMGLongwave(), timedelta(seconds=1000))
diff --git a/tests/test_initialization.py b/tests/test_initialization.py
index 1bba7544..02b7ba2d 100644
--- a/tests/test_initialization.py
+++ b/tests/test_initialization.py
@@ -1,7 +1,7 @@
 from climt import (
     get_default_state, Frierson06LongwaveOpticalDepth, GrayLongwaveRadiation, HeldSuarez,
     GridScaleCondensation, BergerSolarInsolation, SimplePhysics, RRTMGLongwave, RRTMGShortwave,
-    EmanuelConvection, SlabSurface, GFSDynamicalCore, DcmipInitialConditions, IceSheet,
+    EmanuelConvection, SlabSurface, DcmipInitialConditions, IceSheet,
     Instellation, get_grid
 )
 import random
@@ -205,12 +205,6 @@ class ComponentQuantityInitializationTests(unittest.TestCase):
         func = create_3d_grid_test_for(cls)
         locals()[func.__name__] = func
 
-    def test_GFSDynamicalCore(self):
-        grid = get_grid(nx=12, ny=16, nz=28)
-        component = GFSDynamicalCore()
-        state = get_default_state([component], grid_state=grid)
-        call_component(component, state)
-
     def test_component_pairs(self):
         random.seed(0)
         for _ in range(self.pair_tests):
@@ -236,40 +230,6 @@ def test_component_triplets(self):
             call_component(component3, state)
 
 
-class TestFullMoistGFSDycoreWithPhysics(unittest.TestCase):
-
-    def get_component_instance(self):
-        # Create Radiation Prognostic
-        radiation = RRTMGLongwave()
-        # Create Convection Prognostic
-        convection = EmanuelConvection()
-        # Create a SimplePhysics Prognostic
-        boundary_layer = TimeDifferencingWrapper(
-            SimplePhysics()
-        )
-        return GFSDynamicalCore(
-            [radiation, convection, boundary_layer]
-        )
-
-    def test_component_3d_grid(self):
-        grid = get_grid(nx=16, ny=16, nz=16)
-        component = self.get_component_instance()
-        state = get_default_state([component], grid_state=grid)
-        call_component(component, state)
-
-
-class TestGFSDycoreWith32VerticalLevels(unittest.TestCase):
-
-    def get_component_instance(self):
-        return GFSDynamicalCore()
-
-    def test_component_3d_grid(self):
-        grid = get_grid(nx=16, ny=16, nz=32)
-        component = self.get_component_instance()
-        state = get_default_state([component], grid_state=grid)
-        call_component(component, state)
-
-
 def test_3d_initialization_is_full_based_on_wildcard():
     grid_state = get_grid(nx=10, ny=10, nz=20)
     rrtmg_shortwave = RRTMGShortwave()