LCOV - code coverage report
Current view: top level - dynamics/fv - dynamics_vars.F90 (source / functions) Hit Total Coverage
Test: coverage.info Lines: 282 403 70.0 %
Date: 2025-03-14 01:21:06 Functions: 2 8 25.0 %

          Line data    Source code
       1             : module dynamics_vars
       2             : 
       3             : !-----------------------------------------------------------------------
       4             : ! CAM fvcore internal variables
       5             : !
       6             : ! !REVISION HISTORY:
       7             : !   01.06.06   Sawyer     Consolidated from various code snippets
       8             : !   03.06.25   Sawyer     Cleaned up, used ParPatternCopy (Create)
       9             : !   03.08.05   Sawyer     Removed rayf_init and hswf_init, related vars
      10             : !   03.10.22   Sawyer     pmgrid removed (now spmd_dyn)
      11             : !   03.11.18   Sawyer     Removed set_eta (ak, bk, now read from restart)
      12             : !   03.12.04   Sawyer     Moved T_FVDYCORE_GRID here (removed some vars)
      13             : !   04.08.25   Sawyer     Removed all module data members, now GRID only
      14             : !   04.10.06   Sawyer     Added spmd_dyn vars here; ESMF transpose vars
      15             : !   05.04.12   Sawyer     Added support for r4/r8 tracers
      16             : !   05.05.24   Sawyer     CAM/GEOS5 merge (removed GEOS_mod dependencies)
      17             : !   05.06.10   Sawyer     Scaled down version for CAM (no ESMF)
      18             : !   05.11.10   Sawyer     Removed dyn_interface (now in dyn_comp)
      19             : !   06.03.01   Sawyer     Removed m_ttrans, q_to_qxy, qxy_to_q, etc.
      20             : !   06.05.09   Sawyer     Added CONSV to dyn_state (conserve energy)
      21             : !   06.08.27   Sawyer     Removed unused ESMF code for RouteHandle
      22             : !-----------------------------------------------------------------------
      23             : 
      24             : use shr_kind_mod,       only: r8=>shr_kind_r8
      25             : use pmgrid,             only: plon, plat, plev
      26             : 
      27             : use decompmodule,       only: decomptype
      28             : use ghostmodule,        only: ghosttype
      29             : 
      30             : use cam_logfile,        only: iulog
      31             : use cam_abortutils,     only: endrun
      32             : 
      33             : #if defined(SPMD)
      34             : use parutilitiesmodule, only: parpatterntype, REAL4, INT4
      35             : #endif
      36             : 
      37             : implicit none
      38             : private
      39             : save
      40             : 
      41             : public :: &
      42             :    t_fvdycore_vars,      &
      43             :    t_fvdycore_grid,      &
      44             :    t_fvdycore_constants, &
      45             :    t_fvdycore_state,     &
      46             :    grid_vars_init,       &
      47             :    dynamics_clean
      48             : 
      49             : #ifdef SPMD
      50             : public :: spmd_vars_init
      51             : #endif
      52             : 
      53             : ! T_FVDYCORE_VARS contains the prognostic variables for FVdycore
      54             : 
      55             : type T_FVDYCORE_VARS
      56             :    real(r8), dimension(:,:,:  ), pointer     :: U      ! U winds (D-grid)
      57             :    real(r8), dimension(:,:,:  ), pointer     :: V      ! V winds (D-grid)
      58             :    real(r8), dimension(:,:,:  ), pointer     :: PT     ! scaled virtual pot. temp.
      59             :    real(r8), dimension(:,:,:  ), pointer     :: PE     ! Pressure at layer edges
      60             :    real(r8), dimension(:,:,:  ), pointer     :: PKZ    ! P^kappa mean
      61             :    real(r8), dimension(:,:,:,:), pointer     :: tracer ! Tracers
      62             : end type T_FVDYCORE_VARS
      63             : 
      64             : ! T_FVDYCORE_GRID contains information about the horizontal and vertical
      65             : ! discretization and decompositions.
      66             :  
      67             : type T_FVDYCORE_GRID
      68             : 
      69             :    ! PILGRIM communication information
      70             : 
      71             :    integer :: twod_decomp = 0  ! 1 for multi-2D decompositions, 0 otherwise
      72             :               ! To assure that the latitudinal decomposition operates
      73             :               ! as efficiently as before, a separate parameter "twod_decomp" has
      74             :               ! been defined; a value of 1 refers to the multi-2D decomposition with
      75             :               ! transposes; a value of 0 means that the decomposition is effectively
      76             :               ! one-dimensional, thereby enabling the transpose logic to be skipped;
      77             :               ! there is an option to force computation of transposes even for case
      78             :               ! where decomposition is effectively 1-D.
      79             : 
      80             :    integer :: npes_xy= 1    ! number of PEs for XY decomposition
      81             :    integer :: npes_yz= 1    ! number of PEs for YZ decomposition
      82             :    integer :: myid_y = 0    ! subdomain index (0-based) in latitude (y)
      83             :    integer :: myid_z = 0    ! subdomain index (0 based) in level (z)
      84             :    integer :: npr_y  = 1    ! number of subdomains in y
      85             :    integer :: npr_z  = 1    ! number of subdomains in z
      86             : 
      87             :    integer :: myidxy_x = 0  ! subdomain index (0-based) in longitude (x) (second. decomp.)
      88             :    integer :: myidxy_y = 0  ! subdomain index (0 based) in latitude (y) (second. decomp.)
      89             :    integer :: nprxy_x = 1   ! number of subdomains in x (second. decomp.)
      90             :    integer :: nprxy_y = 1   ! number of subdomains in y (second. decomp.)
      91             :    integer :: iam = 0       ! 
      92             : 
      93             :    integer :: mod_method = 0  ! 1 for mpi derived types with transposes, 0 for contiguous buffers
      94             :    integer :: mod_geopk = 0   ! 1 for mpi derived types with transposes, 0 for contiguous buffers
      95             :    integer :: mod_gatscat = 0 ! 1 for mpi derived types with transposes, 0 for contiguous buffers
      96             : 
      97             :    type(decomptype) :: strip2d, strip2dx, strip3dxyz, strip3dxzy,          &
      98             :                        strip3dxyzp, strip3zaty, strip3dxzyp,               &
      99             :                        strip3yatz, strip3yatzp, strip3zatypt,              &
     100             :                        strip3kxyz, strip3kxzy, strip3kxyzp, strip3kxzyp,   &
     101             :                        strip3dyz, checker3kxy
     102             : 
     103             :    integer :: commdyn           ! communicator for all dynamics
     104             :    integer :: commxy            ! communicator for XY decomposition
     105             :    integer :: commyz            ! communicator for YZ decomposition
     106             :    integer :: commnyz           ! communicator for multiple YZ decomposition
     107             : 
     108             :    integer :: comm_y            ! communicator in latitude
     109             :    integer :: comm_z            ! communicator in vertical
     110             :    integer :: commxy_x          ! communicator in longitude (xy second. decomp.)
     111             :    integer :: commxy_y          ! communicator in latitude (xy second. decomp.)
     112             :    logical :: geopkdist         ! use distributed method for geopotential calculation 
     113             :                                 !  with 2D decomp.
     114             :    logical :: geopk16byte       ! use Z-parallel distributed method for geopotential 
     115             :                                 !  calculation with 2D decomp.; otherwise use Z-serial
     116             :                                 !  pipeline algorithm when using distributed algoritm
     117             :    integer :: geopkblocks       ! number of stages to use in Z-serial pipeline 
     118             :                                 !  (non-transpose) geopotential algorithm
     119             :    integer :: modc_dynrun(4)    ! 1: mod_comm irregular underlying communication method for dyn_run/misc
     120             :                                 ! 2: mod_comm irregular communication handshaking for dyn_run/misc
     121             :                                 ! 3: mod_comm irregular communication send protocol for dyn_run/misc
     122             :                                 ! 4: mod_comm irregular communication nonblocking request throttle for dyn_run/misc
     123             :    integer :: modc_cdcore(4)    ! 1: mod_comm irregular underlying communication method for cd_core/geopk
     124             :                                 ! 2: mod_comm irregular communication handshaking for cd_core/geopk
     125             :                                 ! 3: geopk_d and mod_comm irregular communication send protocol for cd_core/geopk
     126             :                                 ! 4: mod_comm irregular communication nonblocking request throttle for cd_core/geopk
     127             :    integer :: modc_gather(4)    ! 1: mod_comm irregular underlying communication method for gather
     128             :                                 ! 2: mod_comm irregular communication handshaking for gather
     129             :                                 ! 3: mod_comm irregular communication send protocol for gather
     130             :                                 ! 4: mod_comm irregular communication nonblocking request throttle for gather
     131             :    integer :: modc_scatter(4)   ! 1: mod_comm irregular underlying communication method for scatter
     132             :                                 ! 2: mod_comm irregular communication handshaking for scatter
     133             :                                 ! 3: mod_comm irregular communication send protocol for scatter
     134             :                                 ! 4: mod_comm irregular communication nonblocking request throttle for scatter 
     135             :    integer :: modc_tracer(4)    ! 1: mod_comm irregular underlying communication method for multiple tracers
     136             :                                 ! 2: mod_comm irregular communication handshaking for multiple tracers
     137             :                                 ! 3: mod_comm irregular communication send protocol for multiple tracers
     138             :                                 ! 4: mod_comm irregular communication nonblocking request throttle for multiple tracers 
     139             :    integer :: modc_onetwo       ! one or two simultaneous mod_comm irregular communications (excl. tracers)
     140             :    integer :: modc_tracers      ! max number of tracers for simultaneous mod_comm irregular communications
     141             : 
     142             : #if defined(SPMD)
     143             :    type (ghosttype)        :: ghostu_yz, ghostv_yz, ghostpt_yz,                   &
     144             :                               ghostpe_yz, ghostpkc_yz
     145             :    type (parpatterntype)   :: u_to_uxy, uxy_to_u, v_to_vxy, vxy_to_v,             &
     146             :                               ikj_yz_to_xy, ikj_xy_to_yz,                         &
     147             :                               ijk_yz_to_xy, ijk_xy_to_yz,                         &
     148             :                               pe_to_pexy, pexy_to_pe,                             &
     149             :                               pt_to_ptxy, ptxy_to_pt, pkxy_to_pkc,                &
     150             :                               r4_xy_to_yz, r4_yz_to_xy, q3_to_qxy3, qxy3_to_q3,   &
     151             :                               xy2d_to_yz2d, yz2d_to_xy2d, scatter_3d, gather_3d,  &
     152             :                               g_2dxy_r8, g_2dxy_r4, g_2dxy_i4,                    &
     153             :                               s_2dxy_r8, s_2dxy_r4, s_2dxy_i4,                    &
     154             :                               g_3dxyz_r8, g_3dxyz_r4, g_3dxyzp_r8, g_3dxyzp_r4,   &
     155             :                               s_3dxyz_r8, s_3dxyz_r4, s_3dxyzp_r8, s_3dxyzp_r4
     156             : #endif
     157             : 
     158             :    ! END PILGRIM communication information
     159             : 
     160             : 
     161             :    integer  :: JFIRST = 1       ! Start latitude (exclusive)
     162             :    integer  :: JLAST  = plat    ! End latitude (exclusive)
     163             :  
     164             :    integer  :: ng_c = 0         ! Ccore ghosting
     165             :    integer  :: ng_d = 0         ! Dcore ghosting
     166             :    integer  :: ng_s = 0         ! Staggered grid ghosting for
     167             :                                 ! certain arrays, max(ng_c+1,ng_d)
     168             :    ! For 2D decomposition
     169             : 
     170             :    integer  :: IFIRSTXY = 1     ! Start longitude (exclusive)
     171             :    integer  :: ILASTXY  = plon  ! End longitude (exclusive)
     172             :    integer  :: JFIRSTXY = 1     ! Start latitude (exclusive)
     173             :    integer  :: JLASTXY  = plat  ! End latitude (exclusive)
     174             : 
     175             :    integer  :: IM               ! Full longitude dim
     176             :    integer  :: JM               ! Full latitude dim (including poles)
     177             : 
     178             :    real(r8) :: DL
     179             :    real(r8) :: DP
     180             :    real(r8) :: ACAP
     181             :    real(r8) :: RCAP
     182             : 
     183             :    real(r8), dimension(:), pointer :: COSP             ! Cosine of lat angle -- volume mean
     184             :    real(r8), dimension(:), pointer :: SINP             ! Sine of lat angle -- volume mean
     185             :    real(r8), dimension(:), pointer :: COSE             ! Cosine at finite volume edge
     186             :    real(r8), dimension(:), pointer :: SINE             ! Sine at finite volume edge
     187             :    real(r8), dimension(:), pointer :: ACOSP            ! Reciprocal of cosine of lat angle
     188             : 
     189             :    real(r8), dimension(:), pointer :: ACOSU            ! Reciprocal of cosine of lat angle (staggered)
     190             : 
     191             :    real(r8), dimension(:), pointer :: COSLON           ! Cosine of longitudes - volume center
     192             :    real(r8), dimension(:), pointer :: SINLON           ! Sine of longitudes - volume center
     193             :    real(r8), dimension(:), pointer :: COSL5            ! Cosine of longitudes - volume center
     194             :    real(r8), dimension(:), pointer :: SINL5            ! Sine of longitudes - volume center
     195             :  
     196             :    ! Variables which are used repeatedly in CD_CORE
     197             : 
     198             :    integer ::       js2g0
     199             :    integer ::       jn2g0
     200             :    integer ::       jn1g1
     201             : 
     202             :    real(r8), pointer :: trigs(:)
     203             :    real(r8), pointer :: fc(:), f0(:)
     204             :    real(r8), pointer :: dc(:,:), de(:,:), sc(:), se(:)
     205             :    real(r8), pointer :: cdx(:,:), cdy(:,:)
     206             :    real(r8), pointer :: cdx4(:,:), cdy4(:,:)                           ! for div4 damping
     207             :    real(r8), pointer :: cdxde(:,:), cdxdp(:,:),cdyde(:,:),cdydp(:,:)   ! for del2 damping
     208             :    real(r8), pointer :: cdxdiv(:,:), cdydiv(:,:), cdtau4(:,:)          ! for del2 damping
     209             : 
     210             :    real(r8), pointer :: dcdiv4(:,:), dediv4(:,:), scdiv4(:), sediv4(:) ! for div4 damping
     211             : 
     212             :    real(r8), pointer :: dtdx(:), dtdxe(:), txe5(:), dtxe5(:)
     213             :    real(r8), pointer :: dyce(:),   dx(:) ,  rdx(:),    cy(:)
     214             :    real(r8), pointer :: dtdx2(:), dtdx4(:),  dxdt(:), dxe(:)
     215             :    real(r8), pointer :: cye(:),    dycp(:),  rdxe(:)
     216             : 
     217             :    real(r8) :: rdy, dtdy, dydt, dtdy5, tdy5
     218             :    real(r8) :: dt0 = 0
     219             : 
     220             :    integer  :: ifax(13)
     221             : 
     222             :    real(r8) ::  zt_c
     223             :    real(r8) ::  zt_d
     224             : 
     225             :    ! This part refers to the vertical grid
     226             : 
     227             :    integer                         :: KM              ! Numer of levels
     228             :    integer                         :: KMAX            ! KM+1 (?)
     229             : 
     230             :    ! For 2D decomposition
     231             : 
     232             :    integer                         :: KFIRST = 1      ! Start level (exclusive)
     233             :    integer                         :: KLAST  = plev   ! End level (exclusive)
     234             :    integer                         :: KLASTP = plev+1 ! klast+1, except km+1 when klastp=km+1
     235             : 
     236             :    integer                         :: KORD            ! monotonicity order for mapping (te_map)
     237             :    integer                         :: KS              ! Number of true pressure levels (out of KM+1)
     238             :    real(r8)                        :: PTOP            ! pressure at top (ak(1))
     239             :    real(r8)                        :: PINT            ! initial pressure (ak(km+1))
     240             :    real(r8), dimension(:), pointer :: AK              ! Sigma mapping
     241             :    real(r8), dimension(:), pointer :: BK              ! Sigma mapping
     242             : 
     243             :    ! Tracers
     244             : 
     245             :    integer                         :: NQ              ! Number of advected tracers
     246             :    integer                         :: NTOTQ           ! Total number of tracers (NQ <= NC)
     247             : 
     248             :    integer  :: ct_overlap    ! nonzero for overlap of cd_core and trac2d, 0 otherwise
     249             :    integer  :: trac_decomp   ! size of tracer domain decomposition for trac2d
     250             : 
     251             :    ! Extra subdomain bounds for cd_core/trac2d overlap and trac2d decomposition
     252             :    ! Relevant for secondary yz decomposition only; refers back to primary yz decomposition
     253             :    integer                         :: JFIRSTCT          ! jfirst
     254             :    integer                         :: JLASTCT           ! jlast
     255             :    integer                         :: KFIRSTCT          ! kfirst
     256             :    integer                         :: KLASTCT           ! klast
     257             : 
     258             :    ! Bounds for tracer decomposition
     259             :    integer, dimension(:), pointer  :: ktloa             ! lower tracer index (global map)
     260             :    integer, dimension(:), pointer  :: kthia             ! upper tracer index (global map)
     261             :    integer                         :: ktlo              ! lower tracer index (local)
     262             :    integer                         :: kthi              ! upper tracer index (local)
     263             : 
     264             :    logical :: high_alt  ! high-altitude physics parameters switch
     265             : 
     266             : end type T_FVDYCORE_GRID
     267             : 
     268             : ! Constants used by fvcore
     269             : type T_FVDYCORE_CONSTANTS
     270             :    real(r8)                             :: pi
     271             :    real(r8)                             :: omega    ! angular velocity of earth's rotation  
     272             :    real(r8)                             :: cp       ! heat capacity of air at constant pressure
     273             :    real(r8)                             :: ae       ! radius of the earth (m)
     274             :    real(r8)                             :: rair     ! Gas constant of the air
     275             :    real(r8)                             :: cappa    ! Cappa?
     276             :    real(r8)                             :: zvir     ! RWV/RAIR-1
     277             : end type T_FVDYCORE_CONSTANTS
     278             : 
     279             : type t_fvdycore_state
     280             :    type (t_fvdycore_vars)      :: vars
     281             :    type (t_fvdycore_grid )     :: grid
     282             :    type (t_fvdycore_constants) :: constants
     283             :    real(r8) :: DT            ! Large time step
     284             :    real(r8) :: CHECK_DT      ! Time step to check maxmin
     285             :    integer  :: ICD, JCD      ! Algorithm orders (C Grid)
     286             :    integer  :: IORD, JORD    ! Algorithm orders (D Grid)
     287             :    integer  :: KORD          ! Vertical order
     288             :    integer  :: TE_METHOD     ! method for total energy mapping (te_map)
     289             :    logical  :: CONSV         ! dycore conserves tot. en.
     290             :    integer  :: NSPLIT
     291             :    integer  :: NSPLTRAC
     292             :    integer  :: NSPLTVRM
     293             :    integer  :: FILTCW        ! filter c-grid winds if positive
     294             :    integer  :: fft_flt       ! 0 => FFT/algebraic filter; 1 => FFT filter
     295             :    integer  :: div24del2flag ! 2 for 2nd order div damping, 4 for 4th order div damping,
     296             :                              ! 42 for 4th order div damping plus 2nd order velocity damping
     297             :    real(r8) :: del2coef      ! strength of 2nd order velocity damping
     298             :    logical  :: high_order_top! use normal 4-order PPM calculation near the model top
     299             :    logical  :: am_geom_crrct ! apply correction for angular momentum (AM) conservation in geometry
     300             :    logical  :: am_correction ! apply correction for angular momentum (AM) conservation in SW eqns
     301             :    logical  :: am_fixer      ! apply global fixer to conserve AM
     302             :    logical  :: am_fix_lbl    ! apply global AM fixer level by level
     303             :    logical  :: am_diag       ! turns on an AM diagnostic calculations
     304             : end type t_fvdycore_state
     305             : 
     306             : !========================================================================================
     307             : contains
     308             : !========================================================================================
     309             : 
     310             : #if defined(SPMD)
     311             : 
     312        1536 : subroutine spmd_vars_init(imxy, jmxy, jmyz, kmyz, grid)
     313             : 
     314             :    ! Initialize SPMD related variables.
     315             :    ! !REVISION HISTORY: 
     316             :    !   02.11.08    Sawyer       Creation
     317             :    !   03.05.07    Sawyer       Use ParPatternCopy for q_to_qxy, etc.
     318             :    !   03.07.23    Sawyer       Removed dependency on constituents module
     319             :    !   03.09.10    Sawyer       Reactivated u_to_uxy, etc, redefined pe2pexy
     320             :    !   03.11.19    Sawyer       Merged in CAM code with mod_method
     321             :    !   04.08.25    Sawyer       Added GRID as argument
     322             :    
     323             :    use decompmodule, only: decompcreate, decompfree
     324             :    use ghostmodule, only : ghostcreate, ghostfree
     325             :    use parutilitiesmodule, only : gid, parpatterncreate, parsplit
     326             :    use mpishorthand, only: mpiint
     327             : 
     328             :    ! Arguments
     329             : 
     330             :    integer, dimension(:), intent(in) :: imxy
     331             :    integer, dimension(:), intent(in) :: jmxy
     332             :    integer, dimension(:), intent(in) :: jmyz
     333             :    integer, dimension(:), intent(in) :: kmyz
     334             : 
     335             :    type( t_fvdycore_grid ), intent(inout) :: grid
     336             : 
     337             :    ! local variables:
     338             : 
     339             :    type(decomptype) :: global2d, local2d
     340             : 
     341             :    integer :: im, jm, km         !  Global dims
     342             :    integer :: nq
     343             : 
     344             :    integer :: nprxy_x    ! XY decomp - Nr in X
     345             :    integer :: nprxy_y    ! XY decomp - Nr in Y
     346             :    integer :: npryz_y    ! YZ decomp - Nr in Y
     347             :    integer :: npryz_z    ! YZ decomp - Nr in Z
     348             :    integer :: npes_xy    ! XY decomp - Total nr
     349             :    integer :: npes_yz    ! YZ decomp - Total nr
     350             : 
     351             :    integer :: commxy     ! Communicator for XY decomp
     352             :    integer :: commyz     ! Communicator for YZ decomp
     353             :    integer :: commnyz    ! Communicator for multiple YZ decomp
     354             : 
     355             :    integer :: jfirstxy, jlastxy
     356             :    integer :: jfirst, jlast
     357             :    integer :: kfirst, klast
     358             : 
     359             :    integer :: ng_s, ng_d   !  Ghost widths
     360             : 
     361             :    integer :: rank_y, rank_z, rankxy_x, rankxy_y  ! Currently not used
     362             :    integer :: size_y, size_z, sizexy_x, sizexy_y  ! Currently not used
     363             : 
     364             :    integer :: xdist(1), ydistk(1), zdist1(1), zdistxy(1) ! non-distributed dims
     365        1536 :    integer, allocatable :: xdist_global(:), ydist_global(:) 
     366        1536 :    integer, allocatable :: zdist(:) ! number of levels per subdomain
     367             : 
     368             :    integer :: ig1, ig2, jg1, jg2
     369             :    integer :: jg1d, jg2d, jg1s, jg2s
     370             :    integer :: kg1, kg2, kg2p
     371             : 
     372             :    integer :: ct_overlap
     373             :    integer :: trac_decomp
     374             :    integer :: ktmod, ml
     375             :    integer :: myidmod
     376             :    integer :: ictstuff(4)
     377             :    integer :: kquot, krem, krun, kt, mlt
     378             :    !---------------------------------------------------------------------------
     379             : 
     380        1536 :    im = grid%im
     381        1536 :    jm = grid%jm
     382        1536 :    km = grid%km
     383        1536 :    nq = grid%nq
     384             : 
     385        1536 :    nprxy_x = grid%nprxy_x
     386        1536 :    nprxy_y = grid%nprxy_y
     387        1536 :    npryz_y = grid%npr_y
     388        1536 :    npryz_z = grid%npr_z
     389        1536 :    npes_xy = grid%npes_xy
     390        1536 :    npes_yz = grid%npes_yz
     391             : 
     392        1536 :    commxy  = grid%commxy
     393        1536 :    commyz  = grid%commyz
     394        1536 :    commnyz = grid%commnyz
     395             : 
     396        1536 :    jfirstxy = grid%jfirstxy
     397        1536 :    jlastxy  = grid%jlastxy
     398             : 
     399        1536 :    jfirst = grid%jfirst
     400        1536 :    jlast  = grid%jlast
     401        1536 :    kfirst = grid%kfirst
     402        1536 :    klast  = grid%klast
     403             : 
     404        1536 :    ng_s   = grid%ng_s
     405        1536 :    ng_d   = grid%ng_d
     406             : 
     407             :    ! Split communicators
     408        1536 :    call parsplit(commyz, grid%myid_z, gid, grid%comm_y, rank_y, size_y)
     409        1536 :    call parsplit(commyz, grid%myid_y, gid, grid%comm_z, rank_z, size_z)
     410        1536 :    call parsplit(commxy, grid%myidxy_y, gid, grid%commxy_x, rankxy_x, sizexy_x)
     411        1536 :    call parsplit(commxy, grid%myidxy_x, gid, grid%commxy_y, rankxy_y, sizexy_y)
     412             : 
     413             : 
     414             :    ! create decompositions for CAM data structures
     415             : 
     416        4608 :    allocate(xdist_global(nprxy_x))
     417        4608 :    allocate(ydist_global(nprxy_y))
     418        4608 :    allocate(zdist (npryz_z))
     419        1536 :    xdist(1) = im
     420             : 
     421             :    ! Create PILGRIM decompositions (see decompmodule)
     422             : 
     423        1536 :    if (gid < npes_xy) then
     424       19968 :       xdist_global = 0
     425       99840 :       ydist_global = 0
     426        1536 :       xdist_global(1) = im
     427        1536 :       ydist_global(1) = jm
     428             :       call decompcreate(nprxy_x, nprxy_y, xdist_global,            &
     429        1536 :                         ydist_global, global2d)
     430        1536 :       call decompcreate(nprxy_x, nprxy_y, imxy, jmxy, local2d )
     431             : 
     432             :       ! Decompositions needed on xy decomposition for parpatterncreate
     433             : 
     434        1536 :       call decompcreate( 1, npryz_y, xdist, jmyz, grid%strip2d )
     435             :       call decompcreate( 1, npryz_y, npryz_z, xdist,                &
     436        1536 :                          jmyz, kmyz, grid%strip3dxyz )
     437             :       call decompcreate( "xzy", 1, npryz_z, grid%npr_y, xdist,      &
     438        1536 :                            kmyz, jmyz, grid%strip3dxzy )
     439             : 
     440             :       ! For y communication within z subdomain (klast version)
     441             :       ! Use myidmod to have valid index for inactive processes
     442             :       !  for smaller yz decomposition
     443        1536 :       myidmod = mod(grid%myid_z, grid%npr_z)  ! = myid_z for active yz process
     444        1536 :       zdist1(1) = kmyz(myidmod+1)
     445             :       call decompcreate( 1, npryz_y, 1, xdist, jmyz, zdist1,        &
     446        1536 :                          grid%strip3yatz )
     447             : 
     448             :       ! For z communication within y subdomain
     449             : 
     450        1536 :       ydistk(1) = jmyz(grid%myid_y+1)
     451             :       call decompcreate( 1, 1, npryz_z, xdist, ydistk, kmyz,        &
     452        1536 :                          grid%strip3zaty )
     453             : 
     454             :       ! Arrays dimensioned plev+1
     455             : 
     456       19968 :       zdist(:) = kmyz(:)
     457        1536 :       zdist(npryz_z) = kmyz(npryz_z) + 1
     458             :       call decompcreate( 1, npryz_y, npryz_z, xdist, jmyz, zdist,&
     459        1536 :                          grid%strip3dxyzp )
     460             :       call decompcreate( "xzy", 1, npryz_z, npryz_y,                &
     461        1536 :                          xdist, zdist, jmyz, grid%strip3dxzyp )
     462             : 
     463             :       ! Arrays dimensioned plev+1, within y subdomain
     464             : 
     465        1536 :       ydistk(1) = jmyz(grid%myid_y+1)
     466             :       call decompcreate( "xzy", 1, npryz_z, 1, xdist, zdist, ydistk,   &
     467        1536 :                          grid%strip3zatypt )
     468             : 
     469             :       ! For y communication within z subdomain (klast+1 version)
     470             :       ! Use myidmod to have valid index for inactive processes
     471             :       !  for smaller yz decomposition
     472        1536 :       myidmod = mod(grid%myid_z, grid%npr_z)  ! = myid_z for active yz process
     473        1536 :       zdist1(1) = kmyz(myidmod+1)
     474             :       call decompcreate( 1, npryz_y, 1, xdist, jmyz, zdist1,        &
     475        1536 :                          grid%strip3yatzp )
     476             : 
     477             :       ! For the 2D XY-YZ data transfer, we need a short 3D array
     478       19968 :       zdist(:) = 1    ! One copy on each z PE set
     479             :       call decompcreate( 1, npryz_y, npryz_z,                       &
     480        1536 :                          xdist, jmyz, zdist, grid%strip3dyz )
     481             :    end if
     482             : 
     483             :    ! Secondary xy decomposition
     484             : 
     485        1536 :    if (grid%twod_decomp == 1) then
     486        1536 :       if (gid < npes_xy) then
     487        1536 :          zdistxy(1) = npryz_z     ! All npr_z copies on 1 PE
     488             :          call decompcreate( nprxy_x, nprxy_y, 1,                     &
     489        1536 :                             imxy, jmxy, zdistxy, grid%checker3kxy )
     490        1536 :          zdistxy(1) = km
     491             :          call decompcreate( nprxy_x, nprxy_y, 1,                     &
     492        1536 :                             imxy, jmxy, zdistxy, grid%strip3kxyz )
     493             :          call decompcreate( "xzy", nprxy_x, 1, nprxy_y,              &
     494        1536 :                             imxy, zdistxy, jmxy, grid%strip3kxzy )
     495             : 
     496        1536 :          zdistxy(1) = zdistxy(1) + 1
     497             :          call decompcreate( nprxy_x, nprxy_y, 1,                     &
     498        1536 :                             imxy, jmxy, zdistxy, grid%strip3kxyzp )
     499             :          call decompcreate( "xzy", nprxy_x, 1, nprxy_y,              &
     500        1536 :                             imxy, zdistxy, jmxy, grid%strip3kxzyp )
     501        1536 :          zdistxy(1) = jlastxy - jfirstxy + 1
     502        1536 :          call decompcreate( nprxy_x, 1, imxy, zdistxy, grid%strip2dx )
     503             :       end if
     504             :    end if
     505             : 
     506        1536 :    deallocate(zdist)
     507        1536 :    deallocate(ydist_global)
     508        1536 :    deallocate(xdist_global)
     509             : 
     510        1536 :    if ( grid%twod_decomp == 1 ) then
     511             : 
     512             :       ! Initialize ghost regions
     513             : 
     514             :       ! Set limits for ghostcreate
     515        1536 :       ig1 = 1
     516        1536 :       ig2 = im
     517        1536 :       jg1 = jfirst
     518        1536 :       jg2 = jlast
     519        1536 :       jg1d = jfirst-ng_d
     520        1536 :       jg1s = jfirst-ng_s
     521        1536 :       jg2d = jlast+ng_d
     522        1536 :       jg2s = jlast+ng_s
     523        1536 :       kg1 = kfirst
     524        1536 :       kg2 = klast
     525        1536 :       kg2p = klast+1
     526             : 
     527             :       ! Call ghostcreate with null ranges for non-yz processes
     528        1536 :       if (gid >= npes_yz) then
     529           0 :          ig1 = im/2
     530           0 :          ig2 = ig1 - 1
     531           0 :          jg1 = (jfirst+jlast)/2
     532           0 :          jg2 = jg1 - 1
     533           0 :          jg1d = jg1
     534           0 :          jg1s = jg1
     535           0 :          jg2d = jg2
     536           0 :          jg2s = jg2
     537           0 :          kg1 = (kfirst+klast)/2
     538           0 :          kg2 = kg1 - 1
     539           0 :          kg2p = kg2
     540             :       end if
     541             : 
     542             :       ! Ghosted decompositions needed on xy decomposition for parpatterncreate
     543        1536 :       if (gid < npes_xy) then
     544             :          call ghostcreate( grid%strip3dxyz, gid, im, ig1, ig2, .true., &
     545             :                            jm, jg1d, jg2s, .false., &
     546        1536 :                            km, kg1, kg2, .false., grid%ghostu_yz )
     547             :          call ghostcreate( grid%strip3dxyz, gid, im, ig1, ig2, .true., &
     548             :                            jm, jg1s, jg2d, .false., &
     549        1536 :                            km, kg1, kg2, .false., grid%ghostv_yz )
     550             :          call ghostcreate( grid%strip3dxyz, gid, im, ig1, ig2, .true., &
     551             :                            jm, jg1d, jg2d, .false., &
     552        1536 :                            km, kg1, kg2, .false., grid%ghostpt_yz )
     553             :          call ghostcreate( grid%strip3dxzyp, gid, im, ig1, ig2, .true., &
     554             :                            km+1, kg1, kg2p, .false., &
     555        1536 :                            jm, jg1, jg2, .false., grid%ghostpe_yz)
     556             :          call ghostcreate( grid%strip3dxyzp, gid, im, ig1, ig2, .true., &
     557             :                            jm, jg1, jg2, .false.,       &
     558        1536 :                            km+1, kg1, kg2p, .false., grid%ghostpkc_yz)
     559             :       end if
     560             : 
     561             :       ! Initialize transposes
     562             : 
     563        1536 :       if (gid < npes_xy) then
     564             :          call parpatterncreate(commxy, grid%ghostu_yz, grid%strip3kxyz, &
     565        1536 :                                grid%u_to_uxy, mod_method=grid%mod_method)
     566             :          call parpatterncreate(commxy, grid%strip3kxyz,grid%ghostu_yz, &
     567        1536 :                                grid%uxy_to_u, mod_method=grid%mod_method)
     568             :          call parpatterncreate(commxy, grid%ghostv_yz, grid%strip3kxyz, &
     569        1536 :                                grid%v_to_vxy, mod_method=grid%mod_method)
     570             :          call parpatterncreate(commxy, grid%strip3kxyz, grid%ghostv_yz, &
     571        1536 :                                grid%vxy_to_v, mod_method=grid%mod_method)
     572             :          call parpatterncreate(commxy, grid%strip3dxyz, grid%strip3kxyz,&
     573        1536 :                                grid%ijk_yz_to_xy, mod_method=grid%mod_method)
     574             :          call parpatterncreate(commxy, grid%strip3kxyz, grid%strip3dxyz,&
     575        1536 :                                grid%ijk_xy_to_yz, mod_method=grid%mod_method)
     576             :          call parpatterncreate(commxy, grid%strip3dxzy, grid%strip3kxzy,&
     577        1536 :                                grid%ikj_yz_to_xy, mod_method=grid%mod_method)
     578             :          call parpatterncreate(commxy, grid%strip3kxzy, grid%strip3dxzy,&
     579        1536 :                                grid%ikj_xy_to_yz, mod_method=grid%mod_method)
     580             : 
     581             :          ! Note PE <-> PEXY has been redefined for PEXY ijk, but PE ikj
     582             : 
     583             :          call parpatterncreate(commxy, grid%ghostpe_yz, grid%strip3kxzyp, &
     584        1536 :                                grid%pe_to_pexy, mod_method=grid%mod_method)
     585             :          call parpatterncreate(commxy, grid%strip3kxzyp, grid%ghostpe_yz, &
     586        1536 :                                grid%pexy_to_pe, mod_method=grid%mod_method)
     587             :          call parpatterncreate(commxy, grid%ghostpt_yz, grid%strip3kxyz,  &
     588        1536 :                                grid%pt_to_ptxy, mod_method=grid%mod_method)
     589             :          call parpatterncreate(commxy, grid%strip3kxyz, grid%ghostpt_yz,  &
     590        1536 :                                grid%ptxy_to_pt, mod_method=grid%mod_method)
     591             :          call parpatterncreate(commxy, grid%strip3dxyz, grid%strip3kxyz,  &
     592             :                                grid%r4_yz_to_xy, mod_method=grid%mod_method,  &
     593        1536 :                                T = REAL4 )
     594             :          call parpatterncreate(commxy, grid%strip3kxyz, grid%strip3dxyz,  &
     595             :                                grid%r4_xy_to_yz, mod_method=grid%mod_method,  &
     596        1536 :                                T = REAL4 )
     597             :          call parpatterncreate(commxy, grid%strip3kxyzp, grid%ghostpkc_yz, &
     598        1536 :                                grid%pkxy_to_pkc, mod_method=grid%mod_method)
     599             : 
     600             :          ! These are for 'transposing' 2D arrays from XY YZ
     601             :          call parpatterncreate(commxy, grid%checker3kxy, grid%strip3dyz, &
     602        1536 :                                grid%xy2d_to_yz2d, mod_method=grid%mod_method)
     603             :          call parpatterncreate(commxy, grid%strip3dyz, grid%checker3kxy, &
     604        1536 :                                grid%yz2d_to_xy2d, mod_method=grid%mod_method)
     605             :       end if
     606             : 
     607             :       ! Free unneeded decompositions
     608             : 
     609        1536 :       call decompfree(grid%strip3dxzyp)
     610        1536 :       call decompfree(grid%strip3dyz)
     611        1536 :       call decompfree(grid%strip3yatz)
     612        1536 :       call decompfree(grid%strip3yatzp)
     613        1536 :       call decompfree(grid%strip3zaty)
     614        1536 :       call decompfree(grid%strip3zatypt)
     615        1536 :       call decompfree(grid%strip3kxyz)
     616        1536 :       call decompfree(grid%strip3kxzy)
     617        1536 :       call decompfree(grid%strip3kxyzp)
     618        1536 :       call decompfree(grid%strip3kxzyp)
     619        1536 :       call decompfree(grid%checker3kxy)
     620             : 
     621        1536 :       call ghostfree(grid%ghostu_yz)
     622        1536 :       call ghostfree(grid%ghostv_yz)
     623        1536 :       call ghostfree(grid%ghostpt_yz)
     624        1536 :       call ghostfree(grid%ghostpe_yz)
     625        1536 :       call ghostfree(grid%ghostpkc_yz)
     626             : 
     627             :    end if
     628             : 
     629             :    ! Define scatter and gather patterns for 2D and 3D unghosted arrays
     630             : 
     631        1536 :    if (gid < npes_xy) then
     632             :       call parpatterncreate( commxy, global2d, local2d, grid%s_2dxy_r8, &
     633        1536 :                              mod_method=grid%mod_gatscat )
     634             :       call parpatterncreate( commxy, local2d, global2d, grid%g_2dxy_r8,  &
     635        1536 :                              mod_method=grid%mod_gatscat )
     636             : 
     637             :       call parpatterncreate( commxy, global2d, local2d, grid%s_2dxy_r4, &
     638        1536 :                              mod_method=grid%mod_gatscat, T = REAL4 )
     639             :       call parpatterncreate( commxy, local2d, global2d, grid%g_2dxy_r4,  &
     640        1536 :                              mod_method=grid%mod_gatscat, T = REAL4 )
     641             : 
     642             :       call parpatterncreate( commxy, global2d, local2d, grid%s_2dxy_i4, &
     643        1536 :                              mod_method=grid%mod_gatscat, T = INT4 )
     644             :       call parpatterncreate( commxy, local2d, global2d, grid%g_2dxy_i4,  &
     645        1536 :                              mod_method=grid%mod_gatscat, T = INT4 )
     646             : 
     647             :       ! 3D XYZ patterns, will replace XZY patterns eventually
     648             : 
     649        1536 :       call parpatterncreate( commxy, grid%s_2dxy_r8, grid%s_3dxyz_r8, km )
     650        1536 :       call parpatterncreate( commxy, grid%g_2dxy_r8, grid%g_3dxyz_r8, km )
     651        1536 :       call parpatterncreate( commxy, grid%s_2dxy_r8, grid%s_3dxyzp_r8, km+1 )
     652        1536 :       call parpatterncreate( commxy, grid%g_2dxy_r8, grid%g_3dxyzp_r8, km+1 )
     653             : 
     654        1536 :       call parpatterncreate( commxy, grid%s_2dxy_r4, grid%s_3dxyz_r4, km )
     655        1536 :       call parpatterncreate( commxy, grid%g_2dxy_r4, grid%g_3dxyz_r4, km )
     656        1536 :       call parpatterncreate( commxy, grid%s_2dxy_r4, grid%s_3dxyzp_r4, km+1 )
     657        1536 :       call parpatterncreate( commxy, grid%g_2dxy_r4, grid%g_3dxyzp_r4, km+1 )
     658             : 
     659        1536 :       call decompfree( global2d )
     660        1536 :       call decompfree( local2d )
     661             :    end if
     662             : 
     663             :    ! Secondary subdomain limits for cd_core/trac2d overlap and trac2d decomposition
     664             : 
     665        1536 :    ct_overlap    = grid%ct_overlap
     666        1536 :    trac_decomp   = grid%trac_decomp
     667        1536 :    grid%jfirstct = grid%jfirst
     668        1536 :    grid%jlastct  = grid%jlast
     669        1536 :    grid%kfirstct = grid%kfirst
     670        1536 :    grid%klastct  = grid%klast
     671        1536 :    if (ct_overlap > 0) then
     672             :       mlt = 2
     673        1536 :    elseif (trac_decomp .gt. 1) then
     674             :       mlt = trac_decomp
     675             :    else
     676             :       mlt = 1
     677             :    end if
     678             : 
     679             :    if (mlt > 1) then
     680           0 :       if (gid < npes_yz) then
     681           0 :          ictstuff(1) = grid%jfirstct
     682           0 :          ictstuff(2) = grid%jlastct
     683           0 :          ictstuff(3) = grid%kfirstct
     684           0 :          ictstuff(4) = grid%klastct
     685           0 :          do ml = 2, mlt
     686           0 :             call mpisend(ictstuff, 4, mpiint, gid+(ml-1)*npes_yz, gid+(ml-1)*npes_yz, commnyz)
     687             :          enddo
     688           0 :       elseif (gid < mlt*npes_yz) then
     689           0 :          ktmod = gid/npes_yz
     690           0 :          call mpirecv(ictstuff, 4, mpiint, gid-ktmod*npes_yz, gid, commnyz)
     691           0 :          grid%jfirstct = ictstuff(1)
     692           0 :          grid%jlastct  = ictstuff(2)
     693           0 :          grid%kfirstct = ictstuff(3)
     694           0 :          grid%klastct  = ictstuff(4)
     695             :       end if
     696             :    end if
     697             : 
     698        1536 :    if (trac_decomp .gt. 1) then
     699           0 :       kquot = nq / trac_decomp
     700           0 :       krem = nq - kquot * trac_decomp
     701           0 :       krun = 0
     702           0 :       do kt = 1, krem
     703           0 :          grid%ktloa(kt) = krun + 1
     704           0 :          krun = krun + kquot + 1
     705           0 :          grid%kthia(kt) = krun
     706             :       enddo
     707           0 :       do kt = krem+1, trac_decomp
     708           0 :          grid%ktloa(kt) = krun + 1
     709           0 :          krun = krun + kquot
     710           0 :          grid%kthia(kt) = krun
     711             :       enddo
     712           0 :       ktmod = gid/npes_yz + 1
     713           0 :       ktmod = min(ktmod, trac_decomp)
     714           0 :       grid%ktlo = grid%ktloa(ktmod)
     715           0 :       grid%kthi = grid%kthia(ktmod)
     716             :    endif
     717             : 
     718        3072 : end subroutine spmd_vars_init
     719             : 
     720             : #endif
     721             : 
     722             : !========================================================================================
     723             : 
     724        1536 : subroutine grid_vars_init(pi, ae, om, dt, fft_flt, &
     725             :                           am_geom_crrct, grid)
     726             : 
     727             :    ! Initialize FV specific GRID vars
     728             :    ! 
     729             :    ! !REVISION HISTORY: 
     730             :    !   00.01.10    Grant        Creation using code from SJ Lin
     731             :    !   01.06.06    Sawyer       Modified for dynamics_vars
     732             :    !   04.08.25    Sawyer       Now updates GRID
     733             :    !   05.06.30    Sawyer       Added initializations from cd_core
     734             :    !   06.09.15    Sawyer       PI now passed as argument
     735             : 
     736             :    use pft_module, only: pftinit, pft_cf
     737             : 
     738             :    ! Arguments
     739             :    real(r8), intent(in) :: pi
     740             :    real(r8), intent(in) :: ae     ! radius of the earth (m)
     741             :    real(r8), intent(in) :: om     ! angular velocity of earth's rotation 
     742             :    real(r8), intent(in) :: dt
     743             :    integer,  intent(in) :: fft_flt
     744             :    logical,  intent(in) :: am_geom_crrct
     745             : 
     746             :    type( T_FVDYCORE_GRID ), intent(inout) :: grid
     747             : 
     748             :    ! local variables:
     749             :    integer :: im
     750             :    integer :: jm
     751             : 
     752             :    real(r8) :: ph5      ! This is to ensure 64-bit for any choice of r8
     753             : 
     754             :    integer  :: i, j, imh
     755             :    real(r8) :: zam5, zamda
     756             :    integer  :: js2g0, jn2g0, jn1g1, js2gc, jn1gc
     757             :    integer  :: js2gs, jn2gd, jn1gs
     758             : 
     759        1536 :    real(r8), pointer :: cosp(:), sinp(:), cose(:), sine(:), acosp(:), acosu(:)
     760        1536 :    real(r8), pointer :: coslon(:), sinlon(:), cosl5(:), sinl5(:)
     761             : 
     762             :    real(r8) :: rat, ycrit
     763             :    real(r8) :: dt5
     764             : 
     765             :    character(len=*), parameter :: sub='grid_vars_init'
     766             :    !---------------------------------------------------------------------------
     767             : 
     768        1536 :    im = grid%im
     769        1536 :    jm = grid%jm
     770             : 
     771        1536 :    grid%dl = (pi+pi)/im
     772        1536 :    grid%dp = pi/(jm-1)
     773             : 
     774        4608 :    allocate(grid%cosp(jm))
     775        3072 :    allocate(grid%sinp(jm))
     776        3072 :    allocate(grid%cose(jm))
     777        3072 :    allocate(grid%sine(jm))
     778        3072 :    allocate(grid%acosp(jm))
     779        3072 :    allocate(grid%acosu(jm))
     780             : 
     781        4608 :    allocate(grid%coslon(im))
     782        3072 :    allocate(grid%sinlon(im))
     783        3072 :    allocate(grid%cosl5(im))
     784        3072 :    allocate(grid%sinl5(im))
     785             : 
     786        1536 :    cosp => grid%cosp
     787        1536 :    sinp => grid%sinp
     788        1536 :    cose => grid%cose
     789        1536 :    sine => grid%sine
     790        1536 :    acosp  => grid%acosp
     791        1536 :    acosu  => grid%acosu
     792             : 
     793        1536 :    coslon => grid%coslon
     794        1536 :    sinlon => grid%sinlon
     795        1536 :    cosl5  => grid%cosl5
     796        1536 :    sinl5  => grid%sinl5
     797             : 
     798             :    ! philosophy below: edge values = local true values; centred values = area averages
     799             : 
     800      294912 :    do j = 2, jm
     801      293376 :       ph5  = -0.5_r8*pi + ((j-1)-0.5_r8)*(pi/(jm-1))
     802      294912 :       sine(j) = sin(ph5)
     803             :    end do
     804             : 
     805             :    ! cos(theta) at cell center distretized as
     806             :    !
     807             :    ! cos(theta) = d(sin(theta))/d(theta)
     808             : 
     809        1536 :    cosp( 1) =  0._r8
     810        1536 :    cosp(jm) =  0._r8
     811      293376 :    do j = 2, jm-1
     812      293376 :       cosp(j) = (sine(j+1)-sine(j)) / grid%dp
     813             :    end do
     814             : 
     815             :    ! Define cosine at edges..
     816             :    
     817        1536 :    if (am_geom_crrct) then 
     818           0 :       do j = 2, jm
     819           0 :          ph5     = -0.5_r8*pi + ((j-1)-0.5_r8)*(pi/(jm-1._r8))
     820           0 :          cose(j) = cos(ph5)
     821             :       end do
     822             :    else
     823      294912 :       do j = 2, jm
     824      294912 :          cose(j) = 0.5_r8 * (cosp(j-1) + cosp(j))  ! dsine/dpe between j+1 and j-1
     825             :       end do
     826             :    end if
     827        1536 :    cose(1) = cose(2)
     828             : 
     829      293376 :    do j = 2, jm-1
     830      293376 :       acosu(j) = 2._r8 / (cose(j) + cose(j+1))
     831             :    end do
     832             : 
     833        1536 :    sinp( 1) = -1._r8
     834        1536 :    sinp(jm) =  1._r8
     835        1536 :    if (am_geom_crrct) then 
     836           0 :       do j = 2, jm-1
     837           0 :          sinp(j) =        (cose(j) - cose(j+1))/grid%dp  ! sqrt(cosp^2+sinp^2)=1
     838             :       end do
     839             :    else
     840      293376 :       do j = 2, jm-1
     841      293376 :          sinp(j) = 0.5_r8 * (sine(j) + sine(j+1))        ! 2*sinp*cosp*dp=d(cose^2)
     842             :       end do
     843             :    end if
     844             : 
     845             :    ! Pole cap area and inverse
     846        1536 :    grid%acap = im*(1._r8+sine(2)) / grid%dp
     847        1536 :    grid%rcap = 1._r8 / grid%acap
     848             :  
     849        1536 :    imh = im/2
     850        1536 :    if (im /= 2*imh) then
     851           0 :       write(iulog,*) sub//': ERROR: im must be an even integer'
     852           0 :       call endrun(sub//': ERROR: im must be an even integer')
     853             :    end if
     854             :  
     855             :    ! Define logitude at the center of the volume
     856             :    ! i=1, Zamda = -pi
     857             :  
     858      222720 :    do i = 1, imh
     859      221184 :       zam5          = ((i-1)-0.5_r8) * grid%dl
     860      221184 :       cosl5(i)      =  cos(zam5)
     861      221184 :       cosl5(i+imh)  = -cosl5(i)
     862      221184 :       sinl5(i)      =  sin(zam5)
     863      221184 :       sinl5(i+imh)  = -sinl5(i)
     864      221184 :       zamda         = (i-1)*grid%dl
     865      221184 :       coslon(i)     =  cos(zamda)
     866      221184 :       coslon(i+imh) = -coslon(i)
     867      221184 :       sinlon(i)     =  sin(zamda)
     868      222720 :       sinlon(i+imh) = -sinlon(i)
     869             :    end do
     870             : 
     871        1536 :    acosp( 1) = grid%rcap * im
     872        1536 :    acosp(jm) = grid%rcap * im
     873      293376 :    do j = 2, jm-1
     874      293376 :       acosp(j) = 1._r8 / cosp(j)
     875             :    enddo
     876             : 
     877             :    ! cd_core initializations
     878             : 
     879        4608 :    allocate(grid%dtdx(jm))
     880        3072 :    allocate(grid%dtdx2(jm))
     881        3072 :    allocate(grid%dtdx4(jm))
     882        3072 :    allocate(grid%dtdxe(jm))
     883        3072 :    allocate(grid%dxdt(jm))
     884        3072 :    allocate(grid%dxe(jm))
     885        3072 :    allocate(grid%cye(jm))
     886        3072 :    allocate(grid%dycp(jm))
     887        3072 :    allocate(grid%rdxe(jm))
     888        3072 :    allocate(grid%txe5(jm))
     889        3072 :    allocate(grid%dtxe5(jm))
     890        3072 :    allocate(grid%dyce(jm))
     891        3072 :    allocate(grid%dx(jm))
     892        3072 :    allocate(grid%rdx(jm))
     893        3072 :    allocate(grid%cy(jm))
     894             : 
     895        1536 :    js2g0  = max(2,grid%jfirst)
     896        1536 :    jn2g0  = min(jm-1,grid%jlast)
     897        1536 :    jn1g1  = min(jm,grid%jlast+1)
     898        1536 :    js2gc  = max(2,grid%jfirst-grid%ng_c) ! NG lats on S (starting at 2)
     899        1536 :    jn1gc  = min(jm,grid%jlast+grid%ng_c) ! ng_c lats on N (ending at jm)
     900             : 
     901        1536 :    grid%js2g0  = js2g0
     902        1536 :    grid%jn2g0  = jn2g0
     903        1536 :    grid%jn1g1  = jn1g1
     904             : 
     905        1536 :    js2gs = max(2,grid%jfirst-grid%ng_s)
     906        1536 :    jn2gd = min(jm-1,grid%jlast+grid%ng_d)
     907        1536 :    jn1gs = min(jm,grid%jlast+grid%ng_s)
     908             : 
     909        4608 :    allocate(grid%sc(js2g0:jn2g0))
     910        4608 :    allocate(grid%se(js2g0:jn1g1))
     911        6144 :    allocate(grid%dc(im,js2g0:jn2g0))
     912        6144 :    allocate(grid%de(im,js2g0:jn1g1))
     913             : 
     914        4608 :    allocate(grid%scdiv4(js2gs:jn2gd))   !for filtering of u and v in div4 damping 
     915        4608 :    allocate(grid%sediv4(js2gs:jn1gs))   !for filtering of u and v in div4 damping 
     916        6144 :    allocate(grid%dcdiv4(im,js2gs:jn2gd))!for filtering of u and v in div4 damping 
     917        6144 :    allocate(grid%dediv4(im,js2gs:jn1gs))!for filtering of u and v in div4 damping 
     918             : 
     919        1536 :    call pftinit(im, fft_flt)
     920             : 
     921             :    ! Determine ycrit such that effective DX >= DY
     922        1536 :    rat = real(im,r8)/real(2*(jm-1),r8)
     923        1536 :    ycrit = acos( min(0.81_r8, rat) ) * (180._r8/pi)
     924             : 
     925             :    call pft_cf(im, jm, js2g0, jn2g0, jn1g1, &
     926             :                grid%sc, grid%se, grid%dc, grid%de,  &
     927        1536 :                grid%cosp, grid%cose, ycrit)
     928             : 
     929             :    !for filtering of u and v in div4 damping 
     930             :    call pft_cf(im, jm, js2gs, jn2gd, jn1gs,                             & 
     931             :                grid%scdiv4, grid%sediv4, grid%dcdiv4, grid%dediv4,      & 
     932        1536 :                grid%cosp, grid%cose, ycrit)                               
     933             : 
     934        6144 :    allocate( grid%cdx   (js2g0:jn1g1,grid%kfirst:grid%klast) )
     935        4608 :    allocate( grid%cdy   (js2g0:jn1g1,grid%kfirst:grid%klast) )
     936             : 
     937        4608 :    allocate( grid%cdx4  (js2g0:jn1g1,grid%kfirst:grid%klast) )!for div4 damping
     938        4608 :    allocate( grid%cdy4  (js2g0:jn1g1,grid%kfirst:grid%klast) )!for div4 damping
     939             : 
     940        4608 :    allocate( grid%cdxde (js2g0:jn1g1,grid%kfirst:grid%klast) )!for del2 damping
     941        4608 :    allocate( grid%cdxdp (js2g0:jn1g1,grid%kfirst:grid%klast) )!for del2 damping
     942        4608 :    allocate( grid%cdyde (js2g0:jn1g1,grid%kfirst:grid%klast) )!for del2 damping
     943        4608 :    allocate( grid%cdydp (js2g0:jn1g1,grid%kfirst:grid%klast) )!for del2 damping
     944             : 
     945        6144 :    allocate( grid%cdxdiv(jm,grid%kfirst:grid%klast) )!for div4 damping
     946        4608 :    allocate( grid%cdydiv(jm,grid%kfirst:grid%klast) )!for div4 damping
     947        4608 :    allocate( grid%cdtau4(js2g0:jn1g1,grid%kfirst:grid%klast) )!for div4 damping
     948             : 
     949        4608 :    allocate( grid%f0(grid%jfirst-grid%ng_s-1:grid%jlast+grid%ng_d) )
     950        4608 :    allocate( grid%fc(js2gc:jn1gc) )
     951             : 
     952       16704 :    do j = max(1,grid%jfirst-grid%ng_s-1), min(jm,grid%jlast+grid%ng_d)
     953       16704 :       grid%f0(j) = (om+om)*grid%sinp(j)
     954             :    end do
     955             : 
     956             :    ! Compute coriolis parameter at cell corners.
     957             : 
     958        1536 :    if (am_geom_crrct) then
     959           0 :       do j = js2gc, jn1gc
     960           0 :          grid%fc(j) = (om+om)*grid%sine(j)
     961             :       end do
     962             :    else
     963       12168 :       do j = js2gc, jn1gc                    ! Not the issue with ng_c = ng_d
     964       12168 :          grid%fc(j) = 0.5_r8*(grid%f0(j) + grid%f0(j-1))
     965             :       end do
     966             :    end if
     967             : 
     968        1536 :    grid%dt0 = 0._r8
     969        1536 :    dt5      = 0.5_r8*dt
     970             : 
     971        1536 :    grid%rdy   = 1._r8/(ae*grid%dp)
     972        1536 :    grid%dtdy  = dt *grid%rdy
     973        1536 :    grid%dtdy5 = dt5*grid%rdy
     974        1536 :    grid%dydt  = (ae*grid%dp) / dt
     975        1536 :    grid%tdy5  = 0.5_r8/grid%dtdy
     976             : 
     977        1536 : end subroutine grid_vars_init
     978             : 
     979             : !========================================================================================
     980             : 
     981           0 : subroutine dynamics_clean(grid)
     982             : 
     983             :    ! Arguments
     984             :    type(t_fvdycore_grid), intent(inout)  :: grid
     985             : 
     986             :    ! Temporary data structures
     987             : 
     988           0 :     if(associated(GRID%SINLON          )) deallocate(GRID%SINLON)
     989           0 :     if(associated(GRID%COSLON          )) deallocate(GRID%COSLON)
     990           0 :     if(associated(GRID%SINL5           )) deallocate(GRID%SINL5)
     991           0 :     if(associated(GRID%COSL5           )) deallocate(GRID%COSL5)
     992             : 
     993           0 :     if(associated(GRID%ACOSP           )) deallocate(GRID%ACOSP)
     994           0 :     if(associated(GRID%ACOSU           )) deallocate(GRID%ACOSU)
     995           0 :     if(associated(GRID%SINP            )) deallocate(GRID%SINP)
     996           0 :     if(associated(GRID%COSP            )) deallocate(GRID%COSP)
     997           0 :     if(associated(GRID%SINE            )) deallocate(GRID%SINE)
     998           0 :     if(associated(GRID%COSE            )) deallocate(GRID%COSE)
     999           0 :     if(associated(GRID%AK              )) deallocate(GRID%AK)
    1000           0 :     if(associated(GRID%BK              )) deallocate(GRID%BK)
    1001             : 
    1002             :     ! cd_core variables
    1003             : 
    1004           0 :     if(associated( grid%dtdx  )) deallocate(grid%dtdx)
    1005           0 :     if(associated( grid%dtdx2 )) deallocate(grid%dtdx2)
    1006           0 :     if(associated( grid%dtdx4 )) deallocate(grid%dtdx4)
    1007           0 :     if(associated( grid%dtdxe )) deallocate(grid%dtdxe)
    1008           0 :     if(associated( grid%dxdt  )) deallocate(grid%dxdt)
    1009           0 :     if(associated( grid%dxe   )) deallocate(grid%dxe)
    1010           0 :     if(associated( grid%cye   )) deallocate(grid%cye)
    1011           0 :     if(associated( grid%dycp  )) deallocate(grid%dycp)
    1012           0 :     if(associated( grid%rdxe  )) deallocate(grid%rdxe)
    1013           0 :     if(associated( grid%txe5  )) deallocate(grid%txe5)
    1014           0 :     if(associated( grid%dtxe5 )) deallocate(grid%dtxe5)
    1015           0 :     if(associated( grid%dyce  )) deallocate(grid%dyce)
    1016           0 :     if(associated( grid%dx    )) deallocate(grid%dx)
    1017           0 :     if(associated( grid%rdx   )) deallocate(grid%rdx)
    1018           0 :     if(associated( grid%cy    )) deallocate(grid%cy)
    1019             : 
    1020           0 :     if(associated( grid%sc    )) deallocate(grid%sc)
    1021           0 :     if(associated( grid%se    )) deallocate(grid%se)
    1022           0 :     if(associated( grid%dc    )) deallocate(grid%dc)
    1023           0 :     if(associated( grid%de    )) deallocate(grid%de)
    1024             : 
    1025           0 :     if(associated( grid%cdx   )) deallocate(grid%cdx)
    1026           0 :     if(associated( grid%cdy   )) deallocate(grid%cdy)
    1027           0 :     if(associated( grid%cdx4  )) deallocate(grid%cdx4)  
    1028           0 :     if(associated( grid%cdy4  )) deallocate(grid%cdy4)  
    1029           0 :     if(associated( grid%cdxde )) deallocate(grid%cdxde) 
    1030           0 :     if(associated( grid%cdxdp )) deallocate(grid%cdxdp) 
    1031           0 :     if(associated( grid%cdydp )) deallocate(grid%cdydp) 
    1032           0 :     if(associated( grid%cdyde )) deallocate(grid%cdyde) 
    1033           0 :     if(associated( grid%cdxdiv)) deallocate(grid%cdxdiv)
    1034           0 :     if(associated( grid%cdydiv)) deallocate(grid%cdydiv)
    1035           0 :     if(associated( grid%cdtau4)) deallocate(grid%cdtau4)
    1036             : 
    1037           0 :     if(associated( grid%scdiv4)) deallocate(grid%scdiv4)
    1038           0 :     if(associated( grid%sediv4)) deallocate(grid%sediv4)
    1039           0 :     if(associated( grid%dcdiv4)) deallocate(grid%dcdiv4)
    1040           0 :     if(associated( grid%dediv4)) deallocate(grid%dediv4)
    1041             : 
    1042           0 :     if(associated( grid%f0    )) deallocate(grid%f0)
    1043           0 :     if(associated( grid%fc    )) deallocate(grid%fc)
    1044             : 
    1045             : #if defined(SPMD)
    1046           0 :    call spmd_vars_clean(grid)
    1047             : #endif
    1048             : 
    1049           0 : end subroutine dynamics_clean
    1050             : 
    1051             : !========================================================================================
    1052             : 
    1053             : #if defined(SPMD)
    1054           0 : subroutine spmd_vars_clean(grid)
    1055             : 
    1056             :    use parutilitiesmodule, only : parpatternfree
    1057             : 
    1058             :    ! args
    1059             :    type (T_FVDYCORE_GRID), intent(inout) :: grid
    1060             :    !-----------------------------------------------------------------------
    1061             : 
    1062           0 :    if ( grid%twod_decomp == 1 ) then
    1063             : 
    1064             :       ! Clean transposes
    1065             : 
    1066           0 :       call parpatternfree(grid%commxy, grid%u_to_uxy)
    1067           0 :       call parpatternfree(grid%commxy, grid%uxy_to_u)
    1068           0 :       call parpatternfree(grid%commxy, grid%v_to_vxy)
    1069           0 :       call parpatternfree(grid%commxy, grid%vxy_to_v)
    1070           0 :       call parpatternfree(grid%commxy, grid%ijk_yz_to_xy)
    1071           0 :       call parpatternfree(grid%commxy, grid%ijk_xy_to_yz)
    1072           0 :       call parpatternfree(grid%commxy, grid%ikj_xy_to_yz)
    1073           0 :       call parpatternfree(grid%commxy, grid%ikj_yz_to_xy)
    1074           0 :       call parpatternfree(grid%commxy, grid%pe_to_pexy)
    1075           0 :       call parpatternfree(grid%commxy, grid%pexy_to_pe)
    1076           0 :       call parpatternfree(grid%commxy, grid%pt_to_ptxy)
    1077           0 :       call parpatternfree(grid%commxy, grid%ptxy_to_pt)
    1078           0 :       call parpatternfree(grid%commxy, grid%r4_xy_to_yz)
    1079           0 :       call parpatternfree(grid%commxy, grid%r4_yz_to_xy)
    1080           0 :       call parpatternfree(grid%commxy, grid%pkxy_to_pkc)
    1081           0 :       call parpatternfree(grid%commxy, grid%xy2d_to_yz2d)
    1082           0 :       call parpatternfree(grid%commxy, grid%yz2d_to_xy2d)
    1083             :    endif
    1084             : 
    1085           0 : end subroutine spmd_vars_clean
    1086             : #endif
    1087             : 
    1088             : !========================================================================================
    1089             : 
    1090           0 : end module dynamics_vars
    1091             : 

Generated by: LCOV version 1.14