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 22:39:59 |  | 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