LCOV - code coverage report
Current view: | top level - utils - zonal_mean_mod.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 11 | 681 | 1.6 %
Date: | 2025-01-13 21:54:50 |  | Functions: | 2 | 28 | 7.1 %

              Line data    Source code
           1             : module zonal_mean_mod
           2             : !======================================================================
           3             : !
           4             : ! Purpose: Compute and make use of Zonal Mean values on physgrid
           5             : !
           6             : !    This module implements 3 data structures for the spectral analysis
           7             : !    and synthesis of zonal mean values based on m=0 spherical harmonics.
           8             : !
           9             : !    ZonalMean_t:    For the analysis/synthesis of zonal mean values
          10             : !                    on a 2D grid of points distributed over the
          11             : !                    surface of a sphere.
          12             : !    ZonalProfile_t: For the analysis/synthesis of zonal mean values
          13             : !                    on a meridional grid that spans the latitudes
          14             : !                    from SP to NP
          15             : !    ZonalAverage_t: To calculate zonal mean values via a simple
          16             : !                    area weighted bin-averaging of 2D grid points
          17             : !                    assigned to each latitude band.
          18             : !
          19             : !   NOTE: The weighting of the Zonal Profiles values is scaled such
          20             : !         that ZonalMean_t amplitudes can be used to evaluate values
          21             : !         on the ZonalProfile_t grid and vice-versa.
          22             : !
          23             : !         The ZonalMean_t computes global integrals to compute basis
          24             : !         amplitudes. For distributed environments the cost of these
          25             : !         can be reduced using the The ZonalAverage_t data structures.
          26             : !
          27             : ! USAGE:
          28             : !
          29             : !    (1) Compute Zonal mean amplitudes and synthesize values on 2D/3D physgrid
          30             : !
          31             : !         Usage: type(ZonalMean_t):: ZM
          32             : !         =========================================
          33             : !           call ZM%init(nbas)
          34             : !           ------------------
          35             : !               - Initialize the data structure with 'nbas' basis functions
          36             : !                 for the given physgrid latitudes and areas.
          37             : !
          38             : !               Arguments:
          39             : !                 integer ,intent(in):: nbas     -Number of m=0 spherical harmonics
          40             : !
          41             : !           call ZM%calc_amps(Gdata,Bamp)
          42             : !           -----------------------------
          43             : !               - For the initialized ZonalMean_t; Given Gdata() values on the physgrid,
          44             : !                 compute the zonal mean basis amplitudes Bamp().
          45             : !
          46             : !               Interface: 2D data on the physgrid
          47             : !                 real(r8),intent(in ):: Gdata(pcols,begchunk:endchunk)
          48             : !                 real(r8),intent(out):: Bamp (nbas)
          49             : !
          50             : !               Interface: 3D data on the physgrid
          51             : !                 real(r8),intent(in ):: Gdata(pcols,pver,begchunk:endchunk)
          52             : !                 real(r8),intent(out):: Bamp (nbas,pver)
          53             : !
          54             : !           call ZM%eval_grid(Bamp,Gdata)
          55             : !           -----------------------------
          56             : !               - For the initialized ZonalMean_t; Given Bamp() zonal mean basis
          57             : !                 amplitudes, compute the Gdata() values on the physgrid.
          58             : !
          59             : !               Interface: 2D data on the physgrid
          60             : !                 real(r8),intent(in ):: Bamp (nbas)
          61             : !                 real(r8),intent(out):: Gdata(pcols,begchunk:endchunk)
          62             : !
          63             : !               Interface: 3D data on the physgrid
          64             : !                 real(r8),intent(in ):: Bamp (nbas,pver)
          65             : !                 real(r8),intent(out):: Gdata(pcols,pver,begchunk:endchunk)
          66             : !
          67             : !
          68             : !    (2) Compute Zonal mean amplitudes and synthesize values on Zonal profile grid
          69             : !
          70             : !         Usage: type(ZonalProfile_t):: ZP
          71             : !         =========================================
          72             : !           call ZP%init(lats,area,nlat,nbas,GEN_GAUSSLATS=.true.)
          73             : !           ------------------------------------------------------
          74             : !               - Initialize the data structure for the given number of
          75             : !                 latitudes. Either use the given Latitudes and weights,
          76             : !                 or OPTIONALLY create profile gridpoints and associated
          77             : !                 area weights from SP to NP. Then initialize 'nbas' basis
          78             : !                 functions for the profile gridpoints.
          79             : !                 If the user supplies the lats/area values, the area values must
          80             : !                 be correctly scaled such that the global area adds up to 4PI.
          81             : !                 Otherwise, the ampitudes between ZonalProfile_t and ZonalMean_t
          82             : !                 are not interchangable.
          83             : !
          84             : !               Arguments:
          85             : !                 real(r8),intent(inout):: lats(:) - Latitudes of meridional grid.
          86             : !                 real(r8),intent(inout):: area(:) - Area of each meridional gridpoint.
          87             : !                 integer ,intent(in)   :: nlat    - Number of meridional gridpoints.
          88             : !                 integer ,intent(in)   :: nbas    - Number of m=0 spherical harmonics
          89             : !                 logical ,intent(in),optional:: GEN_GAUSLATS - Flag to generate
          90             : !                                                               lats/areas values.
          91             : !
          92             : !           call ZP%calc_amps(Zdata,Bamp)
          93             : !           -----------------------------
          94             : !               - Given Zdata() on the Zonal profile grid, compute the
          95             : !                 zonal basis amplitudes Bamp().
          96             : !
          97             : !               Interface: 1D data on (nlat) grid
          98             : !                 real(r8),intent(in ):: Zdata(nlat) - Meridional Profile data
          99             : !                 real(r8),intent(out):: Bamp (nbas) - Zonal Basis Amplitudes
         100             : !
         101             : !               Interface: 2D data on (nlat,pver) grid
         102             : !                 real(r8),intent(in ):: Zdata(nlat,pver) - Meridional Profile data
         103             : !                 real(r8),intent(out):: Bamp (nbas,pver) - Zonal Basis Amplitudes
         104             : !
         105             : !           call ZP%eval_grid(Bamp,Zdata)
         106             : !           -----------------------------
         107             : !               - Given Bamp() zonal basis amplitudes, evaluate the Zdata()
         108             : !                 values on the Zonal profile grid.
         109             : !
         110             : !               Interface: 1D data on (nlat) grid
         111             : !                 real(r8),intent(in ):: Bamp (nbas) - Zonal Basis Amplitudes
         112             : !                 real(r8),intent(out):: Zdata(nlat) - Meridional Profile data
         113             : !
         114             : !               Interface: 2D data on (nlat,pver) grid
         115             : !                 real(r8),intent(in ):: Bamp (nbas,pver) - Zonal Basis Amplitudes
         116             : !                 real(r8),intent(out):: Zdata(nlat,pver) - Meridional Profile data
         117             : !
         118             : !    (3) Compute Zonal mean averages (FASTER/LESS-ACCURATE) on Zonal profile grid
         119             : !        (For the created zonal profile, just bin average area weighted
         120             : !         2D/3D physgrid grid values)
         121             : !
         122             : !         Usage: type(ZonalAverage_t):: ZA
         123             : !         =========================================
         124             : !           call ZA%init(lats,area,nlat,GEN_GAUSSLATS=.true.)
         125             : !           --------------------------------------------------
         126             : !               - Given the latitude/area for the nlat meridional gridpoints, initialize
         127             : !                 the ZonalAverage data structure for computing bin-averaging of physgrid
         128             : !                 values. It is assumed that the domain of these gridpoints of the
         129             : !                 profile span latitudes from SP to NP.
         130             : !                 The optional GEN_GAUSSLATS flag allows for the generation of Gaussian
         131             : !                 latitude gridpoints. The generated grid over-writes the given values
         132             : !                 lats and area passed by the user.
         133             : !
         134             : !               Arguments:
         135             : !                 real(r8),intent(inout):: lats(nlat) - Latitudes of meridional grid.
         136             : !                 real(r8),intent(inout):: area(nlat) - Area of meridional gridpoints.
         137             : !                 integer    ,intent(in):: nlat       - Number of meridional gridpoints
         138             : !                 logical,intent(in),optional:: GEN_GAUSLATS - Flag to generate
         139             : !                                                              lats/areas values.
         140             : !
         141             : !           call ZA%binAvg(Gdata,Zdata)
         142             : !           ---------------------------
         143             : !               - For the initialized ZonalAverage_t; Given Gdata() on the physgrid,
         144             : !                 compute bin averages and return Zdata() on the Zonal profile grid.
         145             : !
         146             : !               Interface: 2D data on the physgrid
         147             : !                 real(r8),intent(out):: Gdata(pcols,begchunk:endchunk)
         148             : !                 real(r8),intent(out):: Zdata(nlat)
         149             : !
         150             : !               Interface: 3D data on the physgrid
         151             : !                 real(r8),intent(out):: Gdata(pcols,pver,begchunk:endchunk)
         152             : !                 real(r8),intent(out):: Zdata(nlat,pver)
         153             : !
         154             : !======================================================================
         155             :
         156             :   use shr_kind_mod,    only: r8=>SHR_KIND_R8
         157             :   use phys_grid,       only: get_ncols_p, get_rlat_p, get_wght_all_p, get_nlcols_p
         158             :   use ppgrid,          only: begchunk, endchunk, pcols
         159             :   use shr_reprosum_mod,only: shr_reprosum_calc
         160             :   use cam_abortutils,  only: endrun, handle_allocate_error
         161             :   use spmd_utils,      only: mpicom
         162             :   use physconst,       only: pi
         163             :   use phys_grid,       only: ngcols_p => num_global_phys_cols
         164             :   use cam_logfile,     only: iulog
         165             :
         166             :   implicit none
         167             :   private
         168             :
         169             :   public :: ZonalMean_t
         170             :   public :: ZonalProfile_t
         171             :   public :: ZonalAverage_t
         172             :
         173             :   ! Type definitions
         174             :   !-------------------
         175             :   type ZonalMean_t
         176             :      private
         177             :      integer             :: nbas
         178             :      real(r8),allocatable:: area (:,:)
         179             :      real(r8),allocatable:: basis(:,:,:)
         180             :      real(r8),allocatable:: map  (:,:)
         181             :     contains
         182             :      procedure,pass:: init      => init_ZonalMean
         183             :      generic,public:: calc_amps => calc_ZonalMean_2Damps, &
         184             :                                    calc_ZonalMean_3Damps
         185             :      generic,public:: eval_grid => eval_ZonalMean_2Dgrid, &
         186             :                                    eval_ZonalMean_3Dgrid
         187             :      procedure,private,pass:: calc_ZonalMean_2Damps
         188             :      procedure,private,pass:: calc_ZonalMean_3Damps
         189             :      procedure,private,pass:: eval_ZonalMean_2Dgrid
         190             :      procedure,private,pass:: eval_ZonalMean_3Dgrid
         191             :      procedure, pass :: final => final_ZonalMean
         192             :  end type ZonalMean_t
         193             :
         194             :   type ZonalProfile_t
         195             :      private
         196             :      integer             :: nlat
         197             :      integer             :: nbas
         198             :      real(r8),allocatable:: area (:)
         199             :      real(r8),allocatable:: basis(:,:)
         200             :      real(r8),allocatable:: map  (:,:)
         201             :     contains
         202             :      procedure,pass:: init      => init_ZonalProfile
         203             :      generic,public:: calc_amps => calc_ZonalProfile_1Damps, &
         204             :                                    calc_ZonalProfile_2Damps
         205             :      generic,public:: eval_grid => eval_ZonalProfile_1Dgrid, &
         206             :                                    eval_ZonalProfile_2Dgrid
         207             :      procedure,private,pass:: calc_ZonalProfile_1Damps
         208             :      procedure,private,pass:: calc_ZonalProfile_2Damps
         209             :      procedure,private,pass:: eval_ZonalProfile_1Dgrid
         210             :      procedure,private,pass:: eval_ZonalProfile_2Dgrid
         211             :      procedure, pass :: final => final_ZonalProfile
         212             :   end type ZonalProfile_t
         213             :
         214             :   type ZonalAverage_t
         215             :      private
         216             :      integer             :: nlat
         217             :      real(r8),allocatable:: area   (:)
         218             :      real(r8),allocatable:: a_norm (:)
         219             :      real(r8),allocatable:: area_g (:,:)
         220             :      integer ,allocatable:: idx_map(:,:)
         221             :     contains
         222             :      procedure,pass:: init   => init_ZonalAverage
         223             :      generic,public:: binAvg => calc_ZonalAverage_2DbinAvg, &
         224             :                                 calc_ZonalAverage_3DbinAvg
         225             :      procedure,private,pass:: calc_ZonalAverage_2DbinAvg
         226             :      procedure,private,pass:: calc_ZonalAverage_3DbinAvg
         227             :      procedure, pass :: final => final_ZonalAverage
         228             :   end type ZonalAverage_t
         229             :
         230             :   real(r8), parameter :: halfPI = 0.5_r8*pi
         231             :   real(r8), parameter :: twoPI  = 2.0_r8*pi
         232             :   real(r8), parameter :: fourPI = 4.0_r8*pi
         233             :   real(r8), parameter :: qrtrPI = 0.25_r8*pi
         234             :   real(r8), parameter :: invSqrt4pi = 1._r8/sqrt(fourPI)
         235             :
         236             : contains
         237             :     !=======================================================================
         238           0 :     subroutine init_ZonalMean(this,I_nbas)
         239             :       !
         240             :       ! init_ZonalMean: Initialize the ZonalMean data structures for the
         241             :       !                 physics grid. It is assumed that the domain
         242             :       !                 of these gridpoints spans the surface of the sphere.
         243             :       !                 The representation of basis functions is
         244             :       !                 normalized w.r.t integration over the sphere.
         245             :       !=====================================================================
         246             :       !
         247             :       ! Passed Variables
         248             :       !------------------
         249             :       class(ZonalMean_t) :: this
         250             :       integer ,intent(in):: I_nbas
         251             :       !
         252             :       ! Local Values
         253             :       !--------------
         254           0 :       real(r8),allocatable:: Clats(:,:)
         255           0 :       real(r8),allocatable:: Bcoef(:)
         256           0 :       real(r8),allocatable:: Csum (:,:)
         257           0 :       real(r8),allocatable:: Cvec (:)
         258           0 :       real(r8),allocatable:: Bsum (:,:)
         259           0 :       real(r8),allocatable:: Bnorm(:)
         260           0 :       real(r8),allocatable:: Bcov (:,:)
         261             :       real(r8):: area(pcols),rlat
         262             :
         263             :       integer :: nn,n2,nb,lchnk,ncols,cc
         264             :       integer :: cnum,Cvec_len
         265             :
         266             :       integer :: nlcols, count, astat
         267             :       character(len=*), parameter :: subname = 'init_ZonalMean'
         268             :
         269           0 :       if (I_nbas<1) then
         270           0 :          call endrun('ZonalMean%init: ERROR I_nbas must be greater than 0')
         271             :       end if
         272             :
         273             :       ! Allocate space
         274             :       !-----------------
         275           0 :       if(allocated(this%area )) deallocate(this%area)
         276           0 :       if(allocated(this%basis)) deallocate(this%basis)
         277           0 :       if(allocated(this%map  )) deallocate(this%map)
         278             :
         279           0 :       this%nbas = I_nbas
         280           0 :       allocate(this%area (pcols,begchunk:endchunk), stat=astat)
         281           0 :       call handle_allocate_error(astat, subname, 'this%area')
         282           0 :       allocate(this%basis(pcols,begchunk:endchunk,I_nbas), stat=astat)
         283           0 :       call handle_allocate_error(astat, subname, 'this%basis')
         284           0 :       allocate(this%map  (I_nbas,I_nbas), stat=astat)
         285           0 :       call handle_allocate_error(astat, subname, 'this%map')
         286           0 :       this%area (:,:)   = 0._r8
         287           0 :       this%basis(:,:,:) = 0._r8
         288           0 :       this%map  (:,:)   = 0._r8
         289             :
         290           0 :       Cvec_len = 0
         291           0 :       do nn= 1,this%nbas
         292           0 :          do n2=nn,this%nbas
         293           0 :             Cvec_len = Cvec_len + 1
         294             :          end do
         295             :       end do
         296             :
         297           0 :       nlcols = get_nlcols_p()
         298             :
         299           0 :       allocate(Clats(pcols,begchunk:endchunk), stat=astat)
         300           0 :       call handle_allocate_error(astat, subname, 'Clats')
         301           0 :       allocate(Bcoef(I_nbas/2+1), stat=astat)
         302           0 :       call handle_allocate_error(astat, subname, 'Bcoef')
         303           0 :       allocate(Csum (nlcols, Cvec_len), stat=astat)
         304           0 :       call handle_allocate_error(astat, subname, 'Csum')
         305           0 :       allocate(Cvec (Cvec_len), stat=astat)
         306           0 :       call handle_allocate_error(astat, subname, 'Cvec')
         307           0 :       allocate(Bsum (nlcols, I_nbas), stat=astat)
         308           0 :       call handle_allocate_error(astat, subname, 'Bsum')
         309           0 :       allocate(Bnorm(I_nbas), stat=astat)
         310           0 :       call handle_allocate_error(astat, subname, 'Bnorm')
         311           0 :       allocate(Bcov (I_nbas,I_nbas), stat=astat)
         312           0 :       call handle_allocate_error(astat, subname, 'Bcov')
         313             :
         314           0 :       Bsum(:,:) = 0._r8
         315           0 :       Csum(:,:) = 0._r8
         316             :
         317             :       ! Save a copy of the area weights for each ncol gridpoint
         318             :       ! and convert Latitudes to SP->NP colatitudes in radians
         319             :       !-------------------------------------------------------
         320           0 :       do lchnk=begchunk,endchunk
         321           0 :          ncols = get_ncols_p(lchnk)
         322           0 :          call get_wght_all_p(lchnk, ncols, area)
         323           0 :          do cc = 1,ncols
         324           0 :             rlat=get_rlat_p(lchnk,cc)
         325           0 :             this%area(cc,lchnk) = area(cc)
         326           0 :             Clats    (cc,lchnk) = rlat + halfPI
         327             :          end do
         328             :       end do
         329             :
         330             :       ! Add first basis for the mean values.
         331             :       !------------------------------------------
         332           0 :       this%basis(:,begchunk:endchunk,1) = invSqrt4pi
         333             :
         334             :       ! Loop over the remaining basis functions
         335             :       !---------------------------------------
         336           0 :       do nn=2,this%nbas
         337           0 :          nb = nn-1
         338             :
         339             :          ! Generate coefs for the basis
         340             :          !------------------------------
         341           0 :          call sh_gen_basis_coefs(nb,0,Bcoef)
         342             :
         343             :          ! Create basis for the coefs at each ncol gridpoint
         344             :          !---------------------------------------------------
         345           0 :          do lchnk=begchunk,endchunk
         346           0 :             ncols = get_ncols_p(lchnk)
         347           0 :             do cc = 1,ncols
         348           0 :                call sh_create_basis(nb,0,Clats(cc,lchnk),Bcoef,this%basis(cc,lchnk,nn))
         349             :             end do
         350             :          end do
         351             :       end do ! nn=2,this%nbas
         352             :
         353             :       ! Numerically normalize the basis funnctions
         354             :       !--------------------------------------------------------------
         355           0 :       do nn=1,this%nbas
         356           0 :          count = 0
         357           0 :          do lchnk=begchunk,endchunk
         358           0 :             ncols = get_ncols_p(lchnk)
         359           0 :             do cc = 1,ncols
         360           0 :                count=count+1
         361           0 :                Bsum(count,nn) = this%basis(cc,lchnk,nn)*this%basis(cc,lchnk,nn)*this%area(cc,lchnk)
         362             :             end do
         363             :          end do
         364             :       end do ! nn=1,this%nbas
         365             :
         366           0 :       call shr_reprosum_calc(Bsum, Bnorm, count, nlcols, this%nbas, gbl_count=ngcols_p, commid=mpicom)
         367             :
         368           0 :       do nn=1,this%nbas
         369           0 :          do lchnk=begchunk,endchunk
         370           0 :             ncols = get_ncols_p(lchnk)
         371           0 :             this%basis(:ncols,lchnk,nn) = this%basis(:ncols,lchnk,nn)/sqrt(Bnorm(nn))
         372             :          end do
         373             :       end do ! nn=1,this%nbas
         374             :
         375             :       ! Compute covariance matrix for basis functions
         376             :       ! (Yes, they are theoretically orthonormal, but lets make sure)
         377             :       !---------------------------------------------------------------
         378           0 :       cnum = 0
         379           0 :       do nn= 1,this%nbas
         380           0 :          do n2=nn,this%nbas
         381           0 :             cnum = cnum + 1
         382           0 :             count = 0
         383           0 :             do lchnk=begchunk,endchunk
         384           0 :                ncols = get_ncols_p(lchnk)
         385           0 :                do cc = 1,ncols
         386           0 :                   count=count+1
         387           0 :                   Csum(count,cnum) = this%basis(cc,lchnk,nn)*this%basis(cc,lchnk,n2)*this%area(cc,lchnk)
         388             :                end do
         389             :             end do
         390             :
         391             :          end do
         392             :       end do
         393             :
         394           0 :       call shr_reprosum_calc(Csum, Cvec, count, nlcols, Cvec_len, gbl_count=ngcols_p, commid=mpicom)
         395             :
         396           0 :       cnum = 0
         397           0 :       do nn= 1,this%nbas
         398           0 :          do n2=nn,this%nbas
         399           0 :             cnum = cnum + 1
         400           0 :             Bcov(nn,n2) = Cvec(cnum)
         401           0 :             Bcov(n2,nn) = Cvec(cnum)
         402             :          end do
         403             :       end do
         404             :
         405             :       ! Invert to get the basis amplitude map
         406             :       !--------------------------------------
         407           0 :       call Invert_Matrix(Bcov,this%nbas,this%map)
         408             :
         409             :       ! End Routine
         410             :       !------------
         411           0 :       deallocate(Clats)
         412           0 :       deallocate(Bcoef)
         413           0 :       deallocate(Csum )
         414           0 :       deallocate(Cvec )
         415           0 :       deallocate(Bsum )
         416           0 :       deallocate(Bnorm)
         417           0 :       deallocate(Bcov )
         418             :
         419           0 :     end subroutine init_ZonalMean
         420             :     !=======================================================================
         421             :
         422             :
         423             :     !=======================================================================
         424           0 :     subroutine calc_ZonalMean_2Damps(this,I_Gdata,O_Bamp)
         425             :       !
         426             :       ! calc_ZonalMean_2Damps: Given 2D data values for the ncol gridpoints,
         427             :       !                        compute the zonal mean basis amplitudes.
         428             :       !=====================================================================
         429             :       !
         430             :       ! Passed Variables
         431             :       !------------------
         432             :       class(ZonalMean_t) :: this
         433             :       real(r8),intent(in ) :: I_Gdata(pcols,begchunk:endchunk)
         434             :       real(r8),intent(out) :: O_Bamp(:)
         435             :       !
         436             :       ! Local Values
         437             :       !--------------
         438           0 :       real(r8),allocatable :: Csum(:,:)
         439           0 :       real(r8),allocatable :: Gcov(:)
         440             :       integer :: nn,n2,ncols,lchnk,cc
         441             :       integer :: nlcols, count, astat
         442             :
         443             :       character(len=*), parameter :: subname = 'calc_ZonalMean_2Damps'
         444             :
         445           0 :       nlcols = get_nlcols_p()
         446             :
         447           0 :       allocate(Gcov(this%nbas), stat=astat)
         448           0 :       call handle_allocate_error(astat, subname, 'Gcov')
         449           0 :       allocate(Csum(nlcols, this%nbas), stat=astat)
         450           0 :       call handle_allocate_error(astat, subname, 'Csum')
         451           0 :       Csum(:,:) = 0._r8
         452             :
         453             :       ! Compute Covariance with input data and basis functions
         454             :       !--------------------------------------------------------
         455           0 :       do nn= 1,this%nbas
         456           0 :         count = 0
         457           0 :         do lchnk=begchunk,endchunk
         458           0 :           ncols = get_ncols_p(lchnk)
         459           0 :           do cc = 1,ncols
         460           0 :             count=count+1
         461           0 :             Csum(count,nn) = I_Gdata(cc,lchnk)*this%basis(cc,lchnk,nn)*this%area(cc,lchnk)
         462             :           end do
         463             :         end do
         464             :       end do
         465             :
         466           0 :       call shr_reprosum_calc(Csum, Gcov, count, nlcols, this%nbas, gbl_count=ngcols_p, commid=mpicom)
         467             :
         468             :       ! Multiply by map to get the amplitudes
         469             :       !-------------------------------------------
         470           0 :       do nn=1,this%nbas
         471           0 :         O_Bamp(nn) = 0._r8
         472           0 :         do n2=1,this%nbas
         473           0 :           O_Bamp(nn) = O_Bamp(nn) + this%map(n2,nn)*Gcov(n2)
         474             :         end do
         475             :       end do
         476             :
         477             :       ! End Routine
         478             :       !------------
         479           0 :       deallocate(Csum)
         480           0 :       deallocate(Gcov)
         481             :
         482           0 :     end subroutine calc_ZonalMean_2Damps
         483             :     !=======================================================================
         484             :
         485             :
         486             :     !=======================================================================
         487           0 :     subroutine calc_ZonalMean_3Damps(this,I_Gdata,O_Bamp)
         488             :       !
         489             :       ! calc_ZonalMean_3Damps: Given 3D data values for the ncol,nlev gridpoints,
         490             :       !                        compute the zonal mean basis amplitudes.
         491             :       !=====================================================================
         492             :       !
         493             :       ! Passed Variables
         494             :       !------------------
         495             :       class(ZonalMean_t) :: this
         496             :       real(r8),intent(in ):: I_Gdata(:,:,begchunk:)
         497             :       real(r8),intent(out):: O_Bamp (:,:)
         498             :       !
         499             :       ! Local Values
         500             :       !--------------
         501           0 :       real(r8),allocatable:: Csum (:,:)
         502           0 :       real(r8),allocatable:: Gcov (:)
         503             :       integer:: nn,n2,ncols,lchnk,cc
         504             :       integer:: Nsum,ns,ll
         505             :       integer :: nlcols, count, astat
         506             :
         507             :       integer :: nlev
         508             :       character(len=*), parameter :: subname = 'calc_ZonalMean_3Damps'
         509             :
         510           0 :       nlev = size(I_Gdata,dim=2)
         511             :
         512           0 :       nlcols = get_nlcols_p()
         513           0 :       allocate(Gcov(this%nbas), stat=astat)
         514           0 :       call handle_allocate_error(astat, subname, 'Gcov')
         515           0 :       allocate(Csum(nlcols, this%nbas), stat=astat)
         516           0 :       call handle_allocate_error(astat, subname, 'Csum')
         517             :
         518           0 :       Csum(:,:) = 0._r8
         519           0 :       O_Bamp(:,:) = 0._r8
         520             :
         521             :       ! Compute Covariance with input data and basis functions
         522             :       !--------------------------------------------------------
         523           0 :       do ll= 1,nlev
         524             :
         525           0 :          Csum(:,:) = 0._r8
         526           0 :          Gcov(:) = 0._r8
         527             :
         528           0 :          do nn= 1,this%nbas
         529           0 :             count = 0
         530           0 :             do lchnk=begchunk,endchunk
         531           0 :                ncols = get_ncols_p(lchnk)
         532           0 :                do cc = 1,ncols
         533           0 :                   count=count+1
         534           0 :                   Csum(count,nn) = I_Gdata(cc,ll,lchnk)*this%basis(cc,lchnk,nn)*this%area(cc,lchnk)
         535             :                end do
         536             :             end do
         537             :          end do
         538             :
         539           0 :          call shr_reprosum_calc(Csum, Gcov, count, nlcols, this%nbas, gbl_count=ngcols_p, commid=mpicom)
         540             :
         541             :          ! Multiply by map to get the amplitudes
         542             :          !-------------------------------------------
         543           0 :          do nn=1,this%nbas
         544           0 :             O_Bamp(nn,ll) = 0._r8
         545           0 :             do n2=1,this%nbas
         546           0 :                O_Bamp(nn,ll) = O_Bamp(nn,ll) + this%map(n2,nn)*Gcov(n2)
         547             :             end do
         548             :          end do
         549             :
         550             :       end do
         551             :
         552             :       ! End Routine
         553             :       !------------
         554           0 :       deallocate(Csum)
         555           0 :       deallocate(Gcov)
         556             :
         557           0 :     end subroutine calc_ZonalMean_3Damps
         558             :     !=======================================================================
         559             :
         560             :
         561             :     !=======================================================================
         562           0 :     subroutine eval_ZonalMean_2Dgrid(this,I_Bamp,O_Gdata)
         563             :       !
         564             :       ! eval_ZonalMean_2Dgrid: Given the zonal mean basis amplitudes,
         565             :       !                        compute 2D data values for the ncol gridpoints.
         566             :       !=====================================================================
         567             :       !
         568             :       ! Passed Variables
         569             :       !------------------
         570             :       class(ZonalMean_t) :: this
         571             :       real(r8),intent(in ):: I_Bamp (:)
         572             :       real(r8),intent(out):: O_Gdata(pcols,begchunk:endchunk)
         573             :       !
         574             :       ! Local Values
         575             :       !--------------
         576             :       integer:: nn,ncols,lchnk,cc
         577             :
         578           0 :       O_Gdata(:,:) = 0._r8
         579             :
         580             :       ! Construct grid values from basis amplitudes.
         581             :       !--------------------------------------------------
         582             :
         583           0 :       do nn=1,this%nbas
         584           0 :          do lchnk=begchunk,endchunk
         585           0 :             ncols = get_ncols_p(lchnk)
         586           0 :             do cc = 1,ncols
         587           0 :                O_Gdata(cc,lchnk) = O_Gdata(cc,lchnk) + (I_Bamp(nn)*this%basis(cc,lchnk,nn))
         588             :             end do
         589             :          end do
         590             :       end do
         591             :
         592           0 :     end subroutine eval_ZonalMean_2Dgrid
         593             :     !=======================================================================
         594             :
         595             :
         596             :     !=======================================================================
         597           0 :     subroutine eval_ZonalMean_3Dgrid(this,I_Bamp,O_Gdata)
         598             :       !
         599             :       ! eval_ZonalMean_3Dgrid: Given the zonal mean basis amplitudes,
         600             :       !                      compute 3D data values for the ncol,nlev gridpoints.
         601             :       !=====================================================================
         602             :       !
         603             :       ! Passed Variables
         604             :       !------------------
         605             :       class(ZonalMean_t) :: this
         606             :       real(r8),intent(in ):: I_Bamp (:,:)
         607             :       real(r8),intent(out):: O_Gdata(:,:,begchunk:)
         608             :       !
         609             :       ! Local Values
         610             :       !--------------
         611             :       integer:: nn,ncols,lchnk,cc
         612             :       integer:: ll
         613             :
         614             :       integer :: nlev
         615           0 :       nlev = size(O_Gdata,dim=2)
         616             :
         617           0 :       O_Gdata(:,:,:) = 0._r8
         618             :
         619             :       ! Construct grid values from basis amplitudes.
         620             :       !--------------------------------------------------
         621             :
         622           0 :       do ll = 1,nlev
         623           0 :          do nn=1,this%nbas
         624           0 :             do lchnk=begchunk,endchunk
         625           0 :                ncols = get_ncols_p(lchnk)
         626           0 :                do cc = 1,ncols
         627           0 :                   O_Gdata(cc,ll,lchnk) = O_Gdata(cc,ll,lchnk) + (I_Bamp(nn,ll)*this%basis(cc,lchnk,nn))
         628             :                end do
         629             :             end do
         630             :          end do
         631             :       end do
         632             :
         633           0 :     end subroutine eval_ZonalMean_3Dgrid
         634             :     !=======================================================================
         635             :
         636             :     !=======================================================================
         637        1536 :     subroutine final_ZonalMean(this)
         638             :       class(ZonalMean_t) :: this
         639             :
         640        1536 :       if(allocated(this%area )) deallocate(this%area)
         641        1536 :       if(allocated(this%basis)) deallocate(this%basis)
         642        1536 :       if(allocated(this%map  )) deallocate(this%map)
         643             :
         644        1536 :     end subroutine final_ZonalMean
         645             :     !=======================================================================
         646             :
         647             :     !=======================================================================
         648           0 :     subroutine init_ZonalProfile(this,IO_lats,IO_area,I_nlat,I_nbas,GEN_GAUSSLATS)
         649             :       !
         650             :       ! init_ZonalProfile: Initialize the ZonalProfile data structure for the
         651             :       !                    given nlat gridpoints. It is assumed that the domain
         652             :       !                    of these gridpoints of the profile span latitudes
         653             :       !                    from SP to NP.
         654             :       !                    The representation of basis functions functions is
         655             :       !                    normalized w.r.t integration over the sphere so that
         656             :       !                    when configured for tha same number of basis functions,
         657             :       !                    the calculated amplitudes are interchangable with
         658             :       !                    those for the ZonalMean_t class.
         659             :       !
         660             :       !                    The optional GEN_GAUSSLATS flag allows for the
         661             :       !                    generation of Gaussian latitudes. The generated grid
         662             :       !                    over-writes the values of IO_lats/IO_area passed by
         663             :       !                    the user.
         664             :       !=====================================================================
         665             :       !
         666             :       ! Passed Variables
         667             :       !------------------
         668             :       class(ZonalProfile_t) :: this
         669             :       real(r8)     ,intent(inout):: IO_lats(:)
         670             :       real(r8)     ,intent(inout):: IO_area(:)
         671             :       integer         ,intent(in):: I_nlat
         672             :       integer         ,intent(in):: I_nbas
         673             :       logical,optional,intent(in):: GEN_GAUSSLATS
         674             :       !
         675             :       ! Local Values
         676             :       !--------------
         677           0 :       real(r8),allocatable:: Clats(:)
         678           0 :       real(r8),allocatable:: Bcoef(:)
         679           0 :       real(r8),allocatable:: Bcov (:,:)
         680             :       real(r8):: Bnorm
         681             :       integer :: ii,nn,n2,nb,ierr, astat
         682             :       logical :: generate_lats
         683             :
         684             :       character(len=*), parameter :: subname = 'init_ZonalProfile'
         685             :
         686           0 :       generate_lats = .false.
         687             :
         688           0 :       if (present(GEN_GAUSSLATS)) then
         689           0 :          generate_lats = GEN_GAUSSLATS
         690             :       end if
         691             :
         692             :       ! Allocate space
         693             :       !-----------------
         694           0 :       if(allocated(this%area )) deallocate(this%area)
         695           0 :       if(allocated(this%basis)) deallocate(this%basis)
         696           0 :       if(allocated(this%map  )) deallocate(this%map)
         697             :
         698           0 :       this%nlat = I_nlat
         699           0 :       this%nbas = I_nbas
         700           0 :       allocate(this%area (I_nlat), stat=astat)
         701           0 :       call handle_allocate_error(astat, subname, 'this%area')
         702           0 :       allocate(this%basis(I_nlat,I_nbas), stat=astat)
         703           0 :       call handle_allocate_error(astat, subname, 'this%basis')
         704           0 :       allocate(this%map  (I_nbas,I_nbas), stat=astat)
         705           0 :       call handle_allocate_error(astat, subname, 'this%map')
         706             :
         707           0 :       allocate(Clats(I_nlat), stat=astat)
         708           0 :       call handle_allocate_error(astat, subname, 'Clats')
         709           0 :       allocate(Bcoef(I_nbas/2+1), stat=astat)
         710           0 :       call handle_allocate_error(astat, subname, 'Bcoef')
         711           0 :       allocate(Bcov (I_nbas,I_nbas), stat=astat)
         712           0 :       call handle_allocate_error(astat, subname, 'Bcov')
         713             :
         714             :       ! Optionally create the Latitude Gridpoints
         715             :       ! and their associated area weights. Otherwise
         716             :       ! they need to be supplied by the user.
         717             :       !-----------------------------------------------
         718           0 :       if(generate_lats) then
         719             :
         720             :         ! Create a Gaussian grid from SP to NP
         721             :         !--------------------------------------
         722           0 :         call sh_create_gaus_grid(I_nlat,Clats,IO_area,ierr)
         723           0 :         if (ierr/=0) then
         724           0 :            call endrun('init_ZonalProfile: Error creating Gaussian grid')
         725             :         end if
         726             :
         727             :         ! Convert generated colatitudes SP->NP to Lats and convert
         728             :         ! to degrees and scale the area for global 2D integrals
         729             :         !-----------------------------------------------------------
         730           0 :         do nn=1,I_nlat
         731           0 :           IO_lats(nn) = (45._r8*Clats(nn)/qrtrPI) - 90._r8
         732           0 :           IO_area(nn) = IO_area(nn)*twoPI
         733             :         end do
         734             :       else
         735             :         ! Convert Latitudes to SP->NP colatitudes in radians
         736             :         !----------------------------------------------------
         737           0 :         do nn=1,I_nlat
         738           0 :           Clats(nn) = (IO_lats(nn) + 90._r8)*qrtrPI/45._r8
         739             :         end do
         740             :       endif
         741             :
         742             :       ! Copy the area weights for each nlat
         743             :       ! gridpoint to the data structure
         744             :       !---------------------------------------
         745           0 :       this%area(1:I_nlat) = IO_area(1:I_nlat)
         746             :
         747             :       ! Add first basis for the mean values.
         748             :       !------------------------------------------
         749           0 :       this%basis(:,1) = invSqrt4pi
         750             :       Bnorm = 0._r8
         751           0 :       do ii=1,I_nlat
         752           0 :         Bnorm = Bnorm + (this%basis(ii,1)*this%basis(ii,1)*this%area(ii))
         753             :       end do
         754           0 :       this%basis(:,1) = this%basis(:,1)/sqrt(Bnorm)
         755             :
         756             :       ! Loop over the remaining basis functions
         757             :       !---------------------------------------
         758           0 :       do nn=2,I_nbas
         759           0 :         nb = nn-1
         760             :
         761             :         ! Generate coefs for the basis
         762             :         !------------------------------
         763           0 :         call sh_gen_basis_coefs(nb,0,Bcoef)
         764             :
         765             :         ! Create an un-normalized basis for the
         766             :         ! coefs at each nlat gridpoint
         767             :         !---------------------------------------
         768           0 :         do ii=1,I_nlat
         769           0 :           call sh_create_basis(nb,0,Clats(ii),Bcoef,this%basis(ii,nn))
         770             :         end do
         771             :
         772             :         ! Numerically normalize the basis funnction
         773             :         !--------------------------------------------------------------
         774             :         Bnorm = 0._r8
         775           0 :         do ii=1,I_nlat
         776           0 :           Bnorm = Bnorm + (this%basis(ii,nn)*this%basis(ii,nn)*this%area(ii))
         777             :         end do
         778           0 :         this%basis(:,nn) = this%basis(:,nn)/sqrt(Bnorm)
         779             :
         780             :       end do ! nn=1,I_nbas
         781             :
         782             :       ! Compute covariance matrix for basis functions
         783             :       ! (Yes, they are theoretically orthonormal, but lets make sure)
         784             :       !--------------------------------------------------------------
         785           0 :       do nn=1,I_nbas
         786           0 :         do n2=1,I_nbas
         787           0 :           Bcov(nn,n2) = 0._r8
         788           0 :           do ii=1,I_nlat
         789           0 :             Bcov(nn,n2) = Bcov(nn,n2) + (this%basis(ii,nn)*this%basis(ii,n2)*this%area(ii))
         790             :           end do
         791             :         end do
         792             :       end do
         793             :
         794             :       ! Invert to get the basis amplitude map
         795             :       !--------------------------------------
         796           0 :       call Invert_Matrix(Bcov,I_nbas,this%map)
         797             :
         798             :       ! End Routine
         799             :       !------------
         800           0 :       deallocate(Clats)
         801           0 :       deallocate(Bcoef)
         802           0 :       deallocate(Bcov )
         803             :
         804           0 :     end subroutine init_ZonalProfile
         805             :     !=======================================================================
         806             :
         807             :
         808             :     !=======================================================================
         809           0 :     subroutine calc_ZonalProfile_1Damps(this,I_Zdata,O_Bamp)
         810             :       !
         811             :       ! calc_ZonalProfile_1Damps: Given 1D data values for the nlat zonal
         812             :       !                           profiles gridpoints, compute the zonal
         813             :       !                           profile basis amplitudes.
         814             :       !=====================================================================
         815             :       !
         816             :       ! Passed Variables
         817             :       !------------------
         818             :       class(ZonalProfile_t):: this
         819             :       real(r8),intent(in ):: I_Zdata(:)
         820             :       real(r8),intent(out):: O_Bamp (:)
         821             :       !
         822             :       ! Local Values
         823             :       !--------------
         824           0 :       real(r8),allocatable:: Gcov(:)
         825             :       integer:: ii,nn,n2, astat
         826             :       character(len=*), parameter :: subname = 'calc_ZonalProfile_1Damps'
         827             :
         828             :       ! Compute Covariance with input data and basis functions
         829             :       !--------------------------------------------------------
         830           0 :       allocate(Gcov(this%nbas), stat=astat)
         831           0 :       call handle_allocate_error(astat, subname, 'Gcov')
         832           0 :       do nn=1,this%nbas
         833           0 :         Gcov(nn) = 0._r8
         834           0 :         do ii=1,this%nlat
         835           0 :           Gcov(nn) = Gcov(nn) + (I_Zdata(ii)*this%basis(ii,nn)*this%area(ii))
         836             :         end do
         837             :       end do
         838             :
         839             :       ! Multiply by map to get the amplitudes
         840             :       !-------------------------------------------
         841           0 :       do nn=1,this%nbas
         842           0 :         O_Bamp(nn) = 0._r8
         843           0 :         do n2=1,this%nbas
         844           0 :           O_Bamp(nn) = O_Bamp(nn) + this%map(n2,nn)*Gcov(n2)
         845             :         end do
         846             :       end do
         847             :
         848           0 :       deallocate(Gcov)
         849             :
         850           0 :     end subroutine calc_ZonalProfile_1Damps
         851             :     !=======================================================================
         852             :
         853             :
         854             :     !=======================================================================
         855           0 :     subroutine calc_ZonalProfile_2Damps(this,I_Zdata,O_Bamp)
         856             :       !
         857             :       ! calc_ZonalProfile_2Damps: Given 2D data values for the nlat,nlev zonal
         858             :       !                           profiles gridpoints, compute the zonal
         859             :       !                           profile basis amplitudes.
         860             :       !=====================================================================
         861             :       !
         862             :       ! Passed Variables
         863             :       !------------------
         864             :       class(ZonalProfile_t):: this
         865             :       real(r8),intent(in ):: I_Zdata(:,:)
         866             :       real(r8),intent(out):: O_Bamp (:,:)
         867             :       !
         868             :       ! Local Values
         869             :       !--------------
         870           0 :       real(r8),allocatable:: Gcov(:,:)
         871             :       integer:: ii,nn,n2,ilev
         872             :
         873             :       integer :: nlev, astat
         874             :       character(len=*), parameter :: subname = 'calc_ZonalProfile_2Damps'
         875             :
         876           0 :       nlev = size(I_Zdata,dim=2)
         877             :
         878             :       ! Compute Covariance with input data and basis functions
         879             :       !--------------------------------------------------------
         880           0 :       allocate(Gcov(this%nbas,nlev), stat=astat)
         881           0 :       call handle_allocate_error(astat, subname, 'Gcov')
         882           0 :       do ilev=1,nlev
         883           0 :          do nn=1,this%nbas
         884           0 :             Gcov(nn,ilev) = 0._r8
         885           0 :             do ii=1,this%nlat
         886           0 :                Gcov(nn,ilev) = Gcov(nn,ilev) + (I_Zdata(ii,ilev)*this%basis(ii,nn)*this%area(ii))
         887             :             end do
         888             :          end do
         889             :       end do
         890             :
         891             :       ! Multiply by map to get the amplitudes
         892             :       !-------------------------------------------
         893           0 :       do ilev=1,nlev
         894           0 :          do nn=1,this%nbas
         895           0 :             O_Bamp(nn,ilev) = 0._r8
         896           0 :             do n2=1,this%nbas
         897           0 :                O_Bamp(nn,ilev) = O_Bamp(nn,ilev) + this%map(n2,nn)*Gcov(n2,ilev)
         898             :             end do
         899             :          end do
         900             :       end do
         901           0 :       deallocate(Gcov)
         902             :
         903           0 :     end subroutine calc_ZonalProfile_2Damps
         904             :     !=======================================================================
         905             :
         906             :
         907             :     !=======================================================================
         908           0 :     subroutine eval_ZonalProfile_1Dgrid(this,I_Bamp,O_Zdata)
         909             :       !
         910             :       ! eval_ZonalProfile_1Dgrid: Given the zonal profile basis amplitudes,
         911             :       !                           compute 1D data values for the nlat gridpoints.
         912             :       !=====================================================================
         913             :       !
         914             :       ! Passed Variables
         915             :       !------------------
         916             :       class(ZonalProfile_t):: this
         917             :       real(r8),intent(in ):: I_Bamp (:)
         918             :       real(r8),intent(out):: O_Zdata(:)
         919             :       !
         920             :       ! Local Values
         921             :       !--------------
         922             :       integer:: ii,nn
         923             :
         924             :       ! Construct grid values from basis amplitudes.
         925             :       !--------------------------------------------------
         926           0 :       O_Zdata(1:this%nlat) = 0._r8
         927           0 :       do nn=1,this%nbas
         928           0 :         do ii=1,this%nlat
         929           0 :           O_Zdata(ii) = O_Zdata(ii) + (I_Bamp(nn)*this%basis(ii,nn))
         930             :         end do
         931             :       end do
         932             :
         933           0 :     end subroutine eval_ZonalProfile_1Dgrid
         934             :     !=======================================================================
         935             :
         936             :
         937             :     !=======================================================================
         938           0 :     subroutine eval_ZonalProfile_2Dgrid(this,I_Bamp,O_Zdata)
         939             :       !
         940             :       ! eval_ZonalProfile_2Dgrid: Given the zonal profile basis amplitudes,
         941             :       !                           compute 2D data values for the nlat,nlev gridpoints.
         942             :       !=====================================================================
         943             :       !
         944             :       ! Passed Variables
         945             :       !------------------
         946             :       class(ZonalProfile_t):: this
         947             :       real(r8),intent(in ):: I_Bamp (:,:)
         948             :       real(r8),intent(out):: O_Zdata(:,:)
         949             :       !
         950             :       ! Local Values
         951             :       !--------------
         952             :       integer:: ii,nn,ilev
         953             :
         954             :       integer :: nlev
         955             :
         956           0 :       nlev = size(I_Bamp,dim=2)
         957             :
         958             :       ! Construct grid values from basis amplitudes.
         959             :       !--------------------------------------------------
         960           0 :       O_Zdata(1:this%nlat,1:nlev) = 0._r8
         961           0 :       do nn=1,this%nbas
         962           0 :         do ilev=1,nlev
         963           0 :           do ii=1,this%nlat
         964           0 :             O_Zdata(ii,ilev) = O_Zdata(ii,ilev) + (I_Bamp(nn,ilev)*this%basis(ii,nn))
         965             :           end do
         966             :         end do
         967             :       end do
         968             :
         969           0 :     end subroutine eval_ZonalProfile_2Dgrid
         970             :     !=======================================================================
         971             :
         972             :     !=======================================================================
         973           0 :     subroutine final_ZonalProfile(this)
         974             :       class(ZonalProfile_t) :: this
         975             :
         976           0 :       if(allocated(this%area )) deallocate(this%area)
         977           0 :       if(allocated(this%basis)) deallocate(this%basis)
         978           0 :       if(allocated(this%map  )) deallocate(this%map)
         979             :
         980           0 :     end subroutine final_ZonalProfile
         981             :     !=======================================================================
         982             :
         983             :     !=======================================================================
         984           0 :     subroutine init_ZonalAverage(this,IO_lats,IO_area,I_nlat,GEN_GAUSSLATS)
         985             :       !
         986             :       ! init_ZonalAverage: Initialize the ZonalAverage data structure for the
         987             :       !                    given nlat gridpoints. It is assumed that the domain
         988             :       !                    of these gridpoints of the profile span latitudes
         989             :       !                    from SP to NP.
         990             :       !
         991             :       !                    The optional GEN_GAUSSLATS flag allows for the
         992             :       !                    generation of Gaussian latitudes. The generated grid
         993             :       !                    over-writes the values of IO_lats/IO_area passed by
         994             :       !                    the user.
         995             :       !=====================================================================
         996             :       !
         997             :       ! Passed Variables
         998             :       !------------------
         999             :       class(ZonalAverage_t) :: this
        1000             :       real(r8)     ,intent(inout):: IO_lats(:)
        1001             :       real(r8)     ,intent(inout):: IO_area(:)
        1002             :       integer         ,intent(in):: I_nlat
        1003             :       logical,optional,intent(in):: GEN_GAUSSLATS
        1004             :       !
        1005             :       ! Local Values
        1006             :       !--------------
        1007           0 :       real(r8),allocatable:: Clats (:)
        1008           0 :       real(r8),allocatable:: Glats (:,:)
        1009           0 :       real(r8),allocatable:: BinLat(:)
        1010           0 :       real(r8),allocatable:: Asum  (:,:)
        1011           0 :       real(r8),allocatable:: Anorm (:)
        1012             :       real(r8):: area(pcols),rlat
        1013             :       integer :: nn,jj,ierr, astat
        1014             :       integer :: ncols,lchnk,cc,jlat
        1015             :       integer :: nlcols, count
        1016             :       logical :: generate_lats
        1017             :       character(len=*), parameter :: subname = 'init_ZonalAverage'
        1018             :
        1019           0 :       generate_lats = .false.
        1020             :
        1021           0 :       if (present(GEN_GAUSSLATS)) then
        1022           0 :          generate_lats = GEN_GAUSSLATS
        1023             :       end if
        1024             :
        1025           0 :       nlcols = get_nlcols_p()
        1026             :
        1027             :       ! Allocate space
        1028             :       !-----------------
        1029           0 :       if(allocated(this%area   )) deallocate(this%area)
        1030           0 :       if(allocated(this%a_norm )) deallocate(this%a_norm)
        1031           0 :       if(allocated(this%area_g )) deallocate(this%area_g)
        1032           0 :       if(allocated(this%idx_map)) deallocate(this%idx_map)
        1033             :
        1034           0 :       this%nlat = I_nlat
        1035           0 :       allocate(this%area   (I_nlat), stat=astat)
        1036           0 :       call handle_allocate_error(astat, subname, 'this%area')
        1037           0 :       allocate(this%a_norm (I_nlat), stat=astat)
        1038           0 :       call handle_allocate_error(astat, subname, 'this%a_norm')
        1039           0 :       allocate(this%area_g (pcols,begchunk:endchunk), stat=astat)
        1040           0 :       call handle_allocate_error(astat, subname, 'this%area_g')
        1041           0 :       allocate(this%idx_map(pcols,begchunk:endchunk), stat=astat)
        1042           0 :       call handle_allocate_error(astat, subname, 'this%idx_map')
        1043             :
        1044           0 :       allocate(Clats (I_nlat), stat=astat)
        1045           0 :       call handle_allocate_error(astat, subname, 'Clats')
        1046           0 :       allocate(BinLat(I_nlat+1), stat=astat)
        1047           0 :       call handle_allocate_error(astat, subname, 'BinLat')
        1048           0 :       allocate(Glats (pcols,begchunk:endchunk), stat=astat)
        1049           0 :       call handle_allocate_error(astat, subname, 'Glats')
        1050           0 :       allocate(Asum  (nlcols,I_nlat), stat=astat)
        1051           0 :       call handle_allocate_error(astat, subname, 'Asum')
        1052           0 :       allocate(Anorm (I_nlat), stat=astat)
        1053           0 :       call handle_allocate_error(astat, subname, 'Anorm')
        1054             :
        1055             :       ! Optionally create the Latitude Gridpoints
        1056             :       ! and their associated area weights. Otherwise
        1057             :       ! they need to be supplied by the user.
        1058             :       !-----------------------------------------------
        1059           0 :       if(generate_lats) then
        1060             :
        1061             :         ! Create a Gaussin grid from SP to NP
        1062             :         !--------------------------------------
        1063           0 :         call sh_create_gaus_grid(this%nlat,Clats,IO_area,ierr)
        1064           0 :         if (ierr/=0) then
        1065           0 :            call endrun('init_ZonalAverage: Error creating Gaussian grid')
        1066             :         end if
        1067             :
        1068             :         ! Convert generated colatitudes SP->NP to Lats and convert
        1069             :         ! to degrees and scale the area for global 2D integrals
        1070             :         !-----------------------------------------------------------
        1071           0 :         do nn=1,this%nlat
        1072           0 :           IO_lats(nn) = (45._r8*Clats(nn)/qrtrPI) - 90._r8
        1073           0 :           IO_area(nn) = IO_area(nn)*twoPI
        1074             :         end do
        1075             :       else
        1076             :         ! Convert Latitudes to SP->NP colatitudes in radians
        1077             :         !----------------------------------------------------
        1078           0 :         do nn=1,this%nlat
        1079           0 :           Clats(nn) = (IO_lats(nn) + 90._r8)*qrtrPI/45._r8
        1080             :         end do
        1081             :       endif
        1082             :
        1083             :       ! Copy the Lat grid area weights to the data structure
        1084             :       !-----------------------------------------------------
        1085           0 :       this%area(1:this%nlat) = IO_area(1:this%nlat)
        1086             :
        1087             :       ! Save a copy of the area weights for each 2D gridpoint
        1088             :       ! and convert Latitudes to SP->NP colatitudes in radians
        1089             :       !-------------------------------------------------------
        1090           0 :       do lchnk=begchunk,endchunk
        1091           0 :         ncols = get_ncols_p(lchnk)
        1092           0 :         call get_wght_all_p(lchnk, ncols, area)
        1093           0 :         do cc = 1,ncols
        1094           0 :           rlat=get_rlat_p(lchnk,cc)
        1095           0 :           this%area_g(cc,lchnk) = area(cc)
        1096           0 :           Glats      (cc,lchnk) = rlat + halfPI
        1097             :         end do
        1098             :       end do
        1099             :
        1100             :       ! Set boundaries for Latitude bins
        1101             :       !-----------------------------------
        1102           0 :       BinLat(1)           = 0._r8
        1103           0 :       BinLat(this%nlat+1) = pi
        1104           0 :       do nn=2,this%nlat
        1105           0 :         BinLat(nn) = (Clats(nn-1)+Clats(nn))/2._r8
        1106             :       end do
        1107             :
        1108             :       ! Loop over 2D gridpoints and determine its lat bin index
        1109             :       !---------------------------------------------------------
        1110           0 :       do lchnk=begchunk,endchunk
        1111           0 :         ncols = get_ncols_p(lchnk)
        1112           0 :         do cc = 1,ncols
        1113           0 :           jlat = -1
        1114           0 :           if((Glats(cc,lchnk)<=BinLat(2)).and. &
        1115             :              (Glats(cc,lchnk)>=BinLat(1))      ) then
        1116           0 :             jlat = 1
        1117           0 :           elseif((Glats(cc,lchnk)>=BinLat(this%nlat)  ).and. &
        1118           0 :                  (Glats(cc,lchnk)<=BinLat(this%nlat+1))      ) then
        1119           0 :             jlat = this%nlat
        1120             :           else
        1121           0 :             do jj=2,(this%nlat-1)
        1122           0 :               if((Glats(cc,lchnk)>BinLat(jj  )).and. &
        1123           0 :                  (Glats(cc,lchnk)<=BinLat(jj+1))      ) then
        1124           0 :                 jlat = jj
        1125           0 :                 exit
        1126             :               endif
        1127             :             end do
        1128             :           endif
        1129           0 :           if (jlat<1) then
        1130           0 :             call endrun('ZonalAverage init ERROR: jlat not in range')
        1131             :           endif
        1132           0 :           this%idx_map(cc,lchnk) = jlat
        1133             :         end do
        1134             :       end do
        1135             :
        1136             :       ! Initialize 2D Area sums for each bin
        1137             :       !--------------------------------------
        1138           0 :       Asum(:,:) = 0._r8
        1139           0 :       Anorm(:) = 0._r8
        1140           0 :       count = 0
        1141           0 :       do lchnk=begchunk,endchunk
        1142           0 :         ncols = get_ncols_p(lchnk)
        1143           0 :         do cc = 1,ncols
        1144           0 :           jlat = this%idx_map(cc,lchnk)
        1145           0 :           count=count+1
        1146           0 :           Asum(count,jlat) = this%area_g(cc,lchnk)
        1147             :         end do
        1148             :       end do
        1149             :
        1150           0 :       call shr_reprosum_calc(Asum, Anorm, count, nlcols, I_nlat, gbl_count=ngcols_p, commid=mpicom)
        1151             :
        1152           0 :       this%a_norm = Anorm
        1153             :
        1154           0 :       if (.not.all(Anorm(:)>0._r8)) then
        1155           0 :          write(iulog,*) 'init_ZonalAverage -- ERROR in Anorm values: '
        1156           0 :          do jlat = 1,I_nlat
        1157           0 :             if (.not.Anorm(jlat)>0._r8) then
        1158           0 :                write(iulog,*) ' Anorm(',jlat,'): ', Anorm(jlat)
        1159             :             endif
        1160             :          end do
        1161           0 :          call endrun('init_ZonalAverage -- ERROR in Anorm values')
        1162             :       end if
        1163             :
        1164             :       ! End Routine
        1165             :       !------------
        1166           0 :       deallocate(Clats)
        1167           0 :       deallocate(BinLat)
        1168           0 :       deallocate(Glats)
        1169           0 :       deallocate(Asum)
        1170           0 :       deallocate(Anorm)
        1171             :
        1172           0 :     end subroutine init_ZonalAverage
        1173             :     !=======================================================================
        1174             :
        1175             :
        1176             :     !=======================================================================
        1177           0 :     subroutine calc_ZonalAverage_2DbinAvg(this,I_Gdata,O_Zdata)
        1178             :       !
        1179             :       ! calc_ZonalAverage_2DbinAvg: Given 2D data values for ncol gridpoints,
        1180             :       !                             compute the nlat area weighted binAvg profile
        1181             :       !=====================================================================
        1182             :       !
        1183             :       ! Passed Variables
        1184             :       !------------------
        1185             :       class(ZonalAverage_t):: this
        1186             :       real(r8),intent(in ):: I_Gdata(pcols,begchunk:endchunk)
        1187             :       real(r8),intent(out):: O_Zdata(:)
        1188             :       !
        1189             :       ! Local Values
        1190             :       !--------------
        1191           0 :       real(r8),allocatable:: Asum (:,:)
        1192             :       integer:: nn,ncols,lchnk,cc,jlat
        1193             :       integer :: nlcols, count, astat
        1194             :       character(len=*), parameter :: subname = 'calc_ZonalAverage_2DbinAvg'
        1195             :
        1196           0 :       nlcols = get_nlcols_p()
        1197             :
        1198             :
        1199             :       ! Initialize Zonal profile
        1200             :       !---------------------------
        1201           0 :       allocate(Asum(nlcols,this%nlat), stat=astat)
        1202           0 :       call handle_allocate_error(astat, subname, 'Asum')
        1203           0 :       Asum(:,:) = 0._r8
        1204             :
        1205           0 :       O_Zdata(1:this%nlat) = 0._r8
        1206             :
        1207             :       ! Compute area-weighted sums
        1208             :       !-----------------------------
        1209           0 :       count = 0
        1210           0 :       do lchnk=begchunk,endchunk
        1211           0 :         ncols = get_ncols_p(lchnk)
        1212           0 :         do cc = 1,ncols
        1213           0 :           jlat = this%idx_map(cc,lchnk)
        1214           0 :           count=count+1
        1215           0 :           Asum(count,jlat) = I_Gdata(cc,lchnk)*this%area_g(cc,lchnk)
        1216             :         end do
        1217             :       end do
        1218             :
        1219           0 :       call shr_reprosum_calc(Asum,O_Zdata,count, nlcols, this%nlat,gbl_count=ngcols_p, commid=mpicom)
        1220             :
        1221             :       ! Divide by area norm to get the averages
        1222             :       !-----------------------------------------
        1223           0 :       do nn=1,this%nlat
        1224           0 :         O_Zdata(nn) = O_Zdata(nn)/this%a_norm(nn)
        1225             :       end do
        1226             :
        1227           0 :       deallocate(Asum)
        1228             :
        1229           0 :     end subroutine calc_ZonalAverage_2DbinAvg
        1230             :     !=======================================================================
        1231             :
        1232             :
        1233             :     !=======================================================================
        1234           0 :     subroutine calc_ZonalAverage_3DbinAvg(this,I_Gdata,O_Zdata)
        1235             :       !
        1236             :       ! calc_ZonalAverage_3DbinAvg: Given 3D data values for ncol,nlev gridpoints,
        1237             :       !                             compute the nlat,nlev area weighted binAvg profile
        1238             :       !=====================================================================
        1239             :       !
        1240             :       ! Passed Variables
        1241             :       !------------------
        1242             :       class(ZonalAverage_t):: this
        1243             :       real(r8),intent(in ):: I_Gdata(:,:,begchunk:)
        1244             :       real(r8),intent(out):: O_Zdata(:,:)
        1245             :       !
        1246             :       ! Local Values
        1247             :       !--------------
        1248           0 :       real(r8),allocatable:: Gsum(:)
        1249           0 :       real(r8),allocatable:: Asum(:,:)
        1250             :       integer:: nn,ncols,lchnk,cc,jlat
        1251             :       integer:: Nsum,ilev,ns
        1252             :
        1253             :       integer :: nlev
        1254             :       integer :: nlcols, count, astat
        1255             :       character(len=*), parameter :: subname = 'calc_ZonalAverage_3DbinAvg'
        1256             :
        1257           0 :       nlev = size(I_Gdata,dim=2)
        1258           0 :       nlcols = get_nlcols_p()
        1259             :
        1260             :       ! Initialize Zonal profile
        1261             :       !---------------------------
        1262           0 :       Nsum = this%nlat*nlev
        1263           0 :       allocate(Gsum(Nsum), stat=astat)
        1264           0 :       call handle_allocate_error(astat, subname, 'Gsum')
        1265           0 :       allocate(Asum(nlcols,Nsum), stat=astat)
        1266           0 :       call handle_allocate_error(astat, subname, 'Asum')
        1267           0 :       Asum(:,:) = 0._r8
        1268             :
        1269           0 :       O_Zdata(1:this%nlat,1:nlev) = 0._r8
        1270             :
        1271             :       ! Compute area-weighted sums
        1272             :       !-----------------------------
        1273           0 :       do ilev = 1,nlev
        1274           0 :          count = 0
        1275           0 :          do lchnk=begchunk,endchunk
        1276           0 :             ncols = get_ncols_p(lchnk)
        1277           0 :             do cc = 1,ncols
        1278           0 :                jlat = this%idx_map(cc,lchnk)
        1279           0 :                ns = jlat + (ilev-1)*this%nlat
        1280           0 :                count=count+1
        1281           0 :                Asum(count,ns) = I_Gdata(cc,ilev,lchnk)*this%area_g(cc,lchnk)
        1282             :             end do
        1283             :          end do
        1284             :       end do
        1285             :
        1286           0 :       call shr_reprosum_calc(Asum,Gsum, count, nlcols, Nsum, gbl_count=ngcols_p, commid=mpicom)
        1287             :
        1288             :       ! Divide by area norm to get the averages
        1289             :       !-----------------------------------------
        1290           0 :       do ilev = 1,nlev
        1291           0 :          do nn = 1,this%nlat
        1292           0 :             ns = nn + (ilev-1)*this%nlat
        1293           0 :             O_Zdata(nn,ilev) = Gsum(ns)/this%a_norm(nn)
        1294             :          end do
        1295             :       end do
        1296             :
        1297           0 :       deallocate(Gsum)
        1298           0 :       deallocate(Asum)
        1299             :
        1300           0 :     end subroutine calc_ZonalAverage_3DbinAvg
        1301             :     !=======================================================================
        1302             :
        1303             :     !=======================================================================
        1304        1536 :     subroutine final_ZonalAverage(this)
        1305             :       class(ZonalAverage_t) :: this
        1306             :
        1307        1536 :       if(allocated(this%area   )) deallocate(this%area)
        1308        1536 :       if(allocated(this%a_norm )) deallocate(this%a_norm)
        1309        1536 :       if(allocated(this%area_g )) deallocate(this%area_g)
        1310        1536 :       if(allocated(this%idx_map)) deallocate(this%idx_map)
        1311             :
        1312        1536 :     end subroutine final_ZonalAverage
        1313             :     !=======================================================================
        1314             :
        1315             :
        1316             :     !=======================================================================
        1317           0 :     subroutine Invert_Matrix(I_Mat,Nbas,O_InvMat)
        1318             :       !
        1319             :       ! Invert_Matrix: Given the NbasxNbas matrix, calculate and return
        1320             :       !                the inverse of the matrix.
        1321             :       !
        1322             :       ! Implemented with the LAPACK DGESV routine.
        1323             :       !
        1324             :       !====================================================================
        1325             :       !
        1326             :       ! Passed Variables
        1327             :       !------------------
        1328             :       real(r8), intent(inout) :: I_Mat(:,:)  ! input matrix contains P*L*U
        1329             :                                              ! decomposition on output
        1330             :       integer,  intent(in)    :: Nbas
        1331             :       real(r8), intent(out)   :: O_InvMat(:,:)
        1332             :       !
        1333             :       ! Local Values
        1334             :       !-------------
        1335           0 :       integer, allocatable :: Indx(:)  ! pivot indices
        1336             :       integer :: astat, ii
        1337             :       character(len=*), parameter :: subname = 'Invert_Matrix'
        1338             :       character(len=80) :: msg
        1339             :
        1340             :       external DGESV
        1341             :
        1342             :       ! Allocate work space
        1343             :       !---------------------
        1344           0 :       allocate(Indx(Nbas), stat=astat)
        1345           0 :       call handle_allocate_error(astat, subname, 'Indx')
        1346             :
        1347             :       ! Initialize the inverse array with the identity matrix
        1348             :       !-------------------------------------------------------
        1349           0 :       O_InvMat(:,:) = 0._r8
        1350           0 :       do ii=1,Nbas
        1351           0 :         O_InvMat(ii,ii) = 1._r8
        1352             :       end do
        1353             :
        1354           0 :       call DGESV(Nbas, Nbas, I_Mat, Nbas, Indx, O_InvMat, Nbas, astat)
        1355             :
        1356           0 :       if (astat < 0) then
        1357           0 :          write(msg, '(a, i1, a)') 'argument # ', abs(astat), ' has an illegal value'
        1358           0 :          call endrun(subname//': DGESV error return: '//msg)
        1359           0 :       else if (astat > 0) then
        1360           0 :          call endrun(subname//': DGESV error return: matrix is singular')
        1361             :       end if
        1362             :
        1363           0 :       deallocate(Indx)
        1364             :
        1365           0 :     end subroutine Invert_Matrix
        1366             :     !=======================================================================
        1367             :
        1368             :     !=======================================================================
        1369             :     ! legacy spherepack routines
        1370             :     !=======================================================================
        1371           0 :     subroutine sh_gen_basis_coefs(nn,mm,cp)
        1372             :       !
        1373             :       ! spherepack alfk
        1374             :       !
        1375             :       ! dimension of           real cp(nn/2 + 1)
        1376             :       ! arguments
        1377             :       !
        1378             :       ! purpose                computes fourier coefficients in the trigonometric series
        1379             :       !                        representation of the normalized associated
        1380             :       !                        legendre function pbar(nn,mm,theta) for use by
        1381             :       !                        sh_gen_basis_coefs in calculating pbar(nn,mm,theta).
        1382             :       !
        1383             :       !                        first define the normalized associated
        1384             :       !                        legendre functions
        1385             :       !
        1386             :       !                        pbar(mm,nn,theta) = sqrt((2*nn+1)*factorial(nn-mm)
        1387             :       !                        /(2*factorial(nn+mm)))*sin(theta)**mm/(2**nn*
        1388             :       !                        factorial(nn)) times the (nn+mm)th derivative of
        1389             :       !                        (x**2-1)**nn with respect to x=cos(theta)
        1390             :       !
        1391             :       !                        where theta is colatitude.
        1392             :       !
        1393             :       !                        then subroutine sh_gen_basis_coefs computes the coefficients
        1394             :       !                        cp(k) in the following trigonometric
        1395             :       !                        expansion of pbar(m,n,theta).
        1396             :       !
        1397             :       !                        1) for n even and m even, pbar(mm,nn,theta) =
        1398             :       !                           .5*cp(1) plus the sum from k=1 to k=nn/2
        1399             :       !                           of cp(k+1)*cos(2*k*th)
        1400             :       !
        1401             :       !                        2) for nn even and mm odd, pbar(mm,nn,theta) =
        1402             :       !                           the sum from k=1 to k=nn/2 of
        1403             :       !                           cp(k)*sin(2*k*th)
        1404             :       !
        1405             :       !                        3) for n odd and m even, pbar(mm,nn,theta) =
        1406             :       !                           the sum from k=1 to k=(nn+1)/2 of
        1407             :       !                           cp(k)*cos((2*k-1)*th)
        1408             :       !
        1409             :       !                        4) for nn odd and mm odd,  pbar(mm,nn,theta) =
        1410             :       !                           the sum from k=1 to k=(nn+1)/2 of
        1411             :       !                           cp(k)*sin((2*k-1)*th)
        1412             :       !
        1413             :       ! arguments
        1414             :       !
        1415             :       ! on input               nn
        1416             :       !                          nonnegative integer specifying the degree of
        1417             :       !                          pbar(nn,mm,theta)
        1418             :       !
        1419             :       !                        mm
        1420             :       !                          is the order of pbar(nn,mm,theta). mm can be
        1421             :       !                          any integer however cp is computed such that
        1422             :       !                          pbar(nn,mm,theta) = 0 if abs(m) is greater
        1423             :       !                          than nn and pbar(nn,mm,theta) = (-1)**mm*
        1424             :       !                          pbar(nn,-mm,theta) for negative mm.
        1425             :       !
        1426             :       ! on output              cp
        1427             :       !                          array of length (nn/2)+1
        1428             :       !                          which contains the fourier coefficients in
        1429             :       !                          the trigonometric series representation of
        1430             :       !                          pbar(nn,mm,theta)
        1431             :       !
        1432             :       ! special conditions     none
        1433             :       !
        1434             :       ! algorithm              the highest order coefficient is determined in
        1435             :       !                        closed form and the remainig coefficients are
        1436             :       !                        determined as the solution of a backward
        1437             :       !                        recurrence relation.
        1438             :       !
        1439             :       !=====================================================================
        1440             :       !
        1441             :       ! Passed Variables
        1442             :       !------------------
        1443             :       integer ,intent(in ):: nn
        1444             :       integer ,intent(in ):: mm
        1445             :       real(r8),intent(out):: cp(nn/2+1)
        1446             :       !
        1447             :       ! Local Values
        1448             :       !----------------
        1449             :       real(r8):: fnum,fnmh
        1450             :       real(r8):: pm1
        1451             :       real(r8):: t1,t2
        1452             :       real(r8):: fden
        1453             :       real(r8):: cp2
        1454             :       real(r8):: fnnp1
        1455             :       real(r8):: fnmsq
        1456             :       real(r8):: fk
        1457             :       real(r8):: a1,b1,C1
        1458             :       integer :: ma,nmms2,nex
        1459             :       integer :: ii,jj
        1460             :
        1461             :       real(r8),parameter:: SC10=1024._r8
        1462             :       real(r8),parameter:: SC20=SC10*SC10
        1463             :       real(r8),parameter:: SC40=SC20*SC20
        1464             :
        1465           0 :       cp(1) = 0._r8
        1466           0 :       ma = abs(mm)
        1467           0 :       if(ma>nn) return
        1468             :
        1469           0 :       if((nn-1)<0) then
        1470           0 :         cp(1) = sqrt(2._r8)
        1471           0 :         return
        1472           0 :       elseif((nn-1)==0) then
        1473           0 :         if(ma/=0) then
        1474           0 :           cp(1) = sqrt(.75_r8)
        1475           0 :           if(mm==-1) cp(1) = -cp(1)
        1476             :         else
        1477           0 :           cp(1) = sqrt(1.5_r8)
        1478             :         endif
        1479           0 :         return
        1480             :       else
        1481           0 :         if(mod(nn+ma,2)/=0) then
        1482           0 :           nmms2 = (nn-ma-1)/2
        1483           0 :           fnum  = nn + ma + 2
        1484           0 :           fnmh  = nn - ma + 2
        1485           0 :           pm1   = -1._r8
        1486             :         else
        1487           0 :           nmms2 = (nn-ma)/2
        1488           0 :           fnum  = nn + ma + 1
        1489           0 :           fnmh  = nn - ma + 1
        1490           0 :           pm1   = 1._r8
        1491             :         endif
        1492             :       endif
        1493             :
        1494           0 :       t1   = 1._r8/SC20
        1495           0 :       nex  = 20
        1496           0 :       fden = 2._r8
        1497           0 :       if(nmms2>=1) then
        1498           0 :         do ii = 1,nmms2
        1499           0 :           t1 = fnum*t1/fden
        1500           0 :           if (t1>SC20) then
        1501           0 :             t1  = t1/SC40
        1502           0 :             nex = nex + 40
        1503             :           endif
        1504           0 :           fnum = fnum + 2._r8
        1505           0 :           fden = fden + 2._r8
        1506             :         end do
        1507             :       endif
        1508             :
        1509           0 :       if(mod(ma/2,2)/=0) then
        1510           0 :         t1 = -t1/2._r8**(nn-1-nex)
        1511             :       else
        1512           0 :         t1 =  t1/2._r8**(nn-1-nex)
        1513             :       endif
        1514           0 :       t2 = 1._r8
        1515           0 :       if(ma/=0) then
        1516           0 :         do ii = 1,ma
        1517           0 :           t2   = fnmh*t2/ (fnmh+pm1)
        1518           0 :           fnmh = fnmh + 2._r8
        1519             :         end do
        1520             :       endif
        1521             :
        1522           0 :       cp2   = t1*sqrt((nn+.5_r8)*t2)
        1523           0 :       fnnp1 = nn*(nn+1)
        1524           0 :       fnmsq = fnnp1 - 2._r8*ma*ma
        1525             :
        1526           0 :       if((mod(nn,2)==0).and.(mod(ma,2)==0)) then
        1527           0 :         jj = 1+(nn+1)/2
        1528             :       else
        1529           0 :         jj = (nn+1)/2
        1530             :       endif
        1531             :
        1532           0 :       cp(jj) = cp2
        1533           0 :       if(mm<0) then
        1534           0 :         if(mod(ma,2)/=0) cp(jj) = -cp(jj)
        1535             :       endif
        1536           0 :       if(jj<=1) return
        1537             :
        1538           0 :       fk = nn
        1539           0 :       a1 = (fk-2._r8)*(fk-1._r8) - fnnp1
        1540           0 :       b1 = 2._r8* (fk*fk-fnmsq)
        1541           0 :       cp(jj-1) = b1*cp(jj)/a1
        1542             :
        1543           0 :       jj = jj - 1
        1544           0 :       do while(jj>1)
        1545           0 :         fk = fk - 2._r8
        1546           0 :         a1 = (fk-2._r8)*(fk-1._r8) - fnnp1
        1547           0 :         b1 = -2._r8*(fk*fk-fnmsq)
        1548           0 :         c1 = (fk+1._r8)*(fk+2._r8) - fnnp1
        1549           0 :         cp(jj-1) = -(b1*cp(jj)+c1*cp(jj+1))/a1
        1550           0 :         jj = jj - 1
        1551             :      end do
        1552             :
        1553             :     end subroutine sh_gen_basis_coefs
        1554             :     !=======================================================================
        1555             :
        1556             :     !=======================================================================
        1557           0 :     subroutine sh_create_basis(nn,mm,theta,cp,pb)
        1558             :       !
        1559             :       ! spherepack lfpt
        1560             :       !
        1561             :       ! dimension of
        1562             :       ! arguments
        1563             :       !                        cp((nn/2)+1)
        1564             :       !
        1565             :       ! purpose                routine sh_create_basis uses coefficients computed by
        1566             :       !                        routine sh_gen_basis_coefs to compute the
        1567             :       !                        normalized associated legendre function pbar(nn,mm,theta)
        1568             :       !                        at colatitude theta.
        1569             :       !
        1570             :       ! arguments
        1571             :       !
        1572             :       ! on input               nn
        1573             :       !                          nonnegative integer specifying the degree of
        1574             :       !                          pbar(nn,mm,theta)
        1575             :       !                        mm
        1576             :       !                          is the order of pbar(nn,mm,theta). mm can be
        1577             :       !                          any integer however pbar(nn,mm,theta) = 0
        1578             :       !                          if abs(mm) is greater than nn and
        1579             :       !                          pbar(nn,mm,theta) = (-1)**mm*pbar(nn,-mm,theta)
        1580             :       !                          for negative mm.
        1581             :       !
        1582             :       !                        theta
        1583             :       !                          colatitude in radians
        1584             :       !
        1585             :       !                        cp
        1586             :       !                          array of length (nn/2)+1
        1587             :       !                          containing coefficients computed by routine
        1588             :       !                          sh_gen_basis_coefs
        1589             :       !
        1590             :       ! on output              pb
        1591             :       !                          variable containing pbar(n,m,theta)
        1592             :       !
        1593             :       ! special conditions     calls to routine sh_create_basis must be preceded by an
        1594             :       !                        appropriate call to routine sh_gen_basis_coefs.
        1595             :       !
        1596             :       ! algorithm              the trigonometric series formula used by
        1597             :       !                        routine sh_create_basis to calculate pbar(nn,mm,theta) at
        1598             :       !                        colatitude theta depends on mm and nn as follows:
        1599             :       !
        1600             :       !                           1) for nn even and mm even, the formula is
        1601             :       !                              .5*cp(1) plus the sum from k=1 to k=n/2
        1602             :       !                              of cp(k)*cos(2*k*theta)
        1603             :       !                           2) for nn even and mm odd. the formula is
        1604             :       !                              the sum from k=1 to k=nn/2 of
        1605             :       !                              cp(k)*sin(2*k*theta)
        1606             :       !                           3) for nn odd and mm even, the formula is
        1607             :       !                              the sum from k=1 to k=(nn+1)/2 of
        1608             :       !                              cp(k)*cos((2*k-1)*theta)
        1609             :       !                           4) for nn odd and mm odd, the formula is
        1610             :       !                              the sum from k=1 to k=(nn+1)/2 of
        1611             :       !                              cp(k)*sin((2*k-1)*theta)
        1612             :       !
        1613             :       !=====================================================================
        1614             :       integer, intent(in) :: nn,mm
        1615             :       real(r8), intent(in) :: theta
        1616             :       real(r8), intent(in) :: cp(:)
        1617             :       real(r8), intent(out) :: pb
        1618             :
        1619             :       real(r8) :: cdt
        1620             :       real(r8) :: sdt
        1621             :       real(r8) :: ct
        1622             :       real(r8) :: st
        1623             :       real(r8) :: summ
        1624             :       real(r8) :: cth
        1625             :
        1626             :       integer:: ma,nmod,mmod,kdo
        1627             :       integer:: kp1,kk
        1628             :
        1629           0 :       pb = 0._r8
        1630           0 :       ma = abs(mm)
        1631           0 :       if(ma>nn) return
        1632             :
        1633           0 :       if(nn<=0) then
        1634           0 :         if(ma<=0) then
        1635           0 :           pb = sqrt(.5_r8)
        1636           0 :           return
        1637             :         endif
        1638             :       endif
        1639             :
        1640           0 :       nmod = mod(nn,2)
        1641           0 :       mmod = mod(ma,2)
        1642             :
        1643           0 :       if(nmod<=0) then
        1644           0 :         if(mmod<=0) then
        1645           0 :           kdo = nn/2 + 1
        1646           0 :           cdt = cos(theta+theta)
        1647           0 :           sdt = sin(theta+theta)
        1648           0 :           ct  = 1._r8
        1649           0 :           st  = 0._r8
        1650           0 :           summ = .5_r8*cp(1)
        1651           0 :           do kp1 = 2,kdo
        1652           0 :             cth = cdt*ct - sdt*st
        1653           0 :             st  = sdt*ct + cdt*st
        1654           0 :             ct  = cth
        1655           0 :             summ = summ + cp(kp1)*ct
        1656             :           end do
        1657           0 :           pb = summ
        1658           0 :           return
        1659             :         endif
        1660           0 :         kdo = nn/2
        1661           0 :         cdt = cos(theta+theta)
        1662           0 :         sdt = sin(theta+theta)
        1663           0 :         ct  = 1._r8
        1664           0 :         st  = 0._r8
        1665           0 :         summ = 0._r8
        1666           0 :         do kk = 1,kdo
        1667           0 :           cth = cdt*ct - sdt*st
        1668           0 :           st  = sdt*ct + cdt*st
        1669           0 :           ct  = cth
        1670           0 :           summ = summ + cp(kk)*st
        1671             :         end do
        1672           0 :         pb = summ
        1673           0 :         return
        1674             :       endif
        1675             :
        1676           0 :       kdo = (nn+1)/2
        1677           0 :       if(mmod<=0) then
        1678           0 :         cdt =  cos(theta+theta)
        1679           0 :         sdt =  sin(theta+theta)
        1680           0 :         ct  =  cos(theta)
        1681           0 :         st  = -sin(theta)
        1682           0 :         summ = 0._r8
        1683           0 :         do kk = 1,kdo
        1684           0 :           cth = cdt*ct - sdt*st
        1685           0 :           st  = sdt*ct + cdt*st
        1686           0 :           ct  = cth
        1687           0 :           summ = summ + cp(kk)*ct
        1688             :         end do
        1689           0 :         pb = summ
        1690           0 :         return
        1691             :       endif
        1692             :
        1693           0 :       cdt =  cos(theta+theta)
        1694           0 :       sdt =  sin(theta+theta)
        1695           0 :       ct  =  cos(theta)
        1696           0 :       st  = -sin(theta)
        1697           0 :       summ = 0._r8
        1698           0 :       do kk = 1,kdo
        1699           0 :         cth = cdt*ct - sdt*st
        1700           0 :         st  = sdt*ct + cdt*st
        1701           0 :         ct  = cth
        1702           0 :         summ = summ + cp(kk)*st
        1703             :       end do
        1704           0 :       pb = summ
        1705             :
        1706             :     end subroutine sh_create_basis
        1707             :     !=======================================================================
        1708             :
        1709             :     !=======================================================================
        1710           0 :     subroutine sh_create_gaus_grid(nlat,theta,wts,ierr)
        1711             :       !
        1712             :       ! spherepack gaqd
        1713             :       !  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
        1714             :       !  .                                                             .
        1715             :       !  .                  copyright (c) 2001 by ucar                 .
        1716             :       !  .                                                             .
        1717             :       !  .       university corporation for atmospheric research       .
        1718             :       !  .                                                             .
        1719             :       !  .                      all rights reserved                    .
        1720             :       !  .                                                             .
        1721             :       !  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
        1722             :       !
        1723             :       !                             February 2002
        1724             :       !
        1725             :       !     gauss points and weights are computed using the fourier-newton
        1726             :       !     described in "on computing the points and weights for
        1727             :       !     gauss-legendre quadrature", paul n. swarztrauber, siam journal
        1728             :       !     on scientific computing (DOI 10.1137/S1064827500379690).
        1729             :       !     This routine is faster and more accurate than older program
        1730             :       !     with the same name.
        1731             :       !
        1732             :       !     computes the nlat gaussian colatitudes and weights
        1733             :       !     in double precision. the colatitudes are in radians and lie in the
        1734             :       !     in the interval (0,pi).
        1735             :       !
        1736             :       !     input parameters
        1737             :       !
        1738             :       !     nlat    the number of gaussian colatitudes in the interval (0,pi)
        1739             :       !             (between the two poles).  nlat must be greater than zero.
        1740             :       !
        1741             :       !     output parameters
        1742             :       !
        1743             :       !     theta   a double precision array with length nlat
        1744             :       !             containing the gaussian colatitudes in
        1745             :       !             increasing radians on the interval (0,pi).
        1746             :       !
        1747             :       !     wts     a double precision array with lenght nlat
        1748             :       !             containing the gaussian weights.
        1749             :       !
        1750             :       !     ierror = 0 no errors
        1751             :       !            = 1 if nlat<=0
        1752             :       !
        1753             :       !===================================================================
        1754             :       !
        1755             :       ! Passed variables
        1756             :       !-----------------
        1757             :       integer ,intent(in ) :: nlat
        1758             :       real(r8),intent(out) :: theta(nlat)
        1759             :       real(r8),intent(out) :: wts(nlat)
        1760             :       integer ,intent(out) :: ierr
        1761             :       !
        1762             :       ! Local Values
        1763             :       !-------------
        1764             :       real(r8):: sgnd
        1765             :       real(r8):: xx,dtheta,dthalf
        1766             :       real(r8):: cmax,zprev,zlast,zero,zhold,pb,dpb,dcor,summ,cz
        1767             :       integer :: mnlat,ns2,nhalf,nix,it,ii
        1768             :
        1769             :       real(r8), parameter :: eps = epsilon(1._r8)
        1770             :
        1771             :       ! check work space length
        1772             :       !------------------------
        1773           0 :       if(nlat<=0) then
        1774           0 :         ierr = 1
        1775           0 :         return
        1776             :       endif
        1777           0 :       ierr = 0
        1778             :
        1779             :       ! compute weights and points analytically when nlat=1,2
        1780             :       !-------------------------------------------------------
        1781           0 :       if(nlat==1) then
        1782           0 :         theta(1) = acos(0._r8)
        1783           0 :         wts  (1) = 2._r8
        1784           0 :         return
        1785           0 :       elseif(nlat==2) then
        1786           0 :         xx       = sqrt(1._r8/3._r8)
        1787           0 :         theta(1) = acos( xx)
        1788           0 :         theta(2) = acos(-xx)
        1789           0 :         wts  (1) = 1._r8
        1790           0 :         wts  (2) = 1._r8
        1791           0 :         return
        1792             :       endif
        1793             :
        1794             :       ! Proceed for nlat > 2
        1795             :       !----------------------
        1796           0 :       mnlat = mod(nlat,2)
        1797           0 :       ns2   = nlat/2
        1798           0 :       nhalf = (nlat+1)/2
        1799             :
        1800           0 :       call sh_fourier_coefs_dp(nlat,cz,theta(ns2+1),wts(ns2+1))
        1801             :
        1802           0 :       dtheta = halfPI/nhalf
        1803           0 :       dthalf = dtheta/2._r8
        1804           0 :       cmax   = .2_r8*dtheta
        1805             :
        1806             :       ! estimate first point next to theta = pi/2
        1807             :       !-------------------------------------------
        1808           0 :       if(mnlat/=0) then
        1809           0 :         zero  = halfPI - dtheta
        1810           0 :         zprev = halfPI
        1811           0 :         nix   = nhalf - 1
        1812             :       else
        1813           0 :         zero = halfPI - dthalf
        1814           0 :         nix  = nhalf
        1815             :       endif
        1816             :
        1817           0 :       do while(nix/=0)
        1818             :          dcor = huge(1._r8)
        1819             :          it = 0
        1820           0 :          do while (abs(dcor) > eps*abs(zero))
        1821             :             it = it + 1
        1822             :             ! newton iterations
        1823             :             !-----------------------
        1824           0 :             call sh_legp_dlegp_theta(nlat,zero,cz,theta(ns2+1),wts(ns2+1),pb,dpb)
        1825           0 :             dcor = pb/dpb
        1826           0 :             if(dcor.ne.0._r8) then
        1827           0 :                sgnd = dcor/abs(dcor)
        1828             :             else
        1829             :                sgnd = 1._r8
        1830             :             endif
        1831           0 :             dcor = sgnd*min(abs(dcor),cmax)
        1832           0 :             zero = zero - dcor
        1833             :          end do
        1834             :
        1835           0 :          theta(nix) = zero
        1836           0 :          zhold      = zero
        1837             :
        1838             :          !   wts(nix) = (nlat+nlat+1)/(dpb*dpb)
        1839             :          ! yakimiw's formula permits using old pb and dpb
        1840             :          !--------------------------------------------------
        1841           0 :          wts(nix) = (nlat+nlat+1)/ (dpb+pb*dcos(zlast)/dsin(zlast))**2
        1842           0 :          nix      = nix - 1
        1843           0 :          if(nix==nhalf-1) zero = 3._r8*zero - pi
        1844           0 :          if(nix<nhalf-1) zero = zero + zero - zprev
        1845           0 :          zprev = zhold
        1846             :       end do
        1847             :
        1848             :       ! extend points and weights via symmetries
        1849             :       !-------------------------------------------
        1850           0 :       if(mnlat/=0) then
        1851           0 :         theta(nhalf) = halfPI
        1852           0 :         call sh_legp_dlegp_theta(nlat,halfPI,cz,theta(ns2+1),wts(ns2+1),pb,dpb)
        1853           0 :         wts(nhalf) = (nlat+nlat+1)/ (dpb*dpb)
        1854             :       endif
        1855             :
        1856           0 :       do ii = 1,ns2
        1857           0 :         wts  (nlat-ii+1) = wts(ii)
        1858           0 :         theta(nlat-ii+1) = pi - theta(ii)
        1859             :       end do
        1860             :
        1861             :       summ = 0._r8
        1862           0 :       do ii = 1,nlat
        1863           0 :         summ = summ + wts(ii)
        1864             :       end do
        1865           0 :       do ii = 1,nlat
        1866           0 :         wts(ii) = 2._r8*wts(ii)/summ
        1867             :       end do
        1868             :
        1869             :     end subroutine sh_create_gaus_grid
        1870             :     !=======================================================================
        1871             :
        1872             :     !=======================================================================
        1873           0 :     subroutine sh_fourier_coefs_dp(nn,cz,cp,dcp)
        1874             :       !
        1875             :       ! spherepack cpdp
        1876             :       !
        1877             :       !     computes the fourier coefficients of the legendre
        1878             :       !     polynomial p_n^0 and its derivative.
        1879             :       !     nn is the degree and nn/2 or (nn+1)/2
        1880             :       !     coefficients are returned in cp depending on whether
        1881             :       !     nn is even or odd. The same number of coefficients
        1882             :       !     are returned in dcp. For nn even the constant
        1883             :       !     coefficient is returned in cz.
        1884             :       !=====================================================================
        1885             :       !
        1886             :       ! Passed variables
        1887             :       !-----------------
        1888             :       integer, intent(in) :: nn
        1889             :       real(r8), intent(out) :: cz
        1890             :       real(r8), intent(out) :: cp(nn/2+1)
        1891             :       real(r8), intent(out) :: dcp(nn/2+1)
        1892             :       !
        1893             :       ! Local Values
        1894             :       !--------------
        1895             :       real(r8):: t1,t2,t3,t4
        1896             :       integer :: ncp,jj
        1897             :
        1898           0 :       ncp = (nn+1)/2
        1899           0 :       t1  = -1._r8
        1900           0 :       t2  = nn + 1._r8
        1901           0 :       t3  = 0._r8
        1902           0 :       t4  = nn + nn + 1._r8
        1903           0 :       if(mod(nn,2)==0) then
        1904           0 :         cp(ncp) = 1._r8
        1905           0 :         do jj = ncp,2,-1
        1906           0 :           t1 = t1 + 2._r8
        1907           0 :           t2 = t2 - 1._r8
        1908           0 :           t3 = t3 + 1._r8
        1909           0 :           t4 = t4 - 2._r8
        1910           0 :           cp(jj-1) = (t1*t2)/ (t3*t4)*cp(jj)
        1911             :         end do
        1912           0 :         t1 = t1 + 2._r8
        1913           0 :         t2 = t2 - 1._r8
        1914           0 :         t3 = t3 + 1._r8
        1915           0 :         t4 = t4 - 2._r8
        1916           0 :         cz = (t1*t2)/ (t3*t4)*cp(1)
        1917           0 :         do jj = 1,ncp
        1918           0 :           dcp(jj) = (jj+jj)*cp(jj)
        1919             :         end do
        1920             :       else
        1921           0 :         cp(ncp) = 1._r8
        1922           0 :         do jj = ncp - 1,1,-1
        1923           0 :           t1 = t1 + 2._r8
        1924           0 :           t2 = t2 - 1._r8
        1925           0 :           t3 = t3 + 1._r8
        1926           0 :           t4 = t4 - 2._r8
        1927           0 :           cp(jj) = (t1*t2)/ (t3*t4)*cp(jj+1)
        1928             :         end do
        1929           0 :         do jj = 1,ncp
        1930           0 :           dcp(jj) = (jj+jj-1)*cp(jj)
        1931             :         end do
        1932             :       endif
        1933             :
        1934           0 :     end subroutine sh_fourier_coefs_dp
        1935             :     !=======================================================================
        1936             :
        1937             :     !=======================================================================
        1938           0 :     subroutine sh_legp_dlegp_theta(nn,theta,cz,cp,dcp,pb,dpb)
        1939             :       !
        1940             :       ! spherepack tpdp
        1941             :       !
        1942             :       !     computes pn(theta) and its derivative dpb(theta) with
        1943             :       !     respect to theta
        1944             :       !=====================================================================
        1945             :       !
        1946             :       ! Passed variables
        1947             :       !------------------
        1948             :       integer, intent(in) :: nn
        1949             :       real(r8), intent(in) :: theta
        1950             :       real(r8), intent(in) :: cz
        1951             :       real(r8), intent(in) :: cp (nn/2+1)
        1952             :       real(r8), intent(in) :: dcp(nn/2+1)
        1953             :       real(r8), intent(out) :: pb
        1954             :       real(r8), intent(out) :: dpb
        1955             :       !
        1956             :       ! Local Values
        1957             :       !--------------
        1958             :       real(r8):: cdt,sdt,cth,sth,chh
        1959             :       integer :: kdo,kk
        1960             :
        1961           0 :       cdt = dcos(theta+theta)
        1962           0 :       sdt = dsin(theta+theta)
        1963           0 :       if(mod(nn,2)==0) then
        1964             :         ! n even
        1965             :         !----------
        1966           0 :         kdo = nn/2
        1967           0 :         pb  = .5_r8*cz
        1968           0 :         dpb = 0._r8
        1969           0 :         if(nn>0) then
        1970             :           cth = cdt
        1971             :           sth = sdt
        1972           0 :           do kk = 1,kdo
        1973           0 :             pb  = pb  +  cp(kk)*cth
        1974           0 :             dpb = dpb - dcp(kk)*sth
        1975           0 :             chh = cdt*cth - sdt*sth
        1976           0 :             sth = sdt*cth + cdt*sth
        1977           0 :             cth = chh
        1978             :           end do
        1979             :         endif
        1980             :       else
        1981             :         ! n odd
        1982             :         !-----------
        1983           0 :         kdo = (nn+1)/2
        1984           0 :         pb  = 0._r8
        1985           0 :         dpb = 0._r8
        1986           0 :         cth = dcos(theta)
        1987           0 :         sth = dsin(theta)
        1988           0 :         do kk = 1,kdo
        1989           0 :           pb  = pb  +  cp(kk)*cth
        1990           0 :           dpb = dpb - dcp(kk)*sth
        1991           0 :           chh = cdt*cth - sdt*sth
        1992           0 :           sth = sdt*cth + cdt*sth
        1993           0 :           cth = chh
        1994             :         end do
        1995             :       endif
        1996             :
        1997           0 :     end subroutine sh_legp_dlegp_theta
        1998             :     !=======================================================================
        1999             :
        2000           0 : end module zonal_mean_mod

Generated by: LCOV version 1.14