LCOV - code coverage report
Current view: | top level - dynamics/se/dycore - prim_driver_mod.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 150 | 235 | 63.8 %
Date: | 2025-01-13 21:54:50 |  | Functions: | 5 | 8 | 62.5 %

              Line data    Source code
           1             : !#define _DBG_ print *,"file: ",__FILE__," line: ",__LINE__," ithr: ",hybrid%ithr
           2             : #define _DBG_
           3             : module prim_driver_mod
           4             :   use shr_kind_mod,           only: r8=>shr_kind_r8
           5             :   use cam_logfile,            only: iulog
           6             :   use cam_abortutils,         only: endrun
           7             :   use dimensions_mod,         only: np, nlev, nelem, nelemd, GlobalUniqueCols, qsize, nc,nhc
           8             :   use hybrid_mod,             only: hybrid_t, config_thread_region, PrintHybrid
           9             :   use derivative_mod,         only: derivative_t
          10             :   use fvm_control_volume_mod, only: fvm_struct
          11             :
          12             :   use element_mod,            only: element_t, timelevels, allocate_element_desc
          13             :   use thread_mod ,            only: horz_num_threads, vert_num_threads, tracer_num_threads
          14             :   use thread_mod ,            only: omp_set_nested
          15             :   use perf_mod,               only: t_startf, t_stopf
          16             :   use prim_init,              only: gp, fvm_corners, fvm_points
          17             :
          18             :   implicit none
          19             :   private
          20             :   public :: prim_init2, prim_run_subcycle, prim_finalize
          21             :   public :: prim_set_dry_mass
          22             : contains
          23             :
          24             : !=============================================================================!
          25             :
          26        1536 :   subroutine prim_init2(elem, fvm, hybrid, nets, nete, tl, hvcoord)
          27             :     use dimensions_mod,         only: irecons_tracer, fvm_supercycling
          28             :     use dimensions_mod,         only: fv_nphys, nc
          29             :     use parallel_mod,           only: syncmp
          30             :     use se_dyn_time_mod,        only: timelevel_t, tstep, phys_tscale, nsplit, TimeLevel_Qdp
          31             :     use se_dyn_time_mod,        only: nsplit_baseline,rsplit_baseline
          32             :     use prim_state_mod,         only: prim_printstate
          33             :     use control_mod,            only: runtype, topology, rsplit, qsplit, rk_stage_user,         &
          34             :                                       nu, nu_q, nu_div, hypervis_subcycle, hypervis_subcycle_q, &
          35             :                                       hypervis_subcycle_sponge, variable_nsplit
          36             :     use fvm_mod,                only: fill_halo_fvm,ghostBufQnhc_h
          37             :     use thread_mod,             only: omp_get_thread_num
          38             :     use global_norms_mod,       only: print_cfl
          39             :     use hybvcoord_mod,          only: hvcoord_t
          40             :     use prim_advection_mod,     only: prim_advec_init2,deriv
          41             :     use prim_advance_mod,       only: compute_omega
          42             :     use physconst,              only: rga, cappa, cpair, tref, lapse_rate
          43             :     use cam_thermo,             only: get_dp_ref
          44             :     use physconst,              only: pstd
          45             :
          46             :     type (element_t), intent(inout) :: elem(:)
          47             :     type (fvm_struct), intent(inout)    :: fvm(:)
          48             :     type (hybrid_t), intent(in) :: hybrid
          49             :
          50             :     type (TimeLevel_t), intent(inout)    :: tl              ! time level struct
          51             :     type (hvcoord_t), intent(inout)      :: hvcoord         ! hybrid vertical coordinate struct
          52             :
          53             :     integer, intent(in)                     :: nets  ! starting thread element number (private)
          54             :     integer, intent(in)                     :: nete  ! ending thread element number   (private)
          55             :
          56             :
          57             :     ! ==================================
          58             :     ! Local variables
          59             :     ! ==================================
          60             :
          61             : !   variables used to calculate CFL
          62             :     real (kind=r8) :: dtnu            ! timestep*viscosity parameter
          63             :     real (kind=r8) :: dt_dyn_del2_sponge
          64             :     real (kind=r8) :: dt_tracer_vis      ! viscosity timestep used in tracers
          65             :     real (kind=r8) :: dt_dyn_vis      ! viscosity timestep
          66             :     real (kind=r8) :: dt_remap        ! remapping timestep
          67             :
          68             :     real (kind=r8) :: dp,dp0,T1,T0,pmid_ref(np,np)
          69        3072 :     real (kind=r8) :: ps_ref(np,np,nets:nete)
          70             :
          71             :     integer :: i,j,k,ie,t,q
          72             :     integer :: n0,n0_qdp
          73             :
          74             :
          75       12336 :     do ie=nets,nete
          76    12085200 :       elem(ie)%derived%FM=0.0_r8
          77     5907600 :       elem(ie)%derived%FT=0.0_r8
          78    17735136 :       elem(ie)%derived%FQ=0.0_r8
          79             :     end do
          80             :
          81             :     ! ==========================
          82             :     ! begin executable code
          83             :     ! ==========================
          84             :     !call prim_advance_init(hybrid%par,elem)
          85             :
          86             :     ! compute most restrictive dt*nu for use by variable res viscosity:
          87             :     ! compute timestep seen by viscosity operator:
          88        1536 :     dt_dyn_vis = tstep
          89        1536 :     dt_dyn_del2_sponge = tstep
          90        1536 :     dt_tracer_vis=tstep*qsplit
          91        1536 :     dt_remap=dt_tracer_vis*rsplit
          92             :     ! compute most restrictive condition:
          93             :     ! note: dtnu ignores subcycling
          94        1536 :     dtnu=max(dt_dyn_vis*max(nu,nu_div), dt_tracer_vis*nu_q)
          95             :     ! compute actual viscosity timesteps with subcycling
          96        1536 :     dt_tracer_vis = dt_tracer_vis/hypervis_subcycle_q
          97        1536 :     dt_dyn_vis = dt_dyn_vis/hypervis_subcycle
          98        1536 :     dt_dyn_del2_sponge = dt_dyn_del2_sponge/hypervis_subcycle_sponge
          99        1536 :     if (variable_nsplit) then
         100           0 :        nsplit_baseline=nsplit
         101           0 :        rsplit_baseline=rsplit
         102             :     end if
         103             :     ! ==================================
         104             :     ! Initialize derivative structure
         105             :     ! ==================================
         106        1536 :     call Prim_Advec_Init2(fvm_corners, fvm_points)
         107        1536 :     if (fv_nphys>0.and.nc.ne.fv_nphys) then
         108             :       !
         109             :       ! need to fill halo for dp_coupling for fvm2phys mapping
         110             :       !
         111           0 :       call fill_halo_fvm(ghostBufQnhc_h,elem,fvm,hybrid,nets,nete,nhc,1,nlev,nlev)
         112             :     end if
         113             : !    !$OMP BARRIER
         114             : !    if (hybrid%ithr==0) then
         115             : !       call syncmp(hybrid%par)
         116             : !    end if
         117             : !    !$OMP BARRIER
         118             :
         119        1536 :     if (topology /= "cube") then
         120           0 :        call endrun('Error: only cube topology supported for primaitve equations')
         121             :     endif
         122             :
         123             :     ! CAM has set tstep based on dtime before calling prim_init2(),
         124             :     ! so only now does HOMME learn the timstep.  print them out:
         125             :     call print_cfl(elem,hybrid,nets,nete,dtnu,&
         126             :          !p top and p mid levels
         127             :          hvcoord%hyai(1)*hvcoord%ps0,hvcoord%hyam(:)*hvcoord%ps0+hvcoord%hybm(:)*pstd,&
         128             :          !dt_remap,dt_tracer_fvm,dt_tracer_se
         129             :          tstep*qsplit*rsplit,tstep*qsplit*fvm_supercycling,tstep*qsplit,&
         130             :          !dt_dyn,dt_dyn_visco,dt_tracer_visco, dt_phys
         131       41472 :          tstep,dt_dyn_vis,dt_dyn_del2_sponge,dt_tracer_vis,tstep*nsplit*qsplit*rsplit)
         132             :
         133        1536 :     if (hybrid%masterthread) then
         134           2 :        if (phys_tscale/=0) then
         135           0 :           write(iulog,'(a,2f9.2)') "CAM physics timescale:        ",phys_tscale
         136             :        endif
         137           2 :        write(iulog,'(a,2f9.2)') "CAM dtime (dt_phys):             ",tstep*nsplit*qsplit*rsplit
         138             :
         139           2 :        write(iulog,*) "CAM-SE uses dry-mass vertical coordinates"
         140             :      end if
         141             :
         142        1536 :      n0=tl%n0
         143        1536 :      call TimeLevel_Qdp( tl, qsplit, n0_qdp)
         144        1536 :      call compute_omega(hybrid,n0,n0_qdp,elem,deriv,nets,nete,dt_remap,hvcoord)
         145             :      !
         146             :      ! pre-compute pressure-level thickness reference profile
         147             :      !
         148       12336 :      do ie=nets,nete
         149       10800 :        call get_dp_ref(hvcoord%hyai, hvcoord%hybi, hvcoord%ps0, elem(ie)%state%phis(:,:), &
         150       23136 :             elem(ie)%derived%dp_ref(:,:,:), ps_ref(:,:,ie))
         151             :      end do
         152             :      !
         153             :      ! pre-compute reference temperature profile (Simmons and Jiabin, 1991, QJRMS, Section 2a
         154             :      !                                            doi: https://doi.org/10.1002/qj.49711749703c)
         155             :      !
         156             :      !  Tref = T0+T1*Exner
         157             :      !  T1 = .0065*Tref*Cp/g ! = ~191
         158             :      !  T0 = Tref-T1         ! = ~97
         159             :      !
         160        1536 :      T1 = lapse_rate*Tref*cpair*rga
         161        1536 :      T0 = Tref-T1
         162       12336 :      do ie=nets,nete
         163      293136 :        do k=1,nlev
         164     5896800 :          pmid_ref =hvcoord%hyam(k)*hvcoord%ps0 + hvcoord%hybm(k)*ps_ref(:,:,ie)
         165      280800 :          dp0 = ( hvcoord%hyai(k+1) - hvcoord%hyai(k) )*hvcoord%ps0 + &
         166      280800 :                ( hvcoord%hybi(k+1) - hvcoord%hybi(k) )*hvcoord%ps0
         167      291600 :          if (hvcoord%hybm(k)>0) then
         168     4309200 :            elem(ie)%derived%T_ref(:,:,k)    = T0+T1*(pmid_ref/hvcoord%ps0)**cappa
         169             :            !
         170             :            ! pel@ucar.edu: resolved noise issue over Antartica
         171             :            !
         172     4309200 :            elem(ie)%derived%dp_ref(:,:,k)   = elem(ie)%derived%dp_ref(:,:,k)-dp0
         173             :          else
         174     1587600 :            elem(ie)%derived%T_ref(:,:,k)    = 0.0_r8
         175             :          end if
         176             :        end do
         177             :      end do
         178             :
         179        1536 :      if (hybrid%masterthread) write(iulog,*) "initial state:"
         180        1536 :      call prim_printstate(elem, tl, hybrid,nets,nete, fvm)
         181             :
         182        1536 :   end subroutine prim_init2
         183             :
         184             : !=======================================================================================================!
         185             :
         186             :
         187      738816 :   subroutine prim_run_subcycle(elem, fvm, hybrid,nets,nete, dt, tl, hvcoord,nsubstep, single_column, omega_cn)
         188             : !
         189             : !   advance all variables (u,v,T,ps,Q,C) from time t to t + dt_q
         190             : !
         191             : ! for the RK schemes:
         192             : !   input:
         193             : !       tl%nm1   not used
         194             : !       tl%n0    data at time t
         195             : !       tl%np1   new values at t+dt_q
         196             : !
         197             : !   then we update timelevel pointers:
         198             : !       tl%nm1 = tl%n0
         199             : !       tl%n0  = tl%np1
         200             : !   so that:
         201             : !       tl%nm1   tracers:  t    dynamics:  t+(qsplit-1)*dt
         202             : !       tl%n0    time t + dt_q
         203             : !
         204             : ! for the implicit schemes:
         205             : !
         206             : !   input:
         207             : !       tl%nm1   variables at t-1 level are stored fro BDF2 scheme
         208             : !       tl%n0    data at time t
         209             : !       tl%np1   new values at t+dt_q
         210             : !       generally dt_q = t for BDF2, so its t+1
         211             : !
         212             : !   then we update timelevel pointers:
         213             : !       tl%nm1 = tl%n0
         214             : !       tl%n0  = tl%np1
         215             : !   so that:
         216             : !       tl%nm1   tracers:  t    dynamics:  t+(qsplit-1)*dt
         217             : !       tl%n0    time t + dt_q
         218             : !
         219             : !
         220        1536 :     use hybvcoord_mod, only : hvcoord_t
         221             :     use se_dyn_time_mod,        only: TimeLevel_t, timelevel_update, timelevel_qdp, nsplit
         222             :     use control_mod,            only: statefreq,qsplit, rsplit, variable_nsplit, dribble_in_rsplit_loop
         223             :     use prim_advance_mod,       only: applycamforcing
         224             :     use prim_advance_mod,       only: tot_energy_dyn,compute_omega
         225             :     use prim_state_mod,         only: prim_printstate, adjust_nsplit
         226             :     use prim_advection_mod,     only: vertical_remap, deriv
         227             :     use thread_mod,             only: omp_get_thread_num
         228             :     use perf_mod   ,            only: t_startf, t_stopf
         229             :     use fvm_mod    ,            only: fill_halo_fvm, ghostBufQnhc_h
         230             :     use dimensions_mod,         only: use_cslam,fv_nphys
         231             :     use fvm_mapping,            only: cslam2gll
         232             :     type (element_t) , intent(inout) :: elem(:)
         233             :     type(fvm_struct), intent(inout)  :: fvm(:)
         234             :     type (hybrid_t), intent(in)      :: hybrid  ! distributed parallel structure (shared)
         235             :     type (hvcoord_t), intent(in)     :: hvcoord         ! hybrid vertical coordinate struct
         236             :     integer, intent(in)              :: nets  ! starting thread element number (private)
         237             :     integer, intent(in)              :: nete  ! ending thread element number   (private)
         238             :     real(kind=r8), intent(in)        :: dt  ! "timestep dependent" timestep
         239             :     type (TimeLevel_t), intent(inout):: tl
         240             :     integer, intent(in)              :: nsubstep  ! nsubstep = 1 .. nsplit
         241             :     logical,              intent(in) :: single_column
         242             :     real (kind=r8)    , intent(inout):: omega_cn(2,nets:nete) !min and max of vertical Courant number
         243             :
         244             :     real(kind=r8)   :: dt_q, dt_remap, dt_phys
         245             :     integer         :: ie, q,k,n0_qdp,np1_qdp,r, nstep_end,region_num_threads,i,j
         246             :     real (kind=r8)  :: dp_np1(np,np)
         247     1477632 :     real (kind=r8)  :: dp_start(np,np,nlev+1,nets:nete),dp_end(np,np,nlev,nets:nete)
         248             :     logical         :: compute_diagnostics
         249             :     ! ===================================
         250             :     ! Main timestepping loop
         251             :     ! ===================================
         252      738816 :     dt_q = dt*qsplit
         253      738816 :     nstep_end = tl%nstep + qsplit
         254      738816 :     dt_remap=dt_q*rsplit
         255      738816 :     nstep_end = tl%nstep + qsplit*rsplit  ! nstep at end of this routine
         256      738816 :     dt_phys   = nsplit*dt_remap
         257             :
         258             :     ! compute diagnostics for STDOUT
         259      738816 :     compute_diagnostics=.false.
         260             :
         261      738816 :     if (statefreq>0) then
         262           0 :       if (MODULO(nstep_end,statefreq)==0 .or. nstep_end==tl%nstep0) then
         263           0 :         compute_diagnostics=.true.
         264             :       endif
         265             :     end if
         266             :     !
         267             :     ! initialize variables for computing vertical Courant number
         268             :     !
         269      738816 :     if (variable_nsplit.or.compute_diagnostics) then
         270           0 :       if (nsubstep==1) then
         271           0 :         do ie=nets,nete
         272           0 :           omega_cn(1,ie) = 0.0_r8
         273           0 :           omega_cn(2,ie) = 0.0_r8
         274             :         end do
         275             :       end if
         276           0 :       do ie=nets,nete
         277           0 :         dp_start(:,:,1:nlev,ie) = elem(ie)%state%dp3d(:,:,:,tl%n0)
         278           0 :         dp_start(:,:,nlev+1,ie) = elem(ie)%state%dp3d(:,:,nlev,tl%n0)
         279             :       end do
         280             :     endif
         281             :
         282             :
         283      738816 :     call TimeLevel_Qdp( tl, qsplit, n0_qdp)
         284             :
         285      738816 :     if (dribble_in_rsplit_loop==0) then
         286      738816 :       call tot_energy_dyn(elem,fvm,nets,nete,tl%n0,n0_qdp,'dAF')
         287      738816 :       call ApplyCAMForcing(elem,fvm,tl%n0,n0_qdp,dt_remap,dt_phys,nets,nete,nsubstep)
         288      738816 :       call tot_energy_dyn(elem,fvm,nets,nete,tl%n0,n0_qdp,'dBD')
         289             :     end if
         290     2955264 :     do r=1,rsplit
         291     2216448 :       if (r.ne.1) call TimeLevel_update(tl,"leapfrog")
         292             :       !
         293             :       ! if nsplit==1 and physics time-step is long then there will be noise in the
         294             :       ! pressure field; hence "dripple" in tendencies
         295             :       !
         296     2216448 :       if (dribble_in_rsplit_loop==1) then
         297           0 :          call tot_energy_dyn(elem,fvm,nets,nete,tl%n0,n0_qdp,'dAF')
         298           0 :          call ApplyCAMForcing(elem,fvm,tl%n0,n0_qdp,dt,dt_phys,nets,nete,MAX(nsubstep,r))
         299           0 :          call tot_energy_dyn(elem,fvm,nets,nete,tl%n0,n0_qdp,'dBD')
         300             :       end if
         301             :       !
         302             :       ! right after physics overwrite Qdp with CSLAM values
         303             :       !
         304     2216448 :       if (use_cslam.and.nsubstep==1.and.r==1) then
         305      369408 :          call tot_energy_dyn(elem,fvm,nets,nete,tl%n0,n0_qdp,'dAF')
         306      369408 :          call cslam2gll(elem, fvm, hybrid,nets,nete, tl%n0, n0_qdp)
         307      369408 :          call tot_energy_dyn(elem,fvm,nets,nete,tl%n0,n0_qdp,'dBD')
         308             :       end if
         309     2216448 :       call tot_energy_dyn(elem,fvm,nets,nete,tl%n0,n0_qdp,'dBL')
         310     2216448 :       if (single_column) then
         311             :          ! Single Column Case
         312             :          ! Loop over rsplit vertically lagrangian timesteps
         313           0 :          call prim_step_scm(elem, fvm, hybrid,nets,nete, dt, tl, hvcoord,r)
         314             :       else
         315     2216448 :          call prim_step(elem, fvm, hybrid,nets,nete, dt, tl, hvcoord,r,nsubstep==nsplit,dt_remap)
         316             :       end if
         317     2955264 :       call tot_energy_dyn(elem,fvm,nets,nete,tl%np1,n0_qdp,'dAL')
         318             :     enddo
         319             :
         320             :
         321             :     ! defer final timelevel update until after remap and diagnostics
         322      738816 :     call TimeLevel_Qdp( tl, qsplit, n0_qdp, np1_qdp)
         323             :
         324             :     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
         325             :     !
         326             :     !  apply vertical remap
         327             :     !  always for tracers
         328             :     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
         329             :
         330      738816 :     call tot_energy_dyn(elem,fvm,nets,nete,tl%np1,np1_qdp,'dAD')
         331             :
         332      738816 :     if (variable_nsplit.or.compute_diagnostics) then
         333             :       !
         334             :       ! initialize variables for computing vertical Courant number
         335             :       !
         336           0 :       do ie=nets,nete
         337           0 :         dp_end(:,:,:,ie) = elem(ie)%state%dp3d(:,:,:,tl%np1)
         338             :       end do
         339             :     end if
         340      738816 :     call t_startf('vertical_remap')
         341      738816 :     call vertical_remap(hybrid,elem,fvm,hvcoord,tl%np1,np1_qdp,nets,nete)
         342      738816 :     call t_stopf('vertical_remap')
         343             :
         344             :     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
         345             :     ! time step is complete.
         346             :     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
         347      738816 :     call tot_energy_dyn(elem,fvm,nets,nete,tl%np1,np1_qdp,'dAR')
         348             :
         349      738816 :     if (nsubstep==nsplit.and. .not. single_column) then
         350      369408 :       call compute_omega(hybrid,tl%np1,np1_qdp,elem,deriv,nets,nete,dt_remap,hvcoord)
         351             :     end if
         352             :
         353             :     ! now we have:
         354             :     !   u(nm1)   dynamics at  t+dt_remap - 2*dt
         355             :     !   u(n0)    dynamics at  t+dt_remap - dt
         356             :     !   u(np1)   dynamics at  t+dt_remap
         357             :     !
         358             :     !   Q(1)   Q at t+dt_remap
         359             :
         360             :
         361             :
         362             :     ! =================================
         363             :     ! update dynamics time level pointers
         364             :     ! =================================
         365      738816 :     call TimeLevel_update(tl,"leapfrog")
         366             :     ! note: time level update for fvm tracers takes place in fvm_mod
         367             :
         368             :     ! now we have:
         369             :     !   u(nm1)   dynamics at  t+dt_remap - dt       (Robert-filtered)
         370             :     !   u(n0)    dynamics at  t+dt_remap
         371             :     !   u(np1)   undefined
         372             :
         373             :
         374             :     !
         375             :     ! Compute vertical Courant numbers
         376             :     !
         377      738816 :     if (variable_nsplit.or.compute_diagnostics) then
         378           0 :       do ie=nets,nete
         379           0 :         do k=1,nlev
         380           0 :           do j=1,np
         381           0 :             do i=1,np
         382           0 :               if (dp_end(i,j,k,ie)<dp_start(i,j,k,ie)) then
         383           0 :                 omega_cn(1,ie) = MIN((dp_end(i,j,k,ie)-dp_start(i,j,k,ie))/dp_start(i,j,k,ie),omega_cn(1,ie))
         384           0 :                 omega_cn(2,ie) = MAX((dp_end(i,j,k,ie)-dp_start(i,j,k,ie))/dp_start(i,j,k,ie),omega_cn(2,ie))
         385             :               else
         386           0 :                 omega_cn(1,ie) = MIN((dp_end(i,j,k,ie)-dp_start(i,j,k,ie))/dp_start(i,j,k+1,ie),omega_cn(1,ie))
         387           0 :                 omega_cn(2,ie) = MAX((dp_end(i,j,k,ie)-dp_start(i,j,k,ie))/dp_start(i,j,k+1,ie),omega_cn(2,ie))
         388             :               end if
         389             :             end do
         390             :           end do
         391             :         end do
         392             :       end do
         393           0 :       if (nsubstep==nsplit.and.variable_nsplit) then
         394           0 :          call t_startf('adjust_nsplit')
         395           0 :          call adjust_nsplit(elem, tl, hybrid,nets,nete, fvm, omega_cn)
         396           0 :          call t_stopf('adjust_nsplit')
         397             :       end if
         398             :     end if
         399             :
         400             :     ! ============================================================
         401             :     ! Print some diagnostic information
         402             :     ! ============================================================
         403           0 :     if (compute_diagnostics) then
         404           0 :       call prim_printstate(elem, tl, hybrid,nets,nete, fvm, omega_cn)
         405             :     end if
         406             :
         407      738816 :     if (use_cslam.and.nsubstep==nsplit.and.nc.ne.fv_nphys) then
         408             :       !
         409             :       ! fill the fvm halo for mapping in d_p_coupling if
         410             :       ! physics grid resolution is different than fvm resolution
         411             :       !
         412           0 :       call fill_halo_fvm(ghostBufQnhc_h, elem,fvm,hybrid,nets,nete,nhc,1,nlev,nlev)
         413             :     end if
         414             :
         415      738816 :   end subroutine prim_run_subcycle
         416             :
         417             :
         418     2216448 :   subroutine prim_step(elem, fvm, hybrid,nets,nete, dt, tl, hvcoord, rstep, last_step,dt_remap)
         419             :     !
         420             :     !   Take qsplit dynamics steps and one tracer step
         421             :     !   for vertically lagrangian option, this subroutine does only the horizontal step
         422             :     !
         423             :     !   input:
         424             :     !       tl%nm1   not used
         425             :     !       tl%n0    data at time t
         426             :     !       tl%np1   new values at t+dt_q
         427             :     !
         428             :     !   then we update timelevel pointers:
         429             :     !       tl%nm1 = tl%n0
         430             :     !       tl%n0  = tl%np1
         431             :     !   so that:
         432             :     !       tl%nm1   tracers:  t    dynamics:  t+(qsplit-1)*dt
         433             :     !       tl%n0    time t + dt_q
         434             :     !
         435      738816 :     use hybvcoord_mod,          only: hvcoord_t
         436             :     use se_dyn_time_mod,        only: TimeLevel_t, timelevel_update
         437             :     use control_mod,            only: statefreq, qsplit, nu_p
         438             :     use thread_mod,             only: omp_get_thread_num
         439             :     use prim_advance_mod,       only: prim_advance_exp
         440             :     use prim_advection_mod,     only: prim_advec_tracers_remap, prim_advec_tracers_fvm, deriv
         441             :     use derivative_mod,         only: subcell_integration
         442             :     use hybrid_mod,             only: set_region_num_threads, config_thread_region, get_loop_ranges
         443             :     use dimensions_mod,         only: use_cslam,fvm_supercycling,fvm_supercycling_jet
         444             :     use dimensions_mod,         only: kmin_jet, kmax_jet
         445             :     use fvm_mod,                only: ghostBufQnhc_vh,ghostBufQ1_vh, ghostBufFlux_vh
         446             :     use fvm_mod,                only: ghostBufQ1_h,ghostBufQnhcJet_h, ghostBufFluxJet_h
         447             :     use se_dyn_time_mod,        only: timelevel_qdp
         448             :     use fvm_mapping,            only: cslam2gll
         449             : #ifdef waccm_debug
         450             :   use cam_history, only: outfld
         451             : #endif
         452             :
         453             :
         454             :     type (element_t) ,  intent(inout) :: elem(:)
         455             :     type(fvm_struct),   intent(inout) :: fvm(:)
         456             :     type (hybrid_t),    intent(in)    :: hybrid  ! distributed parallel structure (shared)
         457             :     type (hvcoord_t),   intent(in)    :: hvcoord         ! hybrid vertical coordinate struct
         458             :     integer,            intent(in)    :: nets  ! starting thread element number (private)
         459             :     integer,            intent(in)    :: nete  ! ending thread element number   (private)
         460             :     real(kind=r8),      intent(in)    :: dt  ! "timestep dependent" timestep
         461             :     type (TimeLevel_t), intent(inout) :: tl
         462             :     integer, intent(in)               :: rstep ! vertical remap subcycling step
         463             :     logical, intent(in)               :: last_step! last step before d_p_coupling
         464             :     real(kind=r8), intent(in)         :: dt_remap
         465             :
         466             :     type (hybrid_t):: hybridnew,hybridnew2
         467             :     real(kind=r8)  :: st, st1, dp, dt_q
         468             :     integer        :: ie,t,q,k,i,j,n, n_Q
         469             :     integer        :: ithr
         470             :     integer        :: region_num_threads
         471             :     integer        :: kbeg,kend
         472             :     integer        :: n0_qdp, np1_qdp
         473             :
         474             :     real (kind=r8) :: tempdp3d(np,np), x
         475             :     real (kind=r8) :: tempmass(nc,nc)
         476             :     real (kind=r8) :: tempflux(nc,nc,4)
         477             :
         478             :     real (kind=r8) :: dp_np1(np,np)
         479             :
         480             :
         481     2216448 :     dt_q = dt*qsplit
         482             :     ! ===============
         483             :     ! initialize mean flux accumulation variables and save some variables at n0
         484             :     ! for use by advection
         485             :     ! ===============
         486    17800848 :     do ie=nets,nete
         487 17438943600 :       elem(ie)%derived%vn0=0              ! mean horizontal mass flux
         488  8524666800 :       elem(ie)%derived%omega=0
         489    15584400 :       if (nu_p>0) then
         490  8524666800 :          elem(ie)%derived%dpdiss_ave=0
         491  8524666800 :          elem(ie)%derived%dpdiss_biharmonic=0
         492             :       endif
         493             :
         494             :       ! dp at time t:  use floating lagrangian levels:
         495  8526883248 :       elem(ie)%derived%dp(:,:,:)=elem(ie)%state%dp3d(:,:,:,tl%n0)
         496             :     enddo
         497             :
         498             :     ! ===============
         499             :     ! Dynamical Step
         500             :     ! ===============
         501     2216448 :     n_Q = tl%n0  ! n_Q = timelevel of FV tracers at time t.  need to save this
         502             :                  ! SE tracers only carry 2 timelevels
         503             :
         504     2216448 :     call t_startf('prim_advance_exp')
         505             : !    ithr   = 0 ! omp_get_thread_num()
         506             : !    vybrid = hybrid_create(hybrid%par,ithr)
         507             :
         508             :     call prim_advance_exp(elem, fvm, deriv, hvcoord,   &
         509     2216448 :          hybrid, dt, tl, nets, nete)
         510             :
         511     2216448 :     call t_stopf('prim_advance_exp')
         512             :
         513     2216448 :     do n=2,qsplit
         514           0 :        call TimeLevel_update(tl,"leapfrog")
         515             :
         516           0 :        call t_startf('prim_advance_exp')
         517             :
         518             :        call prim_advance_exp(elem, fvm, deriv, hvcoord,   &
         519           0 :             hybrid, dt, tl, nets, nete)
         520             :
         521     2216448 :     call t_stopf('prim_advance_exp')
         522             :
         523             :        ! defer final timelevel update until after Q update.
         524             :   enddo
         525             : #ifdef HOMME_TEST_SUB_ELEMENT_MASS_FLUX
         526             :     if (use_cslam.and.rstep==1) then
         527             :       do ie=nets,nete
         528             :       do k=1,nlev
         529             :         tempdp3d = elem(ie)%state%dp3d(:,:,k,tl%np1) - &
         530             :                    elem(ie)%derived%dp(:,:,k)
         531             :         call subcell_integration(tempdp3d, np, nc, elem(ie)%metdet,tempmass)
         532             :         tempflux = dt_q*elem(ie)%sub_elem_mass_flux(:,:,:,k)
         533             :         do i=1,nc
         534             :         do j=1,nc
         535             :           x = SUM(tempflux(i,j,:))
         536             :           if (ABS(tempmass(i,j)).lt.1e-11_r8 .and. 1e-11_r8.lt.ABS(x)) then
         537             :             write(iulog,*) __FILE__,__LINE__,"**CSLAM mass-flux ERROR***",ie,k,i,j,tempmass(i,j),x
         538             :             call endrun('**CSLAM mass-flux ERROR***')
         539             :           elseif (1e-5_r8.lt.ABS((tempmass(i,j)-x)/tempmass(i,j))) then
         540             :             write(iulog,*) __FILE__,__LINE__,"**CSLAM mass-flux ERROR**",ie,k,i,j,tempmass(i,j),x,&
         541             :                    ABS((tempmass(i,j)-x)/tempmass(i,j))
         542             :             call endrun('**CSLAM mass-flux ERROR**')
         543             :           endif
         544             :         end do
         545             :         end do
         546             :       end do
         547             :       end do
         548             :     end if
         549             : #endif
         550             :     ! current dynamics state variables:
         551             :     !    derived%dp              =  dp at start of timestep
         552             :     !    derived%vn0             =  mean horiz. flux:   U*dp
         553             :     ! rsplit>0
         554             :     !        state%v(:,:,:,np1)      = velocity on lagrangian levels
         555             :     !        state%dp3d(:,:,:,np1)   = dp3d
         556             :     !
         557             :
         558             :
         559             :     ! ===============
         560             :     ! Tracer Advection.
         561             :     ! in addition, this routine will apply the DSS to:
         562             :     !        derived%omega           =
         563             :     ! Tracers are always vertically lagrangian.
         564             :     ! ===============
         565             :     ! Advect tracers if their count is > 0.
         566             :     ! special case in CAM: if CSLAM tracers are turned on , qsize=1 but this tracer should
         567             :     ! not be advected.  This will be cleaned up when the physgrid is merged into CAM trunk
         568             :     ! Currently advecting all species
         569     2216448 :     if (.not.use_cslam) then
         570           0 :       call t_startf('prim_advec_tracers_remap')
         571           0 :       region_num_threads=tracer_num_threads
         572           0 :       call omp_set_nested(.true.)
         573             :       !$OMP PARALLEL NUM_THREADS(region_num_threads), DEFAULT(SHARED), PRIVATE(hybridnew)
         574           0 :       hybridnew = config_thread_region(hybrid,'tracer')
         575           0 :       call Prim_Advec_Tracers_remap(elem, deriv,hvcoord,hybridnew,dt_q,tl,nets,nete)
         576             :       !$OMP END PARALLEL
         577           0 :       call omp_set_nested(.false.)
         578           0 :       call t_stopf('prim_advec_tracers_remap')
         579             :    else
         580             :       !
         581             :       ! only run fvm transport every fvm_supercycling rstep
         582             :       !
         583             :       ! FVM transport
         584             :       !
         585     2216448 :       if ((mod(rstep,fvm_supercycling) == 0).and.(mod(rstep,fvm_supercycling_jet) == 0)) then
         586             :
         587             : !        call omp_set_nested(.true.)
         588             : !        !$OMP PARALLEL NUM_THREADS(vert_num_threads), DEFAULT(SHARED), PRIVATE(hybridnew2,kbeg,kend)
         589             : !        hybridnew2 = config_thread_region(hybrid,'vertical')
         590             : !        call get_loop_ranges(hybridnew2,kbeg=kbeg,kend=kend)
         591             :         call Prim_Advec_Tracers_fvm(elem,fvm,hvcoord,hybrid,&
         592      738816 :              dt_q,tl,nets,nete,ghostBufQnhc_vh,ghostBufQ1_vh, ghostBufFlux_vh,1,nlev)
         593             : !        !$OMP END PARALLEL
         594             : !        call omp_set_nested(.false.)
         595             :         !
         596             :         ! to avoid accumulation of truncation error overwrite CSLAM surface pressure with SE
         597             :         ! surface pressure
         598             :         !
         599     5933616 :         do ie=nets,nete
         600             :           !
         601             :           ! overwrite PSDRY on CSLAM grid with SE PSDRY integrated over CSLAM control volume
         602             :           !
         603             :           !          call subcell_integration(elem(ie)%state%psdry(:,:), np, nc, elem(ie)%metdet,fvm(ie)%psc)
         604             :           !          fvm(ie)%psc = fvm(ie)%psc*fvm(ie)%inv_se_area_sphere
         605             :           !
         606             :           ! Update CSLAM surface pressure
         607             :           !
         608    21518016 :           do j=1,nc
         609    67532400 :             do i=1,nc
         610  1277920800 :               fvm(ie)%psc(i,j) = sum(fvm(ie)%dp_fvm(i,j,:)) +  hvcoord%hyai(1)*hvcoord%ps0
         611             :             end do
         612             :           end do
         613             :        end do
         614      738816 :        call TimeLevel_Qdp( tl, qsplit, n0_qdp, np1_qdp)
         615      738816 :        if (.not.last_step) call cslam2gll(elem, fvm, hybrid,nets,nete, tl%np1, np1_qdp)
         616     1477632 :       else if ((mod(rstep,fvm_supercycling_jet) == 0)) then
         617             :         !
         618             :         ! shorter fvm time-step in jet region
         619             :         !
         620             :         call Prim_Advec_Tracers_fvm(elem,fvm,hvcoord,hybrid,&
         621           0 :              dt_q,tl,nets,nete,ghostBufQnhcJet_h,ghostBufQ1_h, ghostBufFluxJet_h,kmin_jet,kmax_jet)
         622             :       end if
         623             :
         624             : #ifdef waccm_debug
         625             :       do ie=nets,nete
         626             :         call outfld('CSLAM_gamma', RESHAPE(fvm(ie)%CSLAM_gamma(:,:,:,1), &
         627             :              (/nc*nc,nlev/)), nc*nc, ie)
         628             :       end do
         629             : #endif
         630             :    endif
         631             :
         632     2216448 :    end subroutine prim_step
         633           0 :    subroutine prim_step_scm(elem, fvm, hybrid,nets,nete, dt, tl, hvcoord, rstep)
         634             :      !
         635             :      !   prim_step version for single column model (SCM)
         636             :      !   Here we simply want to compute the floating level tendency
         637             :      !   based on the prescribed large scale vertical velocity
         638             :      !   Take qsplit dynamics steps and one tracer step
         639             :      !   for vertically lagrangian option, this subroutine does only
         640             :      !   the horizontal step
         641             :      !
         642             :      !   input:
         643             :      !       tl%nm1   not used
         644             :      !       tl%n0    data at time t
         645             :      !       tl%np1   new values at t+dt_q
         646             :      !
         647             :      !   then we update timelevel pointers:
         648             :      !       tl%nm1 = tl%n0
         649             :      !       tl%n0  = tl%np1
         650             :      !   so that:
         651             :      !       tl%nm1   tracers:  t    dynamics:  t+(qsplit-1)*dt
         652             :      !       tl%n0    time t + dt_q
         653             :      !
         654     2216448 :     use hybvcoord_mod,          only: hvcoord_t
         655             :     use se_dyn_time_mod,        only: TimeLevel_t, timelevel_update
         656             :     use control_mod,            only: statefreq, qsplit, nu_p
         657             :     use prim_advection_mod,     only: deriv
         658             :     use hybrid_mod,             only: config_thread_region, get_loop_ranges
         659             :
         660             :     type (element_t) ,  intent(inout) :: elem(:)
         661             :     type(fvm_struct),   intent(inout) :: fvm(:)
         662             :     type (hybrid_t),    intent(in)    :: hybrid  ! distributed parallel structure (shared)
         663             :     type (hvcoord_t),   intent(in)    :: hvcoord         ! hybrid vertical coordinate struct
         664             :     integer,            intent(in)    :: nets  ! starting thread element number (private)
         665             :     integer,            intent(in)    :: nete  ! ending thread element number   (private)
         666             :     real(kind=r8),      intent(in)    :: dt  ! "timestep dependent" timestep
         667             :     type (TimeLevel_t), intent(inout) :: tl
         668             :     integer, intent(in)               :: rstep ! vertical remap subcycling step
         669             :
         670             :     integer        :: ie,n
         671             :
         672             :     ! ===============
         673             :     ! initialize mean flux accumulation variables and save some variables at n0
         674             :     ! for use by advection
         675             :     ! ===============
         676           0 :     do ie=nets,nete
         677           0 :       elem(ie)%derived%vn0=0              ! mean horizontal mass flux
         678           0 :       if (nu_p>0) then
         679           0 :          elem(ie)%derived%dpdiss_ave=0
         680           0 :          elem(ie)%derived%dpdiss_biharmonic=0
         681             :       endif
         682           0 :       elem(ie)%derived%dp(:,:,:)=elem(ie)%state%dp3d(:,:,:,tl%n0)
         683             :     enddo
         684             :
         685             :     ! ===============
         686             :     ! Dynamical Step
         687             :     ! ===============
         688             :
         689           0 :     call t_startf('set_prescribed_scm')
         690             :
         691             :     call set_prescribed_scm(elem, fvm, deriv, hvcoord,   &
         692           0 :          hybrid, dt, tl, nets, nete)
         693             :
         694           0 :     call t_stopf('set_prescribed_scm')
         695             :
         696           0 :     do n=2,qsplit
         697           0 :        call TimeLevel_update(tl,"leapfrog")
         698             :
         699           0 :        call t_startf('set_prescribed_scm')
         700             :
         701             :        call set_prescribed_scm(elem, fvm, deriv, hvcoord,   &
         702           0 :             hybrid, dt, tl, nets, nete)
         703             :
         704           0 :        call t_stopf('set_prescribed_scm')
         705             :     enddo
         706             :
         707           0 :   end subroutine prim_step_scm
         708             : !=======================================================================================================!
         709             :
         710             :
         711           0 :   subroutine prim_finalize(hybrid)
         712             :     type (hybrid_t), intent(in)           :: hybrid  ! distributed parallel structure (shared)
         713             :
         714             :     ! ==========================
         715             :     ! end of the hybrid program
         716             :     ! ==========================
         717           0 :   end subroutine prim_finalize
         718             :
         719             : !=========================================================================================
         720             :
         721         768 :     subroutine prim_set_dry_mass(elem, hvcoord,initial_global_ave_dry_ps,q)
         722             :       use element_mod,      only: element_t
         723             :       use hybvcoord_mod ,   only: hvcoord_t
         724             :       use dimensions_mod,   only: nelemd, nlev, np
         725             :       use constituents,     only: cnst_type, qmin, pcnst
         726             :       use cam_logfile,      only: iulog
         727             :       use spmd_utils,       only: masterproc
         728             :
         729             :       type (element_t)     , intent(inout):: elem(:)
         730             :       type (hvcoord_t)     , intent(in)   :: hvcoord
         731             :       real (kind=r8), intent(in)   :: initial_global_ave_dry_ps
         732             :       real (kind=r8), intent(inout):: q(np,np,nlev,nelemd,pcnst)
         733             :
         734             :       ! local
         735             :       real (kind=r8)               :: global_ave_ps_inic,dp_tmp, factor(np,np,nlev)
         736             :       integer                      :: ie, i, j ,k, m_cnst
         737             :
         738         768 :       if (initial_global_ave_dry_ps == 0) return;
         739             :
         740         768 :       call get_global_ave_surface_pressure(elem, global_ave_ps_inic)
         741             :
         742        6168 :       do ie=1,nelemd
         743      113400 :         elem(ie)%state%psdry(:,:)=elem(ie)%state%psdry(:,:)*(initial_global_ave_dry_ps/global_ave_ps_inic)
         744      145800 :         do k=1,nlev
         745      707400 :           do j = 1,np
         746     2948400 :             do i = 1,np
         747     4492800 :               dp_tmp =  ((hvcoord%hyai(k+1) - hvcoord%hyai(k))*hvcoord%ps0)+&
         748     6739200 :                    ((hvcoord%hybi(k+1) - hvcoord%hybi(k))*elem(ie)%state%psdry(i,j))
         749     2246400 :               factor(i,j,k) = elem(ie)%state%dp3d(i,j,k,1)/dp_tmp
         750     9547200 :               elem(ie)%state%dp3d(i,j,k,:) = dp_tmp
         751             :             end do
         752             :           end do
         753             :         end do
         754             :         !
         755             :         ! conserve initial condition mass of 'wet' tracers (following dryairm.F90 for FV dycore)
         756             :         ! and conserve mixing ratio (not mass) of 'dry' tracers
         757             :         !
         758       22368 :         do  m_cnst=1,pcnst
         759       21600 :           if (cnst_type(m_cnst).ne.'dry') then
         760      437400 :             do k=1,nlev
         761     2122200 :               do j = 1,np
         762     8845200 :                 do i = 1,np
         763     6739200 :                   q(i,j,k,ie,m_cnst) = q(i,j,k,ie,m_cnst)*factor(i,j,k)
         764     8424000 :                   q(i,j,k,ie,m_cnst) = max(qmin(m_cnst),q(i,j,k,ie,m_cnst))
         765             :                 end do
         766             :               end do
         767             :             end do
         768             :           end if
         769             :         end do
         770             :       end do
         771         768 :       if (masterproc) then
         772           1 :         write (iulog,*) "------ info from prim_set_dry_mass -----------------------------------------------------------"
         773           1 :         write (iulog,*) "Scaling dry surface pressure to global average of = ",&
         774           2 :              initial_global_ave_dry_ps/100.0_r8,"hPa"
         775           1 :         write (iulog,*) "Average dry surface pressure in initial condition = ",global_ave_ps_inic/100.0_r8,"hPa"
         776           1 :         write (iulog,*) "Average dry surface pressure change               = ",&
         777           2 :              initial_global_ave_dry_ps-global_ave_ps_inic,"Pa"
         778           1 :         write (iulog,*) "Mixing ratios that are wet have been scaled so that total mass of tracer is conserved"
         779           1 :         write (iulog,*) "Mixing ratios that are dry have not been changed (mass not conserved in scaling process)"
         780           1 :         write (iulog,*) "------ end info from prim_set_dry_mass -------------------------------------------------------"
         781             :       endif
         782             :     end subroutine prim_set_dry_mass
         783             :
         784         768 :     subroutine get_global_ave_surface_pressure(elem, global_ave_ps_inic)
         785             :       use element_mod       , only : element_t
         786             :       use dimensions_mod    , only : np
         787             :       use global_norms_mod  , only : global_integral
         788             :       use hybrid_mod        , only : config_thread_region, get_loop_ranges, hybrid_t
         789             :       use parallel_mod      , only : par
         790             :
         791             :       type (element_t)     , intent(in)   :: elem(:)
         792             :       real (kind=r8), intent(out)  :: global_ave_ps_inic
         793             :
         794             :       ! local
         795         768 :       real (kind=r8), allocatable  :: tmp(:,:,:)
         796             :       type (hybrid_t)                     :: hybrid
         797             :       integer                             :: ie, nets, nete
         798             :
         799             :       !JMD $OMP PARALLEL NUM_THREADS(horz_num_threads), DEFAULT(SHARED), PRIVATE(hybrid,nets,nete,n)
         800             :       !JMD        hybrid = config_thread_region(par,'horizontal')
         801         768 :       hybrid = config_thread_region(par,'serial')
         802         768 :       call get_loop_ranges(hybrid,ibeg=nets,iend=nete)
         803        2304 :       allocate(tmp(np,np,nets:nete))
         804             :
         805        6168 :       do ie=nets,nete
         806      114168 :         tmp(:,:,ie)=elem(ie)%state%psdry(:,:)
         807             :       enddo
         808         768 :       global_ave_ps_inic = global_integral(elem, tmp(:,:,nets:nete),hybrid,np,nets,nete)
         809         768 :       deallocate(tmp)
         810         768 :     end subroutine get_global_ave_surface_pressure
         811             :
         812           0 :     subroutine set_prescribed_scm(elem, fvm, deriv, hvcoord,   &
         813             :          hybrid, dt, tl, nets, nete)
         814             :       use control_mod,       only: tstep_type, qsplit
         815             :       use derivative_mod,    only: derivative_t
         816             :       use dimensions_mod,    only: np, nlev
         817             :       use element_mod,       only: element_t
         818             :       use hybvcoord_mod,     only: hvcoord_t
         819             :       use hybrid_mod,        only: hybrid_t
         820             :       use se_dyn_time_mod,   only: TimeLevel_t,  timelevel_qdp
         821             :       use fvm_control_volume_mod, only: fvm_struct
         822             :       implicit none
         823             :
         824             :       type (element_t), intent(inout), target   :: elem(:)
         825             :       type(fvm_struct)     , intent(inout) :: fvm(:)
         826             :       type (derivative_t)  , intent(in) :: deriv
         827             :       type (hvcoord_t)                  :: hvcoord
         828             :       type (hybrid_t)      , intent(in) :: hybrid
         829             :       real (kind=r8), intent(in) :: dt
         830             :       type (TimeLevel_t)   , intent(in) :: tl
         831             :       integer              , intent(in) :: nets
         832             :       integer              , intent(in) :: nete
         833             :
         834             :       ! Local
         835             :       integer        :: ie,nm1,n0,np1,k,qn0,qnp1,p
         836             :       real(kind=r8)  :: eta_dot_dpdn(np,np,nlev+1)
         837             :
         838             :
         839           0 :       nm1   = tl%nm1
         840           0 :       n0    = tl%n0
         841           0 :       np1   = tl%np1
         842             :
         843           0 :       call TimeLevel_Qdp(tl, qsplit, qn0, qnp1)  ! compute current Qdp() timelevel
         844             :
         845           0 :       do ie=nets,nete
         846           0 :          do k=1,nlev
         847           0 :             eta_dot_dpdn(:,:,k)=elem(ie)%derived%omega(:,:,k)
         848             :          enddo
         849           0 :          eta_dot_dpdn(:,:,nlev+1) = eta_dot_dpdn(:,:,nlev)
         850             :
         851           0 :          do k=1,nlev
         852           0 :             elem(ie)%state%dp3d(:,:,k,np1) = elem(ie)%state%dp3d(:,:,k,n0) &
         853           0 :                  + dt*(eta_dot_dpdn(:,:,k+1) - eta_dot_dpdn(:,:,k))
         854             :          enddo
         855             :
         856           0 :          do k=1,nlev
         857           0 :             elem(ie)%state%T(:,:,k,np1) = elem(ie)%state%T(:,:,k,n0)
         858             :          enddo
         859             :
         860           0 :          do p=1,qsize
         861           0 :             do k=1,nlev
         862           0 :                elem(ie)%state%Qdp(:,:,k,p,qnp1) = elem(ie)%state%Qdp(:,:,k,p,qn0) &
         863           0 :                     + elem(ie)%state%Qdp(:,:,k,p,qn0)/elem(ie)%state%dp3d(:,:,k,n0) * &
         864           0 :                     dt*(eta_dot_dpdn(:,:,k+1) - eta_dot_dpdn(:,:,k))
         865             :             enddo
         866             :          enddo
         867             :       enddo
         868           0 :     end subroutine set_prescribed_scm
         869             :
         870             : end module prim_driver_mod

Generated by: LCOV version 1.14