LCOV - code coverage report
Current view: | top level - dynamics/se - dyn_comp.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 672 | 1012 | 66.4 %
Date: | 2024-12-17 17:57:11 |  | Functions: | 12 | 17 | 70.6 %

              Line data    Source code
           1             : module dyn_comp
           2             :
           3             : ! CAM interfaces to the SE Dynamical Core
           4             :
           5             : use shr_kind_mod,           only: r8=>shr_kind_r8, shr_kind_cl
           6             : use physconst,              only: pi
           7             : use spmd_utils,             only: iam, masterproc
           8             : use constituents,           only: pcnst, cnst_get_ind, cnst_name, cnst_longname, &
           9             :                                   cnst_read_iv, qmin, cnst_type, tottnam,        &
          10             :                                   cnst_is_a_water_species
          11             : use cam_control_mod,        only: initial_run
          12             : use cam_initfiles,          only: initial_file_get_id, topo_file_get_id, pertlim
          13             : use phys_control,           only: use_gw_front, use_gw_front_igw
          14             : use dyn_grid,               only: ini_grid_name, timelevel, hvcoord, edgebuf, &
          15             :                                   ini_grid_hdim_name
          16             :
          17             : use cam_grid_support,       only: cam_grid_id, cam_grid_get_gcid, &
          18             :                                   cam_grid_dimensions, &
          19             :                                   cam_grid_get_latvals, cam_grid_get_lonvals,  &
          20             :                                   max_hcoordname_len
          21             : use cam_map_utils,          only: iMap
          22             :
          23             : use inic_analytic,          only: analytic_ic_active, analytic_ic_set_ic
          24             : use dyn_tests_utils,        only: vcoord=>vc_dry_pressure
          25             :
          26             : use cam_history,            only: outfld, hist_fld_active, fieldname_len
          27             : use cam_history_support,    only: max_fieldname_len
          28             : use time_manager,           only: get_step_size
          29             :
          30             : use ncdio_atm,              only: infld
          31             : use pio,                    only: file_desc_t, pio_seterrorhandling, PIO_BCAST_ERROR, &
          32             :                                   pio_inq_dimid, pio_inq_dimlen, PIO_NOERR
          33             :
          34             : use infnan,                 only: isnan
          35             : use cam_logfile,            only: iulog
          36             : use cam_abortutils,         only: endrun
          37             : use shr_sys_mod,            only: shr_sys_flush
          38             :
          39             : use parallel_mod,           only: par
          40             : use hybrid_mod,             only: hybrid_t
          41             : use dimensions_mod,         only: nelemd, nlev, np, npsq, ntrac, nc, fv_nphys
          42             : use dimensions_mod,         only: qsize, use_cslam
          43             : use element_mod,            only: element_t, elem_state_t
          44             : use fvm_control_volume_mod, only: fvm_struct
          45             : use se_dyn_time_mod,        only: nsplit
          46             : use edge_mod,               only: initEdgeBuffer, edgeVpack, edgeVunpack, FreeEdgeBuffer
          47             : use edgetype_mod,           only: EdgeBuffer_t
          48             : use bndry_mod,              only: bndry_exchange
          49             : use se_single_column_mod,   only: scm_setinitial
          50             : use scamMod,                only: single_column, readiopdata, use_iop, setiopupdate_init
          51             : use hycoef,                 only: hyai, hybi, ps0
          52             :
          53             : implicit none
          54             : private
          55             : save
          56             :
          57             : public ::          &
          58             :      dyn_import_t, &
          59             :      dyn_export_t, &
          60             :      dyn_readnl,   &
          61             :      dyn_register, &
          62             :      dyn_init,     &
          63             :      dyn_run,      &
          64             :      dyn_final
          65             :
          66             : type dyn_import_t
          67             :   type (element_t),  pointer :: elem(:) => null()
          68             :   type (fvm_struct), pointer :: fvm(:) => null()
          69             : end type dyn_import_t
          70             :
          71             : type dyn_export_t
          72             :   type (element_t),  pointer :: elem(:) => null()
          73             :   type (fvm_struct), pointer :: fvm(:) => null()
          74             : end type dyn_export_t
          75             :
          76             : ! Namelist
          77             : logical, public, protected :: write_restart_unstruct
          78             :
          79             : ! Frontogenesis indices
          80             : integer, public    :: frontgf_idx      = -1
          81             : integer, public    :: frontga_idx      = -1
          82             :
          83             : interface read_dyn_var
          84             :   module procedure read_dyn_field_2d
          85             :   module procedure read_dyn_field_3d
          86             : end interface read_dyn_var
          87             :
          88             : real(r8), parameter :: rad2deg = 180.0_r8 / pi
          89             : real(r8), parameter :: deg2rad = pi / 180.0_r8
          90             : real(r8), parameter :: rarea_sphere = 1.0_r8 / (4.0_r8*PI)
          91             :
          92             : !===============================================================================
          93             : contains
          94             : !===============================================================================
          95             :
          96        1536 : subroutine dyn_readnl(NLFileName)
          97             :    use air_composition,only: thermodynamic_active_species_num
          98             :    use namelist_utils, only: find_group_name
          99             :    use namelist_mod,   only: homme_set_defaults, homme_postprocess_namelist
         100             :    use units,          only: getunit, freeunit
         101             :    use spmd_utils,     only: masterproc, masterprocid, mpicom, npes
         102             :    use spmd_utils,     only: mpi_real8, mpi_integer, mpi_character, mpi_logical
         103             :    use dyn_grid,       only: se_write_grid_file, se_grid_filename, se_write_gll_corners
         104             :    use dp_mapping,     only: nphys_pts
         105             :    use native_mapping, only: native_mapping_readnl
         106             :
         107             :    use control_mod,    only: hypervis_subcycle, hypervis_subcycle_sponge
         108             :    use control_mod,    only: hypervis_subcycle_q, statefreq, runtype
         109             :    use control_mod,    only: nu, nu_div, nu_p, nu_q, nu_top, qsplit, rsplit
         110             :    use control_mod,    only: vert_remap_uvTq_alg, vert_remap_tracer_alg
         111             :    use control_mod,    only: tstep_type, rk_stage_user
         112             :    use control_mod,    only: ftype, limiter_option, partmethod
         113             :    use control_mod,    only: topology, variable_nsplit
         114             :    use control_mod,    only: fine_ne, hypervis_power, hypervis_scaling
         115             :    use control_mod,    only: max_hypervis_courant, statediag_numtrac,refined_mesh
         116             :    use control_mod,    only: molecular_diff, pgf_formulation, dribble_in_rsplit_loop
         117             :    use control_mod,    only: sponge_del4_nu_div_fac, sponge_del4_nu_fac, sponge_del4_lev
         118             :    use dimensions_mod, only: ne, npart
         119             :    use dimensions_mod, only: large_Courant_incr
         120             :    use dimensions_mod, only: fvm_supercycling, fvm_supercycling_jet
         121             :    use dimensions_mod, only: kmin_jet, kmax_jet
         122             :    use params_mod,     only: SFCURVE
         123             :    use parallel_mod,   only: initmpi
         124             :    use thread_mod,     only: initomp, max_num_threads
         125             :    use thread_mod,     only: horz_num_threads, vert_num_threads, tracer_num_threads
         126             :    ! Dummy argument
         127             :    character(len=*), intent(in) :: NLFileName
         128             :
         129             :    ! Local variables
         130             :    integer                      :: unitn, ierr,k
         131             :
         132             :    ! SE Namelist variables
         133             :    integer                      :: se_fine_ne
         134             :    integer                      :: se_ftype
         135             :    integer                      :: se_statediag_numtrac
         136             :    integer                      :: se_fv_nphys
         137             :    real(r8)                     :: se_hypervis_power
         138             :    real(r8)                     :: se_hypervis_scaling
         139             :    integer                      :: se_hypervis_subcycle
         140             :    integer                      :: se_hypervis_subcycle_sponge
         141             :    integer                      :: se_hypervis_subcycle_q
         142             :    integer                      :: se_limiter_option
         143             :    real(r8)                     :: se_max_hypervis_courant
         144             :    character(len=SHR_KIND_CL)   :: se_mesh_file
         145             :    integer                      :: se_ne
         146             :    integer                      :: se_npes
         147             :    integer                      :: se_nsplit
         148             :    real(r8)                     :: se_nu
         149             :    real(r8)                     :: se_nu_div
         150             :    real(r8)                     :: se_nu_p
         151             :    real(r8)                     :: se_nu_top
         152             :    real(r8)                     :: se_sponge_del4_nu_fac
         153             :    real(r8)                     :: se_sponge_del4_nu_div_fac
         154             :    integer                      :: se_sponge_del4_lev
         155             :    integer                      :: se_qsplit
         156             :    logical                      :: se_refined_mesh
         157             :    integer                      :: se_rsplit
         158             :    integer                      :: se_statefreq
         159             :    integer                      :: se_tstep_type
         160             :    character(len=32)            :: se_vert_remap_T
         161             :    character(len=32)            :: se_vert_remap_uvTq_alg
         162             :    character(len=32)            :: se_vert_remap_tracer_alg
         163             :    integer                      :: se_horz_num_threads
         164             :    integer                      :: se_vert_num_threads
         165             :    integer                      :: se_tracer_num_threads
         166             :    logical                      :: se_write_restart_unstruct
         167             :    logical                      :: se_large_Courant_incr
         168             :    integer                      :: se_fvm_supercycling
         169             :    integer                      :: se_fvm_supercycling_jet
         170             :    integer                      :: se_kmin_jet
         171             :    integer                      :: se_kmax_jet
         172             :    real(r8)                     :: se_molecular_diff
         173             :    integer                      :: se_pgf_formulation
         174             :    integer                      :: se_dribble_in_rsplit_loop
         175             :    namelist /dyn_se_inparm/        &
         176             :       se_fine_ne,                  & ! For refined meshes
         177             :       se_ftype,                    & ! forcing type
         178             :       se_statediag_numtrac,        &
         179             :       se_fv_nphys,                 &
         180             :       se_hypervis_power,           &
         181             :       se_hypervis_scaling,         &
         182             :       se_hypervis_subcycle,        &
         183             :       se_hypervis_subcycle_sponge, &
         184             :       se_hypervis_subcycle_q,      &
         185             :       se_limiter_option,           &
         186             :       se_max_hypervis_courant,     &
         187             :       se_mesh_file,                & ! Refined mesh definition file
         188             :       se_ne,                       &
         189             :       se_npes,                     &
         190             :       se_nsplit,                   & ! # of dyn steps per physics timestep
         191             :       se_nu,                       &
         192             :       se_nu_div,                   &
         193             :       se_nu_p,                     &
         194             :       se_nu_top,                   &
         195             :       se_sponge_del4_nu_fac,       &
         196             :       se_sponge_del4_nu_div_fac,   &
         197             :       se_sponge_del4_lev,          &
         198             :       se_qsplit,                   &
         199             :       se_refined_mesh,             &
         200             :       se_rsplit,                   &
         201             :       se_statefreq,                & ! number of steps per printstate call
         202             :       se_tstep_type,               &
         203             :       se_vert_remap_T,             &
         204             :       se_vert_remap_uvTq_alg,      &
         205             :       se_vert_remap_tracer_alg,    &
         206             :       se_write_grid_file,          &
         207             :       se_grid_filename,            &
         208             :       se_write_gll_corners,        &
         209             :       se_horz_num_threads,         &
         210             :       se_vert_num_threads,         &
         211             :       se_tracer_num_threads,       &
         212             :       se_write_restart_unstruct,   &
         213             :       se_large_Courant_incr,       &
         214             :       se_fvm_supercycling,         &
         215             :       se_fvm_supercycling_jet,     &
         216             :       se_kmin_jet,                 &
         217             :       se_kmax_jet,                 &
         218             :       se_molecular_diff,           &
         219             :       se_pgf_formulation,          &
         220             :       se_dribble_in_rsplit_loop
         221             :    !--------------------------------------------------------------------------
         222             :
         223             :    ! defaults for variables not set by build-namelist
         224        1536 :    se_fine_ne                  = -1
         225        1536 :    se_hypervis_power           = 0
         226        1536 :    se_hypervis_scaling         = 0
         227        1536 :    se_max_hypervis_courant     = 1.0e99_r8
         228        1536 :    se_mesh_file                = ''
         229        1536 :    se_npes                     = npes
         230        1536 :    se_write_restart_unstruct   = .false.
         231             :
         232             :    ! Read the namelist (dyn_se_inparm)
         233        1536 :    call MPI_barrier(mpicom, ierr)
         234        1536 :    if (masterproc) then
         235           2 :       write(iulog, *) "dyn_readnl: reading dyn_se_inparm namelist..."
         236           2 :       unitn = getunit()
         237           2 :       open( unitn, file=trim(NLFileName), status='old' )
         238           2 :       call find_group_name(unitn, 'dyn_se_inparm', status=ierr)
         239           2 :       if (ierr == 0) then
         240           2 :          read(unitn, dyn_se_inparm, iostat=ierr)
         241           2 :          if (ierr /= 0) then
         242           0 :             call endrun('dyn_readnl: ERROR reading dyn_se_inparm namelist')
         243             :          end if
         244             :       end if
         245           2 :       close(unitn)
         246           2 :       call freeunit(unitn)
         247             :    end if
         248             :
         249             :    ! Broadcast namelist values to all PEs
         250        1536 :    call MPI_bcast(se_fine_ne, 1, mpi_integer, masterprocid, mpicom, ierr)
         251        1536 :    call MPI_bcast(se_ftype, 1, mpi_integer, masterprocid, mpicom, ierr)
         252        1536 :    call MPI_bcast(se_statediag_numtrac, 1, mpi_integer, masterprocid, mpicom, ierr)
         253        1536 :    call MPI_bcast(se_hypervis_power, 1, mpi_real8, masterprocid, mpicom, ierr)
         254        1536 :    call MPI_bcast(se_hypervis_scaling, 1, mpi_real8, masterprocid, mpicom, ierr)
         255        1536 :    call MPI_bcast(se_hypervis_subcycle, 1, mpi_integer, masterprocid, mpicom, ierr)
         256        1536 :    call MPI_bcast(se_hypervis_subcycle_sponge, 1, mpi_integer, masterprocid, mpicom, ierr)
         257        1536 :    call MPI_bcast(se_hypervis_subcycle_q, 1, mpi_integer, masterprocid, mpicom, ierr)
         258        1536 :    call MPI_bcast(se_limiter_option, 1, mpi_integer, masterprocid, mpicom, ierr)
         259        1536 :    call MPI_bcast(se_max_hypervis_courant, 1, mpi_real8, masterprocid, mpicom, ierr)
         260        1536 :    call MPI_bcast(se_mesh_file, SHR_KIND_CL,  mpi_character, masterprocid, mpicom, ierr)
         261        1536 :    call MPI_bcast(se_ne, 1, mpi_integer, masterprocid, mpicom, ierr)
         262        1536 :    call MPI_bcast(se_npes, 1, mpi_integer, masterprocid, mpicom, ierr)
         263        1536 :    call MPI_bcast(se_nsplit, 1, mpi_integer, masterprocid, mpicom, ierr)
         264        1536 :    call MPI_bcast(se_nu, 1, mpi_real8, masterprocid, mpicom, ierr)
         265        1536 :    call MPI_bcast(se_nu_div, 1, mpi_real8, masterprocid, mpicom, ierr)
         266        1536 :    call MPI_bcast(se_nu_p, 1, mpi_real8, masterprocid, mpicom, ierr)
         267        1536 :    call MPI_bcast(se_nu_top, 1, mpi_real8, masterprocid, mpicom, ierr)
         268        1536 :    call MPI_bcast(se_sponge_del4_nu_fac, 1, mpi_real8, masterprocid, mpicom, ierr)
         269        1536 :    call MPI_bcast(se_sponge_del4_nu_div_fac, 1, mpi_real8, masterprocid, mpicom, ierr)
         270        1536 :    call MPI_bcast(se_sponge_del4_lev, 1, mpi_integer, masterprocid, mpicom, ierr)
         271        1536 :    call MPI_bcast(se_qsplit, 1, mpi_integer, masterprocid, mpicom, ierr)
         272        1536 :    call MPI_bcast(se_refined_mesh, 1, mpi_logical, masterprocid, mpicom, ierr)
         273        1536 :    call MPI_bcast(se_rsplit, 1, mpi_integer, masterprocid, mpicom, ierr)
         274        1536 :    call MPI_bcast(se_statefreq, 1, mpi_integer, masterprocid, mpicom, ierr)
         275        1536 :    call MPI_bcast(se_tstep_type, 1, mpi_integer, masterprocid, mpicom, ierr)
         276        1536 :    call MPI_bcast(se_vert_remap_T, 32, mpi_character, masterprocid, mpicom, ierr)
         277        1536 :    call MPI_bcast(se_vert_remap_uvTq_alg, 32, mpi_character, masterprocid, mpicom, ierr)
         278        1536 :    call MPI_bcast(se_vert_remap_tracer_alg, 32, mpi_character, masterprocid, mpicom, ierr)
         279        1536 :    call MPI_bcast(se_fv_nphys, 1, mpi_integer, masterprocid, mpicom, ierr)
         280        1536 :    call MPI_bcast(se_write_grid_file, 16,  mpi_character, masterprocid, mpicom, ierr)
         281        1536 :    call MPI_bcast(se_grid_filename, shr_kind_cl, mpi_character, masterprocid, mpicom, ierr)
         282        1536 :    call MPI_bcast(se_write_gll_corners, 1,  mpi_logical, masterprocid, mpicom, ierr)
         283        1536 :    call MPI_bcast(se_horz_num_threads, 1, MPI_integer, masterprocid, mpicom,ierr)
         284        1536 :    call MPI_bcast(se_vert_num_threads, 1, MPI_integer, masterprocid, mpicom,ierr)
         285        1536 :    call MPI_bcast(se_tracer_num_threads, 1, MPI_integer, masterprocid, mpicom,ierr)
         286        1536 :    call MPI_bcast(se_write_restart_unstruct, 1, mpi_logical, masterprocid, mpicom, ierr)
         287        1536 :    call MPI_bcast(se_large_Courant_incr, 1, mpi_logical, masterprocid, mpicom, ierr)
         288        1536 :    call MPI_bcast(se_fvm_supercycling, 1, mpi_integer, masterprocid, mpicom, ierr)
         289        1536 :    call MPI_bcast(se_fvm_supercycling_jet, 1, mpi_integer, masterprocid, mpicom, ierr)
         290        1536 :    call MPI_bcast(se_kmin_jet, 1, mpi_integer, masterprocid, mpicom, ierr)
         291        1536 :    call MPI_bcast(se_kmax_jet, 1, mpi_integer, masterprocid, mpicom, ierr)
         292        1536 :    call MPI_bcast(se_molecular_diff, 1, mpi_real8, masterprocid, mpicom, ierr)
         293        1536 :    call MPI_bcast(se_pgf_formulation, 1, mpi_integer, masterprocid, mpicom, ierr)
         294        1536 :    call MPI_bcast(se_dribble_in_rsplit_loop, 1, mpi_integer, masterprocid, mpicom, ierr)
         295        1536 :    if (se_npes <= 0) then
         296           0 :       call endrun('dyn_readnl: ERROR: se_npes must be > 0')
         297             :    end if
         298             :
         299             :    ! Initialize the SE structure that holds the MPI decomposition information
         300        1536 :    par = initmpi(se_npes)
         301        1536 :    call initomp()
         302             :
         303             :
         304        1536 :    if (se_fvm_supercycling < 0) se_fvm_supercycling = se_rsplit
         305        1536 :    if (se_fvm_supercycling_jet < 0) se_fvm_supercycling_jet = se_rsplit
         306             :
         307             :    ! Go ahead and enforce ne = 0 for refined mesh runs
         308        1536 :    if (se_refined_mesh) then
         309           0 :       se_ne = 0
         310             :    end if
         311             :
         312             :    ! Set HOMME defaults
         313        1536 :    call homme_set_defaults()
         314             :    ! Set HOMME variables not in CAM's namelist but with different CAM defaults
         315        1536 :    partmethod               = SFCURVE
         316        1536 :    npart                    = se_npes
         317             :    ! CAM requires forward-in-time, subcycled dynamics
         318             :    ! RK2 3 stage tracers, sign-preserving conservative
         319        1536 :    rk_stage_user            = 3
         320        1536 :    topology                 = "cube"
         321             :    ! Finally, set the HOMME variables which have different names
         322        1536 :    fine_ne                  = se_fine_ne
         323        1536 :    ftype                    = se_ftype
         324        1536 :    statediag_numtrac        = MIN(se_statediag_numtrac,pcnst)
         325        1536 :    hypervis_power           = se_hypervis_power
         326        1536 :    hypervis_scaling         = se_hypervis_scaling
         327        1536 :    hypervis_subcycle        = se_hypervis_subcycle
         328        1536 :    if (hypervis_subcycle_sponge<0) then
         329           0 :      hypervis_subcycle_sponge = hypervis_subcycle
         330             :    else
         331        1536 :      hypervis_subcycle_sponge = se_hypervis_subcycle_sponge
         332             :    end if
         333        1536 :    hypervis_subcycle_q      = se_hypervis_subcycle_q
         334        1536 :    limiter_option           = se_limiter_option
         335        1536 :    max_hypervis_courant     = se_max_hypervis_courant
         336        1536 :    refined_mesh             = se_refined_mesh
         337        1536 :    ne                       = se_ne
         338        1536 :    nsplit                   = se_nsplit
         339        1536 :    nu                       = se_nu
         340        1536 :    nu_div                   = se_nu_div
         341        1536 :    nu_p                     = se_nu_p
         342        1536 :    nu_q                     = se_nu_p !for tracer-wind consistency nu_q must me equal to nu_p
         343        1536 :    nu_top                   = se_nu_top
         344        1536 :    sponge_del4_nu_fac       = se_sponge_del4_nu_fac
         345        1536 :    sponge_del4_nu_div_fac   = se_sponge_del4_nu_div_fac
         346        1536 :    sponge_del4_lev          = se_sponge_del4_lev
         347        1536 :    qsplit                   = se_qsplit
         348        1536 :    rsplit                   = se_rsplit
         349        1536 :    statefreq                = se_statefreq
         350        1536 :    tstep_type               = se_tstep_type
         351        1536 :    vert_remap_uvTq_alg      = set_vert_remap(se_vert_remap_T, se_vert_remap_uvTq_alg)
         352        1536 :    vert_remap_tracer_alg    = set_vert_remap(se_vert_remap_T, se_vert_remap_tracer_alg)
         353        1536 :    fv_nphys                 = se_fv_nphys
         354        1536 :    large_Courant_incr       = se_large_Courant_incr
         355        1536 :    fvm_supercycling         = se_fvm_supercycling
         356        1536 :    fvm_supercycling_jet     = se_fvm_supercycling_jet
         357        1536 :    kmin_jet                 = se_kmin_jet
         358        1536 :    kmax_jet                 = se_kmax_jet
         359        1536 :    variable_nsplit          = .false.
         360        1536 :    molecular_diff           = se_molecular_diff
         361        1536 :    pgf_formulation          = se_pgf_formulation
         362        1536 :    dribble_in_rsplit_loop   = se_dribble_in_rsplit_loop
         363        1536 :    if (fv_nphys > 0) then
         364             :       ! Use finite volume physics grid and CSLAM for tracer advection
         365        1536 :       nphys_pts = fv_nphys*fv_nphys
         366        1536 :       qsize = thermodynamic_active_species_num ! number tracers advected by GLL
         367        1536 :       ntrac = pcnst                            ! number tracers advected by CSLAM
         368        1536 :       use_cslam = .true.
         369             :    else
         370             :       ! Use GLL grid for physics and tracer advection
         371           0 :       nphys_pts = npsq
         372           0 :       qsize = pcnst
         373           0 :       ntrac = 0
         374           0 :       use_cslam = .false.
         375             :    end if
         376             :
         377        1536 :    if (rsplit < 1) then
         378           0 :       call endrun('dyn_readnl: rsplit must be > 0')
         379             :    end if
         380             :
         381             :    ! if restart or branch run
         382        1536 :    if (.not. initial_run) then
         383         768 :       runtype = 1
         384             :    end if
         385             :
         386             :    ! HOMME wants 'none' to indicate no mesh file
         387        1536 :    if (len_trim(se_mesh_file) == 0) then
         388        1536 :       se_mesh_file = 'none'
         389        1536 :       if (se_refined_mesh) then
         390           0 :          call endrun('dyn_readnl ERROR: se_refined_mesh=.true. but no se_mesh_file')
         391             :       end if
         392             :    end if
         393        1536 :    call homme_postprocess_namelist(se_mesh_file, par)
         394             :
         395             :    ! Set threading numbers to reasonable values
         396        1536 :    if ((se_horz_num_threads == 0) .and. (se_vert_num_threads == 0) .and. (se_tracer_num_threads == 0)) then
         397             :       ! The user has not set any threading values, choose defaults
         398        1536 :       se_horz_num_threads = 1
         399        1536 :       se_vert_num_threads = max_num_threads
         400        1536 :       se_tracer_num_threads = se_vert_num_threads
         401             :    end if
         402        1536 :    if (se_horz_num_threads < 1) then
         403           0 :       se_horz_num_threads = 1
         404             :    end if
         405        1536 :    if (se_vert_num_threads < 1) then
         406           0 :       se_vert_num_threads = 1
         407             :    end if
         408        1536 :    if (se_tracer_num_threads < 1) then
         409           0 :       se_tracer_num_threads = 1
         410             :    end if
         411        1536 :    horz_num_threads = se_horz_num_threads
         412        1536 :    vert_num_threads = se_vert_num_threads
         413        1536 :    tracer_num_threads = se_tracer_num_threads
         414             :
         415        1536 :    write_restart_unstruct = se_write_restart_unstruct
         416             :
         417        1536 :    if (se_kmin_jet<0            ) kmin_jet             = 1
         418        1536 :    if (se_kmax_jet<0            ) kmax_jet             = nlev
         419             :
         420        1536 :    if (masterproc) then
         421           2 :       write(iulog, '(a,i0)')   'dyn_readnl: se_ftype                    = ',ftype
         422           2 :       write(iulog, '(a,i0)')   'dyn_readnl: se_statediag_numtrac        = ',statediag_numtrac
         423           2 :       write(iulog, '(a,i0)')   'dyn_readnl: se_hypervis_subcycle        = ',se_hypervis_subcycle
         424           2 :       write(iulog, '(a,i0)')   'dyn_readnl: se_hypervis_subcycle_sponge = ',se_hypervis_subcycle_sponge
         425           2 :       write(iulog, '(a,i0)')   'dyn_readnl: se_hypervis_subcycle_q      = ',se_hypervis_subcycle_q
         426           2 :       write(iulog, '(a,l4)')   'dyn_readnl: se_large_Courant_incr       = ',se_large_Courant_incr
         427           2 :       write(iulog, '(a,i0)')   'dyn_readnl: se_limiter_option           = ',se_limiter_option
         428           2 :       if (.not. se_refined_mesh) then
         429           2 :          write(iulog, '(a,i0)')'dyn_readnl: se_ne                       = ',se_ne
         430             :       end if
         431           2 :       write(iulog, '(a,i0)')   'dyn_readnl: se_npes                     = ',se_npes
         432           2 :       write(iulog, '(a,i0)')   'dyn_readnl: se_nsplit                   = ',se_nsplit
         433             :       !
         434             :       ! se_nu<0 then coefficients are set automatically in module global_norms_mod
         435             :       !
         436           2 :       if (se_nu_div>0) &
         437           0 :            write(iulog, '(a,e9.2)') 'dyn_readnl: se_nu                       = ',se_nu
         438           2 :       if (se_nu_div>0) &
         439           0 :            write(iulog, '(a,e9.2)') 'dyn_readnl: se_nu_div                   = ',se_nu_div
         440           2 :       if (se_nu_p>0) then
         441           0 :         write(iulog, '(a,e9.2)') 'dyn_readnl: se_nu_p                     = ',se_nu_p
         442           0 :         write(iulog, '(a)') 'Note that nu_q must be the same as nu_p for  mass / tracer inconsistency'
         443             :       end if
         444           2 :       write(iulog, '(a,e9.2)') 'dyn_readnl: se_nu_top                     = ',se_nu_top
         445           2 :       write(iulog, '(a,i0)')   'dyn_readnl: se_qsplit                     = ',se_qsplit
         446           2 :       write(iulog, '(a,i0)')   'dyn_readnl: se_rsplit                     = ',se_rsplit
         447           2 :       write(iulog, '(a,i0)')   'dyn_readnl: se_statefreq                  = ',se_statefreq
         448           2 :       write(iulog, '(a,i0)')   'dyn_readnl: se_tstep_type                 = ',se_tstep_type
         449           2 :       write(iulog, '(a,a)')    'dyn_readnl: se_vert_remap_T               = ',trim(se_vert_remap_T)
         450           2 :       write(iulog, '(a,a)')    'dyn_readnl: se_vert_remap_uvTq_alg        = ',trim(se_vert_remap_uvTq_alg)
         451           2 :       write(iulog, '(a,a)')    'dyn_readnl: se_vert_remap_tracer_alg      = ',trim(se_vert_remap_tracer_alg)
         452           2 :       write(iulog, '(a,i0)')   'dyn_readnl: se_fvm_supercycling           = ',fvm_supercycling
         453           2 :       write(iulog, '(a,i0)')   'dyn_readnl: se_fvm_supercycling_jet       = ',fvm_supercycling_jet
         454           2 :       write(iulog, '(a,i0)')   'dyn_readnl: se_kmin_jet                   = ',kmin_jet
         455           2 :       write(iulog, '(a,i0)')   'dyn_readnl: se_kmax_jet                   = ',kmax_jet
         456             :
         457           2 :       write(iulog, *)   'dyn_readnl: se_sponge_del4_nu_fac         = ',se_sponge_del4_nu_fac
         458           2 :       if (se_sponge_del4_nu_fac < 0) write(iulog, '(a)')   ' (automatically set based on model top location)'
         459           2 :       write(iulog, *)   'dyn_readnl: se_sponge_del4_nu_div_fac     = ',se_sponge_del4_nu_div_fac
         460           2 :       if (se_sponge_del4_nu_div_fac < 0)  write(iulog, '(a)')   ' (automatically set based on model top location)'
         461           2 :       write(iulog, *)   'dyn_readnl: se_sponge_del4_lev            = ',se_sponge_del4_lev
         462           2 :       if (se_sponge_del4_lev < 0)  write(iulog, '(a)')   ' (automatically set based on model top location)'
         463             :
         464           2 :       if (se_refined_mesh) then
         465           0 :          write(iulog, '(a)') 'dyn_readnl: Refined mesh simulation'
         466           0 :          write(iulog, '(a)') 'dyn_readnl: se_mesh_file = ',trim(se_mesh_file)
         467           0 :          if (hypervis_power /= 0) then
         468           0 :            write(iulog, '(a)') 'Using scalar viscosity (Zarzycki et al 2014 JClim)'
         469           0 :            write(iulog, '(a,e11.4)') 'dyn_readnl: se_hypervis_power = ',se_hypervis_power, ', (tensor hyperviscosity)'
         470           0 :            write(iulog, '(a,e11.4)') 'dyn_readnl: se_max_hypervis_courant = ',se_max_hypervis_courant
         471             :          end if
         472           0 :          if (hypervis_scaling /= 0) then
         473           0 :            write(iulog, '(a)') 'Using tensor viscosity (Guba et al., 2014)'
         474           0 :            write(iulog, '(a,e11.4)') 'dyn_readnl: se_hypervis_scaling = ',se_hypervis_scaling
         475             :          end if
         476             :       end if
         477             :
         478           2 :       if (use_cslam) then
         479           2 :          write(iulog, '(a)') 'dyn_readnl: physics will run on FVM points; advection by CSLAM'
         480           2 :          write(iulog,'(a,i0)') 'dyn_readnl: se_fv_nphys = ', fv_nphys
         481             :       else
         482           0 :          write(iulog, '(a)') 'dyn_readnl: physics will run on SE GLL points'
         483             :       end if
         484           2 :       write(iulog, '(a,i0)') 'dyn_readnl: se_horz_num_threads = ',horz_num_threads
         485           2 :       write(iulog, '(a,i0)') 'dyn_readnl: se_vert_num_threads = ',vert_num_threads
         486           2 :       write(iulog, '(a,i0)') 'dyn_readnl: se_tracer_num_threads = ',tracer_num_threads
         487           2 :       if (trim(se_write_grid_file) == 'SCRIP') then
         488           0 :          write(iulog,'(2a)') "dyn_readnl: write SCRIP grid file = ", trim(se_grid_filename)
         489             :       else
         490           2 :          write(iulog,'(a)') "dyn_readnl: do not write grid file"
         491             :       end if
         492           2 :       write(iulog,'(a,l1)') 'dyn_readnl: write gll corners to SEMapping.nc = ', &
         493           4 :                             se_write_gll_corners
         494           2 :       write(iulog,'(a,l1)') 'dyn_readnl: write restart data on unstructured grid = ', &
         495           4 :                             se_write_restart_unstruct
         496             :
         497           2 :       write(iulog, '(a,e9.2)') 'dyn_readnl: se_molecular_diff  = ', molecular_diff
         498             :    end if
         499             :
         500        1536 :    call native_mapping_readnl(NLFileName)
         501             :
         502             :    !---------------------------------------------------------------------------
         503             :    contains
         504             :    !---------------------------------------------------------------------------
         505             :
         506        3072 :       integer function set_vert_remap( remap_T, remap_alg )
         507             :
         508             :          ! Convert namelist input strings to the internally used integers.
         509             :
         510             :          character(len=*), intent(in) :: remap_T    ! scheme for remapping temperature
         511             :          character(len=*), intent(in) :: remap_alg  ! remapping algorithm
         512             :
         513             :          ! check valid remap_T values:
         514        3072 :          if (remap_T /= 'thermal_energy_over_P' .and. remap_T /= 'Tv_over_logP') then
         515           0 :             write(iulog,*)'set_vert_remap: invalid remap_T= ',trim(remap_T)
         516           0 :             call endrun('set_vert_remap: invalid remap_T')
         517             :          end if
         518             :
         519             :          select case (remap_alg)
         520             :          case ('PPM_bc_mirror')
         521           0 :             set_vert_remap = 1
         522             :          case ('PPM_bc_PCoM')
         523           0 :             set_vert_remap = 2
         524             :          case ('PPM_bc_linear_extrapolation')
         525        1536 :             set_vert_remap = 10
         526             :          case ('FV3_PPM')
         527           0 :             if (remap_T == 'thermal_energy_over_P') then
         528             :                set_vert_remap = -4
         529             :             else
         530           0 :                set_vert_remap = -40
         531             :             end if
         532             :          case ('FV3_CS')
         533        1536 :             if (remap_T == 'thermal_energy_over_P') then
         534             :                set_vert_remap = -9
         535             :             else
         536           0 :                set_vert_remap = -90
         537             :             end if
         538             :          case ('FV3_CS_2dz_filter')
         539           0 :             if (remap_T == 'thermal_energy_over_P') then
         540             :                set_vert_remap = -10
         541             :             else
         542           0 :                set_vert_remap = -100
         543             :             end if
         544             :          case ('FV3_non_monotone_CS_2dz_filter')
         545           0 :             if (remap_T == 'thermal_energy_over_P') then
         546             :                set_vert_remap = -11
         547             :             else
         548           0 :                set_vert_remap = -110
         549             :             end if
         550             :          case default
         551           0 :             write(iulog,*)'set_vert_remap: invalid remap_alg= ',trim(remap_alg)
         552        3072 :             call endrun('set_vert_remap: invalid remap_alg')
         553             :          end select
         554             :
         555        1536 :       end function set_vert_remap
         556             :
         557             : end subroutine dyn_readnl
         558             :
         559             : !=========================================================================================
         560             :
         561        1536 : subroutine dyn_register()
         562             :
         563             :    use physics_buffer,  only: pbuf_add_field, dtype_r8
         564             :    use ppgrid,          only: pcols, pver
         565             :
         566             :    ! These fields are computed by the dycore and passed to the physics via the
         567             :    ! physics buffer.
         568             :
         569        1536 :    if (use_gw_front .or. use_gw_front_igw) then
         570             :       call pbuf_add_field("FRONTGF", "global", dtype_r8, (/pcols,pver/),       &
         571        1536 :          frontgf_idx)
         572             :       call pbuf_add_field("FRONTGA", "global", dtype_r8, (/pcols,pver/),       &
         573        1536 :          frontga_idx)
         574             :    end if
         575             :
         576        1536 : end subroutine dyn_register
         577             :
         578             : !=========================================================================================
         579             :
         580        3072 : subroutine dyn_init(dyn_in, dyn_out)
         581        1536 :    use prim_advance_mod,   only: prim_advance_init
         582             :    use dyn_grid,           only: elem, fvm
         583             :    use cam_pio_utils,      only: clean_iodesc_list
         584             :    use physconst,          only: cpair, pstd
         585             :    use air_composition,    only: thermodynamic_active_species_num, thermodynamic_active_species_idx
         586             :    use air_composition,    only: thermodynamic_active_species_idx_dycore
         587             :    use air_composition,    only: thermodynamic_active_species_liq_idx,thermodynamic_active_species_ice_idx
         588             :    use air_composition,    only: thermodynamic_active_species_liq_idx_dycore,thermodynamic_active_species_ice_idx_dycore
         589             :    use air_composition,    only: thermodynamic_active_species_liq_num, thermodynamic_active_species_ice_num
         590             :    use cam_history,        only: addfld, add_default, horiz_only, register_vector_field
         591             :    use gravity_waves_sources, only: gws_init
         592             :
         593             :    use thread_mod,         only: horz_num_threads
         594             :    use hybrid_mod,         only: get_loop_ranges, config_thread_region
         595             :    use dimensions_mod,     only: nu_scale_top
         596             :    use dimensions_mod,     only: ksponge_end, kmvis_ref, kmcnd_ref,rho_ref,km_sponge_factor
         597             :    use dimensions_mod,     only: cnst_name_gll, cnst_longname_gll
         598             :    use dimensions_mod,     only: irecons_tracer_lev, irecons_tracer, kord_tr, kord_tr_cslam
         599             :    use prim_driver_mod,    only: prim_init2
         600             :    use control_mod,        only: molecular_diff, nu_top
         601             :    use test_fvm_mapping,   only: test_mapping_addfld
         602             :    use phys_control,       only: phys_getopts
         603             :    use cam_thermo,         only: get_molecular_diff_coef_reference
         604             :    use control_mod,        only: vert_remap_uvTq_alg, vert_remap_tracer_alg
         605             :    use std_atm_profile,    only: std_atm_height
         606             :    use dyn_tests_utils,    only: vc_dycore, vc_dry_pressure, string_vc, vc_str_lgth
         607             :    use cam_budget,         only: cam_budget_em_snapshot, cam_budget_em_register, thermo_budget_history
         608             :
         609             :    ! Dummy arguments:
         610             :    type(dyn_import_t), intent(out) :: dyn_in
         611             :    type(dyn_export_t), intent(out) :: dyn_out
         612             :
         613             :    ! Local variables
         614             :    integer             :: nets, nete, ie, k, kmol_end, mfound
         615             :    real(r8), parameter :: Tinit = 300.0_r8
         616             :    real(r8)            :: press(1), ptop, tref,z(1)
         617             :
         618             :    type(hybrid_t)      :: hybrid
         619             :
         620             :    integer :: ixcldice, ixcldliq
         621             :    integer :: m_cnst, m
         622             :
         623             :    ! variables for initializing energy and axial angular momentum diagnostics
         624             :    integer, parameter                         :: num_stages = 14
         625             :    character (len = 4), dimension(num_stages) :: stage         = (/"dED","dAF","dBD","dBL","dAL","dAD","dAR","dBF","dBH","dCH","dAH","dBS","dAS","p2d"/)
         626             :    character (len = 70),dimension(num_stages) :: stage_txt = (/&
         627             :       " end of previous dynamics                           ",& !dED
         628             :       " from previous remapping or state passed to dynamics",& !dAF - state in beginning of nsplit loop
         629             :       " state after applying CAM forcing                   ",& !dBD - state after applyCAMforcing
         630             :       " before floating dynamics                           ",& !dBL
         631             :       " after floating dynamics                            ",& !dAL
         632             :       " before vertical remapping                          ",& !dAD - state before vertical remapping
         633             :       " after vertical remapping                           ",& !dAR - state at end of nsplit loop
         634             :       " state passed to parameterizations                  ",& !dBF
         635             :       " state before hypervis                              ",& !dBH
         636             :       " state after hypervis but before adding heating term",& !dCH
         637             :       " state after hypervis                               ",& !dAH
         638             :       " state before sponge layer diffusion                ",& !dBS - state before sponge del2
         639             :       " state after sponge layer diffusion                 ",& !dAS - state after sponge del2
         640             :       " phys2dyn mapping errors (requires ftype-1)         " & !p2d - for assessing phys2dyn mapping errors
         641             :       /)
         642             :
         643             :    integer :: istage
         644             :    character (len=vc_str_lgth) :: vc_str
         645             :
         646             :    logical :: history_budget        ! output tendencies and state variables for budgets
         647             :    integer :: budget_hfile_num
         648             :
         649             :    character(len=*), parameter :: sub = 'dyn_init'
         650             :
         651             :    real(r8) :: km_sponge_factor_local(nlev+1)
         652             :    !----------------------------------------------------------------------------
         653        1536 :    vc_dycore = vc_dry_pressure
         654        1536 :    if (masterproc) then
         655           2 :      call string_vc(vc_dycore,vc_str)
         656           2 :      write(iulog,*) sub//': vertical coordinate dycore   : ',trim(vc_str)
         657             :    end if
         658             :    ! Now allocate and set condenstate vars
         659        4608 :    allocate(cnst_name_gll(qsize))     ! constituent names for gll tracers
         660        4608 :    allocate(cnst_longname_gll(qsize)) ! long name of constituents for gll tracers
         661             :
         662        4608 :    allocate(kord_tr(qsize))
         663       10752 :    kord_tr(:) = vert_remap_tracer_alg
         664        1536 :    if (use_cslam) then
         665        4608 :      allocate(kord_tr_cslam(ntrac))
         666       18432 :      kord_tr_cslam(:) = vert_remap_tracer_alg
         667             :    end if
         668             :
         669             :
         670       10752 :    do m=1,qsize
         671             :      !
         672             :      ! The "_gll" index variables below are used to keep track of condensate-loading tracers
         673             :      ! since they are not necessarily indexed contiguously and not necessarily in the same
         674             :      ! order (physics is in charge of the order)
         675             :      !
         676             :      ! if running with CSLAM then the SE (gll) condensate-loading water tracers are always
         677             :      ! indexed contiguously (q,cldliq,cldice,rain,snow,graupel) - see above
         678             :      !
         679             :      ! CSLAM tracers are always indexed as in physics
         680             :      ! of no CSLAM then SE tracers are always indexed as in physics
         681             :      !
         682       10752 :      if (use_cslam) then
         683             :        !
         684             :        ! note that in this case qsize = thermodynamic_active_species_num
         685             :        !
         686        9216 :        thermodynamic_active_species_idx_dycore(m) = m
         687        9216 :        kord_tr_cslam(thermodynamic_active_species_idx(m)) = vert_remap_uvTq_alg
         688        9216 :        kord_tr(m)                                 = vert_remap_uvTq_alg
         689        9216 :        cnst_name_gll    (m)                       = cnst_name    (thermodynamic_active_species_idx(m))
         690        9216 :        cnst_longname_gll(m)                       = cnst_longname(thermodynamic_active_species_idx(m))
         691             :      else
         692             :        !
         693             :        ! if not running with CSLAM then the condensate-loading water tracers are not necessarily
         694             :        ! indexed contiguously (are indexed as in physics)
         695             :        !
         696           0 :        if (m.le.thermodynamic_active_species_num) then
         697           0 :          thermodynamic_active_species_idx_dycore(m) = thermodynamic_active_species_idx(m)
         698           0 :          kord_tr(thermodynamic_active_species_idx_dycore(m)) = vert_remap_uvTq_alg
         699             :        end if
         700           0 :        cnst_name_gll    (m)                = cnst_name    (m)
         701           0 :        cnst_longname_gll(m)                = cnst_longname(m)
         702             :      end if
         703             :    end do
         704             :
         705        4608 :    do m=1,thermodynamic_active_species_liq_num
         706        3072 :      if (use_cslam) then
         707       21504 :        do mfound=1,qsize
         708       21504 :          if (TRIM(cnst_name(thermodynamic_active_species_liq_idx(m)))==TRIM(cnst_name_gll(mfound))) then
         709        3072 :            thermodynamic_active_species_liq_idx_dycore(m) = mfound
         710             :          end if
         711             :        end do
         712             :      else
         713           0 :        thermodynamic_active_species_liq_idx_dycore(m) = thermodynamic_active_species_liq_idx(m)
         714             :      end if
         715        4608 :      if (masterproc) then
         716           4 :        write(iulog,*) sub//": m,thermodynamic_active_species_idx_liq_dycore: ",m,thermodynamic_active_species_liq_idx_dycore(m)
         717             :      end if
         718             :    end do
         719        6144 :    do m=1,thermodynamic_active_species_ice_num
         720        4608 :      if (use_cslam) then
         721       32256 :        do mfound=1,qsize
         722       32256 :          if (TRIM(cnst_name(thermodynamic_active_species_ice_idx(m)))==TRIM(cnst_name_gll(mfound))) then
         723        4608 :            thermodynamic_active_species_ice_idx_dycore(m) = mfound
         724             :          end if
         725             :        end do
         726             :      else
         727           0 :        thermodynamic_active_species_ice_idx_dycore(m) = thermodynamic_active_species_ice_idx(m)
         728             :      end if
         729        6144 :      if (masterproc) then
         730           6 :        write(iulog,*) sub//": m,thermodynamic_active_species_idx_ice_dycore: ",m,thermodynamic_active_species_ice_idx_dycore(m)
         731             :      end if
         732             :    end do
         733             :
         734             :    !
         735             :    ! Initialize the import/export objects
         736             :    !
         737        1536 :    if(iam < par%nprocs) then
         738        1536 :      dyn_in%elem  => elem
         739        1536 :      dyn_in%fvm   => fvm
         740             :
         741        1536 :      dyn_out%elem => elem
         742        1536 :      dyn_out%fvm  => fvm
         743             :    else
         744           0 :      nullify(dyn_in%elem)
         745           0 :      nullify(dyn_in%fvm)
         746           0 :      nullify(dyn_out%elem)
         747           0 :      nullify(dyn_out%fvm)
         748             :    end if
         749             :
         750        1536 :    call set_phis(dyn_in)
         751             :
         752        1536 :    if (initial_run) then
         753         768 :       call read_inidat(dyn_in)
         754         768 :       if (use_iop .and. masterproc) then
         755           0 :          call setiopupdate_init()
         756           0 :          call readiopdata( hvcoord%hyam, hvcoord%hybm, hvcoord%hyai, hvcoord%hybi, hvcoord%ps0 )
         757           0 :          call scm_setinitial(dyn_in%elem)
         758             :       end if
         759         768 :       call clean_iodesc_list()
         760             :    end if
         761             :    !
         762             :    ! initialize diffusion in dycore
         763             :    !
         764        1536 :    kmol_end = 0
         765        1536 :    if (molecular_diff>0) then
         766             :      !
         767             :      ! molecular diffusion and thermal conductivity reference values
         768             :      !
         769           0 :      if (masterproc) write(iulog,*) sub//": initialize molecular diffusion reference profiles"
         770           0 :      tref = 1000._r8     !mean value at model top for solar max
         771           0 :      km_sponge_factor = molecular_diff
         772           0 :      km_sponge_factor_local = molecular_diff
         773             :      !
         774             :      ! get rho, kmvis and kmcnd at mid-levels
         775             :      !
         776             :      call get_molecular_diff_coef_reference(tref,&
         777             :           (hvcoord%hyam(:)+hvcoord%hybm(:))*hvcoord%ps0,km_sponge_factor,&
         778           0 :           kmvis_ref,kmcnd_ref,rho_ref)
         779             :
         780           0 :      if (masterproc) then
         781           0 :         write(iulog,*) "Molecular viscosity and thermal conductivity reference profile"
         782           0 :         write(iulog,*) "k, p, z, km_sponge_factor, kmvis_ref/rho_ref, kmcnd_ref/(cp*rho_ref):"
         783             :      end if
         784           0 :      do k=1,nlev
         785             :        ! only apply molecular viscosity where viscosity is > 1000 m/s^2
         786           0 :        if (MIN(kmvis_ref(k)/rho_ref(k),kmcnd_ref(k)/(cpair*rho_ref(k)))>1000.0_r8) then
         787           0 :          if (masterproc) then
         788           0 :            press = hvcoord%hyam(k)*hvcoord%ps0+hvcoord%hybm(k)*pstd
         789           0 :            call std_atm_height(press,z)
         790           0 :            write(iulog,'(i3,5e11.4)') k,press, z,km_sponge_factor(k),kmvis_ref(k)/rho_ref(k),kmcnd_ref(k)/(cpair*rho_ref(k))
         791             :          end if
         792           0 :          kmol_end = k
         793             :        else
         794           0 :          kmvis_ref(k) = 1.0_r8
         795           0 :          kmcnd_ref(k) = 1.0_r8
         796             :        end if
         797             :      end do
         798             :    else
         799             :      ! -1.0E6 is an arbitrary unrealistic value.  But it is used in the calculation
         800             :      ! of a diagnostic quantity in global_norms_mod so can't be set to huge or nan.
         801      144384 :      kmvis_ref(:) = -1.0E6_r8
         802      144384 :      kmcnd_ref(:) = -1.0E6_r8
         803      144384 :      rho_ref(:)   = -1.0E6_r8
         804             :    end if
         805             :    !
         806      144384 :    irecons_tracer_lev(:) = irecons_tracer !use high-order CSLAM in all layers
         807             :    !
         808             :    ! compute scaling of traditional sponge layer damping (following cd_core.F90 in CAM-FV)
         809             :    !
         810        1536 :    nu_scale_top(:) = 0.0_r8
         811        1536 :    if (nu_top>0) then
         812        1536 :       ptop  = hvcoord%hyai(1)*hvcoord%ps0
         813        1536 :       if (ptop>300.0_r8) then
         814             :          !
         815             :          ! for low tops the tanh formulae below makes the sponge excessively deep
         816             :          !
         817           0 :          nu_scale_top(1) = 4.0_r8
         818           0 :          nu_scale_top(2) = 2.0_r8
         819           0 :          nu_scale_top(3) = 1.0_r8
         820           0 :          ksponge_end = 3
         821        1536 :       else if (ptop>100.0_r8) then
         822             :          !
         823             :          ! CAM6 top (~225 Pa) or CAM7 low top
         824             :          !
         825             :          ! For backwards compatibility numbers below match tanh profile
         826             :          ! used in FV
         827             :          !
         828           0 :          nu_scale_top(1) = 4.4_r8
         829           0 :          nu_scale_top(2) = 1.3_r8
         830           0 :          nu_scale_top(3) = 3.9_r8
         831           0 :          ksponge_end = 3
         832        1536 :       else if (ptop>1e-1_r8) then
         833             :          !
         834             :          ! CAM7 FMT
         835             :          !
         836        1536 :          nu_scale_top(1) = 3.0_r8
         837        1536 :          nu_scale_top(2) = 1.0_r8
         838        1536 :          nu_scale_top(3) = 0.1_r8
         839        1536 :          nu_scale_top(4) = 0.05_r8
         840        1536 :          ksponge_end = 4
         841           0 :       else if (ptop>1e-4_r8) then
         842             :          !
         843             :          ! WACCM and WACCM-x
         844             :          !
         845           0 :          nu_scale_top(1) = 5.0_r8
         846           0 :          nu_scale_top(2) = 5.0_r8
         847           0 :          nu_scale_top(3) = 5.0_r8
         848           0 :          nu_scale_top(4) = 2.0_r8
         849           0 :          nu_scale_top(5) = 1.0_r8
         850           0 :          nu_scale_top(6) = 0.1_r8
         851           0 :          ksponge_end = 6
         852             :       end if
         853             :    else
         854           0 :       ksponge_end = 0
         855             :    end if
         856        1536 :    ksponge_end = MAX(MAX(ksponge_end,1),kmol_end)
         857        1536 :    if (masterproc) then
         858           2 :      write(iulog,*) sub//": ksponge_end = ",ksponge_end
         859           2 :      write(iulog,*) sub//": sponge layer Laplacian damping"
         860           2 :      write(iulog,*) "k, p, z, nu_scale_top, nu (actual Laplacian damping coefficient)"
         861           2 :      if (nu_top>0) then
         862          12 :        do k=1,ksponge_end+1
         863          20 :          press = (hvcoord%hyam(k)+hvcoord%hybm(k))*hvcoord%ps0
         864          10 :          call std_atm_height(press,z)
         865          20 :          write(iulog,'(i3,4e11.4)') k,press,z,&
         866          32 :               nu_scale_top(k),nu_scale_top(k)*nu_top
         867             :        end do
         868             :      end if
         869             :    end if
         870             :
         871        1536 :    if (iam < par%nprocs) then
         872        1536 :       call prim_advance_init(par,elem)
         873             :       !$OMP PARALLEL NUM_THREADS(horz_num_threads), DEFAULT(SHARED), PRIVATE(hybrid,nets,nete)
         874        1536 :       hybrid = config_thread_region(par,'horizontal')
         875        1536 :       call get_loop_ranges(hybrid, ibeg=nets, iend=nete)
         876        1536 :       call prim_init2(elem, fvm, hybrid, nets, nete, TimeLevel, hvcoord)
         877             :       !$OMP END PARALLEL
         878             :
         879        1536 :       if (use_gw_front .or. use_gw_front_igw) call gws_init(elem)
         880             :    end if  ! iam < par%nprocs
         881             :
         882        3072 :    call addfld ('nu_kmvis',   (/ 'lev' /), 'A', '', 'Molecular viscosity Laplacian coefficient'            , gridname='GLL')
         883        3072 :    call addfld ('nu_kmcnd',   (/ 'lev' /), 'A', '', 'Thermal conductivity Laplacian coefficient'           , gridname='GLL')
         884        3072 :    call addfld ('nu_kmcnd_dp',(/ 'lev' /), 'A', '', 'Thermal conductivity like Laplacian coefficient on dp', gridname='GLL')
         885             :
         886             :
         887             :    ! Forcing from physics on the GLL grid
         888        3072 :    call addfld ('FU',  (/ 'lev' /), 'A', 'm/s2', 'Zonal wind forcing term on GLL grid',     gridname='GLL')
         889        3072 :    call addfld ('FV',  (/ 'lev' /), 'A', 'm/s2', 'Meridional wind forcing term on GLL grid',gridname='GLL')
         890        1536 :    call register_vector_field('FU', 'FV')
         891        3072 :    call addfld ('FT',  (/ 'lev' /), 'A', 'K/s', 'Temperature forcing term on GLL grid',gridname='GLL')
         892             :
         893             :    ! Tracer forcing on fvm (CSLAM) grid and internal CSLAM pressure fields
         894        1536 :    if (use_cslam) then
         895       18432 :       do m = 1, ntrac
         896       16896 :          call addfld (trim(cnst_name(m))//'_fvm',  (/ 'lev' /), 'I', 'kg/kg',   &
         897       50688 :             trim(cnst_longname(m)), gridname='FVM')
         898             :
         899       16896 :          call addfld ('F'//trim(cnst_name(m))//'_fvm',  (/ 'lev' /), 'I', 'kg/kg/s',   &
         900       16896 :             trim(cnst_longname(m))//' mixing ratio forcing term (q_new-q_old) on fvm grid', &
         901       69120 :             gridname='FVM')
         902             :       end do
         903             :
         904        3072 :       call addfld ('dp_fvm' ,(/ 'lev' /), 'I', 'Pa','CSLAM Pressure level thickness', gridname='FVM')
         905        1536 :       call addfld ('PSDRY_fvm',horiz_only, 'I','Pa','CSLAM dry surface pressure'    , gridname='FVM')
         906             :    end if
         907             :
         908       10752 :    do m_cnst = 1, qsize
         909           0 :      call addfld ('F'//trim(cnst_name_gll(m_cnst))//'_gll',  (/ 'lev' /), 'I', 'kg/kg/s',   &
         910       19968 :           trim(cnst_longname_gll(m_cnst))//' mixing ratio forcing term (q_new-q_old) on GLL grid', gridname='GLL')
         911             :    end do
         912             :
         913             :    ! Energy diagnostics and axial angular momentum diagnostics
         914        1536 :    call addfld ('ABS_dPSdt',  horiz_only, 'A', 'Pa/s', 'Absolute surface pressure tendency',gridname='GLL')
         915             :
         916        1536 :    if (use_cslam) then
         917             : #ifdef waccm_debug
         918             :      call addfld ('CSLAM_gamma',  (/ 'lev' /), 'A', '', 'Courant number from CSLAM',     gridname='FVM')
         919             : #endif
         920        1536 :      call addfld ('WV_PDC',   horiz_only, 'A', 'kg/m2','Total column water vapor lost in physics-dynamics coupling',gridname='FVM')
         921        1536 :      call addfld ('WL_PDC',   horiz_only, 'A', 'kg/m2','Total column cloud water lost in physics-dynamics coupling',gridname='FVM')
         922        1536 :      call addfld ('WI_PDC',   horiz_only, 'A', 'kg/m2','Total column cloud ice lost in physics-dynamics coupling'  ,gridname='FVM')
         923        1536 :      call addfld ('TT_PDC',   horiz_only, 'A', 'kg/m2','Total column test tracer lost in physics-dynamics coupling',gridname='FVM')
         924             :    else
         925           0 :      call addfld ('WV_PDC',   horiz_only, 'A', 'kg/m2','Total column water vapor lost in physics-dynamics coupling',gridname='GLL')
         926           0 :      call addfld ('WL_PDC',   horiz_only, 'A', 'kg/m2','Total column cloud water lost in physics-dynamics coupling',gridname='GLL')
         927           0 :      call addfld ('WI_PDC',   horiz_only, 'A', 'kg/m2','Total column cloud ice lost in physics-dynamics coupling'  ,gridname='GLL')
         928           0 :      call addfld ('TT_PDC',   horiz_only, 'A', 'kg/m2','Total column test tracer lost in physics-dynamics coupling',gridname='GLL')
         929             :    end if
         930             :
         931        1536 :    if (thermo_budget_history) then
         932             :       ! Register stages for budgets
         933           0 :       do istage = 1, num_stages
         934           0 :          call cam_budget_em_snapshot(TRIM(ADJUSTL(stage(istage))), 'dyn', &
         935           0 :               longname=TRIM(ADJUSTL(stage_txt(istage))))
         936             :       end do
         937             :       !
         938             :       ! Register tendency (difference) budgets
         939             :       !
         940             :       call cam_budget_em_register('dEdt_floating_dyn'   ,'dAL','dBL','dyn','dif', &
         941           0 :                       longname="dE/dt floating dynamics (dAL-dBL)"               )
         942             :       call cam_budget_em_register('dEdt_vert_remap'     ,'dAR','dAD','dyn','dif', &
         943           0 :                       longname="dE/dt vertical remapping (dAR-dAD)"              )
         944             :       call cam_budget_em_register('dEdt_phys_tot_in_dyn','dBD','dAF','dyn','dif', &
         945           0 :                       longname="dE/dt physics tendency in dynamics (dBD-dAF)"    )
         946             :       call cam_budget_em_register('dEdt_del4'          ,'dCH','dBH','dyn','dif', &
         947           0 :                       longname="dE/dt del4 (dCH-dBH)"                            )
         948             :       call cam_budget_em_register('dEdt_del4_fric_heat','dAH','dCH','dyn','dif', &
         949           0 :                       longname="dE/dt del4 frictional heating (dAH-dCH)"         )
         950             :       call cam_budget_em_register('dEdt_del4_tot'      ,'dAH','dBH','dyn','dif', &
         951           0 :                       longname="dE/dt del4 + del4 frictional heating (dAH-dBH)"  )
         952             :       call cam_budget_em_register('dEdt_del2_sponge'   ,'dAS','dBS','dyn','dif', &
         953           0 :                       longname="dE/dt del2 sponge (dAS-dBS)"                     )
         954             :       !
         955             :       ! Register derived budgets
         956             :       !
         957             :       call cam_budget_em_register('dEdt_dycore'        ,'dEdt_floating_dyn','dEdt_vert_remap'   ,'dyn','sum', &
         958           0 :                       longname="dE/dt adiabatic dynamics"                              )
         959             :       call cam_budget_em_register('dEdt_del2_del4_tot' ,'dEdt_del4_tot'    ,'dEdt_del2_sponge'  ,'dyn','sum', &
         960           0 :                       longname="dE/dt explicit diffusion total"                        )
         961             :       call cam_budget_em_register('dEdt_residual'      ,'dEdt_floating_dyn','dEdt_del2_del4_tot','dyn','dif',&
         962           0 :                       longname="dE/dt residual (dEdt_floating_dyn-dEdt_del2_del4_tot)" )
         963             :    end if
         964             :    !
         965             :    ! add dynamical core tracer tendency output
         966             :    !
         967        1536 :    if (use_cslam) then
         968       18432 :      do m = 1, pcnst
         969       33792 :        call addfld(tottnam(m),(/ 'lev' /),'A','kg/kg/s',trim(cnst_name(m))//' horz + vert',  &
         970       69120 :             gridname='FVM')
         971             :      end do
         972             :    else
         973           0 :      do m = 1, pcnst
         974           0 :        call addfld(tottnam(m),(/ 'lev' /),'A','kg/kg/s',trim(cnst_name(m))//' horz + vert',  &
         975           0 :             gridname='GLL')
         976             :      end do
         977             :    end if
         978        1536 :    call phys_getopts(history_budget_out=history_budget, history_budget_histfile_num_out=budget_hfile_num)
         979        1536 :    if ( history_budget ) then
         980           0 :       call cnst_get_ind('CLDLIQ', ixcldliq)
         981           0 :       call cnst_get_ind('CLDICE', ixcldice)
         982           0 :       call add_default(tottnam(       1), budget_hfile_num, ' ')
         983           0 :       call add_default(tottnam(ixcldliq), budget_hfile_num, ' ')
         984           0 :       call add_default(tottnam(ixcldice), budget_hfile_num, ' ')
         985             :    end if
         986             :
         987        1536 :    call test_mapping_addfld
         988        3072 : end subroutine dyn_init
         989             :
         990             : !=========================================================================================
         991             :
         992      369408 : subroutine dyn_run(dyn_state)
         993        1536 :    use air_composition,  only: thermodynamic_active_species_num, dry_air_species_num
         994             :    use air_composition,  only: thermodynamic_active_species_idx_dycore
         995             :    use prim_driver_mod,  only: prim_run_subcycle
         996             :    use dimensions_mod,   only: cnst_name_gll
         997             :    use se_dyn_time_mod,  only: tstep, nsplit, timelevel_qdp, tevolve
         998             :    use hybrid_mod,       only: config_thread_region, get_loop_ranges
         999             :    use control_mod,      only: qsplit, rsplit, ftype_conserve
        1000             :    use thread_mod,       only: horz_num_threads
        1001             :    use scamMod,          only: single_column, use_3dfrc
        1002             :    use se_single_column_mod, only: apply_SC_forcing,ie_scm
        1003             :
        1004             :    type(dyn_export_t), intent(inout) :: dyn_state
        1005             :
        1006             :    type(hybrid_t) :: hybrid
        1007             :    integer        :: tl_f
        1008             :    integer        :: n
        1009             :    integer        :: nets, nete
        1010             :    integer        :: i, ie, j, k, m, nq, m_cnst
        1011             :    integer        :: n0_qdp, nsplit_local
        1012             :    logical        :: ldiag
        1013             :
        1014             :    real(r8) :: ftmp(npsq,nlev,3)
        1015             :    real(r8) :: dtime
        1016             :    real(r8) :: rec2dt, pdel
        1017             :
        1018      369408 :    real(r8), allocatable, dimension(:,:,:) :: ps_before
        1019      369408 :    real(r8), allocatable, dimension(:,:,:) :: abs_ps_tend
        1020      738816 :    real (kind=r8)                          :: omega_cn(2,nelemd) !min and max of vertical Courant number
        1021             :    integer                                 :: nets_in,nete_in
        1022             :    !----------------------------------------------------------------------------
        1023             :
        1024             : #ifdef debug_coupling
        1025             :    return
        1026             : #endif
        1027             :
        1028      369408 :    nsplit_local = nsplit
        1029      369408 :    tevolve = 0._r8
        1030             :
        1031      369408 :    if (iam >= par%nprocs) return
        1032             :
        1033      369408 :    if (.not. use_3dfrc ) then
        1034      369408 :    ldiag = hist_fld_active('ABS_dPSdt')
        1035      369408 :    if (ldiag) then
        1036           0 :       allocate(ps_before(np,np,nelemd))
        1037           0 :       allocate(abs_ps_tend(np,np,nelemd))
        1038             :
        1039             :    end if
        1040             :
        1041             :    !$OMP PARALLEL NUM_THREADS(horz_num_threads), DEFAULT(SHARED), PRIVATE(hybrid,nets,nete,n,ie,m,i,j,k,ftmp)
        1042      369408 :    hybrid = config_thread_region(par,'horizontal')
        1043      369408 :    call get_loop_ranges(hybrid, ibeg=nets, iend=nete)
        1044             :
        1045      369408 :    dtime = get_step_size()
        1046      369408 :    rec2dt = 1._r8/dtime
        1047             :
        1048      369408 :    tl_f = TimeLevel%n0   ! timelevel which was adjusted by physics
        1049      369408 :    call TimeLevel_Qdp(TimeLevel, qsplit, n0_qdp)!get n0_qdp for diagnostics call
        1050             :
        1051             :    ! output physics forcing
        1052      369408 :    if (hist_fld_active('FU') .or. hist_fld_active('FV') .or.hist_fld_active('FT')) then
        1053           0 :       do ie = nets, nete
        1054           0 :          do k = 1, nlev
        1055           0 :             do j = 1, np
        1056           0 :                do i = 1, np
        1057           0 :                   ftmp(i+(j-1)*np,k,1) = dyn_state%elem(ie)%derived%FM(i,j,1,k)
        1058           0 :                   ftmp(i+(j-1)*np,k,2) = dyn_state%elem(ie)%derived%FM(i,j,2,k)
        1059           0 :                   ftmp(i+(j-1)*np,k,3) = dyn_state%elem(ie)%derived%FT(i,j,k)
        1060             :                end do
        1061             :             end do
        1062             :          end do
        1063             :
        1064           0 :          call outfld('FU', ftmp(:,:,1), npsq, ie)
        1065           0 :          call outfld('FV', ftmp(:,:,2), npsq, ie)
        1066           0 :          call outfld('FT', ftmp(:,:,3), npsq, ie)
        1067             :       end do
        1068             :    end if
        1069             :
        1070     2585856 :    do m = 1, qsize
        1071     2585856 :      if (hist_fld_active('F'//trim(cnst_name_gll(m))//'_gll')) then
        1072           0 :        do ie = nets, nete
        1073           0 :          call outfld('F'//trim(cnst_name_gll(m))//'_gll',&
        1074           0 :               RESHAPE(dyn_state%elem(ie)%derived%FQ(:,:,:,m), (/np*np,nlev/)), npsq, ie)
        1075             :        end do
        1076             :      end if
        1077             :    end do
        1078             :
        1079             :    ! convert elem(ie)%derived%fq to mass tendency
        1080      369408 :    if (.not.use_cslam) then
        1081           0 :      do ie = nets, nete
        1082           0 :        do m = 1, qsize
        1083           0 :          do k = 1, nlev
        1084           0 :            do j = 1, np
        1085           0 :              do i = 1, np
        1086           0 :                dyn_state%elem(ie)%derived%FQ(i,j,k,m) = dyn_state%elem(ie)%derived%FQ(i,j,k,m)* &
        1087           0 :                     rec2dt*dyn_state%elem(ie)%state%dp3d(i,j,k,tl_f)
        1088             :              end do
        1089             :            end do
        1090             :          end do
        1091             :        end do
        1092             :      end do
        1093             :    end if
        1094             :
        1095      369408 :    if (ftype_conserve>0.and..not.use_cslam) then
        1096           0 :      do ie = nets, nete
        1097           0 :        do k=1,nlev
        1098           0 :          do j=1,np
        1099           0 :            do i = 1, np
        1100           0 :              pdel     = dyn_state%elem(ie)%state%dp3d(i,j,k,tl_f)
        1101           0 :              do nq=dry_air_species_num+1,thermodynamic_active_species_num
        1102           0 :                m_cnst = thermodynamic_active_species_idx_dycore(nq)
        1103           0 :                pdel = pdel + (dyn_state%elem(ie)%state%qdp(i,j,k,m_cnst,n0_qdp)+dyn_state%elem(ie)%derived%FQ(i,j,k,m_cnst)*dtime)
        1104             :              end do
        1105           0 :              dyn_state%elem(ie)%derived%FDP(i,j,k) = pdel
        1106             :            end do
        1107             :          end do
        1108             :        end do
        1109             :      end do
        1110             :    end if
        1111             :
        1112      369408 :    if (use_cslam) then
        1113     2966808 :       do ie = nets, nete
        1114    31538208 :          do m = 1, ntrac
        1115  2688309000 :             do k = 1, nlev
        1116 10657132200 :                do j = 1, nc
        1117 34542822600 :                   do i = 1, nc
        1118           0 :                      dyn_state%fvm(ie)%fc(i,j,k,m) = dyn_state%fvm(ie)%fc(i,j,k,m)* &
        1119 31885682400 :                         rec2dt!*dyn_state%fvm(ie)%dp_fvm(i,j,k)
        1120             :                   end do
        1121             :                end do
        1122             :             end do
        1123             :          end do
        1124             :       end do
        1125             :    end if
        1126             :
        1127      369408 :    if (ldiag) then
        1128           0 :       abs_ps_tend(:,:,nets:nete) = 0.0_r8
        1129             :    endif
        1130             :
        1131     1108224 :    do n = 1, nsplit_local
        1132             :
        1133      738816 :       if (ldiag) then
        1134           0 :          do ie = nets, nete
        1135           0 :             ps_before(:,:,ie) = dyn_state%elem(ie)%state%psdry(:,:)
        1136             :          end do
        1137             :       end if
        1138             :
        1139             :       ! forward-in-time RK, with subcycling
        1140      738816 :       if (single_column) then
        1141           0 :          nets_in=ie_scm
        1142           0 :          nete_in=ie_scm
        1143             :       else
        1144      738816 :          nets_in=nets
        1145      738816 :          nete_in=nete
        1146             :       end if
        1147             :       call prim_run_subcycle(dyn_state%elem, dyn_state%fvm, hybrid, nets_in, nete_in, &
        1148      738816 :                              tstep, TimeLevel, hvcoord, n, single_column, omega_cn)
        1149             :
        1150     1108224 :       if (ldiag) then
        1151           0 :          do ie = nets, nete
        1152           0 :             abs_ps_tend(:,:,ie) = abs_ps_tend(:,:,ie) +                                &
        1153           0 :                ABS(ps_before(:,:,ie)-dyn_state%elem(ie)%state%psdry(:,:)) &
        1154           0 :                /(tstep*qsplit*rsplit)
        1155             :          end do
        1156             :       end if
        1157             :
        1158             :    end do
        1159             :
        1160      369408 :    if (ldiag) then
        1161           0 :       do ie=nets,nete
        1162           0 :          abs_ps_tend(:,:,ie)=abs_ps_tend(:,:,ie)/DBLE(nsplit)
        1163           0 :          call outfld('ABS_dPSdt',RESHAPE(abs_ps_tend(:,:,ie),(/npsq/)),npsq,ie)
        1164             :       end do
        1165             :    end if
        1166             :
        1167             :    !$OMP END PARALLEL
        1168             :
        1169             :    if (ldiag) then
        1170           0 :       deallocate(ps_before,abs_ps_tend)
        1171             :    endif
        1172             :
        1173             :    end if ! not use_3dfrc
        1174             :
        1175      369408 :    if (single_column) then
        1176           0 :       call apply_SC_forcing(dyn_state%elem,hvcoord,TimeLevel,3,.false.)
        1177             :    end if
        1178             :
        1179             :    ! output vars on CSLAM fvm grid
        1180      369408 :    call write_dyn_vars(dyn_state)
        1181             :
        1182      369408 : end subroutine dyn_run
        1183             :
        1184             : !===============================================================================
        1185             :
        1186           0 : subroutine dyn_final(DYN_STATE, RESTART_FILE)
        1187             :
        1188             :    type (elem_state_t), target     :: DYN_STATE
        1189             :    character(LEN=*)   , intent(IN) :: RESTART_FILE
        1190             :
        1191      369408 : end subroutine dyn_final
        1192             :
        1193             : !===============================================================================
        1194             :
        1195         768 : subroutine read_inidat(dyn_in)
        1196             :    use air_composition,     only: thermodynamic_active_species_num, dry_air_species_num
        1197             :    use shr_sys_mod,         only: shr_sys_flush
        1198             :    use hycoef,              only: hyai, hybi, ps0
        1199             :    use const_init,          only: cnst_init_default
        1200             :
        1201             :    use element_mod,         only: timelevels
        1202             :    use fvm_mapping,         only: dyn2fvm_mass_vars
        1203             :    use control_mod,         only: runtype
        1204             :    use prim_driver_mod,     only: prim_set_dry_mass
        1205             :    use air_composition,     only: thermodynamic_active_species_idx
        1206             :    use cam_initfiles,       only: scale_dry_air_mass
        1207             :    ! Arguments
        1208             :    type (dyn_import_t), target, intent(inout) :: dyn_in   ! dynamics import
        1209             :
        1210             :    ! Local variables
        1211             :
        1212         768 :    integer(iMap), pointer           :: ldof(:) ! Basic (2D) grid dof
        1213             :
        1214             :    type(file_desc_t), pointer       :: fh_ini, fh_topo
        1215             :
        1216         768 :    type(element_t), pointer         :: elem(:)
        1217             :
        1218         768 :    real(r8), allocatable            :: qtmp(:,:,:,:,:)    ! (np,np,nlev,nelemd,n)
        1219         768 :    real(r8), allocatable            :: dbuf2(:,:)         ! (npsq,nelemd)
        1220         768 :    real(r8), allocatable            :: dbuf3(:,:,:)       ! (npsq,nlev,nelemd)
        1221         768 :    real(r8), allocatable            :: phis_tmp(:,:)      ! (npsp,nelemd)
        1222         768 :    real(r8), allocatable            :: factor_array(:,:,:,:) ! (np,np,nlev,nelemd)
        1223         768 :    logical,  allocatable            :: pmask(:)           ! (npsq*nelemd) unique grid vals
        1224             :
        1225             :    character(len=max_hcoordname_len):: grid_name
        1226         768 :    real(r8), allocatable            :: latvals(:)
        1227         768 :    real(r8), allocatable            :: lonvals(:)
        1228         768 :    real(r8), pointer                :: latvals_deg(:)
        1229         768 :    real(r8), pointer                :: lonvals_deg(:)
        1230             :
        1231             :    integer                          :: ie, k, t
        1232             :    character(len=max_fieldname_len) :: fieldname, fieldname2
        1233             :    logical                          :: found
        1234             :    logical                          :: inic_wet           ! true if initial condition is based on
        1235             :                                                           ! wet pressure and water species
        1236             :    integer                          :: kptr, m_cnst
        1237         768 :    type(EdgeBuffer_t)               :: edge
        1238             :
        1239             :    integer                          :: rndm_seed_sz
        1240         768 :    integer, allocatable             :: rndm_seed(:)
        1241             :    integer                          :: dims(2)
        1242             :    integer                          :: pio_errtype
        1243             :    real(r8)                         :: pertval
        1244             :    integer                          :: i, j, indx, nq
        1245             :    integer                          :: dyn_cols
        1246             :    character(len=128)               :: errmsg
        1247             :    character(len=*), parameter      :: sub='READ_INIDAT'
        1248             :
        1249             :    real(r8)                         :: dp_tmp, pstmp(np,np)
        1250             :
        1251             :    ! Variables for analytic initial conditions
        1252         768 :    integer,  allocatable            :: glob_ind(:)
        1253         768 :    integer,  allocatable            :: m_ind(:)
        1254         768 :    real(r8), allocatable            :: dbuf4(:,:,:,:)
        1255             :    !----------------------------------------------------------------------------
        1256             :
        1257        1536 :    fh_ini  => initial_file_get_id()
        1258         768 :    fh_topo => topo_file_get_id()
        1259             :
        1260         768 :    if (iam < par%nprocs) then
        1261         768 :       elem => dyn_in%elem
        1262             :    else
        1263           0 :       nullify(elem)
        1264             :    end if
        1265             :
        1266        3072 :    allocate(qtmp(np,np,nlev,nelemd,pcnst))
        1267   116076816 :    qtmp = 0._r8
        1268             :
        1269             :    ! Set mask to indicate which columns are active
        1270         768 :    nullify(ldof)
        1271         768 :    call cam_grid_get_gcid(cam_grid_id(ini_grid_name), ldof)
        1272        2304 :    allocate(pmask(npsq*nelemd))
        1273       87168 :    pmask(:) = (ldof /= 0)
        1274             :
        1275             :    ! lat/lon needed in radians
        1276         768 :    latvals_deg => cam_grid_get_latvals(cam_grid_id(ini_grid_name))
        1277         768 :    lonvals_deg => cam_grid_get_lonvals(cam_grid_id(ini_grid_name))
        1278        2304 :    allocate(latvals(np*np*nelemd))
        1279        1536 :    allocate(lonvals(np*np*nelemd))
        1280       87168 :    latvals(:) = latvals_deg(:)*deg2rad
        1281       87168 :    lonvals(:) = lonvals_deg(:)*deg2rad
        1282             :
        1283             :    ! Set PIO to return error codes when reading data from IC file.
        1284         768 :    call pio_seterrorhandling(fh_ini, PIO_BCAST_ERROR, pio_errtype)
        1285             :
        1286             :    ! Get the number of columns in the global GLL grid.
        1287         768 :    call cam_grid_dimensions(ini_grid_name, dims)
        1288         768 :    dyn_cols = dims(1)
        1289             :
        1290             :    ! Set ICs.  Either from analytic expressions or read from file.
        1291             :
        1292         768 :    if (analytic_ic_active() .and. (iam < par%nprocs)) then
        1293             :
        1294             :       ! PHIS has already been set by set_phis.  Get local copy for
        1295             :       ! possible use in setting T and PS in the analytic IC code.
        1296           0 :       allocate(phis_tmp(npsq,nelemd))
        1297           0 :       do ie = 1, nelemd
        1298             :          k = 1
        1299           0 :          do j = 1, np
        1300           0 :             do i = 1, np
        1301           0 :                phis_tmp(k,ie) = elem(ie)%state%phis(i,j)
        1302           0 :                k = k + 1
        1303             :             end do
        1304             :          end do
        1305             :       end do
        1306             :
        1307           0 :       inic_wet = .false.
        1308           0 :       allocate(glob_ind(npsq * nelemd))
        1309           0 :       j = 1
        1310           0 :       do ie = 1, nelemd
        1311           0 :          do i = 1, npsq
        1312             :             ! Create a global(ish) column index
        1313           0 :             glob_ind(j) = elem(ie)%GlobalId
        1314           0 :             j = j + 1
        1315             :          end do
        1316             :       end do
        1317             :
        1318             :       ! First, initialize all the variables, then assign
        1319           0 :       allocate(dbuf4(npsq, nlev, nelemd, (qsize + 4)))
        1320           0 :       dbuf4 = 0.0_r8
        1321           0 :       allocate(m_ind(qsize))
        1322           0 :       do m_cnst = 1, qsize
        1323           0 :          m_ind(m_cnst) = m_cnst
        1324             :       end do
        1325             :
        1326             :       ! Init tracers on the GLL grid.  Note that analytic_ic_set_ic makes
        1327             :       ! use of cnst_init_default for the tracers except water vapor.
        1328             :
        1329             :       call analytic_ic_set_ic(vcoord, latvals, lonvals, glob_ind,  &
        1330             :          PS=dbuf4(:,1,:,(qsize+1)), U=dbuf4(:,:,:,(qsize+2)),      &
        1331             :          V=dbuf4(:,:,:,(qsize+3)), T=dbuf4(:,:,:,(qsize+4)),       &
        1332           0 :          Q=dbuf4(:,:,:,1:qsize), m_cnst=m_ind, mask=pmask(:),      &
        1333           0 :          PHIS_IN=PHIS_tmp)
        1334           0 :       deallocate(m_ind)
        1335           0 :       deallocate(glob_ind)
        1336           0 :       deallocate(phis_tmp)
        1337           0 :       do ie = 1, nelemd
        1338             :          indx = 1
        1339           0 :          do j = 1, np
        1340           0 :             do i = 1, np
        1341             :                ! PS
        1342           0 :                elem(ie)%state%psdry(i,j) = dbuf4(indx, 1, ie, (qsize+1))
        1343             :                ! U
        1344           0 :                elem(ie)%state%v(i,j,1,:,1) = dbuf4(indx, :, ie, (qsize+2))
        1345             :                ! V
        1346           0 :                elem(ie)%state%v(i,j,2,:,1) = dbuf4(indx, :, ie, (qsize+3))
        1347             :                ! T
        1348           0 :                elem(ie)%state%T(i,j,:,1) = dbuf4(indx, :, ie, (qsize+4))
        1349           0 :                indx = indx + 1
        1350             :             end do
        1351             :          end do
        1352             :       end do
        1353             :
        1354             :       ! Tracers to be advected on GLL grid.
        1355             :       ! Note that fvm tracers are initialized below.
        1356           0 :       do m_cnst = 1, qsize
        1357           0 :          do ie = 1, nelemd
        1358           0 :             qtmp(:,:,:,ie,m_cnst) = 0.0_r8
        1359             :             indx = 1
        1360           0 :             do j = 1, np
        1361           0 :                do i = 1, np
        1362             :                   ! Set qtmp at the unique columns only
        1363           0 :                   if (pmask(((ie - 1) * npsq) + indx)) then
        1364           0 :                      qtmp(i,j,:,ie,m_cnst) = dbuf4(indx, :, ie, m_cnst)
        1365             :                   end if
        1366           0 :                   indx = indx + 1
        1367             :                end do
        1368             :             end do
        1369             :          end do
        1370             :       end do
        1371           0 :       deallocate(dbuf4)
        1372             :
        1373             :
        1374             :    else
        1375             :
        1376             :       ! Read ICs from file.  Assume all fields in the initial file are on the GLL grid.
        1377             :
        1378        2304 :       allocate(dbuf2(npsq,nelemd))
        1379        2304 :       allocate(dbuf3(npsq,nlev,nelemd))
        1380             :
        1381             :       ! Check that columns in IC file match grid definition.
        1382         768 :       if (.not. single_column) then
        1383         768 :          call check_file_layout(fh_ini, elem, dyn_cols, 'ncdata', .true.)
        1384             :       end if
        1385             :       ! Read 2-D field
        1386             :
        1387         768 :       fieldname  = 'PS'
        1388         768 :       fieldname2 = 'PSDRY'
        1389         768 :       if (dyn_field_exists(fh_ini, trim(fieldname), required=.false.)) then
        1390         768 :          inic_wet = .true.
        1391         768 :          call read_dyn_var(trim(fieldname), fh_ini, ini_grid_hdim_name, dbuf2)
        1392           0 :       elseif (dyn_field_exists(fh_ini, trim(fieldname2), required=.false.)) then
        1393           0 :          inic_wet = .false.
        1394           0 :          call read_dyn_var(trim(fieldname2), fh_ini, ini_grid_hdim_name, dbuf2)
        1395             :       else
        1396           0 :          call endrun(trim(sub)//': PS or PSDRY must be on GLL grid')
        1397             :       end if
        1398             : #ifndef planet_mars
        1399         768 :       if (iam < par%nprocs) then
        1400       94872 :          if (minval(dbuf2, mask=reshape(pmask, (/npsq,nelemd/))) < 10000._r8) then
        1401         768 :             call endrun(trim(sub)//': Problem reading ps or psdry field -- bad values')
        1402             :          end if
        1403             :       end if
        1404             : #endif
        1405        6168 :       do ie = 1, nelemd
        1406             :          indx = 1
        1407       27768 :          do j = 1, np
        1408      113400 :             do i = 1, np
        1409       86400 :                elem(ie)%state%psdry(i,j) = dbuf2(indx,ie) ! can be either wet or dry ps
        1410      108000 :                indx = indx + 1
        1411             :             end do
        1412             :          end do
        1413             :       end do
        1414             :
        1415             :       ! Read in 3-D fields
        1416             :
        1417         768 :       if (dyn_field_exists(fh_ini, 'U')) then
        1418         768 :          call read_dyn_var('U', fh_ini, ini_grid_hdim_name, dbuf3)
        1419             :       else
        1420           0 :          call endrun(trim(sub)//': U not found')
        1421             :       end if
        1422        6168 :       do ie = 1, nelemd
        1423    64805400 :          elem(ie)%state%v = 0.0_r8
        1424             :          indx = 1
        1425       27768 :          do j = 1, np
        1426      113400 :             do i = 1, np
        1427     8121600 :                elem(ie)%state%v(i,j,1,:,1) = dbuf3(indx,:,ie)
        1428      108000 :                indx = indx + 1
        1429             :             end do
        1430             :          end do
        1431             :       end do
        1432             :
        1433         768 :       if (dyn_field_exists(fh_ini, 'V')) then
        1434         768 :          call read_dyn_var('V', fh_ini, ini_grid_hdim_name, dbuf3)
        1435             :       else
        1436           0 :          call endrun(trim(sub)//': V not found')
        1437             :       end if
        1438        6168 :       do ie = 1, nelemd
        1439             :          indx = 1
        1440       27768 :          do j = 1, np
        1441      113400 :             do i = 1, np
        1442     8121600 :                elem(ie)%state%v(i,j,2,:,1) = dbuf3(indx,:,ie)
        1443      108000 :                indx = indx + 1
        1444             :             end do
        1445             :          end do
        1446             :       end do
        1447             :
        1448         768 :       if (dyn_field_exists(fh_ini, 'T')) then
        1449         768 :          call read_dyn_var('T', fh_ini, ini_grid_hdim_name, dbuf3)
        1450             :       else
        1451           0 :          call endrun(trim(sub)//': T not found')
        1452             :       end if
        1453        6168 :       do ie=1,nelemd
        1454    31660200 :          elem(ie)%state%T = 0.0_r8
        1455             :          indx = 1
        1456       27768 :          do j = 1, np
        1457      113400 :             do i = 1, np
        1458     8121600 :                elem(ie)%state%T(i,j,:,1) = dbuf3(indx,:,ie)
        1459      108000 :                indx = indx + 1
        1460             :             end do
        1461             :          end do
        1462             :       end do
        1463             :
        1464         768 :       if (pertlim .ne. 0.0_r8) then
        1465           0 :          if (masterproc) then
        1466           0 :             write(iulog,*) trim(sub), ': Adding random perturbation bounded', &
        1467           0 :                'by +/- ', pertlim, ' to initial temperature field'
        1468             :          end if
        1469             :
        1470           0 :          call random_seed(size=rndm_seed_sz)
        1471           0 :          allocate(rndm_seed(rndm_seed_sz))
        1472             :
        1473           0 :          do ie = 1, nelemd
        1474             :             ! seed random number generator based on element ID
        1475             :             ! (possibly include a flag to allow clock-based random seeding)
        1476           0 :             rndm_seed = elem(ie)%GlobalId
        1477           0 :             call random_seed(put=rndm_seed)
        1478           0 :             do i = 1, np
        1479           0 :                do j = 1, np
        1480           0 :                   do k = 1, nlev
        1481           0 :                      call random_number(pertval)
        1482           0 :                      pertval = 2.0_r8*pertlim*(0.5_r8 - pertval)
        1483           0 :                      elem(ie)%state%T(i,j,k,1) = elem(ie)%state%T(i,j,k,1)*(1.0_r8 + pertval)
        1484             :                   end do
        1485             :                end do
        1486             :             end do
        1487             :          end do
        1488             :
        1489           0 :          deallocate(rndm_seed)
        1490             :       end if
        1491             :
        1492             :       ! Cleanup
        1493         768 :       deallocate(dbuf2)
        1494         768 :       deallocate(dbuf3)
        1495             :
        1496             :    end if ! analytic_ic_active
        1497             :
        1498             :    ! Read in or cold-initialize all the tracer fields.
        1499             :    ! Data is read in on the GLL grid.
        1500             :    ! Both GLL and FVM tracer fields are initialized based on the
        1501             :    ! dimension qsize or ntrac for GLL or FVM tracers respectively.
        1502             :    ! Data is only read in on GLL so if FVM tracers are active,
        1503             :    ! interpolation is performed.
        1504             :    !
        1505             :    ! If analytic ICs are being used, we allow constituents in an initial
        1506             :    ! file to overwrite mixing ratios set by the default constituent initialization
        1507             :    ! except for the water species.
        1508             :
        1509         768 :    if (ntrac > qsize) then
        1510         768 :       if (ntrac < pcnst) then
        1511           0 :          write(errmsg, '(a,3(i0,a))') ': ntrac (',ntrac,') > qsize (',qsize, &
        1512           0 :             ') but < pcnst (',pcnst,')'
        1513           0 :          call endrun(trim(sub)//errmsg)
        1514             :       end if
        1515           0 :    else if (qsize < pcnst) then
        1516           0 :       write(errmsg, '(a,2(i0,a))') ': qsize (',qsize,') < pcnst (',pcnst,')'
        1517           0 :       call endrun(trim(sub)//errmsg)
        1518             :    end if
        1519             :
        1520             :    ! If using analytic ICs the initial file only needs the horizonal grid
        1521             :    ! dimension checked in the case that the file contains constituent mixing
        1522             :    ! ratios.
        1523        3072 :    do m_cnst = 1, pcnst
        1524        3072 :       if (cnst_read_iv(m_cnst) .and. .not. cnst_is_a_water_species(cnst_name(m_cnst))) then
        1525         768 :          if (dyn_field_exists(fh_ini, trim(cnst_name(m_cnst)), required=.false.)) then
        1526         768 :             call check_file_layout(fh_ini, elem, dyn_cols, 'ncdata', .true.)
        1527         768 :             exit
        1528             :          end if
        1529             :       end if
        1530             :    end do
        1531             :
        1532        2304 :    allocate(dbuf3(npsq,nlev,nelemd))
        1533             :
        1534        9216 :    do m_cnst = 1, pcnst
        1535             :
        1536        8448 :       if (analytic_ic_active() .and. cnst_is_a_water_species(cnst_name(m_cnst))) cycle
        1537             :
        1538        8448 :       found = .false.
        1539        8448 :       if (cnst_read_iv(m_cnst)) then
        1540        8448 :          found = dyn_field_exists(fh_ini, trim(cnst_name(m_cnst)), required=.false.)
        1541             :       end if
        1542             :
        1543        8448 :       if (found) then
        1544        8448 :          call read_dyn_var(trim(cnst_name(m_cnst)), fh_ini, ini_grid_hdim_name, dbuf3)
        1545             :       else
        1546           0 :          call cnst_init_default(m_cnst, latvals, lonvals, dbuf3, pmask)
        1547             :       end if
        1548             :
        1549       68616 :       do ie = 1, nelemd
        1550             :          ! Copy tracers defined on GLL grid into Eulerian array
        1551             :          ! Make sure tracers have at least minimum value
        1552     5592048 :          do k=1, nlev
        1553             :             indx = 1
        1554    27680400 :             do j = 1, np
        1555   116008200 :                do i = 1, np
        1556             :                   ! Set qtmp at the unique columns only: zero non-unique columns
        1557    88387200 :                   if (pmask(((ie - 1) * npsq) + indx)) then
        1558    49719846 :                      qtmp(i,j, k, ie, m_cnst) = max(qmin(m_cnst),dbuf3(indx,k,ie))
        1559             :                   else
        1560    38667354 :                      qtmp(i,j, k, ie, m_cnst) = 0.0_r8
        1561             :                   end if
        1562   110484000 :                   indx = indx + 1
        1563             :                end do
        1564             :             end do
        1565             :          end do
        1566             :       end do
        1567             :
        1568             :    end do ! pcnst
        1569             :
        1570             :    ! Cleanup
        1571         768 :    deallocate(dbuf3)
        1572             :
        1573             :    ! Put the error handling back the way it was
        1574         768 :    call pio_seterrorhandling(fh_ini, pio_errtype)
        1575             :
        1576             :    ! Cleanup
        1577         768 :    deallocate(pmask)
        1578         768 :    deallocate(latvals)
        1579         768 :    deallocate(lonvals)
        1580             :
        1581         768 :    if (associated(ldof)) then
        1582         768 :       deallocate(ldof)
        1583             :       nullify(ldof)
        1584             :    end if
        1585             :
        1586             :    ! once we've read or initialized all the fields we do a boundary exchange to
        1587             :    ! update the redundent columns in the dynamics
        1588         768 :    if(iam < par%nprocs) then
        1589         768 :       call initEdgeBuffer(par, edge, elem, (3+pcnst)*nlev + 2 )
        1590             :    end if
        1591        6168 :    do ie = 1, nelemd
        1592        5400 :       kptr = 0
        1593        5400 :       call edgeVpack(edge, elem(ie)%state%psdry,1,kptr,ie)
        1594        5400 :       kptr = kptr + 1
        1595        5400 :       call edgeVpack(edge, elem(ie)%state%v(:,:,:,:,1),2*nlev,kptr,ie)
        1596        5400 :       kptr = kptr + (2 * nlev)
        1597        5400 :       call edgeVpack(edge, elem(ie)%state%T(:,:,:,1),nlev,kptr,ie)
        1598        5400 :       kptr = kptr + nlev
        1599   116073768 :       call edgeVpack(edge, qtmp(:,:,:,ie,:),nlev*pcnst,kptr,ie)
        1600             :    end do
        1601         768 :    if(iam < par%nprocs) then
        1602         768 :       call bndry_exchange(par,edge,location='read_inidat')
        1603             :    end if
        1604        6168 :    do ie = 1, nelemd
        1605        5400 :       kptr = 0
        1606        5400 :       call edgeVunpack(edge, elem(ie)%state%psdry,1,kptr,ie)
        1607        5400 :       kptr = kptr + 1
        1608        5400 :       call edgeVunpack(edge, elem(ie)%state%v(:,:,:,:,1),2*nlev,kptr,ie)
        1609        5400 :       kptr = kptr + (2 * nlev)
        1610        5400 :       call edgeVunpack(edge, elem(ie)%state%T(:,:,:,1),nlev,kptr,ie)
        1611        5400 :       kptr = kptr + nlev
        1612   232141368 :       call edgeVunpack(edge, qtmp(:,:,:,ie,:),nlev*pcnst,kptr,ie)
        1613             :    end do
        1614             :
        1615         768 :    if (inic_wet) then
        1616             :       !
        1617             :       ! convert to dry
        1618             :       !
        1619             :       ! (this has to be done after edge-exchange since shared points between elements are only
        1620             :       !  initialized in one element and not the other!)
        1621             :       !
        1622         768 :       if (par%masterproc) then
        1623           1 :          write(iulog,*) 'Convert specific/wet mixing ratios to dry'
        1624             :       end if
        1625             :
        1626        2304 :       allocate(factor_array(np,np,nlev,nelemd))
        1627             :       !
        1628             :       ! compute: factor_array = 1/(1-sum(q))
        1629             :       !
        1630    10552368 :       factor_array(:,:,:,:) = 1.0_r8
        1631        6168 :       do ie = 1, nelemd
        1632       38568 :          do k = dry_air_species_num+1, thermodynamic_active_species_num
        1633       32400 :             m_cnst = thermodynamic_active_species_idx(k)
        1634    63315000 :             factor_array(:,:,:,ie) = factor_array(:,:,:,ie) - qtmp(:,:,:,ie,m_cnst)
        1635             :          end do
        1636             :       end do
        1637    10552368 :       factor_array(:,:,:,:) = 1.0_r8/factor_array(:,:,:,:)
        1638             :
        1639        9216 :       do m_cnst = 1, pcnst
        1640        9216 :          if (cnst_type(m_cnst) == 'wet') then
        1641       67848 :             do ie = 1, nelemd
        1642     5592048 :                do k = 1, nlev
        1643    27680400 :                   do j = 1, np
        1644   116008200 :                      do i = 1, np
        1645             :
        1646             :                         ! convert wet mixing ratio to dry
        1647    88387200 :                         qtmp(i,j,k,ie,m_cnst) = qtmp(i,j,k,ie,m_cnst) * factor_array(i,j,k,ie)
        1648             :
        1649             :                         ! truncate negative values if they were not analytically specified
        1650   110484000 :                         if (.not. analytic_ic_active()) then
        1651    88387200 :                            qtmp(i,j,k,ie,m_cnst) = max(qmin(m_cnst), qtmp(i,j,k,ie,m_cnst))
        1652             :                         end if
        1653             :                      end do
        1654             :                   end do
        1655             :                end do
        1656             :             end do
        1657             :          end if
        1658             :       end do
        1659             :
        1660             :       ! initialize dp3d and qdp
        1661             :       !
        1662             :       ! compute: factor_array = 1/(1+sum(q))
        1663             :
        1664    10552368 :       factor_array(:,:,:,:) = 1.0_r8
        1665        6168 :       do ie = 1, nelemd
        1666       38568 :          do k = dry_air_species_num+1, thermodynamic_active_species_num
        1667       32400 :             m_cnst = thermodynamic_active_species_idx(k)
        1668    63315000 :             factor_array(:,:,:,ie) = factor_array(:,:,:,ie) + qtmp(:,:,:,ie,m_cnst)
        1669             :          end do
        1670             :       end do
        1671    10552368 :       factor_array(:,:,:,:) = 1.0_r8/factor_array(:,:,:,:)
        1672        6168 :       do ie = 1, nelemd
        1673             :          ! pstmp is the wet ps
        1674      113400 :          pstmp = elem(ie)%state%psdry(:,:)
        1675             :          ! start accumulating the dry air pressure differences across each layer
        1676      113400 :          elem(ie)%state%psdry(:,:) = hyai(1)*ps0
        1677      508368 :          do k=1,nlev
        1678     2516400 :             do j = 1,np
        1679    10546200 :                do i = 1,np
        1680    16070400 :                   dp_tmp = ((hyai(k+1) - hyai(k))*ps0) + &
        1681    24105600 :                            ((hybi(k+1) - hybi(k))*pstmp(i,j))
        1682     8035200 :                   if (.not. analytic_ic_active()) then
        1683             :
        1684             :                      ! if analytic_ic then the surface pressure is already dry
        1685             :                      ! (note that it is not correct to convert to moist pressure
        1686             :                      ! in analytic_ic and not have the #ifndef statement here
        1687             :                      ! since the dry levels are in a different location than
        1688             :                      ! what is obtained from algorithm below)
        1689             :
        1690             :                      ! convert dp_tmp to dry
        1691     8035200 :                      dp_tmp = dp_tmp*factor_array(i,j,k,ie)
        1692             :                   end if
        1693             :
        1694    32140800 :                   elem(ie)%state%dp3d(i,j,k,:) = dp_tmp
        1695             :
        1696             :                   ! compute dry surface pressure; note that at this point
        1697             :                   !
        1698             :                   ! dp3d .NE. (hyai(k+1) - hyai(k))*ps0 + (hybi(k+1) - hybi(k))*ps(i,j)
        1699             :
        1700    10044000 :                   elem(ie)%state%psdry(i,j) = elem(ie)%state%psdry(i,j)+elem(ie)%state%dp3d(i,j,k,1)
        1701             :                end do
        1702             :             end do
        1703             :          end do
        1704             :       end do
        1705             :
        1706         768 :       deallocate(factor_array)
        1707             :
        1708             :    else
        1709             :
        1710             :       ! initial condition is based on dry surface pressure and constituents
        1711             :       !
        1712             :       ! we only need to initialize state%dp3d
        1713             :
        1714           0 :       do ie = 1, nelemd
        1715           0 :          do k = 1, nlev
        1716           0 :             do j = 1, np
        1717           0 :                do i = 1, np
        1718           0 :                   elem(ie)%state%dp3d(i,j,k,:) = (hyai(k+1) - hyai(k))*ps0 + &
        1719           0 :                                                  (hybi(k+1) - hybi(k))*elem(ie)%state%psdry(i,j)
        1720             :                end do
        1721             :             end do
        1722             :          end do
        1723             :       end do
        1724             :    end if
        1725             :
        1726             :    ! If scale_dry_air_mass > 0.0 then scale dry air mass to scale_dry_air_mass global average dry pressure
        1727         768 :    if (runtype == 0) then
        1728         768 :      if (scale_dry_air_mass > 0.0_r8) then
        1729         768 :        if (iam < par%nprocs) then
        1730         768 :          call prim_set_dry_mass(elem, hvcoord, scale_dry_air_mass, qtmp)
        1731             :        end if
        1732             :      end if
        1733             :    end if
        1734             :
        1735             :    ! store Q values:
        1736             :    !
        1737             :    ! if CSLAM is NOT active then state%Qdp for all constituents
        1738             :    ! if CSLAM active then we only advect water vapor and condensate
        1739             :    ! loading tracers in state%qdp
        1740             :
        1741         768 :    if (use_cslam) then
        1742        6168 :       do ie = 1, nelemd
        1743       38568 :          do nq = 1, thermodynamic_active_species_num
        1744       32400 :             m_cnst = thermodynamic_active_species_idx(nq)
        1745     3051000 :             do k = 1, nlev
        1746    15098400 :                do j = 1, np
        1747    63277200 :                   do i = 1, np
        1748           0 :                      elem(ie)%state%Qdp(i,j,k,nq,:) = &
        1749   108475200 :                                  elem(ie)%state%dp3d(i,j,k,1)*qtmp(i,j,k,ie,m_cnst)
        1750             :                   end do
        1751             :                end do
        1752             :             end do
        1753             :          end do
        1754             :       end do
        1755             :    else
        1756           0 :       do ie = 1, nelemd
        1757           0 :          do m_cnst = 1, qsize
        1758           0 :             do k = 1, nlev
        1759           0 :                do j = 1, np
        1760           0 :                   do i = 1, np
        1761           0 :                      elem(ie)%state%Qdp(i,j,k,m_cnst,:)=&
        1762           0 :                         elem(ie)%state%dp3d(i,j,k,1)*qtmp(i,j,k,ie,m_cnst)
        1763             :                   end do
        1764             :                end do
        1765             :             end do
        1766             :          end do
        1767             :       end do
        1768             :    end if
        1769             :
        1770             :    ! interpolate fvm tracers and fvm pressure variables
        1771             :
        1772         768 :    if (use_cslam) then
        1773         768 :       if (par%masterproc) then
        1774           1 :          write(iulog,*) 'Initializing dp_fvm from spectral element dp'
        1775             :       end if
        1776             :
        1777        6168 :       do ie = 1, nelemd
        1778             :
        1779             :          ! note that the area over fvm cells as computed from subcell_integration is up to 1.0E-6
        1780             :          ! different than the areas (exact) computed by CSLAM
        1781             :          !
        1782             :          ! Map the constituents which are also to be transported by dycore
        1783           0 :          call dyn2fvm_mass_vars(elem(ie)%state%dp3d(:,:,:,1),elem(ie)%state%psdry(:,:),&
        1784        5400 :             qtmp(:,:,:,ie,1:ntrac),&
        1785           0 :             dyn_in%fvm(ie)%dp_fvm(1:nc,1:nc,:),dyn_in%fvm(ie)%psC(1:nc,1:nc),&
        1786           0 :             dyn_in%fvm(ie)%c(1:nc,1:nc,:,1:ntrac),&
        1787   272884368 :             ntrac,elem(ie)%metdet,dyn_in%fvm(ie)%inv_se_area_sphere(1:nc,1:nc))
        1788             :       end do
        1789             :
        1790         768 :       if(par%masterproc) then
        1791           1 :          write(iulog,*) 'FVM tracers, FVM pressure variables and se_area_sphere initialized.'
        1792             :       end if
        1793             :
        1794             :    end if    ! (use_cslam)
        1795             :
        1796             :    ! Cleanup
        1797         768 :    deallocate(qtmp)
        1798             :
        1799        6168 :    do ie = 1, nelemd
        1800       16968 :       do t = 2, timelevels
        1801    43200000 :          elem(ie)%state%v(:,:,:,:,t) = elem(ie)%state%v(:,:,:,:,1)
        1802    21108600 :          elem(ie)%state%T(:,:,:,t)   = elem(ie)%state%T(:,:,:,1)
        1803             :       end do
        1804             :    end do
        1805             :
        1806         768 :    if(iam < par%nprocs) then
        1807         768 :       call FreeEdgeBuffer(edge)
        1808             :    end if
        1809             :
        1810        1536 : end subroutine read_inidat
        1811             :
        1812             :
        1813             : !========================================================================================
        1814             :
        1815        1536 : subroutine set_phis(dyn_in)
        1816             :
        1817             :    ! Set PHIS according to the following rules.
        1818             :    !
        1819             :    ! 1) If a topo file is specified use it.  This option has highest precedence.
        1820             :    ! 2) If not using topo file, but analytic_ic option is on, use analytic phis.
        1821             :    ! 3) Set phis = 0.0.
        1822             :    !
        1823             :    ! If using the physics grid then the topo file will be on that grid since its
        1824             :    ! contents are primarily for the physics parameterizations, and the values of
        1825             :    ! PHIS should be consistent with the values of sub-grid variability (e.g., SGH)
        1826             :    ! which are computed on the physics grid.  In this case phis on the physics grid
        1827             :    ! will be interpolated to the GLL grid.
        1828             :
        1829             :
        1830             :    ! Arguments
        1831             :    type (dyn_import_t), target, intent(inout) :: dyn_in   ! dynamics import
        1832             :
        1833             :    ! local variables
        1834             :    type(file_desc_t), pointer       :: fh_topo
        1835             :
        1836        1536 :    type(element_t), pointer         :: elem(:)
        1837             :
        1838        1536 :    real(r8), allocatable            :: phis_tmp(:,:)      ! (npsp,nelemd)
        1839        1536 :    real(r8), allocatable            :: phis_phys_tmp(:,:) ! (fv_nphys**2,nelemd)
        1840             :
        1841             :    integer                          :: i, ie, indx, j, kptr
        1842             :    integer                          :: ierr, pio_errtype
        1843             :
        1844             :    character(len=max_fieldname_len) :: fieldname
        1845             :    character(len=max_fieldname_len) :: fieldname_gll
        1846             :    character(len=max_hcoordname_len):: grid_name
        1847             :    integer                          :: dims(2)
        1848             :    integer                          :: dyn_cols
        1849             :    integer                          :: ncol_did
        1850             :    integer                          :: ncol_size
        1851             :
        1852        1536 :    integer(iMap), pointer           :: ldof(:)            ! Basic (2D) grid dof
        1853        1536 :    logical,  allocatable            :: pmask(:)           ! (npsq*nelemd) unique columns
        1854             :
        1855             :    ! Variables for analytic initial conditions
        1856        1536 :    integer,  allocatable            :: glob_ind(:)
        1857             :    logical,  allocatable            :: pmask_phys(:)
        1858        1536 :    real(r8), pointer                :: latvals_deg(:)
        1859        1536 :    real(r8), pointer                :: lonvals_deg(:)
        1860        1536 :    real(r8), allocatable            :: latvals(:)
        1861        1536 :    real(r8), allocatable            :: lonvals(:)
        1862             :    real(r8), allocatable            :: latvals_phys(:)
        1863             :    real(r8), allocatable            :: lonvals_phys(:)
        1864             :
        1865             :    character(len=*), parameter      :: sub='set_phis'
        1866             :    !----------------------------------------------------------------------------
        1867             :
        1868        3072 :    fh_topo => topo_file_get_id()
        1869             :
        1870        1536 :    if (iam < par%nprocs) then
        1871        1536 :       elem => dyn_in%elem
        1872             :    else
        1873           0 :       nullify(elem)
        1874             :    end if
        1875             :
        1876        4608 :    allocate(phis_tmp(npsq,nelemd))
        1877      185136 :    phis_tmp = 0.0_r8
        1878             :
        1879        1536 :    if (use_cslam) then
        1880        6144 :       allocate(phis_phys_tmp(fv_nphys**2,nelemd))
        1881      109536 :       phis_phys_tmp = 0.0_r8
        1882       12336 :       do ie=1,nelemd
        1883    53245536 :         elem(ie)%sub_elem_mass_flux=0.0_r8
        1884             : #ifdef waccm_debug
        1885             :         dyn_in%fvm(ie)%CSLAM_gamma = 0.0_r8
        1886             : #endif
        1887             :       end do
        1888             :    end if
        1889             :
        1890             :    ! Set mask to indicate which columns are active in GLL grid.
        1891        1536 :    nullify(ldof)
        1892        1536 :    call cam_grid_get_gcid(cam_grid_id('GLL'), ldof)
        1893        4608 :    allocate(pmask(npsq*nelemd))
        1894      174336 :    pmask(:) = (ldof /= 0)
        1895        1536 :    deallocate(ldof)
        1896             :
        1897        1536 :    if (associated(fh_topo)) then
        1898             :
        1899             :       ! Set PIO to return error flags.
        1900        1536 :       call pio_seterrorhandling(fh_topo, PIO_BCAST_ERROR, pio_errtype)
        1901             :
        1902             :       ! Set name of grid object which will be used to read data from file
        1903             :       ! into internal data structure via PIO.
        1904        1536 :       if (single_column) then
        1905           0 :          grid_name = 'SCM'
        1906             :       else
        1907        1536 :          if (fv_nphys == 0) then
        1908           0 :             grid_name = 'GLL'
        1909             :          else
        1910        1536 :             grid_name = 'physgrid_d'
        1911             :          end if
        1912             :       end if
        1913             :
        1914             :       ! Get number of global columns from the grid object and check that
        1915             :       ! it matches the file data.
        1916        1536 :       call cam_grid_dimensions(grid_name, dims)
        1917        1536 :       dyn_cols = dims(1)
        1918             :
        1919             :       ! The dimension of the unstructured grid in the TOPO file is 'ncol'.
        1920        1536 :       ierr = pio_inq_dimid(fh_topo, 'ncol', ncol_did)
        1921        1536 :       if (ierr /= PIO_NOERR) then
        1922           0 :          call endrun(sub//': dimension ncol not found in bnd_topo file')
        1923             :       end if
        1924        1536 :       ierr = pio_inq_dimlen(fh_topo, ncol_did, ncol_size)
        1925        1536 :       if (ncol_size /= dyn_cols .and. .not. single_column) then
        1926           0 :          if (masterproc) then
        1927           0 :             write(iulog,*) sub//': ncol_size=', ncol_size, ' : dyn_cols=', dyn_cols
        1928             :          end if
        1929           0 :          call endrun(sub//': ncol size in bnd_topo file does not match grid definition')
        1930             :       end if
        1931             :
        1932        1536 :       fieldname = 'PHIS'
        1933        1536 :       fieldname_gll = 'PHIS_gll'
        1934        1536 :       if (use_cslam.and.dyn_field_exists(fh_topo, trim(fieldname_gll),required=.false.)) then
        1935             :          !
        1936             :          ! If physgrid it is recommended to read in PHIS on the GLL grid and then
        1937             :          ! map to the physgrid in d_p_coupling
        1938             :          !
        1939             :          ! This requires a topo file with PHIS_gll on it ...
        1940             :          !
        1941        1536 :          if (masterproc) then
        1942           2 :             write(iulog, *) "Reading in PHIS on GLL grid (mapped to physgrid in d_p_coupling)"
        1943             :          end if
        1944        1536 :          call read_dyn_var(fieldname_gll, fh_topo, 'ncol_gll', phis_tmp)
        1945           0 :       else if (dyn_field_exists(fh_topo, trim(fieldname))) then
        1946           0 :          if (.not.use_cslam) then
        1947           0 :             if (masterproc) then
        1948           0 :                write(iulog, *) "Reading in PHIS"
        1949             :             end if
        1950           0 :             call read_dyn_var(fieldname, fh_topo, 'ncol', phis_tmp)
        1951             :          else
        1952             :             !
        1953             :             ! For backwards compatibility we allow reading in PHIS on the physgrid
        1954             :             ! which is then mapped to the GLL grid and back to the physgrid in d_p_coupling
        1955             :             ! (the latter is to avoid noise in derived quantities such as PSL)
        1956             :             !
        1957           0 :             if (masterproc) then
        1958           0 :                write(iulog, *) "Reading in PHIS on physgrid"
        1959           0 :                write(iulog, *) "Recommended to read in PHIS on GLL grid"
        1960             :             end if
        1961           0 :             call read_phys_field_2d(fieldname, fh_topo, 'ncol', phis_phys_tmp)
        1962             :             call map_phis_from_physgrid_to_gll(dyn_in%fvm, elem, phis_phys_tmp, &
        1963           0 :                  phis_tmp, pmask)
        1964           0 :             deallocate(phis_phys_tmp)
        1965             :          end if
        1966             :       else
        1967           0 :          call endrun(sub//': Could not find PHIS field on input datafile')
        1968             :       end if
        1969             :
        1970             :       ! Put the error handling back the way it was
        1971        1536 :       call pio_seterrorhandling(fh_topo, pio_errtype)
        1972             :
        1973           0 :    else if (analytic_ic_active() .and. (iam < par%nprocs)) then
        1974             :
        1975             :       ! lat/lon needed in radians
        1976           0 :       latvals_deg => cam_grid_get_latvals(cam_grid_id('GLL'))
        1977           0 :       lonvals_deg => cam_grid_get_lonvals(cam_grid_id('GLL'))
        1978           0 :       allocate(latvals(np*np*nelemd))
        1979           0 :       allocate(lonvals(np*np*nelemd))
        1980           0 :       latvals(:) = latvals_deg(:)*deg2rad
        1981           0 :       lonvals(:) = lonvals_deg(:)*deg2rad
        1982             :
        1983           0 :       allocate(glob_ind(npsq*nelemd))
        1984           0 :       j = 1
        1985           0 :       do ie = 1, nelemd
        1986           0 :          do i = 1, npsq
        1987             :             ! Create a global(ish) column index
        1988           0 :             glob_ind(j) = elem(ie)%GlobalId
        1989           0 :             j = j + 1
        1990             :          end do
        1991             :       end do
        1992             :       call analytic_ic_set_ic(vcoord, latvals, lonvals, glob_ind, &
        1993           0 :                               PHIS_OUT=phis_tmp, mask=pmask(:))
        1994           0 :       deallocate(glob_ind)
        1995             :
        1996             :    end if
        1997             :
        1998        1536 :    deallocate(pmask)
        1999             :
        2000             :    ! Set PHIS in element objects
        2001       12336 :    do ie = 1, nelemd
        2002      226800 :       elem(ie)%state%phis = 0.0_r8
        2003             :       indx = 1
        2004       55536 :       do j = 1, np
        2005      226800 :          do i = 1, np
        2006      172800 :             elem(ie)%state%phis(i,j) = phis_tmp(indx, ie)
        2007      216000 :             indx = indx + 1
        2008             :          end do
        2009             :       end do
        2010             :    end do
        2011        1536 :    deallocate(phis_tmp)
        2012             :
        2013             :    ! boundary exchange to update the redundent columns in the element objects
        2014       12336 :    do ie = 1, nelemd
        2015       10800 :       kptr = 0
        2016       12336 :       call edgeVpack(edgebuf, elem(ie)%state%phis, 1, kptr, ie)
        2017             :    end do
        2018        1536 :    if(iam < par%nprocs) then
        2019        1536 :       call bndry_exchange(par, edgebuf, location=sub)
        2020             :    end if
        2021       12336 :    do ie = 1, nelemd
        2022       10800 :       kptr = 0
        2023       12336 :       call edgeVunpack(edgebuf, elem(ie)%state%phis,1,kptr,ie)
        2024             :    end do
        2025             :
        2026        3840 : end subroutine set_phis
        2027             :
        2028             : !========================================================================================
        2029             :
        2030        1536 : subroutine check_file_layout(file, elem, dyn_cols, file_desc, dyn_ok)
        2031             :
        2032             :    ! This routine is only called when data will be read from the initial file.  It is not
        2033             :    ! called when the initial file is only supplying vertical coordinate info.
        2034             :
        2035             :    type(file_desc_t), pointer       :: file
        2036             :    type(element_t),   pointer       :: elem(:)
        2037             :    integer,           intent(in)    :: dyn_cols
        2038             :    character(len=*),  intent(in)    :: file_desc
        2039             :    logical,           intent(in)    :: dyn_ok ! .true. iff ncol_d is okay
        2040             :
        2041             :    integer                          :: ncol_did, ncol_size
        2042             :    integer                          :: ierr
        2043             :    integer                          :: ie, i, j
        2044             :    integer                          :: indx
        2045        3072 :    real(r8)                         :: dbuf2(npsq, nelemd)
        2046             :    logical                          :: found
        2047             :    character(len=max_fieldname_len) :: coordname
        2048             :
        2049             :    character(len=*), parameter      :: sub = 'check_file_layout'
        2050             :    !----------------------------------------------------------------------------
        2051             :
        2052             :    ! Check that number of columns in IC file matches grid definition.
        2053        1536 :    if (trim(ini_grid_hdim_name) == 'none') then
        2054             :       call endrun(sub//': ERROR: no horizontal dimension in initial data file. &
        2055           0 :          &Cannot read data from file')
        2056             :    end if
        2057             :
        2058        1536 :    ierr = pio_inq_dimid(file, trim(ini_grid_hdim_name), ncol_did)
        2059        1536 :    if (ierr /= PIO_NOERR) then
        2060             :       call endrun(sub//': ERROR: '//trim(ini_grid_hdim_name)//' dimension not found in ' &
        2061           0 :          //trim(file_desc)//' file')
        2062             :    end if
        2063             :
        2064        1536 :    ierr = pio_inq_dimlen(file, ncol_did, ncol_size)
        2065        1536 :    if (ncol_size /= dyn_cols) then
        2066           0 :       if (masterproc) then
        2067           0 :          write(iulog, '(a,2(a,i0))') trim(sub), ': ncol_size=', ncol_size, &
        2068           0 :              ' : dyn_cols=', dyn_cols
        2069             :       end if
        2070           0 :       call endrun(sub//': ERROR: dimension '//trim(ini_grid_hdim_name)//' size not same as in ncdata file')
        2071             :    end if
        2072             :
        2073             :    ! Set coordinate name associated with ini_grid_hdim_name.
        2074        1536 :    if (trim(ini_grid_hdim_name) == 'ncol') then
        2075           0 :       coordname = 'lat'
        2076             :    else
        2077        1536 :       coordname = 'lat_d'
        2078             :    end if
        2079             :
        2080             :    !! Check to make sure file is in correct order
        2081        1536 :    call read_dyn_var(coordname, file, ini_grid_hdim_name, dbuf2)
        2082        1536 :    found = .true.
        2083       12336 :    do ie = 1, nelemd
        2084       10800 :       indx = 1
        2085       55536 :       do j = 1, np
        2086      226800 :          do i = 1, np
        2087      172800 :             if (abs(dbuf2(indx,ie)) > 1.e-12_r8) then
        2088       96484 :                if (abs((elem(ie)%spherep(i,j)%lat*rad2deg - dbuf2(indx,ie)) / &
        2089             :                     dbuf2(indx,ie)) > 1.0e-10_r8) then
        2090             :                   write(iulog, '(2a,4(i0,a),f11.5,a,f11.5)')                  &
        2091           0 :                        "ncdata file latitudes not in correct column order",   &
        2092           0 :                        ' on task ', iam, ': elem(', ie, ')%spherep(', i,      &
        2093           0 :                        ', ', j, ')%lat = ', elem(ie)%spherep(i,j)%lat,        &
        2094           0 :                        ' /= ', dbuf2(indx, ie)*deg2rad
        2095           0 :                   call shr_sys_flush(iulog)
        2096           0 :                   found = .false.
        2097             :                end if
        2098             :             end if
        2099      216000 :             indx = indx + 1
        2100             :          end do
        2101             :       end do
        2102             :    end do
        2103        1536 :    if (.not. found) then
        2104           0 :       call endrun("ncdata file latitudes not in correct column order")
        2105             :    end if
        2106             :
        2107        1536 :    if (trim(ini_grid_hdim_name) == 'ncol') then
        2108           0 :       coordname = 'lon'
        2109             :    else
        2110        1536 :       coordname = 'lon_d'
        2111             :    end if
        2112             :
        2113        1536 :    call read_dyn_var(coordname, file, ini_grid_hdim_name, dbuf2)
        2114       12336 :    do ie = 1, nelemd
        2115       10800 :       indx = 1
        2116       55536 :       do j = 1, np
        2117      226800 :          do i = 1, np
        2118      172800 :             if (abs(dbuf2(indx,ie)) > 1.e-12_r8) then
        2119       97200 :                if (abs((elem(ie)%spherep(i,j)%lon*rad2deg - dbuf2(indx,ie)) / &
        2120             :                     dbuf2(indx,ie)) > 1.0e-10_r8) then
        2121             :                   write(iulog, '(2a,4(i0,a),f11.5,a,f11.5)')                  &
        2122           0 :                        "ncdata file longitudes not in correct column order",  &
        2123           0 :                        ' on task ', iam, ': elem(', ie, ')%spherep(', i,      &
        2124           0 :                        ', ', j, ')%lon = ', elem(ie)%spherep(i,j)%lon,        &
        2125           0 :                        ' /= ', dbuf2(indx, ie)*deg2rad
        2126           0 :                   call shr_sys_flush(iulog)
        2127           0 :                   found = .false.
        2128             :                end if
        2129             :             end if
        2130      216000 :             indx = indx + 1
        2131             :          end do
        2132             :       end do
        2133             :    end do
        2134        1536 :    if (.not. found) then
        2135           0 :       call endrun("ncdata file longitudes not in correct column order")
        2136             :    end if
        2137        1536 : end subroutine check_file_layout
        2138             :
        2139             : !========================================================================================
        2140             :
        2141       13824 : logical function dyn_field_exists(fh, fieldname, required)
        2142             :
        2143             :    use pio,            only: var_desc_t, PIO_inq_varid
        2144             :    use pio,            only: PIO_NOERR
        2145             :
        2146             :    type(file_desc_t), intent(in) :: fh
        2147             :    character(len=*),  intent(in) :: fieldname
        2148             :    logical, optional, intent(in) :: required
        2149             :
        2150             :    ! Local variables
        2151             :    logical                  :: found
        2152             :    logical                  :: field_required
        2153             :    integer                  :: ret
        2154             :    type(var_desc_t)         :: varid
        2155             :    character(len=128)       :: errormsg
        2156             :    !--------------------------------------------------------------------------
        2157             :
        2158       13824 :    if (present(required)) then
        2159       11520 :       field_required = required
        2160             :    else
        2161             :       field_required = .true.
        2162             :    end if
        2163             :
        2164       13824 :    ret = PIO_inq_varid(fh, trim(fieldname), varid)
        2165       13824 :    found = (ret == PIO_NOERR)
        2166       13824 :    if (.not. found) then
        2167           0 :       if (field_required) then
        2168           0 :          write(errormsg, *) trim(fieldname),' was not present in the input file.'
        2169           0 :          call endrun('DYN_FIELD_EXISTS: '//errormsg)
        2170             :       end if
        2171             :    end if
        2172             :
        2173       13824 :    dyn_field_exists = found
        2174             :
        2175       13824 : end function dyn_field_exists
        2176             :
        2177             : !========================================================================================
        2178             :
        2179        5376 : subroutine read_dyn_field_2d(fieldname, fh, dimname, buffer)
        2180             :
        2181             :    ! Dummy arguments
        2182             :    character(len=*),  intent(in)    :: fieldname
        2183             :    type(file_desc_t), intent(inout) :: fh
        2184             :    character(len=*),  intent(in)    :: dimname
        2185             :    real(r8),          intent(inout) :: buffer(:, :)
        2186             :
        2187             :    ! Local variables
        2188             :    logical                  :: found
        2189             :    real(r8)                 :: fillvalue
        2190             :    !----------------------------------------------------------------------------
        2191             :
        2192      647976 :    buffer = 0.0_r8
        2193             :    call infld(trim(fieldname), fh, dimname, 1, npsq, 1, nelemd, buffer,    &
        2194        5376 :         found, gridname=ini_grid_name, fillvalue=fillvalue)
        2195        5376 :    if(.not. found) then
        2196           0 :       call endrun('READ_DYN_FIELD_2D: Could not find '//trim(fieldname)//' field on input datafile')
        2197             :    end if
        2198             :
        2199             :    ! This code allows use of compiler option to set uninitialized values
        2200             :    ! to NaN.  In that case infld can return NaNs where the element GLL points
        2201             :    ! are not "unique columns"
        2202             :    ! Set NaNs or fillvalue points to zero
        2203      647976 :    where (isnan(buffer) .or. (buffer==fillvalue)) buffer = 0.0_r8
        2204             :
        2205        5376 : end subroutine read_dyn_field_2d
        2206             :
        2207             : !========================================================================================
        2208             :
        2209       10752 : subroutine read_dyn_field_3d(fieldname, fh, dimname, buffer)
        2210             :
        2211             :    ! Dummy arguments
        2212             :    character(len=*),  intent(in)    :: fieldname
        2213             :    type(file_desc_t), intent(inout) :: fh
        2214             :    character(len=*),  intent(in)    :: dimname
        2215             :    real(r8),          intent(inout) :: buffer(:,:,:)
        2216             :
        2217             :    ! Local variables
        2218             :    logical                  :: found
        2219             :    real(r8)                 :: fillvalue
        2220             :    !----------------------------------------------------------------------------
        2221             :
        2222   119609952 :    buffer = 0.0_r8
        2223             :    call infld(trim(fieldname), fh, dimname, 'lev',  1, npsq, 1, nlev,         &
        2224       10752 :         1, nelemd, buffer, found, gridname=ini_grid_name, fillvalue=fillvalue)
        2225       10752 :    if(.not. found) then
        2226           0 :       call endrun('READ_DYN_FIELD_3D: Could not find '//trim(fieldname)//' field on input datafile')
        2227             :    end if
        2228             :
        2229             :    ! This code allows use of compiler option to set uninitialized values
        2230             :    ! to NaN.  In that case infld can return NaNs where the element GLL points
        2231             :    ! are not "unique columns"
        2232             :    ! Set NaNs or fillvalue points to zero
        2233   119609952 :    where (isnan(buffer) .or. (buffer == fillvalue)) buffer = 0.0_r8
        2234             :
        2235       10752 : end subroutine read_dyn_field_3d
        2236             :
        2237             : !========================================================================================
        2238             :
        2239           0 : subroutine read_phys_field_2d(fieldname, fh, dimname, buffer)
        2240             :
        2241             :    ! Dummy arguments
        2242             :    character(len=*),  intent(in)    :: fieldname
        2243             :    type(file_desc_t), intent(inout) :: fh
        2244             :    character(len=*),  intent(in)    :: dimname
        2245             :    real(r8),          intent(inout) :: buffer(:, :)
        2246             :
        2247             :    ! Local variables
        2248             :    logical                  :: found
        2249             :    !----------------------------------------------------------------------------
        2250             :
        2251             :    call infld(trim(fieldname), fh, dimname, 1, fv_nphys**2, 1, nelemd, buffer,    &
        2252           0 :       found, gridname='physgrid_d')
        2253           0 :    if(.not. found) then
        2254           0 :       call endrun('READ_PHYS_FIELD_2D: Could not find '//trim(fieldname)//' field on input datafile')
        2255             :    end if
        2256             :
        2257           0 : end subroutine read_phys_field_2d
        2258             :
        2259             : !========================================================================================
        2260             :
        2261           0 : subroutine map_phis_from_physgrid_to_gll(fvm,elem,phis_phys_tmp,phis_tmp,pmask)
        2262             :
        2263             :    use hybrid_mod,         only: get_loop_ranges, config_thread_region
        2264             :    use dimensions_mod,     only: nhc_phys
        2265             :    use fvm_mapping,        only: phys2dyn
        2266             :    use thread_mod,         only: horz_num_threads
        2267             :
        2268             :    type(element_t),   intent(inout) :: elem(:)
        2269             :    type (fvm_struct), intent(in)    :: fvm(:)
        2270             :    real(r8)         , intent(in)    :: phis_phys_tmp(fv_nphys**2,nelemd) !physgrid phis
        2271             :    real(r8)         , intent(inout) :: phis_tmp(npsq,nelemd)             !gll phis
        2272             :    logical          , intent(in)    :: pmask(npsq*nelemd)
        2273             :
        2274             :    type(hybrid_t)                   :: hybrid
        2275             :    integer                          :: nets, nete, ie,i,j,indx
        2276           0 :    real(r8),            allocatable :: fld_phys(:,:,:,:,:),fld_gll(:,:,:,:,:)
        2277             :    logical                          :: llimiter(1)
        2278             :    !----------------------------------------------------------------------------
        2279             :
        2280             :    !!$OMP PARALLEL NUM_THREADS(horz_num_threads), DEFAULT(SHARED), PRIVATE(hybrid,nets,nete,ie)
        2281             :    !hybrid = config_thread_region(par,'horizontal')
        2282           0 :    hybrid = config_thread_region(par,'serial')
        2283             :
        2284           0 :    call get_loop_ranges(hybrid, ibeg=nets, iend=nete)
        2285             :
        2286           0 :    allocate(fld_phys(1-nhc_phys:fv_nphys+nhc_phys,1-nhc_phys:fv_nphys+nhc_phys,1,1,nets:nete))
        2287           0 :    allocate(fld_gll(np,np,1,1,nets:nete))
        2288           0 :    fld_phys = 0.0_r8
        2289           0 :    do ie = nets, nete
        2290           0 :       fld_phys(1:fv_nphys,1:fv_nphys,1,1,ie) = RESHAPE(phis_phys_tmp(:,ie),(/fv_nphys,fv_nphys/))
        2291             :    end do
        2292           0 :    llimiter = .true.
        2293           0 :    call phys2dyn(hybrid,elem,fld_phys,fld_gll,nets,nete,1,1,fvm,llimiter,halo_filled=.false.)
        2294           0 :    do ie = nets,nete
        2295             :       indx = 1
        2296           0 :       do j = 1, np
        2297           0 :          do i = 1, np
        2298           0 :             if (pmask(((ie - 1) * npsq) + indx)) then
        2299           0 :                phis_tmp(indx,ie) = fld_gll(i,j,1,1,ie)
        2300             :             else
        2301           0 :                phis_tmp(indx,ie) = 0.0_r8
        2302             :             end if
        2303           0 :             indx = indx + 1
        2304             :          end do
        2305             :       end do
        2306             :    end do
        2307           0 :    deallocate(fld_phys)
        2308           0 :    deallocate(fld_gll)
        2309             :    !!$OMP END PARALLEL
        2310           0 : end subroutine map_phis_from_physgrid_to_gll
        2311             :
        2312             : !========================================================================================
        2313             :
        2314      369408 : subroutine write_dyn_vars(dyn_out)
        2315             :
        2316             :    type (dyn_export_t), intent(inout) :: dyn_out ! Dynamics export container
        2317             :
        2318             :    character(len=fieldname_len) :: tfname
        2319             :    integer                      :: ie, m
        2320             :    !----------------------------------------------------------------------------
        2321             :
        2322      369408 :    if (use_cslam) then
        2323     2966808 :       do ie = 1, nelemd
        2324           0 :          call outfld('dp_fvm', RESHAPE(dyn_out%fvm(ie)%dp_fvm(1:nc,1:nc,:), &
        2325     2597400 :                                        (/nc*nc,nlev/)), nc*nc, ie)
        2326           0 :          call outfld('PSDRY_fvm', RESHAPE(dyn_out%fvm(ie)%psc(1:nc,1:nc), &
        2327     2597400 :                                           (/nc*nc/)), nc*nc, ie)
        2328    31538208 :          do m = 1, ntrac
        2329    28571400 :             tfname = trim(cnst_name(m))//'_fvm'
        2330           0 :             call outfld(tfname, RESHAPE(dyn_out%fvm(ie)%c(1:nc,1:nc,:,m), &
        2331    28571400 :                                         (/nc*nc,nlev/)), nc*nc, ie)
        2332             :
        2333    28571400 :             tfname = 'F'//trim(cnst_name(m))//'_fvm'
        2334    57142800 :             call outfld(tfname, RESHAPE(dyn_out%fvm(ie)%fc(1:nc,1:nc,:,m),&
        2335    88311600 :                                         (/nc*nc,nlev/)), nc*nc, ie)
        2336             :          end do
        2337             :       end do
        2338             :    end if
        2339             :
        2340           0 : end subroutine write_dyn_vars
        2341             :
        2342             : !=========================================================================================
        2343           0 : end module dyn_comp

Generated by: LCOV version 1.14