LCOV - code coverage report
Current view: top level - physics/clubb/src/CLUBB_core - grid_class.F90 (source / functions) Hit Total Coverage
Test: coverage.info Lines: 115 349 33.0 %
Date: 2024-12-17 17:57:11 Functions: 10 23 43.5 %

          Line data    Source code
       1             : !------------------------------------------------------------------------
       2             : ! $Id$
       3             : !=============================================================================== 
       4             : module grid_class
       5             : 
       6             :   ! Description:
       7             :   !
       8             :   ! Definition of a grid class and associated functions
       9             :   !
      10             :   ! The grid specification is as follows:
      11             :   !
      12             :   !     +                ================== zm(nz) =========GP=========
      13             :   !     |
      14             :   !     |
      15             :   ! 1/dzt(nz) +          ------------------ zt(nz) ---------GP---------
      16             :   !     |        |
      17             :   !     |        |
      18             :   !     + 1/dzm(nz-1)    ================== zm(nz-1) ==================
      19             :   !              |
      20             :   !              |
      21             :   !              +       ------------------ zt(nz-1) ------------------
      22             :   !
      23             :   !                                           .
      24             :   !                                           .
      25             :   !                                           .
      26             :   !                                           .
      27             :   !
      28             :   !                      ================== zm(k+1) ===================
      29             :   !
      30             :   !
      31             :   !              +       ------------------ zt(k+1) -------------------
      32             :   !              |
      33             :   !              |
      34             :   !     +    1/dzm(k)    ================== zm(k) =====================
      35             :   !     |        |
      36             :   !     |        |
      37             :   ! 1/dzt(k)     +       ------------------ zt(k) ---------------------
      38             :   !     |
      39             :   !     |
      40             :   !     +                ================== zm(k-1) ===================
      41             :   !
      42             :   !
      43             :   !                      ------------------ zt(k-1) -------------------
      44             :   !
      45             :   !                                           .
      46             :   !                                           .
      47             :   !                                           .
      48             :   !                                           .
      49             :   !
      50             :   !     +                ================== zm(2) =====================
      51             :   !     |
      52             :   !     |
      53             :   ! 1/dzt(2)     +       ------------------ zt(2) ---------------------
      54             :   !     |        |
      55             :   !     |        |
      56             :   !     +    1/dzm(1)    ================== zm(1) ============GP=======  zm_init
      57             :   !              |       //////////////////////////////////////////////  surface
      58             :   !              |
      59             :   !              +       ------------------ zt(1) ------------GP-------
      60             :   !
      61             :   !
      62             :   ! The variable zm(k) stands for the momentum level altitude at momentum
      63             :   ! level k; the variable zt(k) stands for the thermodynamic level altitude at
      64             :   ! thermodynamic level k; the variable invrs_dzt(k) is the inverse distance
      65             :   ! between momentum levels (over a central thermodynamic level k); and the
      66             :   ! variable invrs_dzm(k) is the inverse distance between thermodynamic levels
      67             :   ! (over a central momentum level k).  Please note that in the above diagram,
      68             :   ! "invrs_dzt" is denoted "dzt", and "invrs_dzm" is denoted "dzm", such that
      69             :   ! 1/dzt is the distance between successive momentum levels k-1 and k (over a
      70             :   ! central thermodynamic level k), and 1/dzm is the distance between successive
      71             :   ! thermodynamic levels k and k+1 (over a central momentum level k).
      72             :   !
      73             :   ! The grid setup is compatible with a stretched (unevely-spaced) grid.  Thus,
      74             :   ! the distance between successive grid levels may not always be constant.
      75             :   !
      76             :   ! The following diagram is an example of a stretched grid that is defined on
      77             :   ! momentum levels.  The thermodynamic levels are placed exactly halfway
      78             :   ! between the momentum levels.  However, the momentum levels do not fall
      79             :   ! halfway between the thermodynamic levels.
      80             :   !
      81             :   !        =============== zm(k+1) ===============
      82             :   !
      83             :   !
      84             :   !
      85             :   !        --------------- zt(k+1) ---------------
      86             :   !
      87             :   !
      88             :   !
      89             :   !        ===============  zm(k)  ===============
      90             :   !
      91             :   !        ---------------  zt(k)  ---------------
      92             :   !
      93             :   !        =============== zm(k-1) ===============
      94             :   !
      95             :   ! The following diagram is an example of a stretched grid that is defined on
      96             :   ! thermodynamic levels.  The momentum levels are placed exactly halfway
      97             :   ! between the thermodynamic levels.  However, the thermodynamic levels do not
      98             :   ! fall halfway between the momentum levels.
      99             :   !
     100             :   !        --------------- zt(k+1) ---------------
     101             :   !
     102             :   !
     103             :   !
     104             :   !        ===============  zm(k)  ===============
     105             :   !
     106             :   !
     107             :   !
     108             :   !        ---------------  zt(k)  ---------------
     109             :   !
     110             :   !        =============== zm(k-1) ===============
     111             :   !
     112             :   !        --------------- zt(k-1) ---------------
     113             :   !
     114             :   ! NOTE:  Any future code written for use in the CLUBB parameterization should
     115             :   !        use interpolation formulas consistent with a stretched grid.  The
     116             :   !        simplest way to do so is to call the appropriate interpolation
     117             :   !        function from this module.  Interpolations should *not* be handled in
     118             :   !        the form of:  ( var_zm(k) + var_zm(k-1) ) / 2; *nor* in the form of:
     119             :   !        0.5_core_rknd*( var_zt(k+1) + var_zt(k) ).  Rather, all explicit interpolations
     120             :   !        should call zt2zm or zm2zt; while interpolations for a variable being
     121             :   !        solved for implicitly in the code should use gr%weights_zt2zm (which
     122             :   !        refers to interp_weights_zt2zm_imp), or gr%weights_zm2zt (which
     123             :   !        refers to interp_weights_zm2zt_imp).
     124             :   !
     125             :   ! Momentum level 1 is placed at altitude zm_init, which is usually at the
     126             :   ! surface.  However, in general, zm_init can be at any altitude defined by the
     127             :   ! user.
     128             :   !
     129             :   ! GP indicates ghost points. Variables located at those levels are not
     130             :   ! prognosed, but only used for boundary conditions.
     131             :   !
     132             :   ! Chris Golaz, 7/17/99
     133             :   ! modified 9/10/99
     134             :   ! schemena, modified 6/11/2014 - Restructered code to add cubic/linear flag
     135             : 
     136             :   !  References:
     137             :   !
     138             :   !  https://arxiv.org/pdf/1711.03675v1.pdf#nameddest=url:clubb_grid
     139             :   !
     140             :   !  Section 3c, p. 3548 /Numerical discretization/ of:
     141             :   !   ``A PDF-Based Model for Boundary Layer Clouds. Part I:
     142             :   !     Method and Model Description'' Golaz, et al. (2002)
     143             :   !     JAS, Vol. 59, pp. 3540--3551.
     144             :   !-----------------------------------------------------------------------
     145             : 
     146             :   use clubb_precision, only: &
     147             :       core_rknd ! Variable(s)
     148             : 
     149             :   implicit none
     150             : 
     151             :   public :: grid, zt2zm, zm2zt, zt2zm2zt, zm2zt2zm, & 
     152             :             ddzm, ddzt, & 
     153             :             setup_grid, cleanup_grid, setup_grid_heights, &
     154             :             read_grid_heights, flip
     155             : 
     156             :   private :: t_above, t_below, m_above, m_below
     157             : 
     158             :   private ! Default Scoping
     159             : 
     160             :   ! Constant parameters
     161             :   integer, parameter :: & 
     162             :     t_above = 1, & ! Upper thermodynamic level index (gr%weights_zt2zm).
     163             :     t_below = 2, & ! Lower thermodynamic level index (gr%weights_zt2zm).
     164             :     m_above = 1, & ! Upper momentum level index (gr%weights_zm2zt).
     165             :     m_below = 2    ! Lower momentum level index (gr%weights_zm2zt).
     166             :   
     167             :   type grid
     168             : 
     169             :     integer :: nz ! Number of points in the grid
     170             :     !   Note: Fortran 90/95 prevents an allocatable array from appearing
     171             :     !   within a derived type.  However, Fortran 2003 does not!!!!!!!!
     172             :     real( kind = core_rknd ), allocatable, dimension(:,:) :: &
     173             :       zm, & ! Momentum grid
     174             :       zt    ! Thermo grid
     175             :       
     176             :     real( kind = core_rknd ), allocatable, dimension(:,:) :: &
     177             :       invrs_dzm, & ! The inverse spacing between thermodynamic grid
     178             :                    ! levels; centered over momentum grid levels.
     179             :       invrs_dzt    ! The inverse spacing between momentum grid levels;
     180             :                    ! centered over thermodynamic grid levels.
     181             : 
     182             :     real( kind = core_rknd ), allocatable, dimension(:,:) :: &
     183             :       dzm, &  ! Spacing between thermodynamic grid levels; centered over
     184             :               ! momentum grid levels
     185             :       dzt     ! Spcaing between momentum grid levels; centered over
     186             :               ! thermodynamic grid levels
     187             : 
     188             :     ! These weights are normally used in situations
     189             :     ! where a momentum level variable is being
     190             :     ! solved for implicitly in an equation and
     191             :     ! needs to be interpolated to the thermodynamic grid levels.
     192             :     real( kind = core_rknd ), allocatable, dimension(:,:,:) :: weights_zm2zt, & 
     193             :     ! These weights are normally used in situations where a
     194             :     ! thermodynamic level variable is being solved for implicitly in an equation
     195             :     ! and needs to be interpolated to the momentum grid levels.
     196             :                                      weights_zt2zm
     197             : 
     198             :   end type grid
     199             : 
     200             :   !   The grid is defined here so that it is common throughout the module.
     201             :   !   The implication is that only one grid can be defined !
     202             : 
     203             : 
     204             : !   Modification for using CLUBB in a host model (i.e. one grid per column)
     205             : 
     206             : 
     207             :   ! Interfaces provided for function overloading
     208             : 
     209             :   interface zt2zm
     210             :     ! For l_cubic_interp = .true.
     211             :     ! This version uses cublic spline interpolation of Stefen (1990).
     212             :     !
     213             :     ! For l_cubic_interp = .false.
     214             :     ! This performs a linear extension at the highest grid level and therefore
     215             :     ! does not guarantee, for positive definite quantities (e.g. wp2), that the
     216             :     ! extended point is indeed positive definite.  Positive definiteness can be
     217             :     ! ensured with a max statement.
     218             :     ! In the future, we could add a flag (lposdef) and, when needed, apply the
     219             :     ! max statement directly within interpolated_azm and interpolated_azmk.
     220             :     module procedure redirect_interpolated_azm_k    ! Works over a single vertical level
     221             :     module procedure redirect_interpolated_azm_1D   ! Works over all vertical levels 
     222             :     module procedure redirect_interpolated_azm_2D   ! Works over all vertical levels and columns
     223             :   end interface
     224             : 
     225             :   interface zm2zt
     226             :     ! For l_cubic_interp = .true.
     227             :     ! This version uses cublic spline interpolation of Stefen (1990).
     228             :     !
     229             :     ! For l_cubic_interp = .false.
     230             :     ! This performs a linear extension at the lowest grid level and therefore
     231             :     ! does not guarantee, for positive definite quantities (e.g. wp2), that the
     232             :     ! extended point is indeed positive definite.  Positive definiteness can be
     233             :     ! ensured with a max statement.
     234             :     ! In the future, we could add a flag (lposdef) and, when needed, apply the
     235             :     ! max statement directly within interpolated_azt and interpolated_aztk.
     236             :     module procedure redirect_interpolated_azt_k    ! Works over a single vertical level
     237             :     module procedure redirect_interpolated_azt_1D   ! Works over all vertical levels 
     238             :     module procedure redirect_interpolated_azt_2D   ! Works over all vertical levels and columns
     239             :   end interface
     240             : 
     241             :   ! Vertical derivative functions
     242             :   interface ddzm
     243             :     module procedure gradzm_1D    ! Works over all vertical levels 
     244             :     module procedure gradzm_2D    ! Works over all vertical levels and columns
     245             :   end interface
     246             : 
     247             :   interface ddzt
     248             :     module procedure gradzt_1D    ! Works over all vertical levels 
     249             :     module procedure gradzt_2D    ! Works over all vertical levels and columns
     250             :   end interface
     251             : 
     252             :   contains
     253             : 
     254             :   !=============================================================================
     255     4467528 :   subroutine setup_grid( nzmax, ngrdcol, sfc_elevation, l_implemented,  &
     256     4467528 :                          grid_type, deltaz, zm_init, zm_top,            &
     257     4467528 :                          momentum_heights, thermodynamic_heights,       &
     258             :                          gr )
     259             : 
     260             :     ! Description:
     261             :     !   Grid Constructor
     262             :     !
     263             :     !   This subroutine sets up the CLUBB vertical grid.
     264             :     !
     265             :     ! References:
     266             :     !   ``Equations for CLUBB'',  Sec. 8,  Grid Configuration.
     267             :     !-----------------------------------------------------------------------
     268             : 
     269             :     use constants_clubb, only:  & 
     270             :         fstderr ! Variable(s)
     271             : 
     272             :     use error_code, only: &
     273             :         clubb_at_least_debug_level, &   ! Procedure
     274             :         err_code, &                     ! Error indicator
     275             :         clubb_fatal_error, &            ! Constant
     276             :         err_header                      ! String
     277             : 
     278             :     use clubb_precision, only: &
     279             :         core_rknd ! Variable(s)
     280             : 
     281             :     implicit none
     282             : 
     283             :     ! Constant parameters
     284             :     integer, parameter :: & 
     285             :       NWARNING = 250 ! Issue a warning if nzmax exceeds this number.
     286             : 
     287             :     ! Input Variables
     288             :     integer, intent(in) ::  & 
     289             :       nzmax, &  ! Number of vertical levels in grid      [#]
     290             :       ngrdcol
     291             : 
     292             :     type(grid), target, intent(inout) :: gr
     293             : 
     294             :     real( kind = core_rknd ), dimension(ngrdcol), intent(in) ::  &
     295             :       sfc_elevation  ! Elevation of ground level    [m AMSL]
     296             : 
     297             :     ! Flag to see if CLUBB is running on it's own,
     298             :     ! or if it's implemented as part of a host model.
     299             :     logical, intent(in) :: l_implemented
     300             : 
     301             :     ! If CLUBB is running on it's own, this option determines if it is using:
     302             :     ! 1) an evenly-spaced grid;
     303             :     ! 2) a stretched (unevenly-spaced) grid entered on the thermodynamic grid
     304             :     !    levels (with momentum levels set halfway between thermodynamic levels);
     305             :     !    or
     306             :     ! 3) a stretched (unevenly-spaced) grid entered on the momentum grid levels
     307             :     !    (with thermodynamic levels set halfway between momentum levels).
     308             :     integer, intent(in) :: grid_type
     309             : 
     310             :     ! If the CLUBB model is running by itself, and is using an evenly-spaced
     311             :     ! grid (grid_type = 1), it needs the vertical grid spacing and
     312             :     ! momentum-level starting altitude as input.
     313             :     real( kind = core_rknd ), dimension(ngrdcol), intent(in) ::  & 
     314             :       deltaz,   & ! Vertical grid spacing                  [m]
     315             :       zm_init,  & ! Initial grid altitude (momentum level) [m]
     316             :       zm_top      ! Maximum grid altitude (momentum level) [m]
     317             : 
     318             :     ! If the CLUBB parameterization is implemented in a host model, it needs to
     319             :     ! use the host model's momentum level altitudes and thermodynamic level
     320             :     ! altitudes.
     321             :     ! If the CLUBB model is running by itself, but is using a stretched grid
     322             :     ! entered on thermodynamic levels (grid_type = 2), it needs to use the
     323             :     ! thermodynamic level altitudes as input.
     324             :     ! If the CLUBB model is running by itself, but is using a stretched grid
     325             :     ! entered on momentum levels (grid_type = 3), it needs to use the momentum
     326             :     ! level altitudes as input.
     327             :     real( kind = core_rknd ), intent(in), dimension(ngrdcol,nzmax) ::  & 
     328             :       momentum_heights,   & ! Momentum level altitudes (input)      [m]
     329             :       thermodynamic_heights ! Thermodynamic level altitudes (input) [m]
     330             : 
     331             :     ! Local Variables
     332             :     integer :: &
     333             :       begin_height, &  ! Lower bound for *_heights arrays [-]
     334             :       end_height       ! Upper bound for *_heights arrays [-]
     335             : 
     336             :     integer :: ierr, & ! Allocation stat
     337             :                i, k    ! Loop index
     338             : 
     339             : 
     340             :     ! ---- Begin Code ----
     341             : 
     342             :     ! Define the grid size
     343             : 
     344     4467528 :     if ( nzmax > NWARNING .and. clubb_at_least_debug_level( 1 ) ) then
     345             :       write(fstderr,*) "Warning:  running with vertical grid "// & 
     346           0 :                        "which is larger than", NWARNING, "levels."
     347           0 :       write(fstderr,*) "This may take a lot of CPU time and memory."
     348             :     end if
     349             : 
     350     4467528 :     gr%nz = nzmax
     351             : 
     352             :     ! Default bounds
     353     4467528 :     begin_height = 1
     354     4467528 :     end_height   = gr%nz
     355             : 
     356             :     !---------------------------------------------------
     357     4467528 :     if ( .not. l_implemented ) then
     358             :       
     359             :       ! Since l_implemented=.false., we have to calculate the number of vertical levels
     360             :       ! depending on the grid type. Currently, clubb can only have multiple columns
     361             :       ! if implemented in a host model, but if in the future we have a standalone 
     362             :       ! version with multiple columns, we may end up with columns that have a
     363             :       ! different number of vertical levels. CLUBB cannot handle this however, so
     364             :       ! we use only the first column to determine the number of vertical levels.
     365             : 
     366           0 :       if ( grid_type == 1 ) then
     367             : 
     368             :         ! Determine the number of grid points given the spacing
     369             :         ! to fit within the bounds without going over.
     370           0 :         gr%nz = floor( ( zm_top(1) - zm_init(1) + deltaz(1) ) / deltaz(1) )
     371             : 
     372           0 :       else if( grid_type == 2 ) then! Thermo
     373             : 
     374             :         ! Find begin_height (lower bound)
     375             :         k = gr%nz
     376             : 
     377           0 :         do while( thermodynamic_heights(1,k) >= zm_init(1) .and. k > 1 )
     378             : 
     379           0 :           k = k - 1
     380             : 
     381             :         end do
     382             : 
     383           0 :         if( thermodynamic_heights(1,k) >= zm_init(1) ) then
     384             : 
     385           0 :           write(fstderr,*) err_header, "Stretched zt grid cannot fulfill zm_init requirement"
     386           0 :           err_code = clubb_fatal_error
     387           0 :           return
     388             : 
     389             :         else
     390             : 
     391           0 :           begin_height = k
     392             : 
     393             :         end if
     394             : 
     395             :         ! Find end_height (upper bound)
     396             :         k = gr%nz
     397             : 
     398           0 :         do while( thermodynamic_heights(1,k) > zm_top(1) .and. k > 1 )
     399             : 
     400           0 :           k = k - 1
     401             : 
     402             :         end do
     403             : 
     404           0 :         if( zm_top(1) < thermodynamic_heights(1,k) ) then
     405             : 
     406           0 :           write(fstderr,*) err_header, "Stretched zt grid cannot fulfill zm_top requirement"
     407           0 :           err_code = clubb_fatal_error
     408           0 :           return
     409             : 
     410             :         else
     411             : 
     412           0 :           end_height = k
     413             : 
     414           0 :           gr%nz = size( thermodynamic_heights(1,begin_height:end_height ) )
     415             : 
     416             :         end if
     417             : 
     418           0 :       else if( grid_type == 3 ) then ! Momentum
     419             : 
     420             :         ! Find begin_height (lower bound)
     421             :         k = 1
     422             : 
     423           0 :         do while( momentum_heights(1,k) < zm_init(1) .and. k < gr%nz )
     424             : 
     425           0 :           k = k + 1
     426             : 
     427             :         end do
     428             : 
     429           0 :         if( momentum_heights(1,k) < zm_init(1) ) then
     430             : 
     431           0 :           write(fstderr,*) err_header, "Stretched zm grid cannot fulfill zm_init requirement"
     432           0 :           err_code = clubb_fatal_error
     433           0 :           return
     434             : 
     435             :         else
     436             : 
     437           0 :           begin_height = k
     438             : 
     439             :         end if
     440             : 
     441             :         ! Find end_height (lower bound)
     442             :         k = gr%nz
     443             : 
     444           0 :         do while( momentum_heights(1,k) > zm_top(1) .and. k > 1 )
     445             : 
     446           0 :           k = k - 1
     447             : 
     448             :         end do
     449             : 
     450           0 :         if( momentum_heights(1,k) > zm_top(1) ) then
     451             : 
     452           0 :           write(fstderr,*) err_header, "Stretched zm grid cannot fulfill zm_top requirement"
     453           0 :           err_code = clubb_fatal_error
     454           0 :           return
     455             : 
     456             :         else
     457             : 
     458           0 :           end_height = k
     459             : 
     460           0 :           gr%nz = size( momentum_heights(1,begin_height:end_height ) )
     461             : 
     462             :         end if
     463             : 
     464             :       endif ! grid_type
     465             : 
     466             :     endif ! .not. l_implemented
     467             : 
     468             :     !---------------------------------------------------
     469             : 
     470             :     ! Allocate memory for the grid levels
     471             :     allocate( gr%zm(ngrdcol,gr%nz), gr%zt(ngrdcol,gr%nz), & 
     472             :               gr%dzm(ngrdcol,gr%nz), gr%dzt(ngrdcol,gr%nz), &
     473             :               gr%invrs_dzm(ngrdcol,gr%nz), gr%invrs_dzt(ngrdcol,gr%nz),  & 
     474             :               gr%weights_zm2zt(ngrdcol,gr%nz,m_above:m_below), & 
     475             :               gr%weights_zt2zm(ngrdcol,gr%nz,t_above:t_below), & 
     476    89350560 :               stat=ierr )
     477             : 
     478     4467528 :     if ( ierr /= 0 ) then
     479           0 :       write(fstderr,*) err_header, "In setup_grid: allocation of grid variables failed."
     480           0 :       err_code = clubb_fatal_error
     481           0 :       return
     482             :     end if
     483             : 
     484             :     ! Set the values for the derived types used for heights, derivatives, and
     485             :     ! interpolation from the momentum/thermodynamic grid
     486             :     call setup_grid_heights( &
     487             :                   gr%nz, ngrdcol, & ! intent(in)
     488             :                   l_implemented, grid_type,  & ! intent(in)
     489             :                   deltaz, zm_init,  & ! intent(in)
     490           0 :                   momentum_heights(:,begin_height:end_height),  & ! intent(in) 
     491             :                   thermodynamic_heights(:,begin_height:end_height), & ! intent(in)
     492     4467528 :                   gr ) ! intent(inout)
     493             : 
     494    74597328 :     do i = 1, ngrdcol
     495    74597328 :       if ( sfc_elevation(i) > gr%zm(i,1) ) then
     496             :         write(fstderr,*) "The altitude of the lowest momentum level, "        &
     497             :                          // "gr%zm(1,1), must be at or above the altitude of "  &
     498             :                          // "the surface, sfc_elevation.  The lowest model "  &
     499           0 :                          // "momentum level cannot be below the surface."
     500           0 :         write(fstderr,*) "Altitude of lowest momentum level =", gr%zm(i,1)
     501           0 :         write(fstderr,*) "Altitude of the surface =", sfc_elevation(i)
     502           0 :         err_code = clubb_fatal_error
     503           0 :         return
     504             :       endif
     505             :     end do
     506             : 
     507             :     return
     508             : 
     509             :   end subroutine setup_grid
     510             : 
     511             :   !=============================================================================
     512           0 :   subroutine cleanup_grid( gr )
     513             : 
     514             :     ! Description:
     515             :     !   De-allocates the memory for the grid
     516             :     !
     517             :     ! References:
     518             :     !   None
     519             :     !------------------------------------------------------------------------------
     520             :     use constants_clubb, only: &
     521             :         fstderr  ! Constant(s)
     522             : 
     523             :     implicit none
     524             :   
     525             :     type(grid), target, intent(inout) :: gr
     526             : 
     527             :     ! Local Variable(s)
     528             :     integer :: ierr
     529             : 
     530             :     ! ----- Begin Code -----
     531             : 
     532             :     ! Allocate memory for grid levels
     533             :     deallocate( gr%zm, gr%zt, & 
     534             :                 gr%dzm, gr%dzt, &
     535             :                 gr%invrs_dzm, gr%invrs_dzt,  & 
     536             :                 gr%weights_zm2zt, gr%weights_zt2zm, & 
     537           0 :                 stat=ierr )
     538             : 
     539           0 :     if ( ierr /= 0 ) then
     540           0 :       write(fstderr,*) "Grid deallocation failed."
     541             :     end if
     542             : 
     543             : 
     544           0 :     return
     545             : 
     546             :   end subroutine cleanup_grid
     547             : 
     548             :   !=============================================================================
     549     4467528 :   subroutine setup_grid_heights( nz, ngrdcol, &
     550             :                                  l_implemented, grid_type,  & 
     551     4467528 :                                  deltaz, zm_init, momentum_heights,  & 
     552     4467528 :                                  thermodynamic_heights, &
     553             :                                  gr )
     554             : 
     555             :     ! Description:
     556             :     !   Sets the heights and interpolation weights of the column.
     557             :     !   This is seperated from setup_grid for those host models that have heights
     558             :     !   that vary with time.
     559             :     ! References:
     560             :     !   None
     561             :     !------------------------------------------------------------------------------
     562             : 
     563             :     use constants_clubb, only: &
     564             :         fstderr  ! Constant(s)
     565             : 
     566             :     use clubb_precision, only: &
     567             :         core_rknd  ! Variable(s)
     568             : 
     569             :     use error_code, only: &
     570             :         err_code, &                     ! Error indicator
     571             :         clubb_fatal_error, &            ! Constant
     572             :         err_header                      ! String
     573             : 
     574             :     implicit none
     575             : 
     576             :     ! Input Variables
     577             :     integer, intent(in) :: &
     578             :       nz, &
     579             :       ngrdcol
     580             : 
     581             :     type (grid), target, intent(inout) :: gr
     582             : 
     583             :     ! Flag to see if CLUBB is running on it's own,
     584             :     ! or if it's implemented as part of a host model.
     585             :     logical, intent(in) :: l_implemented
     586             : 
     587             :     ! If CLUBB is running on it's own, this option determines if it is using:
     588             :     ! 1) an evenly-spaced grid;
     589             :     ! 2) a stretched (unevenly-spaced) grid entered on the thermodynamic grid
     590             :     !    levels (with momentum levels set halfway between thermodynamic levels);
     591             :     !    or
     592             :     ! 3) a stretched (unevenly-spaced) grid entered on the momentum grid levels
     593             :     !    (with thermodynamic levels set halfway between momentum levels).
     594             :     integer, intent(in) :: grid_type
     595             : 
     596             :     ! If the CLUBB model is running by itself, and is using an evenly-spaced
     597             :     ! grid (grid_type = 1), it needs the vertical grid spacing and
     598             :     ! momentum-level starting altitude as input.
     599             :     real( kind = core_rknd ), dimension(ngrdcol), intent(in) ::  & 
     600             :       deltaz,   & ! Vertical grid spacing                  [m]
     601             :       zm_init     ! Initial grid altitude (momentum level) [m]
     602             : 
     603             : 
     604             :     ! If the CLUBB parameterization is implemented in a host model, it needs to
     605             :     ! use the host model's momentum level altitudes and thermodynamic level
     606             :     ! altitudes.
     607             :     ! If the CLUBB model is running by itself, but is using a stretched grid
     608             :     ! entered on thermodynamic levels (grid_type = 2), it needs to use the
     609             :     ! thermodynamic level altitudes as input.
     610             :     ! If the CLUBB model is running by itself, but is using a stretched grid
     611             :     ! entered on momentum levels (grid_type = 3), it needs to use the momentum
     612             :     ! level altitudes as input.
     613             :     real( kind = core_rknd ), intent(in), dimension(ngrdcol,nz) ::  & 
     614             :       momentum_heights,   & ! Momentum level altitudes (input)      [m]
     615             :       thermodynamic_heights ! Thermodynamic level altitudes (input) [m]
     616             : 
     617             :     integer :: k, i
     618             : 
     619             :     ! ---- Begin Code ----
     620             : 
     621     4467528 :     if ( .not. l_implemented ) then
     622             : 
     623             : 
     624           0 :       if ( grid_type == 1 ) then
     625             : 
     626             :         ! Evenly-spaced grid.
     627             :         ! Momentum level altitudes are defined based on the grid starting
     628             :         ! altitude, zm_init, the constant grid-spacing, deltaz, and the number
     629             :         ! of grid levels, gr%nz.
     630             : 
     631             :         ! Define momentum level altitudes. The first momentum level is at
     632             :         ! altitude zm_init.
     633           0 :         do k = 1, nz, 1
     634           0 :           do i = 1, ngrdcol
     635           0 :             gr%zm(i,k) = zm_init(i) + real( k-1, kind = core_rknd ) * deltaz(i)
     636             :           end do
     637             :         end do
     638             : 
     639             :         ! Define thermodynamic level altitudes.  Thermodynamic level altitudes
     640             :         ! are located at the central altitude levels, exactly halfway between
     641             :         ! momentum level altitudes.  The lowermost thermodynamic level is
     642             :         ! found by taking 1/2 the altitude difference between the bottom two
     643             :         ! momentum levels and subtracting that value from the bottom momentum
     644             :         ! level.  The first thermodynamic level is below zm_init.
     645           0 :         do i = 1, ngrdcol
     646           0 :           gr%zt(i,1) = zm_init(i) - ( 0.5_core_rknd * deltaz(i) )
     647             :         end do
     648             :         
     649           0 :         do k = 2, nz, 1
     650           0 :           do i = 1, ngrdcol
     651           0 :             gr%zt(i,k) = 0.5_core_rknd * ( gr%zm(i,k) + gr%zm(i,k-1) )
     652             :           end do
     653             :         enddo
     654             : 
     655             : 
     656           0 :       elseif ( grid_type == 2 ) then
     657             : 
     658             :         ! Stretched (unevenly-spaced) grid:  stretched thermodynamic levels.
     659             :         ! Thermodynamic levels are defined according to a stretched grid that
     660             :         ! is entered through the use of an input file.  This is similar to a
     661             :         ! SAM-style stretched grid.
     662             : 
     663             :         ! Define thermodynamic level altitudes.
     664           0 :         do k = 1, nz, 1
     665           0 :           do i = 1, ngrdcol
     666           0 :             gr%zt(i,k) = thermodynamic_heights(i,k)
     667             :           end do
     668             :         enddo
     669             : 
     670             :         ! Define momentum level altitudes.  Momentum level altitudes are
     671             :         ! located at the central altitude levels, exactly halfway between
     672             :         ! thermodynamic level altitudes.  The uppermost momentum level
     673             :         ! altitude is found by taking 1/2 the altitude difference between the
     674             :         ! top two thermodynamic levels and adding that value to the top
     675             :         ! thermodynamic level.
     676           0 :         do k = 1, nz-1, 1
     677           0 :           do i = 1, ngrdcol
     678           0 :             gr%zm(i,k) = 0.5_core_rknd * ( gr%zt(i,k+1) + gr%zt(i,k) )
     679             :           end do
     680             :         enddo
     681             :         
     682           0 :         do i = 1, ngrdcol
     683           0 :           gr%zm(i,nz) = gr%zt(i,nz) +  & 
     684           0 :                0.5_core_rknd * ( gr%zt(i,nz) - gr%zt(i,nz-1) )
     685             :         end do
     686             : 
     687           0 :       elseif ( grid_type == 3 ) then
     688             : 
     689             :         ! Stretched (unevenly-spaced) grid:  stretched momentum levels.
     690             :         ! Momentum levels are defined according to a stretched grid that is
     691             :         ! entered through the use of an input file.  This is similar to a
     692             :         ! WRF-style stretched grid.
     693             : 
     694             :         ! Define momentum level altitudes.
     695           0 :         do k = 1, nz, 1
     696           0 :           do i = 1, ngrdcol
     697           0 :             gr%zm(i,k) = momentum_heights(i,k)
     698             :           end do
     699             :         enddo
     700             : 
     701             :         ! Define thermodynamic level altitudes.  Thermodynamic level altitudes
     702             :         ! are located at the central altitude levels, exactly halfway between
     703             :         ! momentum level altitudes.  The lowermost thermodynamic level
     704             :         ! altitude is found by taking 1/2 the altitude difference between the
     705             :         ! bottom two momentum levels and subtracting that value from the
     706             :         ! bottom momentum level.
     707           0 :         do i = 1, ngrdcol
     708           0 :           gr%zt(i,1) = gr%zm(i,1) - 0.5_core_rknd * ( gr%zm(i,2) - gr%zm(i,1) )
     709             :         end do
     710             :         
     711           0 :         do k = 2, nz, 1
     712           0 :           do i = 1, ngrdcol
     713           0 :             gr%zt(i,k) = 0.5_core_rknd * ( gr%zm(i,k) + gr%zm(i,k-1) )
     714             :           end do
     715             :         enddo
     716             : 
     717             : 
     718             :       else
     719             : 
     720             :         ! Invalid grid type.
     721           0 :         write(fstderr,*) err_header, "Invalid grid type: ", grid_type, & 
     722           0 :                          ".  Valid options are 1, 2, or 3."
     723           0 :         err_code = clubb_fatal_error
     724           0 :         return
     725             : 
     726             : 
     727             :       endif
     728             : 
     729             : 
     730             :     else
     731             : 
     732             :       ! The CLUBB parameterization is implemented in a host model.
     733             :       ! Use the host model's momentum level altitudes and thermodynamic level
     734             :       ! altitudes to set up the CLUBB grid.
     735             : 
     736             :       ! Momentum level altitudes from host model.
     737   384207408 :       do k = 1, nz, 1
     738  6345240408 :         do i = 1, ngrdcol
     739  6340772880 :           gr%zm(i,k) = momentum_heights(i,k)
     740             :         end do
     741             :       enddo
     742             : 
     743             :       ! Thermodynamic level altitudes from host model after possible grid-index
     744             :       ! adjustment for CLUBB interface.
     745   384207408 :       do k = 1, nz, 1
     746  6345240408 :         do i = 1, ngrdcol
     747  6340772880 :           gr%zt(i,k) = thermodynamic_heights(i,k)
     748             :         end do
     749             :       enddo
     750             : 
     751             : 
     752             :     endif ! not l_implemented
     753             : 
     754             : 
     755             :     ! Define dzm, the spacing between thermodynamic grid levels; centered over
     756             :     ! momentum grid levels
     757   379739880 :     do k=1,nz-1
     758  6270643080 :       do i = 1, ngrdcol
     759  6266175552 :         gr%dzm(i,k) = gr%zt(i,k+1) - gr%zt(i,k)
     760             :       end do
     761             :     enddo
     762             :     
     763    74597328 :     do i = 1, ngrdcol
     764    74597328 :       gr%dzm(i,nz) = gr%dzm(i,nz-1)
     765             :     end do
     766             : 
     767             :     ! Define dzt, the spacing between momentum grid levels; centered over
     768             :     ! thermodynamic grid levels
     769   379739880 :     do k=2,nz
     770  6270643080 :       do i = 1, ngrdcol
     771  6266175552 :         gr%dzt(i,k) = gr%zm(i,k) - gr%zm(i,k-1)
     772             :       end do
     773             :     enddo
     774             :     
     775    74597328 :     do i = 1, ngrdcol
     776    74597328 :       gr%dzt(i,1) = gr%dzt(i,2)
     777             :     end do
     778             : 
     779             :     ! Define invrs_dzm, which is the inverse spacing between thermodynamic grid
     780             :     ! levels; centered over momentum grid levels.
     781   379739880 :     do k=1,nz-1
     782  6270643080 :       do i = 1, ngrdcol
     783  6266175552 :         gr%invrs_dzm(i,k) = 1._core_rknd / ( gr%zt(i,k+1) - gr%zt(i,k) )
     784             :       end do
     785             :     enddo
     786             :     
     787    74597328 :     do i = 1, ngrdcol
     788    74597328 :       gr%invrs_dzm(i,nz) = gr%invrs_dzm(i,nz-1)
     789             :     end do
     790             : 
     791             : 
     792             :     ! Define invrs_dzt, which is the inverse spacing between momentum grid
     793             :     ! levels; centered over thermodynamic grid levels.
     794   379739880 :     do k=2,nz
     795  6270643080 :       do i = 1, ngrdcol
     796  6266175552 :         gr%invrs_dzt(i,k) = 1._core_rknd / ( gr%zm(i,k) - gr%zm(i,k-1) )
     797             :       end do
     798             :     enddo
     799             :     
     800    74597328 :     do i = 1, ngrdcol
     801    74597328 :       gr%invrs_dzt(i,1) = gr%invrs_dzt(i,2)
     802             :     end do
     803             : 
     804             : 
     805             :     ! Interpolation Weights: zm grid to zt grid.
     806             :     ! The grid index (k) is the index of the level on the thermodynamic (zt)
     807             :     ! grid.  The result is the weights of the upper and lower momentum levels
     808             :     ! (that sandwich the thermodynamic level) applied to that thermodynamic
     809             :     ! level.  These weights are normally used in situations where a momentum
     810             :     ! level variable is being solved for implicitly in an equation, and the
     811             :     ! aforementioned variable needs to be interpolated from three successive
     812             :     ! momentum levels (the central momentum level, as well as one momentum level
     813             :     ! above and below the central momentum level) to the intermediate
     814             :     ! thermodynamic grid levels that sandwich the central momentum level.
     815             :     ! For more information, see the comments in function interpolated_aztk_imp.
     816             :     call calc_zm2zt_weights( nz, ngrdcol, & ! In
     817     4467528 :                              gr )           ! InOut
     818             : 
     819             : 
     820             :     ! Interpolation Weights: zt grid to zm grid.
     821             :     ! The grid index (k) is the index of the level on the momentum (zm) grid.
     822             :     ! The result is the weights of the upper and lower thermodynamic levels
     823             :     ! (that sandwich the momentum level) applied to that momentum level.  These
     824             :     ! weights are normally used in situations where a thermodynamic level
     825             :     ! variable is being solved for implicitly in an equation, and the
     826             :     ! aforementioned variable needs to be interpolated from three successive
     827             :     ! thermodynamic levels (the central thermodynamic level, as well as one
     828             :     ! thermodynamic level above and below the central thermodynamic level) to
     829             :     ! the intermediate momentum grid levels that sandwich the central
     830             :     ! thermodynamic level.
     831             :     ! For more information, see the comments in function interpolated_azmk_imp.
     832             :     call calc_zt2zm_weights( nz, ngrdcol, & ! In
     833     4467528 :                              gr )           ! InOut
     834             : 
     835     4467528 :     return
     836             :   end subroutine setup_grid_heights
     837             : 
     838             :   !=============================================================================
     839           0 :   subroutine read_grid_heights( nzmax, grid_type,  & 
     840             :                                 zm_grid_fname, zt_grid_fname, & 
     841             :                                 file_unit, &
     842           0 :                                 momentum_heights, & 
     843           0 :                                 thermodynamic_heights )
     844             : 
     845             :     ! Description:
     846             :     ! This subroutine is used foremost in cases where the grid_type corresponds
     847             :     ! with the stretched (unevenly-spaced) grid options (either grid_type = 2 or
     848             :     ! grid_type = 3).  This subroutine reads in the values of the stretched grid
     849             :     ! altitude levels for either the thermodynamic level grid or the momentum
     850             :     ! level grid.  This subroutine also handles basic error checking for all
     851             :     ! three grid types.
     852             :     !------------------------------------------------------------------------
     853             : 
     854             :     use constants_clubb, only:  & 
     855             :         fstderr ! Variable(s)
     856             : 
     857             :     use file_functions, only:  & 
     858             :         file_read_1d ! Procedure(s)
     859             : 
     860             :     use clubb_precision, only: &
     861             :         core_rknd ! Variable(s)
     862             : 
     863             :     use error_code, only: &
     864             :         err_code, &                     ! Error indicator
     865             :         clubb_fatal_error, &            ! Constant
     866             :         err_header                      ! String
     867             : 
     868             :     implicit none
     869             : 
     870             :     ! Input Variables.
     871             : 
     872             :     ! Declared number of vertical levels.
     873             :     integer, intent(in) :: & 
     874             :       nzmax
     875             : 
     876             :     ! If CLUBB is running on it's own, this option determines if it is using:
     877             :     ! 1) an evenly-spaced grid;
     878             :     ! 2) a stretched (unevenly-spaced) grid entered on the thermodynamic grid
     879             :     !    levels (with momentum levels set halfway between thermodynamic levels);
     880             :     !    or
     881             :     ! 3) a stretched (unevenly-spaced) grid entered on the momentum grid levels
     882             :     !    (with thermodynamic levels set halfway between momentum levels).
     883             :     integer, intent(in) :: & 
     884             :       grid_type
     885             : 
     886             :     character(len=*), intent(in) :: & 
     887             :       zm_grid_fname,  & ! Path and filename of file for momentum level altitudes
     888             :       zt_grid_fname     ! Path and filename of file for thermodynamic level altitudes
     889             : 
     890             :     integer, intent(in) :: &
     891             :       file_unit   ! Unit number for zt_grid_fname & zm_grid_fname (based on the OpenMP thread)
     892             : 
     893             :     ! Output Variables.
     894             : 
     895             :     ! If the CLUBB model is running by itself, but is using a stretched grid
     896             :     ! entered on thermodynamic levels (grid_type = 2), it needs to use the
     897             :     ! thermodynamic level altitudes as input.
     898             :     ! If the CLUBB model is running by itself, but is using a stretched grid
     899             :     ! entered on momentum levels (grid_type = 3), it needs to use the momentum
     900             :     ! level altitudes as input.
     901             :     real( kind = core_rknd ), dimension(nzmax), intent(out) :: & 
     902             :       momentum_heights,   & ! Momentum level altitudes (file input)      [m]
     903             :       thermodynamic_heights ! Thermodynamic level altitudes (file input) [m]
     904             : 
     905             :     ! Local Variables.
     906             : 
     907             :     integer :: &
     908             :       zt_level_count,  & ! Number of altitudes found in zt_grid_fname
     909             :       zm_level_count     ! Number of altitudes found in zm_grid_fname
     910             : 
     911             :     integer :: input_status   ! Status of file being read:
     912             :     ! > 0 ==> error reading file.
     913             :     ! = 0 ==> no error and more file to be read.
     914             :     ! < 0 ==> end of file indicator.
     915             : 
     916             :     ! Generic variable for storing file data while counting the number
     917             :     ! of file entries.
     918             :     real( kind = core_rknd ) :: generic_input_item
     919             : 
     920             :     integer :: k   ! Loop index
     921             : 
     922             :     ! ---- Begin Code ----
     923             : 
     924             :     ! Declare the momentum level altitude array and the thermodynamic level
     925             :     ! altitude array to be 0 until overwritten.
     926           0 :     momentum_heights(1:nzmax) = 0.0_core_rknd
     927           0 :     thermodynamic_heights(1:nzmax) = 0.0_core_rknd
     928             : 
     929             :     ! Avoid uninitialized memory
     930           0 :     generic_input_item = 0.0_core_rknd
     931             : 
     932             : 
     933           0 :     if ( grid_type == 1 ) then
     934             : 
     935             :       ! Evenly-spaced grid.
     936             :       ! Grid level altitudes are based on a constant distance between them and
     937             :       ! a starting point for the bottom of the grid.
     938             : 
     939             :       ! As a way of error checking, make sure that there isn't any file entry
     940             :       ! for either momentum level altitudes or thermodynamic level altitudes.
     941           0 :       if ( zm_grid_fname /= '' ) then
     942           0 :         write(fstderr,*) err_header, & 
     943             :            "An evenly-spaced grid has been selected. " & 
     944           0 :            // " Please reset zm_grid_fname to ''."
     945           0 :         err_code = clubb_fatal_error
     946           0 :         return
     947             :       endif
     948           0 :       if ( zt_grid_fname /= '' ) then
     949           0 :         write(fstderr,*) err_header, & 
     950             :            "An evenly-spaced grid has been selected. " & 
     951           0 :            // " Please reset zt_grid_fname to ''."
     952           0 :         err_code = clubb_fatal_error
     953           0 :         return
     954             :       endif
     955             : 
     956             : 
     957           0 :     elseif ( grid_type == 2 ) then
     958             : 
     959             :       ! Stretched (unevenly-spaced) grid:  stretched thermodynamic levels.
     960             :       ! Thermodynamic levels are defined according to a stretched grid that is
     961             :       ! entered through the use of an input file.  Momentum levels are set
     962             :       ! halfway between thermodynamic levels.  This is similar to a SAM-style
     963             :       ! stretched grid.
     964             : 
     965             :       ! As a way of error checking, make sure that there isn't any file entry
     966             :       ! for momentum level altitudes.
     967           0 :       if ( zm_grid_fname /= '' ) then
     968           0 :         write(fstderr,*) err_header, & 
     969             :            "Thermodynamic level altitudes have been selected " & 
     970             :            // "for use in a stretched (unevenly-spaced) grid. " & 
     971           0 :            // " Please reset zm_grid_fname to ''."
     972           0 :         err_code = clubb_fatal_error
     973           0 :         return
     974             :       endif
     975             : 
     976             : !$omp critical
     977             :       ! Open the file zt_grid_fname.
     978             :       open( unit=file_unit, file=zt_grid_fname,  &
     979           0 :             status='old', action='read' )
     980             : 
     981             :       ! Find the number of thermodynamic level altitudes listed
     982             :       ! in file zt_grid_fname.
     983           0 :       zt_level_count = 0
     984           0 :       do
     985             :         read( unit=file_unit, fmt=*, iostat=input_status )  & 
     986           0 :            generic_input_item
     987           0 :         if ( input_status < 0 ) exit   ! end of file indicator
     988           0 :         if ( input_status > 0 ) then
     989           0 :           write(fstderr,*) err_header, &  ! error reading input
     990           0 :              "Error reading thermodynamic level input file."
     991           0 :           err_code = clubb_fatal_error
     992           0 :           exit
     993             :         end if
     994           0 :         zt_level_count = zt_level_count + 1
     995             :       enddo
     996             : 
     997             :       ! Close the file zt_grid_fname.
     998           0 :       close( unit=file_unit )
     999             : !$omp end critical
    1000             :       
    1001           0 :       if ( err_code == clubb_fatal_error ) return
    1002             : 
    1003             :       ! Check that the number of thermodynamic grid altitudes in the input file
    1004             :       ! matches the declared number of CLUBB grid levels (nzmax).
    1005           0 :       if ( zt_level_count /= nzmax ) then
    1006             :         write(fstderr,*)  & 
    1007             :            "The number of thermodynamic grid altitudes " & 
    1008             :            // "listed in file " // trim(zt_grid_fname)  & 
    1009             :            // " does not match the number of CLUBB grid " & 
    1010           0 :            // "levels specified in the model.in file."
    1011             :         write(fstderr,*) & 
    1012           0 :            "Number of thermodynamic grid altitudes listed:  ", & 
    1013           0 :            zt_level_count
    1014             :         write(fstderr,*) & 
    1015           0 :            "Number of CLUBB grid levels specified:  ", nzmax
    1016           0 :         err_code = clubb_fatal_error
    1017           0 :         return
    1018             :       endif
    1019             : 
    1020             :       ! Read the thermodynamic level altitudes from zt_grid_fname.
    1021             :       call file_read_1d( file_unit, zt_grid_fname, nzmax, 1,  & ! intent(in)
    1022           0 :                          thermodynamic_heights ) ! intent(out)
    1023             : 
    1024             :       ! Check that each thermodynamic level altitude increases
    1025             :       ! in height as the thermodynamic level grid index increases.
    1026           0 :       do k = 2, nzmax, 1
    1027           0 :         if ( thermodynamic_heights(k)  & 
    1028           0 :              <= thermodynamic_heights(k-1) ) then
    1029             :           write(fstderr,*)  & 
    1030             :              "The declared thermodynamic level grid " & 
    1031             :              // "altitudes are not increasing in height " & 
    1032           0 :              // "as grid level index increases."
    1033             :           write(fstderr,*) & 
    1034           0 :              "Grid index:  ", k-1, ";", & 
    1035           0 :              "  Thermodynamic level altitude:  ", & 
    1036           0 :              thermodynamic_heights(k-1)
    1037             :           write(fstderr,*) & 
    1038           0 :              "Grid index:  ", k, ";", & 
    1039           0 :              "  Thermodynamic level altitude:  ", & 
    1040           0 :              thermodynamic_heights(k)
    1041           0 :           err_code = clubb_fatal_error
    1042           0 :           return
    1043             :         endif
    1044             :       enddo
    1045             : 
    1046             : 
    1047           0 :     elseif ( grid_type == 3 ) then
    1048             : 
    1049             :       ! Stretched (unevenly-spaced) grid:  stretched momentum levels.
    1050             :       ! Momentum levels are defined according to a stretched grid that is
    1051             :       ! entered through the use of an input file.  Thermodynamic levels are set
    1052             :       ! halfway between momentum levels.  This is similar to a WRF-style
    1053             :       ! stretched grid.
    1054             : 
    1055             :       ! As a way of error checking, make sure that there isn't any file entry
    1056             :       ! for thermodynamic level altitudes.
    1057           0 :       if ( zt_grid_fname /= '' ) then
    1058             :         write(fstderr,*) & 
    1059             :            "Momentum level altitudes have been selected " & 
    1060             :            // "for use in a stretched (unevenly-spaced) grid. " & 
    1061           0 :            // " Please reset zt_grid_fname to ''."
    1062           0 :         err_code = clubb_fatal_error
    1063           0 :         return
    1064             :       endif
    1065             : 
    1066             :       ! Open the file zm_grid_fname.
    1067             :       open( unit=file_unit, file=zm_grid_fname,  & 
    1068           0 :             status='old', action='read' )
    1069             : 
    1070             :       ! Find the number of momentum level altitudes
    1071             :       ! listed in file zm_grid_fname.
    1072           0 :       zm_level_count = 0
    1073           0 :       do
    1074             :         read( unit=file_unit, fmt=*, iostat=input_status ) & 
    1075           0 :            generic_input_item
    1076           0 :         if ( input_status < 0 ) exit   ! end of file indicator
    1077           0 :         if ( input_status > 0 ) then
    1078             : 
    1079           0 :           write(fstderr,*) err_header, &! error reading input
    1080           0 :                            "Error reading momentum level input file."
    1081           0 :           err_code = clubb_fatal_error
    1082           0 :           return
    1083             :         end if
    1084           0 :         zm_level_count = zm_level_count + 1
    1085             :       enddo
    1086             : 
    1087             :       ! Close the file zm_grid_fname.
    1088           0 :       close( unit=file_unit )
    1089             : 
    1090             :       ! Check that the number of momentum grid altitudes in the input file
    1091             :       ! matches the declared number of CLUBB grid levels (nzmax).
    1092           0 :       if ( zm_level_count /= nzmax ) then
    1093             :         write(fstderr,*) & 
    1094             :            "The number of momentum grid altitudes " & 
    1095             :            // "listed in file " // trim(zm_grid_fname) & 
    1096             :            // " does not match the number of CLUBB grid " & 
    1097           0 :            // "levels specified in the model.in file."
    1098             :         write(fstderr,*) & 
    1099           0 :            "Number of momentum grid altitudes listed:  ", & 
    1100           0 :            zm_level_count
    1101             :         write(fstderr,*) & 
    1102           0 :            "Number of CLUBB grid levels specified:  ", nzmax
    1103           0 :         err_code = clubb_fatal_error
    1104           0 :         return
    1105             :       endif
    1106             : 
    1107             :       ! Read the momentum level altitudes from zm_grid_fname.
    1108             :       call file_read_1d( file_unit, zm_grid_fname, nzmax, 1, & ! intent(in) 
    1109           0 :                          momentum_heights ) ! intent(out)
    1110             : 
    1111             :       ! Check that each momentum level altitude increases in height as the
    1112             :       ! momentum level grid index increases.
    1113           0 :       do k = 2, nzmax, 1
    1114           0 :         if ( momentum_heights(k)  & 
    1115           0 :              <= momentum_heights(k-1) ) then
    1116             :           write(fstderr,*)  & 
    1117             :              "The declared momentum level grid " & 
    1118             :              // "altitudes are not increasing in height " & 
    1119           0 :              // "as grid level index increases."
    1120             :           write(fstderr,*) & 
    1121           0 :              "Grid index:  ", k-1, ";", & 
    1122           0 :              "  Momentum level altitude:  ", & 
    1123           0 :              momentum_heights(k-1)
    1124             :           write(fstderr,*) & 
    1125           0 :              "Grid index:  ", k, ";", & 
    1126           0 :              "  Momentum level altitude:  ", & 
    1127           0 :              momentum_heights(k)
    1128           0 :           err_code = clubb_fatal_error
    1129           0 :           return
    1130             :         endif
    1131             :       enddo
    1132             : 
    1133             : 
    1134             :     endif
    1135             : 
    1136             : 
    1137             :     ! The purpose of this if statement is to avoid a compiler warning.
    1138             :     if ( generic_input_item > 0.0_core_rknd ) then
    1139             :       ! Do nothing
    1140             :     endif
    1141             :     ! Joshua Fasching  June 2008
    1142             : 
    1143             :     return
    1144             : 
    1145             :   end subroutine read_grid_heights
    1146             :   
    1147             :   !=============================================================================
    1148             :   ! Wrapped in interface zt2zm
    1149           0 :   function redirect_interpolated_azm_k( gr, azt, k, zm_min )
    1150             : 
    1151             :     ! Description:
    1152             :     ! Calls the appropriate corresponding function based on l_cubic_temp
    1153             :     !-----------------------------------------------------------------------
    1154             : 
    1155             :     use clubb_precision, only: &
    1156             :         core_rknd  ! Variable(s)
    1157             : 
    1158             :     use model_flags, only: &
    1159             :         l_cubic_interp, & ! Variable(s)
    1160             :         l_quintic_poly_interp
    1161             : 
    1162             :     use constants_clubb, only: &
    1163             :         fstdout ! Variable
    1164             : 
    1165             :     implicit none
    1166             : 
    1167             :     !---------------------------- Input Variables ----------------------------
    1168             :     type (grid), target, intent(in) :: gr
    1169             : 
    1170             :     real( kind = core_rknd ), intent(in), dimension(gr%nz) :: &
    1171             :       azt    ! Variable on thermodynamic grid levels    [units vary]
    1172             :       
    1173             :     integer, intent(in) :: &
    1174             :       k
    1175             : 
    1176             :     !---------------------------- Return Variable ----------------------------
    1177             :     real( kind = core_rknd ) :: &
    1178             :       redirect_interpolated_azm_k    ! Variable when interp. to momentum levels
    1179             : 
    1180             :     !---------------------------- Optional Input Variable ----------------------------
    1181             :     real( kind = core_rknd ), optional, intent(in) :: &
    1182             :       zm_min
    1183             : 
    1184             :     !---------------------------- Local Variables ----------------------------
    1185             :     real( kind = core_rknd ), dimension(1,gr%nz) :: &
    1186           0 :       redirect_interpolated_azm_col    ! Variable when interp. to momentum levels
    1187             :       
    1188             :     !---------------------------- Begin Code ----------------------------
    1189             : 
    1190           0 :     redirect_interpolated_azm_col = redirect_interpolated_azm_2D( gr%nz, 1, gr, azt, zm_min )
    1191             : 
    1192           0 :     redirect_interpolated_azm_k = redirect_interpolated_azm_col(1,k)
    1193             : 
    1194             :     return
    1195             :   end function redirect_interpolated_azm_k
    1196             :   
    1197             :   !=============================================================================
    1198             :   ! Wrapped in interface zt2zm
    1199           0 :   function redirect_interpolated_azm_1D( gr, azt, zm_min )
    1200             : 
    1201             :     ! Description:
    1202             :     ! Calls the appropriate corresponding function based on l_cubic_temp
    1203             :     !-----------------------------------------------------------------------
    1204             : 
    1205             :     use clubb_precision, only: &
    1206             :         core_rknd  ! Variable(s)
    1207             : 
    1208             :     use model_flags, only: &
    1209             :         l_cubic_interp, & ! Variable(s)
    1210             :         l_quintic_poly_interp
    1211             : 
    1212             :     use constants_clubb, only: &
    1213             :         fstdout ! Variable
    1214             : 
    1215             :     implicit none
    1216             : 
    1217             :     !---------------------------- Input Variables ----------------------------
    1218             :     type (grid), target, intent(in) :: gr
    1219             : 
    1220             :     real( kind = core_rknd ), intent(in), dimension(gr%nz) :: &
    1221             :       azt    ! Variable on thermodynamic grid levels    [units vary]
    1222             : 
    1223             :     !---------------------------- Return Variable ----------------------------
    1224             :     real( kind = core_rknd ), dimension(gr%nz) :: &
    1225             :       redirect_interpolated_azm_1D    ! Variable when interp. to momentum levels
    1226             : 
    1227             :     !---------------------------- Optional Input Variable ----------------------------
    1228             :     real( kind = core_rknd ), optional, intent(in) :: &
    1229             :       zm_min
    1230             : 
    1231             :     !---------------------------- Local Variables ----------------------------
    1232             :     real( kind = core_rknd ), dimension(1,gr%nz) :: &
    1233           0 :       azt_col    ! Variable on thermodynamic grid levels    [units vary]
    1234             : 
    1235             :     real( kind = core_rknd ), dimension(1,gr%nz) :: &
    1236           0 :       redirect_interpolated_azm_col    ! Variable when interp. to momentum levels
    1237             : 
    1238             :     !---------------------------- Begin Code ----------------------------
    1239             : 
    1240           0 :     azt_col(1,:) = azt
    1241             : 
    1242           0 :     redirect_interpolated_azm_col = redirect_interpolated_azm_2D( gr%nz, 1, gr, azt_col, zm_min )
    1243             : 
    1244           0 :     redirect_interpolated_azm_1D = redirect_interpolated_azm_col(1,:)
    1245             : 
    1246           0 :     return
    1247             :   end function redirect_interpolated_azm_1D
    1248             : 
    1249             :   !=============================================================================
    1250             :   ! Wrapped in interface zt2zm
    1251   393142464 :   function redirect_interpolated_azm_2D( nz, ngrdcol, gr, azt, zm_min )
    1252             : 
    1253             :     ! Description:
    1254             :     ! Calls the appropriate corresponding function based on l_cubic_temp
    1255             :     !-----------------------------------------------------------------------
    1256             : 
    1257             :     use clubb_precision, only: &
    1258             :         core_rknd  ! Variable(s)
    1259             : 
    1260             :     use model_flags, only: &
    1261             :         l_cubic_interp, & ! Variable(s)
    1262             :         l_quintic_poly_interp
    1263             : 
    1264             :     use constants_clubb, only: &
    1265             :         fstdout ! Variable
    1266             : 
    1267             :     implicit none
    1268             :     
    1269             :     !---------------------------- Input Variables ----------------------------
    1270             :     integer, intent(in) :: &
    1271             :       nz, &
    1272             :       ngrdcol
    1273             : 
    1274             :     type (grid), target, intent(in) :: gr
    1275             : 
    1276             :     real( kind = core_rknd ), intent(in), dimension(ngrdcol,nz) :: &
    1277             :       azt    ! Variable on thermodynamic grid levels    [units vary]
    1278             : 
    1279             :     !---------------------------- Return Variable ----------------------------
    1280             :     real( kind = core_rknd ), dimension(ngrdcol,nz) :: &
    1281             :       redirect_interpolated_azm_2D    ! Variable when interp. to momentum levels
    1282             :       
    1283             :     !---------------------------- Optional Input Variable ----------------------------
    1284             :     real( kind = core_rknd ), optional, intent(in) :: &
    1285             :       zm_min
    1286             : 
    1287             :     !---------------------------- Begin Code ----------------------------
    1288             : 
    1289             :     ! Sanity Check
    1290             :     if (l_quintic_poly_interp) then
    1291             :       if (.not. l_cubic_interp) then
    1292             :         write (fstdout, *) "Error: Model flag l_quintic_poly_interp should not be true if "&
    1293             :                          //"l_cubic_interp is false."
    1294             :         error stop
    1295             :       end if
    1296             :     end if
    1297             : 
    1298             :     ! Redirect
    1299             :     if ( l_cubic_interp .and. nz >= 3 ) then
    1300             :       redirect_interpolated_azm_2D = cubic_interpolated_azm_2D( nz, ngrdcol, gr, azt )
    1301             :     else
    1302             :       call linear_interpolated_azm_2D( nz, ngrdcol, gr, azt, &
    1303             :                                        redirect_interpolated_azm_2D, &
    1304   393142464 :                                        zm_min )
    1305             :     end if
    1306             : 
    1307   393142464 :     return
    1308             :   end function redirect_interpolated_azm_2D
    1309             :   
    1310             :   !=============================================================================
    1311             :   ! Wrapped in interface zm2zt
    1312           0 :   function redirect_interpolated_azt_k( gr, azm, k, zt_min )
    1313             : 
    1314             :     ! Description:
    1315             :     ! Calls the appropriate corresponding function based on l_cubic_temp
    1316             :     !-----------------------------------------------------------------------
    1317             : 
    1318             :     use clubb_precision, only: &
    1319             :         core_rknd  ! Variable(s)
    1320             : 
    1321             :     use model_flags, only: &
    1322             :         l_cubic_interp, & ! Variable(s)
    1323             :         l_quintic_poly_interp
    1324             : 
    1325             :     use constants_clubb, only: &
    1326             :         fstdout ! Variable
    1327             : 
    1328             :     implicit none
    1329             : 
    1330             :     type (grid), target, intent(in) :: gr
    1331             : 
    1332             :     !---------------------------- Input Variables ----------------------------
    1333             :     real( kind = core_rknd ), intent(in), dimension(gr%nz) :: &
    1334             :       azm    ! Variable on momentum grid levels    [units vary]
    1335             :       
    1336             :     integer, intent(in) :: &
    1337             :       k   ! Vertical level
    1338             : 
    1339             :     !---------------------------- Return Variable ----------------------------
    1340             :     real( kind = core_rknd ) :: &
    1341             :       redirect_interpolated_azt_k    ! Variable when interp. to momentum levels
    1342             : 
    1343             :     !---------------------------- Optional Input Variable ----------------------------
    1344             :     real( kind = core_rknd ), optional, intent(in) :: &
    1345             :       zt_min
    1346             : 
    1347             :     !---------------------------- Local Variables ----------------------------
    1348             : 
    1349             :     real( kind = core_rknd ), dimension(1,gr%nz) :: &
    1350           0 :       redirect_interpolated_azt_col    ! Variable when interp. to momentum levels
    1351             :     
    1352             :     !---------------------------- Begin Code ----------------------------
    1353             : 
    1354           0 :     redirect_interpolated_azt_col = redirect_interpolated_azt_2D( gr%nz, 1, gr, azm, zt_min )
    1355             : 
    1356           0 :     redirect_interpolated_azt_k = redirect_interpolated_azt_col(1,k)
    1357             : 
    1358             :     return
    1359             :   end function redirect_interpolated_azt_k
    1360             :   
    1361             : 
    1362             :   !=============================================================================
    1363             :   ! Wrapped in interface zm2zt
    1364           0 :   function redirect_interpolated_azt_1D( gr, azm, zt_min )
    1365             : 
    1366             :     ! Description:
    1367             :     ! Calls the appropriate corresponding function based on l_cubic_temp
    1368             :     !-----------------------------------------------------------------------
    1369             : 
    1370             :     use clubb_precision, only: &
    1371             :         core_rknd  ! Variable(s)
    1372             : 
    1373             :     use model_flags, only: &
    1374             :         l_cubic_interp, & ! Variable(s)
    1375             :         l_quintic_poly_interp
    1376             : 
    1377             :     use constants_clubb, only: &
    1378             :         fstdout ! Variable
    1379             : 
    1380             :     implicit none
    1381             : 
    1382             :     type (grid), target, intent(in) :: gr
    1383             : 
    1384             :     !---------------------------- Input Variables ----------------------------
    1385             :     real( kind = core_rknd ), intent(in), dimension(gr%nz) :: &
    1386             :       azm    ! Variable on momentum grid levels    [units vary]
    1387             : 
    1388             :     !---------------------------- Return Variable ----------------------------
    1389             :     real( kind = core_rknd ), dimension(gr%nz) :: &
    1390             :       redirect_interpolated_azt_1D    ! Variable when interp. to momentum levels
    1391             : 
    1392             :     !---------------------------- Optional Input Variable ----------------------------
    1393             :     real( kind = core_rknd ), optional, intent(in) :: &
    1394             :       zt_min
    1395             : 
    1396             :     !---------------------------- Local Variables ----------------------------
    1397             :     real( kind = core_rknd ), dimension(1,gr%nz) :: &
    1398           0 :       azm_col    ! Variable on thermodynamic grid levels    [units vary]
    1399             : 
    1400             :     real( kind = core_rknd ), dimension(1,gr%nz) :: &
    1401           0 :       redirect_interpolated_azt_col    ! Variable when interp. to momentum levels
    1402             :     
    1403             :     !---------------------------- Begin Code ----------------------------
    1404             : 
    1405           0 :     azm_col(1,:) = azm
    1406             : 
    1407           0 :     redirect_interpolated_azt_col = redirect_interpolated_azt_2D( gr%nz, 1, gr, azm_col, zt_min )
    1408             : 
    1409           0 :     redirect_interpolated_azt_1D = redirect_interpolated_azt_col(1,:)
    1410             : 
    1411           0 :     return
    1412             :   end function redirect_interpolated_azt_1D
    1413             : 
    1414             :   !=============================================================================
    1415             :   ! Wrapped in interface zm2zt
    1416   424415160 :   function redirect_interpolated_azt_2D( nz, ngrdcol, gr, azm, zt_min )
    1417             : 
    1418             :     ! Description:
    1419             :     ! Calls the appropriate corresponding function based on l_cubic_temp
    1420             :     !-----------------------------------------------------------------------
    1421             : 
    1422             :     use clubb_precision, only: &
    1423             :         core_rknd  ! Variable(s)
    1424             : 
    1425             :     use model_flags, only: &
    1426             :         l_cubic_interp, & ! Variable(s)
    1427             :         l_quintic_poly_interp
    1428             : 
    1429             :     use constants_clubb, only: &
    1430             :         fstdout ! Variable
    1431             : 
    1432             :     implicit none
    1433             :     
    1434             :     !---------------------------- Input Variables ----------------------------
    1435             :     integer, intent(in) :: &
    1436             :       nz, &
    1437             :       ngrdcol
    1438             : 
    1439             :     type (grid), target, intent(in) :: gr
    1440             : 
    1441             :     real( kind = core_rknd ), intent(in), dimension(ngrdcol,nz) :: &
    1442             :       azm    ! Variable on momentum grid levels    [units vary]
    1443             : 
    1444             :     !---------------------------- Return Variable----------------------------
    1445             :     real( kind = core_rknd ), dimension(ngrdcol,nz) :: &
    1446             :       redirect_interpolated_azt_2D    ! Variable when interp. to momentum levels
    1447             : 
    1448             :     !---------------------------- Optional Input Variable ----------------------------
    1449             :     real( kind = core_rknd ), optional, intent(in) :: &
    1450             :       zt_min
    1451             :     
    1452             :     !---------------------------- Begin Code ----------------------------
    1453             : 
    1454             :     ! Sanity Check
    1455             :     if (l_quintic_poly_interp) then
    1456             :       if (.not. l_cubic_interp) then
    1457             :         write (fstdout, *) "Error: Model flag l_quintic_poly_interp should not be true if "&
    1458             :                          //"l_cubic_interp is false."
    1459             :         error stop
    1460             :       end if
    1461             :     end if
    1462             : 
    1463             :     ! Redirect
    1464             :     if (l_cubic_interp .and. nz >= 3 ) then
    1465             :       redirect_interpolated_azt_2D = cubic_interpolated_azt_2D( nz, ngrdcol, gr, azm )
    1466             :     else
    1467             :       call linear_interpolated_azt_2D( nz, ngrdcol, gr, azm, &
    1468             :                                        redirect_interpolated_azt_2D, &
    1469   424415160 :                                        zt_min )
    1470             :     end if
    1471             : 
    1472   424415160 :     return
    1473             :   end function redirect_interpolated_azt_2D
    1474             :   
    1475             :   !=============================================================================
    1476   393142464 :   subroutine linear_interpolated_azm_2D( nz, ngrdcol, gr, azt, &
    1477   393142464 :                                          linear_interpolated_azm, &
    1478             :                                          zm_min )
    1479             : 
    1480             :     ! Description:
    1481             :     ! Function to interpolate a variable located on the thermodynamic grid
    1482             :     ! levels (azt) to the momentum grid levels (azm).  This function inputs the
    1483             :     ! entire azt array and outputs the results as an azm array.  The
    1484             :     ! formulation used is compatible with a stretched (unevenly-spaced) grid.
    1485             :     !-----------------------------------------------------------------------
    1486             : 
    1487             :     use clubb_precision, only: &
    1488             :         core_rknd  ! Variable(s)
    1489             : 
    1490             :     use interpolation, only: &
    1491             :         linear_interp_factor  ! Procedure(s)
    1492             : 
    1493             :     implicit none
    1494             :     
    1495             :     integer, intent(in) :: &
    1496             :       nz, &
    1497             :       ngrdcol
    1498             : 
    1499             :     type (grid), target, intent(in) :: gr
    1500             : 
    1501             :     ! ------------------------------ Input Variable ------------------------------
    1502             :     real( kind = core_rknd ), intent(in), dimension(ngrdcol,nz) :: &
    1503             :       azt    ! Variable on thermodynamic grid levels    [units vary]
    1504             : 
    1505             :     ! ------------------------------ Return Variable ------------------------------
    1506             :     real( kind = core_rknd ), intent(out), dimension(ngrdcol,nz) :: &
    1507             :       linear_interpolated_azm    ! Variable when interp. to momentum levels
    1508             : 
    1509             :     !---------------------------- Optional Input Variable ----------------------------
    1510             :     real( kind = core_rknd ), optional, intent(in) :: &
    1511             :       zm_min
    1512             :       
    1513             :     ! ------------------------------ Local Variable ------------------------------
    1514             :     integer :: i, k  ! Grid level loop index
    1515             : 
    1516             :     ! ------------------------------ Begin Code ------------------------------
    1517             : 
    1518             :     !$acc data copyin( azt, gr, gr%weights_zt2zm, gr%zt, gr%zm ) &
    1519             :     !$acc      copyout( linear_interpolated_azm )
    1520             : 
    1521             :     !$acc parallel loop gang vector default(present)
    1522  6564564864 :     do i = 1, ngrdcol
    1523  6564564864 :       linear_interpolated_azm(i,1) = azt(i,2)
    1524             :     end do
    1525             :     !$acc end parallel loop
    1526             : 
    1527             :     ! Interpolate the value of a thermodynamic-level variable to the central
    1528             :     ! momentum level, k, between two successive thermodynamic levels using
    1529             :     ! linear interpolation.
    1530             :     !$acc parallel loop gang vector collapse(2) default(present)
    1531 33023966976 :     do k = 2, nz-1
    1532 >54525*10^7 :       do i = 1, ngrdcol
    1533 >10244*10^8 :         linear_interpolated_azm(i,k) = gr%weights_zt2zm(i,k,1) &
    1534 >15693*10^8 :                                           * ( azt(i,k+1) - azt(i,k) ) + azt(i,k)
    1535             :       end do
    1536             :     end do
    1537             :     !$acc end parallel loop
    1538             : 
    1539             :     ! Set the value of the thermodynamic-level variable, azt, at the uppermost
    1540             :     ! level of the model, which is a momentum level.  The name of the variable
    1541             :     ! when interpolated/extended to momentum levels is azm.
    1542             :     ! Use a linear extension based on the values of azt at levels gr%nz and
    1543             :     ! gr%nz-1 to find the value of azm at level gr%nz (the uppermost level
    1544             :     ! in the model).
    1545             :     !$acc parallel loop gang vector default(present)
    1546  6564564864 :     do i = 1, ngrdcol
    1547 12342844800 :       linear_interpolated_azm(i,nz) &
    1548 12342844800 :         = ( ( azt(i,nz) - azt(i,nz-1) ) / ( gr%zt(i,nz) - gr%zt(i,nz-1) ) ) & 
    1549 25078832064 :           * ( gr%zm(i,nz) - gr%zt(i,nz) ) + azt(i,nz)
    1550             :     end do
    1551             :     !$acc end parallel loop
    1552             : 
    1553   393142464 :     if ( present(zm_min) ) then
    1554             :       !$acc parallel loop gang vector collapse(2) default(present) copyin(zm_min)
    1555 16905125952 :       do k = 1, nz
    1556 >27919*10^7 :         do i = 1, ngrdcol
    1557 >27899*10^7 :           linear_interpolated_azm(i,k) = max( linear_interpolated_azm(i,k), zm_min )
    1558             :         end do
    1559             :       end do
    1560             :       !$acc end parallel loop
    1561             :     end if
    1562             : 
    1563             :     !$acc end data
    1564             : 
    1565   393142464 :     return
    1566             : 
    1567             :   end subroutine linear_interpolated_azm_2D
    1568             : 
    1569             :   !=============================================================================
    1570           0 :   function zt2zm2zt( nz, ngrdcol, gr, azt, zt_min )
    1571             : 
    1572             :     ! Description:
    1573             :     !    Function to interpolate a variable located on the thermodynamic grid
    1574             :     !    levels (azt) to the momentum grid levels (azm), then interpolate back
    1575             :     !    to thermodynamic grid levels (azt).
    1576             :     !
    1577             :     ! Note:
    1578             :     !   This is intended for smoothing variables.
    1579             :     !-----------------------------------------------------------------------------
    1580             : 
    1581             :     use clubb_precision, only: &
    1582             :         core_rknd  ! Variable(s)
    1583             : 
    1584             :     implicit none
    1585             :     
    1586             :     ! ------------------------------ Input Variable ------------------------------
    1587             :     integer, intent(in) :: &
    1588             :       nz, &
    1589             :       ngrdcol
    1590             : 
    1591             :     type (grid), target, intent(in) :: gr
    1592             : 
    1593             :     real( kind = core_rknd ), intent(in), dimension(ngrdcol,nz) :: &
    1594             :       azt    ! Variable on thermodynamic grid levels    [units vary]
    1595             : 
    1596             :     ! ------------------------------ Return Variable ------------------------------
    1597             :     real( kind = core_rknd ), dimension(ngrdcol,nz) :: &
    1598             :       zt2zm2zt    ! Variable when interp. to momentum levels
    1599             : 
    1600             :     !---------------------------- Optional Input Variable ----------------------------
    1601             :     real( kind = core_rknd ), optional, intent(in) :: &
    1602             :       zt_min
    1603             : 
    1604             :     ! ------------------------------ Local Variable ------------------------------
    1605             :     real( kind = core_rknd ), dimension(ngrdcol,nz) :: &
    1606           0 :       azt_zm
    1607             : 
    1608             :     ! ------------------------------ Begin Code ------------------------------
    1609             : 
    1610             :     !$acc data copyin( azt ) &
    1611             :     !$acc     copyout( zt2zm2zt ) &
    1612             :     !$acc      create( azt_zm )
    1613             : 
    1614             :     ! Interpolate azt to momentum levels 
    1615           0 :     azt_zm = zt2zm( nz, ngrdcol, gr, azt )
    1616             : 
    1617             :     ! Interpolate back to termodynamic levels
    1618           0 :     zt2zm2zt = zm2zt( nz, ngrdcol, gr, azt_zm, zt_min )
    1619             : 
    1620             :     !$acc end data
    1621             : 
    1622           0 :     return 
    1623             : 
    1624             :   end function zt2zm2zt
    1625             :   
    1626             :   !=============================================================================
    1627    26805168 :   function zm2zt2zm( nz, ngrdcol, gr, azm, zm_min )
    1628             : 
    1629             :     ! Description:
    1630             :     !    Function to interpolate a variable located on the momentum grid 
    1631             :     !    levels(azm) to thermodynamic grid levels (azt), then interpolate 
    1632             :     !    back to momentum grid levels (azm).
    1633             :     !
    1634             :     ! Note:
    1635             :     !   This is intended for smoothing variables.
    1636             :     !-----------------------------------------------------------------------------
    1637             : 
    1638             :     use clubb_precision, only: &
    1639             :         core_rknd  ! Variable(s)
    1640             : 
    1641             :     implicit none
    1642             :     
    1643             :     ! ------------------------------ Input Variable ------------------------------
    1644             :     integer, intent(in) :: &
    1645             :       nz, &
    1646             :       ngrdcol
    1647             : 
    1648             :     type (grid), target, intent(in) :: gr
    1649             : 
    1650             :     real( kind = core_rknd ), intent(in), dimension(ngrdcol,nz) :: &
    1651             :       azm    ! Variable on momentum grid levels    [units vary]
    1652             : 
    1653             :     ! ------------------------------ Return Variable ------------------------------
    1654             :     real( kind = core_rknd ), dimension(ngrdcol,nz) :: &
    1655             :       zm2zt2zm    ! Variable when interp. to momentum levels
    1656             : 
    1657             :     !---------------------------- Optional Input Variable ----------------------------
    1658             :     real( kind = core_rknd ), optional, intent(in) :: &
    1659             :       zm_min
    1660             : 
    1661             :     ! ------------------------------ Local Variable ------------------------------
    1662             :     real( kind = core_rknd ), dimension(ngrdcol,nz) :: &
    1663    26805168 :       azm_zt
    1664             : 
    1665             :     ! ------------------------------ Begin Code ------------------------------
    1666             : 
    1667             :     !$acc data copyin( azm ) &
    1668             :     !$acc     copyout( zm2zt2zm ) &
    1669             :     !$acc      create( azm_zt )
    1670             : 
    1671             :     ! Interpolate azt to termodynamic levels 
    1672    26805168 :     azm_zt = zm2zt( nz, ngrdcol, gr, azm )
    1673             : 
    1674             :     ! Interpolate back to momentum levels
    1675    26805168 :     zm2zt2zm = zt2zm( nz, ngrdcol, gr, azm_zt, zm_min )
    1676             : 
    1677             :     !$acc end data
    1678             : 
    1679    26805168 :     return 
    1680             : 
    1681             :   end function zm2zt2zm
    1682             : 
    1683             :   !=============================================================================
    1684             :   function cubic_interpolated_azm_2D( nz, ngrdcol, gr, azt )
    1685             : 
    1686             :     ! Description:
    1687             :     !   Function to interpolate a variable located on the momentum grid
    1688             :     !   levels (azt) to the thermodynamic grid levels (azm).  This function outputs the
    1689             :     !   value of azt at a all grid levels using Steffen's monotonic cubic
    1690             :     !   interpolation implemented by Tak Yamaguchi.
    1691             :     ! 
    1692             :     ! References:
    1693             :     !   None
    1694             :     !-----------------------------------------------------------------------
    1695             : 
    1696             :     use interpolation, only: &
    1697             :         mono_cubic_interp  ! Procedure(s)
    1698             : 
    1699             :     use clubb_precision, only: &
    1700             :         core_rknd ! Variable(s)
    1701             : 
    1702             :     implicit none
    1703             : 
    1704             :     ! Input Variables
    1705             :     integer, intent(in) :: &
    1706             :       nz, &
    1707             :       ngrdcol
    1708             :     
    1709             :     type (grid), target, intent(in) :: gr
    1710             :     
    1711             :     real( kind = core_rknd ), intent(in), dimension(ngrdcol,nz) :: &
    1712             :       azt
    1713             : 
    1714             :     ! Return Variable
    1715             :     real( kind = core_rknd ), dimension(ngrdcol,nz) :: &
    1716             :       cubic_interpolated_azm_2D
    1717             : 
    1718             :     ! Local Variable(s)
    1719             :     integer :: &
    1720             :       k, i
    1721             :       
    1722             :     integer :: km1, k00, kp1, kp2
    1723             : 
    1724             :     ! ---- Begin Code ----
    1725             : 
    1726             :     do k = 1, nz 
    1727             :       do i = 1, ngrdcol
    1728             : 
    1729             :         ! k levels are based on Tak's find_indices subroutine -dschanen 24 Oct 2011
    1730             :         if ( k == nz-1 ) then
    1731             :           km1 = nz-2
    1732             :           kp1 = nz
    1733             :           kp2 = nz
    1734             :           k00 = nz-1
    1735             :         else if ( k == nz ) then ! Extrapolation
    1736             :           km1 = nz
    1737             :           kp1 = nz
    1738             :           kp2 = nz
    1739             :           k00 = nz-1
    1740             :         else if ( k == 1 ) then
    1741             :           km1 = 1
    1742             :           kp1 = 2
    1743             :           kp2 = 3
    1744             :           k00 = 1
    1745             :         else
    1746             :           km1 = k-1
    1747             :           kp1 = k+1
    1748             :           kp2 = k+2
    1749             :           k00 = k
    1750             :         end if
    1751             :         
    1752             :         ! Do the actual interpolation.
    1753             :         ! Use a cubic monotonic spline interpolation.
    1754             :         cubic_interpolated_azm_2D(i,k) = &
    1755             :           mono_cubic_interp( gr%zm(i,k), km1, k00, kp1, kp2, &
    1756             :                              gr%zt(i,km1), gr%zt(i,k00), gr%zt(i,kp1), gr%zt(i,kp2), &
    1757             :                              azt(i,km1), azt(i,k00), azt(i,kp1), azt(i,kp2) )
    1758             :       end do
    1759             :     end do
    1760             : 
    1761             :     return
    1762             : 
    1763             :   end function cubic_interpolated_azm_2D
    1764             : 
    1765             :   !=============================================================================
    1766     4467528 :   subroutine calc_zt2zm_weights( nz, ngrdcol, &
    1767             :                                       gr ) 
    1768             : 
    1769             :     ! Description:
    1770             :     ! Function used to help in an interpolation of a variable (var_zt) located
    1771             :     ! on the thermodynamic grid levels (azt) to the momentum grid levels (azm).
    1772             :     ! This function computes a weighting factor for both the upper thermodynamic
    1773             :     ! level (k+1) and the lower thermodynamic level (k) applied to the central
    1774             :     ! momentum level (k).  For the uppermost momentum grid level (k=gr%nz), a
    1775             :     ! weighting factor for both the thermodynamic level at gr%nz and the
    1776             :     ! thermodynamic level at gr%nz-1 are calculated based on the use of a
    1777             :     ! linear extension.  This function outputs the weighting factors at a single
    1778             :     ! momentum grid level (k).  The formulation used is compatible with a
    1779             :     ! stretched (unevenly-spaced) grid.  The weights are defined as follows:
    1780             :     !
    1781             :     ! ---var_zt(k+1)------------------------------------------- t(k+1)
    1782             :     !                       azt_weight(t_above) = factor
    1783             :     ! ===========var_zt(interp)================================ m(k)
    1784             :     !                       azt_weight(t_below) = 1 - factor
    1785             :     ! ---var_zt(k)--------------------------------------------- t(k)
    1786             :     !
    1787             :     ! The vertical indices t(k+1), m(k), and t(k) correspond with altitudes
    1788             :     ! zt(k+1), zm(k), and zt(k), respectively.  The letter "t" is used for
    1789             :     ! thermodynamic levels and the letter "m" is used for momentum levels.
    1790             :     !
    1791             :     ! For all levels k < gr%nz:
    1792             :     !
    1793             :     ! The formula for a linear interpolation is given by:
    1794             :     !
    1795             :     ! var_zt( interp to zm(k) )
    1796             :     ! = [ ( var_zt(k+1) - var_zt(k) ) / ( zt(k+1) - zt(k) ) ]
    1797             :     !     * ( zm(k) - zt(k) ) + var_zt(k);
    1798             :     !
    1799             :     ! which can be rewritten as:
    1800             :     !
    1801             :     ! var_zt( interp to zm(k) )
    1802             :     ! = [ ( zm(k) - zt(k) ) / ( zt(k+1) - zt(k) ) ]
    1803             :     !     * ( var_zt(k+1) - var_zt(k) ) + var_zt(k).
    1804             :     !
    1805             :     ! Furthermore, the formula can be rewritten as:
    1806             :     !
    1807             :     ! var_zt( interp to zm(k) )
    1808             :     ! = factor * var_zt(k+1) + ( 1 - factor ) * var_zt(k);
    1809             :     !
    1810             :     ! where:
    1811             :     !
    1812             :     ! factor = ( zm(k) - zt(k) ) / ( zt(k+1) - zt(k) ).
    1813             :     !
    1814             :     ! One of the important uses of this function is in situations where the
    1815             :     ! variable to be interpolated is being treated IMPLICITLY in an equation.
    1816             :     ! Usually, the variable to be interpolated is involved in a derivative (such
    1817             :     ! as d(var_zt)/dz in the diagram below).  For the term of the equation
    1818             :     ! containing the derivative, grid weights are needed for two interpolations,
    1819             :     ! rather than just one interpolation.  Thus, four grid weights (labeled
    1820             :     ! A(k), B(k), C(k), and D(k) in the diagram below) are needed.
    1821             :     !
    1822             :     ! ---var_zt(k+1)------------------------------------------- t(k+1)
    1823             :     !                                       A(k)
    1824             :     ! ===========var_zt(interp)================================ m(k)
    1825             :     !                                       B(k) = 1 - A(k)
    1826             :     ! ---var_zt(k)-----------d(var_zt)/dz---------------------- t(k)
    1827             :     !                                       C(k)
    1828             :     ! ===========var_zt(interp)================================ m(k-1)
    1829             :     !                                       D(k) = 1 - C(k)
    1830             :     ! ---var_zt(k-1)------------------------------------------- t(k-1)
    1831             :     !
    1832             :     ! The vertical indices t(k+1), m(k), t(k), m(k-1), and t(k-1) correspond
    1833             :     ! with altitudes zt(k+1), zm(k), zt(k), zm(k-1), and zt(k-1), respectively.
    1834             :     ! The letter "t" is used for thermodynamic levels and the letter "m" is used
    1835             :     ! for momentum levels.
    1836             :     !
    1837             :     ! The grid weights, indexed around the central thermodynamic level (k), are
    1838             :     ! defined as follows:
    1839             :     !
    1840             :     ! A(k) = ( zm(k) - zt(k) ) / ( zt(k+1) - zt(k) );
    1841             :     !
    1842             :     ! which is the same as "factor" for the interpolation to momentum
    1843             :     ! level (k).  In the code, this interpolation is referenced as
    1844             :     ! gr%weights_zt2zm(t_above,mk), which can be read as "grid weight in a zt2zm
    1845             :     ! interpolation of the thermodynamic level above momentum level (k) (applied
    1846             :     ! to momentum level (k))".
    1847             :     !
    1848             :     ! B(k) = 1 - [ ( zm(k) - zt(k) ) / ( zt(k+1) - zt(k) ) ]
    1849             :     !      = 1 - A(k);
    1850             :     !
    1851             :     ! which is the same as "1 - factor" for the interpolation to momentum
    1852             :     ! level (k).  In the code, this interpolation is referenced as
    1853             :     ! gr%weights_zt2zm(t_below,mk), which can be read as "grid weight in a zt2zm
    1854             :     ! interpolation of the thermodynamic level below momentum level (k) (applied
    1855             :     ! to momentum level (k))".
    1856             :     !
    1857             :     ! C(k) = ( zm(k-1) - zt(k-1) ) / ( zt(k) - zt(k-1) );
    1858             :     !
    1859             :     ! which is the same as "factor" for the interpolation to momentum
    1860             :     ! level (k-1).  In the code, this interpolation is referenced as
    1861             :     ! gr%weights_zt2zm(t_above,mkm1), which can be read as "grid weight in a
    1862             :     ! zt2zm interpolation of the thermodynamic level above momentum level (k-1)
    1863             :     ! (applied to momentum level (k-1))".
    1864             :     !
    1865             :     ! D(k) = 1 - [ ( zm(k-1) - zt(k-1) ) / ( zt(k) - zt(k-1) ) ]
    1866             :     !      = 1 - C(k);
    1867             :     !
    1868             :     ! which is the same as "1 - factor" for the interpolation to momentum
    1869             :     ! level (k-1).  In the code, this interpolation is referenced as
    1870             :     ! gr%weights_zt2zm(t_below,mkm1), which can be read as "grid weight in a
    1871             :     ! zt2zm interpolation of the thermodynamic level below momentum level (k-1)
    1872             :     ! (applied to momentum level (k-1))".
    1873             :     !
    1874             :     ! Additionally, as long as the central thermodynamic level (k) in the above
    1875             :     ! scenario is not the uppermost thermodynamic level or the lowermost
    1876             :     ! thermodynamic level (k /= gr%nz and k /= 1), the four weighting factors
    1877             :     ! have the following relationships:  A(k) = C(k+1) and B(k) = D(k+1).
    1878             :     !
    1879             :     !
    1880             :     ! Special condition for uppermost grid level, k = gr%nz:
    1881             :     !
    1882             :     ! The uppermost momentum grid level is above the uppermost thermodynamic
    1883             :     ! grid level.  Thus, a linear extension is used at this level.
    1884             :     !
    1885             :     ! For level k = gr%nz:
    1886             :     !
    1887             :     ! The formula for a linear extension is given by:
    1888             :     !
    1889             :     ! var_zt( extend to zm(k) )
    1890             :     ! = [ ( var_zt(k) - var_zt(k-1) ) / ( zt(k) - zt(k-1) ) ]
    1891             :     !     * ( zm(k) - zt(k-1) ) + var_zt(k-1);
    1892             :     !
    1893             :     ! which can be rewritten as:
    1894             :     !
    1895             :     ! var_zt( extend to zm(k) )
    1896             :     ! = [ ( zm(k) - zt(k-1) ) / ( zt(k) - zt(k-1) ) ]
    1897             :     !     * ( var_zt(k) - var_zt(k-1) ) + var_zt(k-1).
    1898             :     !
    1899             :     ! Furthermore, the formula can be rewritten as:
    1900             :     !
    1901             :     ! var_zt( extend to zm(k) )
    1902             :     ! = factor * var_zt(k) + ( 1 - factor ) * var_zt(k-1);
    1903             :     !
    1904             :     ! where:
    1905             :     !
    1906             :     ! factor = ( zm(k) - zt(k-1) ) / ( zt(k) - zt(k-1) ).
    1907             :     !
    1908             :     ! Due to the fact that a linear extension is being used, the value of factor
    1909             :     ! will be greater than 1.  The weight of thermodynamic level k = gr%nz on
    1910             :     ! momentum level k = gr%nz equals the value of factor.  The weight of
    1911             :     ! thermodynamic level k = gr%nz-1 on momentum level k = gr%nz equals
    1912             :     ! 1 - factor, which is less than 0.  However, the sum of the two weights
    1913             :     ! equals 1.
    1914             :     !
    1915             :     !
    1916             :     ! Brian Griffin; September 12, 2008.
    1917             :     !
    1918             :     !-----------------------------------------------------------------------
    1919             : 
    1920             :     use constants_clubb, only: &
    1921             :         one  ! Constant(s)
    1922             : 
    1923             :     use clubb_precision, only: &
    1924             :         core_rknd ! Variable(s)
    1925             : 
    1926             :     implicit none
    1927             : 
    1928             :     ! Constant parameters
    1929             :     integer, parameter :: & 
    1930             :       t_above = 1,  & ! Upper thermodynamic level.
    1931             :       t_below = 2     ! Lower thermodynamic level.
    1932             : 
    1933             :     ! Input Variable
    1934             :     integer, intent(in) :: &
    1935             :       nz, &
    1936             :       ngrdcol
    1937             :     
    1938             :     ! Input/Output Variable
    1939             :     type (grid), target, intent(inout) :: gr
    1940             : 
    1941             :     ! Local Variables
    1942             :     real( kind = core_rknd ) :: factor
    1943             : 
    1944             :     integer :: i, k
    1945             : 
    1946   379739880 :     do k = 1, nz-1
    1947  6270643080 :       do i = 1, ngrdcol
    1948             :         
    1949             :         ! The top model level (gr%nz) is formulated differently because the top
    1950             :         ! momentum level is above the top thermodynamic level.  A linear
    1951             :         ! extension is required, rather than linear interpolation.
    1952             :         ! Note:  Variable "factor" will be greater than 1 in this situation.
    1953           0 :         gr%weights_zt2zm(i,k,t_above) = ( gr%zm(i,k) - gr%zt(i,k) ) &
    1954  5890903200 :                                           / ( gr%zt(i,k+1) - gr%zt(i,k) )
    1955             : 
    1956             :         ! Weight of lower thermodynamic level on momentum level.
    1957  6266175552 :         gr%weights_zt2zm(i,k,t_below) = one - gr%weights_zt2zm(i,k,t_above)
    1958             :         
    1959             :       end do
    1960             :     end do
    1961             :     
    1962             :     ! At most levels, the momentum level is found in-between two
    1963             :     ! thermodynamic levels.  Linear interpolation is used.
    1964    74597328 :     do i = 1, ngrdcol
    1965           0 :       gr%weights_zt2zm(i,nz,t_above) = ( gr%zm(i,nz) - gr%zt(i,nz-1) ) &
    1966    70129800 :                                         / ( gr%zt(i,nz) - gr%zt(i,nz-1) )
    1967             :                                         
    1968    74597328 :       gr%weights_zt2zm(i,nz,t_below) = one - gr%weights_zt2zm(i,nz,t_above)
    1969             :     end do
    1970             : 
    1971     4467528 :     return
    1972             : 
    1973             :   end subroutine calc_zt2zm_weights
    1974             :   
    1975             :   !=============================================================================
    1976   424415160 :   subroutine linear_interpolated_azt_2D( nz, ngrdcol, gr, azm, &
    1977   424415160 :                                          linear_interpolated_azt, &
    1978             :                                          zt_min )
    1979             : 
    1980             :     ! Description:
    1981             :     ! Function to interpolate a variable located on the momentum grid levels
    1982             :     ! (azm) to the thermodynamic grid levels (azt).  This function inputs the
    1983             :     ! entire azm array and outputs the results as an azt array.  The formulation
    1984             :     ! used is compatible with a stretched (unevenly-spaced) grid.
    1985             :     !-----------------------------------------------------------------------
    1986             : 
    1987             :     use clubb_precision, only: &
    1988             :         core_rknd  ! Variable(s)
    1989             : 
    1990             :     use interpolation, only: &
    1991             :         linear_interp_factor  ! Procedure(s)
    1992             : 
    1993             :     implicit none
    1994             :     
    1995             :     integer, intent(in) :: &
    1996             :       nz, &
    1997             :       ngrdcol
    1998             : 
    1999             :     type (grid), target, intent(in) :: gr
    2000             : 
    2001             :     ! ------------------------------ Input Variable ------------------------------
    2002             :     real( kind = core_rknd ), intent(in), dimension(ngrdcol,nz) :: &
    2003             :       azm    ! Variable on momentum grid levels    [units vary]
    2004             : 
    2005             :     ! ------------------------------ Output Variable ------------------------------
    2006             :     real( kind = core_rknd ), intent(out), dimension(ngrdcol,nz) :: &
    2007             :       linear_interpolated_azt    ! Variable when interp. to thermodynamic levels
    2008             : 
    2009             :     !---------------------------- Optional Input Variable ----------------------------
    2010             :     real( kind = core_rknd ), optional, intent(in) :: &
    2011             :       zt_min
    2012             : 
    2013             :     ! ------------------------------ Local Variable ------------------------------
    2014             :     integer :: i, k  ! Grid level loop index
    2015             : 
    2016             :     ! ------------------------------ Begin Code ------------------------------
    2017             : 
    2018             :     !$acc data copyin( azm, gr, gr%weights_zm2zt, gr%zt, gr%zm ) &
    2019             :     !$acc      copyout( linear_interpolated_azt )
    2020             : 
    2021             :     ! Set the value of the momentum-level variable, azm, at the lowermost level
    2022             :     ! of the model (below the model lower boundary), which is a thermodynamic
    2023             :     ! level.  The name of the variable when interpolated/extended to
    2024             :     ! thermodynamic levels is azt.
    2025             :     ! Use a linear extension based on the values of azm at levels 1 and 2 to
    2026             :     ! find the value of azt at level 1 (the lowermost level in the model).
    2027             :     !$acc parallel loop gang vector default(present)
    2028  7086746160 :     do i = 1, ngrdcol
    2029 13324662000 :       linear_interpolated_azt(i,1) &
    2030 13324662000 :         = ( ( azm(i,2) - azm(i,1) ) / ( gr%zm(i,2) - gr%zm(i,1) ) ) & 
    2031 20411408160 :           * ( gr%zt(i,1) - gr%zm(i,1) ) + azm(i,1)
    2032             :     end do
    2033             :     !$acc end parallel loop
    2034             : 
    2035             :     ! Interpolate the value of a momentum-level variable to the central
    2036             :     ! thermodynamic level, k, between two successive momentum levels using
    2037             :     ! linear interpolation.
    2038             :     !$acc parallel loop gang vector collapse(2) default(present)
    2039 36075288600 :     do k = 2, nz
    2040 >59571*10^7 :       do i = 1, ngrdcol
    2041 >11192*10^8 :         linear_interpolated_azt(i,k) = gr%weights_zm2zt(i,k,1) &
    2042 >17145*10^8 :                                        * ( azm(i,k) - azm(i,k-1) ) + azm(i,k-1)
    2043             :       end do
    2044             :     end do
    2045             :     !$acc end parallel loop
    2046             : 
    2047   424415160 :     if ( present(zt_min) ) then
    2048             :       !$acc parallel loop gang vector collapse(2) default(present) copyin(zt_min)
    2049 19210370400 :       do k = 1, nz
    2050 >31726*10^7 :         do i = 1, ngrdcol
    2051 >31703*10^7 :           linear_interpolated_azt(i,k) = max( linear_interpolated_azt(i,k), zt_min )
    2052             :         end do
    2053             :       end do
    2054             :       !$acc end parallel loop
    2055             :     end if
    2056             : 
    2057             :     !$acc end data
    2058             : 
    2059   424415160 :     return
    2060             : 
    2061             :   end subroutine linear_interpolated_azt_2D
    2062             :   
    2063             :   !=============================================================================
    2064             :   function cubic_interpolated_azt_2D( nz, ngrdcol, gr, azm )
    2065             : 
    2066             :     ! Description:
    2067             :     !   Function to interpolate a variable located on the momentum grid
    2068             :     !   levels (azm) to the thermodynamic grid levels (azt).  This function outputs the
    2069             :     !   value of azt at a all grid levels using Steffen's monotonic cubic
    2070             :     !   interpolation implemented by Tak Yamaguchi.
    2071             :     ! 
    2072             :     ! References:
    2073             :     !   None
    2074             :     !-----------------------------------------------------------------------
    2075             : 
    2076             :     use interpolation, only: &
    2077             :         mono_cubic_interp  ! Procedure(s)
    2078             : 
    2079             :     use clubb_precision, only: &
    2080             :         core_rknd ! Variable(s)
    2081             : 
    2082             :     implicit none
    2083             : 
    2084             :     ! Input Variables
    2085             :     integer, intent(in) :: &
    2086             :       nz, &
    2087             :       ngrdcol
    2088             :     
    2089             :     type (grid), target, intent(in) :: gr
    2090             :     
    2091             :     real( kind = core_rknd ), intent(in), dimension(ngrdcol,nz) :: &
    2092             :       azm
    2093             : 
    2094             :     ! Return Variable
    2095             :     real( kind = core_rknd ), dimension(ngrdcol,nz) :: &
    2096             :       cubic_interpolated_azt_2D
    2097             : 
    2098             :     ! Local Variable(s)
    2099             :     integer :: &
    2100             :       k, i
    2101             :       
    2102             :     integer :: km1, k00, kp1, kp2
    2103             : 
    2104             :     ! ---- Begin Code ----
    2105             : 
    2106             :     do k = 1, nz 
    2107             :       do i = 1, ngrdcol
    2108             : 
    2109             :         ! k levels are based on Tak's find_indices subroutine -dschanen 24 Oct 2011
    2110             :         if ( k == nz ) then
    2111             :           km1 = nz-2
    2112             :           kp1 = nz
    2113             :           kp2 = nz
    2114             :           k00 = nz-1
    2115             :         else if ( k == 2 ) then
    2116             :           km1 = 1
    2117             :           kp1 = 2
    2118             :           kp2 = 3
    2119             :           k00 = 1
    2120             :         else if ( k == 1 ) then ! Extrapolation for the ghost point
    2121             :           km1 = nz
    2122             :           k00 = 1
    2123             :           kp1 = 2
    2124             :           kp2 = 3
    2125             :         else
    2126             :           km1 = k-2
    2127             :           kp1 = k
    2128             :           kp2 = k+1
    2129             :           k00 = k-1
    2130             :         end if
    2131             :         
    2132             :         ! Do the actual interpolation.
    2133             :         ! Use a cubic monotonic spline interpolation.
    2134             :         cubic_interpolated_azt_2D(i,k) = &
    2135             :           mono_cubic_interp( gr%zt(i,k), km1, k00, kp1, kp2, &
    2136             :                              gr%zm(i,km1), gr%zm(i,k00), gr%zm(i,kp1), gr%zm(i,kp2), &
    2137             :                              azm(i,km1), azm(i,k00), azm(i,kp1), azm(i,kp2) )
    2138             :       end do
    2139             :     end do
    2140             : 
    2141             :     return
    2142             : 
    2143             :   end function cubic_interpolated_azt_2D
    2144             : 
    2145             :   !=============================================================================
    2146     4467528 :   subroutine calc_zm2zt_weights( nz, ngrdcol, &
    2147             :                                       gr )
    2148             : 
    2149             :     ! Description:
    2150             :     ! Function used to help in an interpolation of a variable (var_zm) located
    2151             :     ! on the momentum grid levels (azm) to the thermodynamic grid levels (azt).
    2152             :     ! This function computes a weighting factor for both the upper momentum
    2153             :     ! level (k) and the lower momentum level (k-1) applied to the central
    2154             :     ! thermodynamic level (k).  For the lowermost thermodynamic grid
    2155             :     ! level (k=1), a weighting factor for both the momentum level at 1 and the
    2156             :     ! momentum level at 2 are calculated based on the use of a linear extension.
    2157             :     ! This function outputs the weighting factors at a single thermodynamic grid
    2158             :     ! level (k).   The formulation used is compatible with a stretched
    2159             :     ! (unevenly-spaced) grid.  The weights are defined as follows:
    2160             :     !
    2161             :     ! ===var_zm(k)============================================= m(k)
    2162             :     !                       azm_weight(m_above) = factor
    2163             :     ! -----------var_zm(interp)-------------------------------- t(k)
    2164             :     !                       azm_weight(m_below) = 1 - factor
    2165             :     ! ===var_zm(k-1)=========================================== m(k-1)
    2166             :     !
    2167             :     ! The vertical indices m(k), t(k), and m(k-1) correspond with altitudes
    2168             :     ! zm(k), zt(k), and zm(k-1), respectively.  The letter "t" is used for
    2169             :     ! thermodynamic levels and the letter "m" is used for momentum levels.
    2170             :     !
    2171             :     ! For all levels k > 1:
    2172             :     !
    2173             :     ! The formula for a linear interpolation is given by:
    2174             :     !
    2175             :     ! var_zm( interp to zt(k) )
    2176             :     ! = [ ( var_zm(k) - var_zm(k-1) ) / ( zm(k) - zm(k-1) ) ]
    2177             :     !     * ( zt(k) - zm(k-1) ) + var_zm(k-1);
    2178             :     !
    2179             :     ! which can be rewritten as:
    2180             :     !
    2181             :     ! var_zm( interp to zt(k) )
    2182             :     ! = [ ( zt(k) - zm(k-1) ) / ( zm(k) - zm(k-1) ) ]
    2183             :     !     * ( var_zm(k) - var_zm(k-1) ) + var_zm(k-1).
    2184             :     !
    2185             :     ! Furthermore, the formula can be rewritten as:
    2186             :     !
    2187             :     ! var_zm( interp to zt(k) )
    2188             :     ! = factor * var_zm(k) + ( 1 - factor ) * var_zm(k-1);
    2189             :     !
    2190             :     ! where:
    2191             :     !
    2192             :     ! factor = ( zt(k) - zm(k-1) ) / ( zm(k) - zm(k-1) ).
    2193             :     !
    2194             :     ! One of the important uses of this function is in situations where the
    2195             :     ! variable to be interpolated is being treated IMPLICITLY in an equation.
    2196             :     ! Usually, the variable to be interpolated is involved in a derivative (such
    2197             :     ! as d(var_zm)/dz in the diagram below).  For the term of the equation
    2198             :     ! containing the derivative, grid weights are needed for two interpolations,
    2199             :     ! rather than just one interpolation.  Thus, four grid weights (labeled
    2200             :     ! A(k), B(k), C(k), and D(k) in the diagram below) are needed.
    2201             :     !
    2202             :     ! ===var_zm(k+1)=========================================== m(k+1)
    2203             :     !                                       A(k)
    2204             :     ! -----------var_zm(interp)-------------------------------- t(k+1)
    2205             :     !                                       B(k) = 1 - A(k)
    2206             :     ! ===var_zm(k)===========d(var_zm)/dz====================== m(k)
    2207             :     !                                       C(k)
    2208             :     ! -----------var_zm(interp)-------------------------------- t(k)
    2209             :     !                                       D(k) = 1 - C(k)
    2210             :     ! ===var_zm(k-1)=========================================== m(k-1)
    2211             :     !
    2212             :     ! The vertical indices m(k+1), t(k+1), m(k), t(k), and m(k-1) correspond
    2213             :     ! with altitudes zm(k+1), zt(k+1), zm(k), zt(k), and zm(k-1), respectively.
    2214             :     ! The letter "t" is used for thermodynamic levels and the letter "m" is used
    2215             :     ! for momentum levels.
    2216             :     !
    2217             :     ! The grid weights, indexed around the central momentum level (k), are
    2218             :     ! defined as follows:
    2219             :     !
    2220             :     ! A(k) = ( zt(k+1) - zm(k) ) / ( zm(k+1) - zm(k) );
    2221             :     !
    2222             :     ! which is the same as "factor" for the interpolation to thermodynamic
    2223             :     ! level (k+1).  In the code, this interpolation is referenced as
    2224             :     ! gr%weights_zm2zt(m_above,tkp1), which can be read as "grid weight in a
    2225             :     ! zm2zt interpolation of the momentum level above thermodynamic
    2226             :     ! level (k+1) (applied to thermodynamic level (k+1))".
    2227             :     !
    2228             :     ! B(k) = 1 - [ ( zt(k+1) - zm(k) ) / ( zm(k+1) - zm(k) ) ]
    2229             :     !      = 1 - A(k);
    2230             :     !
    2231             :     ! which is the same as "1 - factor" for the interpolation to thermodynamic
    2232             :     ! level (k+1).  In the code, this interpolation is referenced as
    2233             :     ! gr%weights_zm2zt(m_below,tkp1), which can be read as "grid weight in a
    2234             :     ! zm2zt interpolation of the momentum level below thermodynamic
    2235             :     ! level (k+1) (applied to thermodynamic level (k+1))".
    2236             :     !
    2237             :     ! C(k) = ( zt(k) - zm(k-1) ) / ( zm(k) - zm(k-1) );
    2238             :     !
    2239             :     ! which is the same as "factor" for the interpolation to thermodynamic
    2240             :     ! level (k).  In the code, this interpolation is referenced as
    2241             :     ! gr%weights_zm2zt(m_above,tk), which can be read as "grid weight in a zm2zt
    2242             :     ! interpolation of the momentum level above thermodynamic level (k) (applied
    2243             :     ! to thermodynamic level (k))".
    2244             :     !
    2245             :     ! D(k) = 1 - [ ( zt(k) - zm(k-1) ) / ( zm(k) - zm(k-1) ) ]
    2246             :     !      = 1 - C(k);
    2247             :     !
    2248             :     ! which is the same as "1 - factor" for the interpolation to thermodynamic
    2249             :     ! level (k).  In the code, this interpolation is referenced as
    2250             :     ! gr%weights_zm2zt(m_below,tk), which can be read as "grid weight in a zm2zt
    2251             :     ! interpolation of the momentum level below thermodynamic level (k) (applied
    2252             :     ! to thermodynamic level (k))".
    2253             :     !
    2254             :     ! Additionally, as long as the central momentum level (k) in the above
    2255             :     ! scenario is not the lowermost momentum level or the uppermost momentum
    2256             :     ! level (k /= 1 and k /= gr%nz), the four weighting factors have the
    2257             :     ! following relationships:  A(k) = C(k+1) and B(k) = D(k+1).
    2258             :     !
    2259             :     !
    2260             :     ! Special condition for lowermost grid level, k = 1:
    2261             :     !
    2262             :     ! The lowermost thermodynamic grid level is below the lowermost momentum
    2263             :     ! grid level.  Thus, a linear extension is used at this level.  It should
    2264             :     ! be noted that the thermodynamic level k = 1 is considered to be below the
    2265             :     ! model lower boundary, which is defined to be at momentum level k = 1.
    2266             :     ! Thus, the values of most variables at thermodynamic level k = 1 are not
    2267             :     ! often needed or referenced.
    2268             :     !
    2269             :     ! For level k = 1:
    2270             :     !
    2271             :     ! The formula for a linear extension is given by:
    2272             :     !
    2273             :     ! var_zm( extend to zt(k) )
    2274             :     ! = [ ( var_zm(k+1) - var_zm(k) ) / ( zm(k+1) - zm(k) ) ]
    2275             :     !     * ( zt(k) - zm(k) ) + var_zm(k);
    2276             :     !
    2277             :     ! which can be rewritten as:
    2278             :     !
    2279             :     ! var_zm( extend to zt(k) )
    2280             :     ! = [ ( zt(k) - zm(k) ) / ( zm(k+1) - zm(k) ) ]
    2281             :     !     * ( var_zm(k+1) - var_zm(k) ) + var_zm(k).
    2282             :     !
    2283             :     ! Furthermore, the formula can be rewritten as:
    2284             :     !
    2285             :     ! var_zm( extend to zt(k) )
    2286             :     ! = factor * var_zm(k+1) + ( 1 - factor ) * var_zm(k);
    2287             :     !
    2288             :     ! where:
    2289             :     !
    2290             :     ! factor = ( zt(k) - zm(k) ) / ( zm(k+1) - zm(k) ).
    2291             :     !
    2292             :     ! Due to the fact that a linear extension is being used, the value of factor
    2293             :     ! will be less than 0.  The weight of the upper momentum level, which is
    2294             :     ! momentum level k = 2, on thermodynamic level k = 1 equals the value of
    2295             :     ! factor.  The weight of the lower momentum level, which is momentum level
    2296             :     ! k = 1, on thermodynamic level k = 1 equals 1 - factor, which is greater
    2297             :     ! than 1.  However, the sum of the weights equals 1.
    2298             :     !
    2299             :     !
    2300             :     ! Brian Griffin; September 12, 2008.
    2301             :     !
    2302             :     !-----------------------------------------------------------------------
    2303             : 
    2304             :     use constants_clubb, only: &
    2305             :         one  ! Constant(s)
    2306             : 
    2307             :     use clubb_precision, only: &
    2308             :         core_rknd  ! Variable(s)
    2309             : 
    2310             :     implicit none
    2311             : 
    2312             :     ! Constant parameters
    2313             :     integer, parameter :: & 
    2314             :       m_above = 1,  & ! Upper momentum level.
    2315             :       m_below = 2     ! Lower momentum level.
    2316             : 
    2317             :     ! Input Variable
    2318             :     integer, intent(in) :: &
    2319             :       nz, &
    2320             :       ngrdcol
    2321             :     
    2322             :     ! Input/Output Variable
    2323             :     type (grid), target, intent(inout) :: gr
    2324             : 
    2325             :     integer :: i, k
    2326             : 
    2327             :     
    2328             :     ! At most levels, the momentum level is found in-between two
    2329             :     ! thermodynamic levels.  Linear interpolation is used.
    2330    74597328 :     do i = 1, ngrdcol
    2331           0 :       gr%weights_zm2zt(i,1,m_above) = ( gr%zt(i,1) - gr%zm(i,1) ) &
    2332    70129800 :                                         / ( gr%zm(i,2) - gr%zm(i,1) )
    2333             :                                         
    2334    74597328 :       gr%weights_zm2zt(i,1,m_below) = one - gr%weights_zm2zt(i,1,m_above)
    2335             :     end do
    2336             : 
    2337   379739880 :     do k = 2, nz
    2338  6270643080 :       do i = 1, ngrdcol
    2339             :         
    2340             :         ! The top model level (gr%nz) is formulated differently because the top
    2341             :         ! momentum level is above the top thermodynamic level.  A linear
    2342             :         ! extension is required, rather than linear interpolation.
    2343             :         ! Note:  Variable "factor" will be greater than 1 in this situation.
    2344           0 :         gr%weights_zm2zt(i,k,m_above) = ( gr%zt(i,k) - gr%zm(i,k-1) ) &
    2345  5890903200 :                                          / ( gr%zm(i,k) - gr%zm(i,k-1) )
    2346             : 
    2347             :         ! Weight of lower thermodynamic level on momentum level.
    2348  6266175552 :         gr%weights_zm2zt(i,k,m_below) = one - gr%weights_zm2zt(i,k,m_above)
    2349             :         
    2350             :       end do
    2351             :     end do
    2352             :     
    2353     4467528 :     return
    2354             : 
    2355             :   end subroutine calc_zm2zt_weights
    2356             :   
    2357             :   !=============================================================================
    2358             :   ! Wrapped in interface ddzm
    2359           0 :   function gradzm_2D( nz, ngrdcol, gr, azm )
    2360             : 
    2361             :     ! Description:
    2362             :     !  2D version of gradzm
    2363             :     !-----------------------------------------------------------------------
    2364             : 
    2365             :     use clubb_precision, only: &
    2366             :         core_rknd ! Variable(s)
    2367             : 
    2368             :     implicit none
    2369             :     
    2370             :     integer, intent(in) :: &
    2371             :       nz, &
    2372             :       ngrdcol
    2373             : 
    2374             :     type (grid), target, intent(in) :: gr
    2375             : 
    2376             :     ! Input Variable
    2377             :     real( kind = core_rknd ), intent(in), dimension(ngrdcol,nz) :: &
    2378             :       azm    ! Variable on momentum grid levels    [units vary]
    2379             : 
    2380             :     ! Return Variable
    2381             :     real( kind = core_rknd ), dimension(ngrdcol,nz) :: &
    2382             :       gradzm_2D    ! Vertical derivative of azm    [units vary / m]
    2383             : 
    2384             :     ! Local Variable
    2385             :     integer :: i, k  ! Grid level loop index
    2386             :     
    2387             :     !$acc data copyin( gr, gr%invrs_dzt, azm ) &
    2388             :     !$acc     copyout( gradzm_2D )
    2389             : 
    2390             :     !$acc parallel loop gang vector default(present)
    2391           0 :     do i = 1, ngrdcol
    2392           0 :       gradzm_2D(i,1) = ( azm(i,2) - azm(i,1) ) * gr%invrs_dzt(i,2)
    2393             :     end do
    2394             :     !$acc end parallel loop
    2395             : 
    2396             :     !$acc parallel loop gang vector collapse(2) default(present)
    2397           0 :     do k = 2, nz
    2398           0 :       do i = 1, ngrdcol
    2399           0 :         gradzm_2D(i,k) = ( azm(i,k) - azm(i,k-1) ) * gr%invrs_dzt(i,k)
    2400             :       end do
    2401             :     end do
    2402             :     !$acc end parallel loop
    2403             : 
    2404             :     !$acc end data
    2405             : 
    2406           0 :     return
    2407             : 
    2408             :   end function gradzm_2D
    2409             :   
    2410             :   !=============================================================================
    2411             :   ! Wrapped in interface ddzm
    2412           0 :   function gradzm_1D( gr, azm )
    2413             : 
    2414             :     ! Description:
    2415             :     !  2D version of gradzm
    2416             :     !-----------------------------------------------------------------------
    2417             : 
    2418             :     use clubb_precision, only: &
    2419             :         core_rknd ! Variable(s)
    2420             : 
    2421             :     implicit none
    2422             : 
    2423             :     type (grid), target, intent(in) :: gr
    2424             : 
    2425             :     ! Input Variable
    2426             :     real( kind = core_rknd ), intent(in), dimension(gr%nz) :: &
    2427             :       azm    ! Variable on momentum grid levels    [units vary]
    2428             : 
    2429             :     ! Return Variable
    2430             :     real( kind = core_rknd ), dimension(gr%nz) :: &
    2431             :       gradzm_1D    ! Vertical derivative of azm    [units vary / m]
    2432             :       
    2433             :     ! Local Variables
    2434             :     real( kind = core_rknd ), dimension(1,gr%nz) :: &
    2435           0 :       azm_col    ! Variable on momentum grid levels    [units vary]
    2436             : 
    2437             :     ! Return Variable
    2438             :     real( kind = core_rknd ), dimension(1,gr%nz) :: &
    2439           0 :       gradzm_1D_col    ! Vertical derivative of azm    [units vary / m]
    2440             : 
    2441           0 :     azm_col(1,:) = azm
    2442             :     
    2443           0 :     gradzm_1D_col = gradzm_2D(gr%nz, 1, gr, azm_col)
    2444             :       
    2445           0 :     gradzm_1D = gradzm_1D_col(1,:)
    2446             : 
    2447           0 :     return
    2448             : 
    2449             :   end function gradzm_1D
    2450             :   
    2451             :   !=============================================================================
    2452             :   ! Wrapped in interface ddzt
    2453   250181568 :   function gradzt_2D( nz, ngrdcol, gr, azt )
    2454             : 
    2455             :     ! Description:
    2456             :     !  2D version of gradzt
    2457             :     !-----------------------------------------------------------------------
    2458             : 
    2459             :     use clubb_precision, only: &
    2460             :         core_rknd ! Variable(s)
    2461             : 
    2462             :     implicit none
    2463             :     
    2464             :     integer, intent(in) :: &
    2465             :       nz, &
    2466             :       ngrdcol
    2467             : 
    2468             :     type (grid), target, intent(in) :: gr
    2469             : 
    2470             :     ! Input Variable
    2471             :     real( kind = core_rknd ), intent(in), dimension(ngrdcol,nz) :: &
    2472             :       azt    ! Variable on thermodynamic grid levels    [units vary]
    2473             : 
    2474             :     ! Output Variable
    2475             :     real( kind = core_rknd ), dimension(ngrdcol,nz) :: &
    2476             :       gradzt_2D    ! Vertical derivative of azt    [units vary / m]
    2477             : 
    2478             :     ! Local Variable
    2479             :     integer :: i, k  ! Grid level loop index
    2480             : 
    2481             :     !$acc data copyin( gr, gr%invrs_dzm, azt ) &
    2482             :     !$acc     copyout( gradzt_2D )
    2483             : 
    2484             :     !$acc parallel loop gang vector collapse(2) default(present)
    2485 21015251712 :     do k = 2, nz-1
    2486 >34697*10^7 :       do i = 1, ngrdcol
    2487 >34672*10^7 :         gradzt_2D(i,k) = ( azt(i,k+1) - azt(i,k) ) * gr%invrs_dzm(i,k)
    2488             :       end do
    2489             :     end do
    2490             :     !$acc end parallel loop
    2491             : 
    2492             :     !$acc parallel loop gang vector default(present)
    2493  4177450368 :     do i = 1, ngrdcol
    2494  3927268800 :       gradzt_2D(i,1) = gradzt_2D(i,2)
    2495  4177450368 :       gradzt_2D(i,nz) = gradzt_2D(i,nz-1)
    2496             :     end do
    2497             :     !$acc end parallel loop
    2498             :     
    2499             :     !$acc end data
    2500             : 
    2501   250181568 :     return
    2502             : 
    2503             :   end function gradzt_2D
    2504             :   
    2505             :   !=============================================================================
    2506             :   ! Wrapped in interface ddzt
    2507           0 :   function gradzt_1D( gr, azt )
    2508             : 
    2509             :     ! Description:
    2510             :     !  2D version of gradzt
    2511             :     !-----------------------------------------------------------------------
    2512             : 
    2513             :     use clubb_precision, only: &
    2514             :         core_rknd ! Variable(s)
    2515             : 
    2516             :     implicit none
    2517             : 
    2518             :     type (grid), target, intent(in) :: gr
    2519             : 
    2520             :     ! Input Variable
    2521             :     real( kind = core_rknd ), intent(in), dimension(gr%nz) :: &
    2522             :       azt    ! Variable on thermodynamic grid levels    [units vary]
    2523             : 
    2524             :     ! Output Variable
    2525             :     real( kind = core_rknd ), dimension(gr%nz) :: &
    2526             :       gradzt_1D    ! Vertical derivative of azt    [units vary / m]
    2527             : 
    2528             :     ! Local Variables
    2529             : 
    2530             :     ! Input Variable
    2531             :     real( kind = core_rknd ), dimension(1,gr%nz) :: &
    2532           0 :       azt_col    ! Variable on thermodynamic grid levels    [units vary]
    2533             : 
    2534             :     real( kind = core_rknd ), dimension(1,gr%nz) :: &
    2535           0 :       gradzt_1D_col    ! Vertical derivative of azt    [units vary / m]
    2536             : 
    2537           0 :     azt_col(1,:) = azt
    2538             :     
    2539           0 :     gradzt_1D_col = gradzt_2D(gr%nz, 1, gr, azt_col)
    2540             :     
    2541           0 :     gradzt_1D = gradzt_1D_col(1,:)
    2542             : 
    2543           0 :     return
    2544             : 
    2545             :   end function gradzt_1D
    2546             : 
    2547             :   !=============================================================================
    2548           0 :   function flip( x, xdim )
    2549             : 
    2550             :     ! Description:
    2551             :     !   Flips a single dimension array (i.e. a vector), so the first element
    2552             :     !   becomes the last and vice versa for the whole column.  This is a
    2553             :     !   necessary part of the code because BUGSrad and CLUBB store altitudes in
    2554             :     !   reverse order.
    2555             :     !
    2556             :     ! References:
    2557             :     !   None
    2558             :     !-------------------------------------------------------------------------
    2559             :     
    2560             :     use clubb_precision, only: &
    2561             :         dp ! double precision
    2562             : 
    2563             :     implicit none
    2564             : 
    2565             :     ! Input Variables
    2566             :     integer, intent(in) :: xdim
    2567             : 
    2568             :     real(kind = dp), dimension(xdim), intent(in) :: x
    2569             : 
    2570             :     ! Output Variables
    2571             :     real(kind = dp), dimension(xdim) :: flip
    2572             : 
    2573             :     ! Local Variables
    2574           0 :     real(kind = dp), dimension(xdim) :: tmp
    2575             : 
    2576             :     integer :: indx
    2577             : 
    2578             : 
    2579             :     ! Get rid of an annoying compiler warning.
    2580             :     indx = 1
    2581             :     indx = indx
    2582             : 
    2583           0 :     forall ( indx = 1 : xdim )
    2584           0 :        tmp(indx) = x((xdim+1) - (indx))
    2585             :     end forall
    2586             : 
    2587           0 :     flip = tmp
    2588             : 
    2589             : 
    2590           0 :     return
    2591             : 
    2592             :   end function flip
    2593             : 
    2594             : !===============================================================================
    2595             : 
    2596           0 : end module grid_class

Generated by: LCOV version 1.14