text:
	SUBROUTINE IN_MRGD  ( mrgdat, mrgflg, ipttyp, iret )
C************************************************************************
C* IN_MRGD								*
C*									*
C* This subroutine breaks the user input for MRGDAT into a flag 	*
C* indicating whether merged data are to be used and a type for		*
C* unmerged data.  The default for MRGFLG is true and for IPTTYP	*
C* is 3.								*
C*									*
C* IN_MRGD  ( MRGDAT, MRGFLG, IPTTYP, IRET )				*
C*									*
C* Input parameters:							*
C*	MRGDAT		CHAR*		User input for MRGDAT		*
C*									*
C* Output parameters:							*
C*	MRGFLG		LOGICAL		Merged file type		*
C*	IPTTYP		INTEGER		Unmerged type			*
C*					  1 = man below 100 mb		*
C*					  2 = man & sig below 100 mb	*
C*					  3 = man & sig below & above	*
C*	IRET		INTEGER		Return code			*
C*					  0 = normal return		*
C**									*
C* Log:									*
C* M. desJardins/GSFC	 4/89						*
C* M. desJardins/GSFC	 4/90	Change default for mrgdat to YES	*
C************************************************************************
	CHARACTER*(*)	mrgdat
	LOGICAL		mrgflg
C*
	CHARACTER	clist (2)*12
C-----------------------------------------------------------------------
	iret = 0
C
C*	Decode input for MRGDAT.
C
	CALL ST_CLST  ( mrgdat, '/', ' ', 2, clist, n, ier )
	CALL ST_LCUC  ( clist (1), clist (1), ier )
	IF  ( clist (1) (1:1) .eq. 'N' )  THEN
	    mrgflg = .false.
	    CALL ST_NUMB  ( clist (2), ipttyp, ier )
	    IF  ( ( ipttyp .lt. 1 ) .or. ( ipttyp .gt. 3 ) )  ipttyp = 3
	  ELSE
	    mrgflg = .true.
	    ipttyp = 0
	END IF
C*
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
program t
  TYPE PERSON
    INTEGER AGE
  END TYPE PERSON
  print *,kind(AGE)
end program t

max_stars_repo_name: maddenp/ppp

max_issues_repo_name: maddenp/ppp

max_forks_repo_name: maddenp/ppp


text:
module update_scal_module

  use bl_types
  use multifab_module
  use define_bc_module
  use ml_layout_module

  implicit none

  private

  public :: update_scal

contains

  subroutine update_scal(mla,nstart,nstop,sold,snew,sflux,scal_force,p0_new,p0_new_cart, &
                         dx,dt,the_bc_level)

    use bl_prof_module
    use bl_constants_module
    use geometry,  only: spherical
    use variables, only: spec_comp, rho_comp
    use network,   only: nspec
    use ml_restrict_fill_module

    type(ml_layout)   , intent(inout) :: mla
    integer           , intent(in   ) :: nstart, nstop
    type(multifab)    , intent(in   ) :: sold(:)
    type(multifab)    , intent(inout) :: snew(:)
    type(multifab)    , intent(in   ) :: sflux(:,:)
    type(multifab)    , intent(in   ) :: scal_force(:)
    real(kind = dp_t) , intent(in   ) :: p0_new(:,0:)
    type(multifab)    , intent(in   ) :: p0_new_cart(:)
    real(kind = dp_t) , intent(in   ) :: dx(:,:),dt
    type(bc_level)    , intent(in   ) :: the_bc_level(:)

    ! local
    real(kind=dp_t), pointer :: sop(:,:,:,:)
    real(kind=dp_t), pointer :: snp(:,:,:,:)
    real(kind=dp_t), pointer :: sfpx(:,:,:,:)
    real(kind=dp_t), pointer :: sfpy(:,:,:,:)
    real(kind=dp_t), pointer :: sfpz(:,:,:,:)
    real(kind=dp_t), pointer :: fp(:,:,:,:)
    real(kind=dp_t), pointer :: p0p(:,:,:,:)

    integer :: lo(mla%dim),hi(mla%dim),dm,nlevs
    integer :: i,n
    integer :: ng_so,ng_sn,ng_sf,ng_f,ng_p

    type(bl_prof_timer), save :: bpt

    call build(bpt, "update_scal")

    dm = mla%dim
    nlevs = mla%nlevel

    ng_so = nghost(sold(1))
    ng_sn = nghost(snew(1))
    ng_sf = nghost(sflux(1,1))
    ng_f  = nghost(scal_force(1))
    ng_p  = nghost(p0_new_cart(1))

    do n=1,nlevs

       do i = 1, nfabs(sold(n))
          sop => dataptr(sold(n),i)
          snp => dataptr(snew(n),i)
          sfpx => dataptr(sflux(n,1),i)
          fp => dataptr(scal_force(n),i)
          lo =  lwb(get_box(sold(n),i))
          hi =  upb(get_box(sold(n),i))
          select case (dm)
          case (1)
             call update_scal_1d(nstart, nstop, &
                                  sop(:,1,1,:), ng_so, &
                                  snp(:,1,1,:), ng_sn, &
                                 sfpx(:,1,1,:), ng_sf, &
                                   fp(:,1,1,:), ng_f, &
                                 p0_new(n,:), lo, hi, dx(n,:), dt)
          case (2)
             sfpy => dataptr(sflux(n,2),i)
             call update_scal_2d(nstart, nstop, &
                                 sop(:,:,1,:), ng_so, snp(:,:,1,:), ng_sn, &
                                 sfpx(:,:,1,:), sfpy(:,:,1,:), ng_sf, &
                                 fp(:,:,1,:), ng_f, &
                                 p0_new(n,:), lo, hi, dx(n,:), dt)
          case (3)
             sfpy => dataptr(sflux(n,2),i)
             sfpz => dataptr(sflux(n,3),i)
             if (spherical .eq. 0) then
                call update_scal_3d_cart(nstart, nstop, &
                                         sop(:,:,:,:), ng_so, snp(:,:,:,:), ng_sn, &
                                         sfpx(:,:,:,:), sfpy(:,:,:,:), sfpz(:,:,:,:), &
                                         ng_sf, fp(:,:,:,:), ng_f, &
                                         p0_new(n,:), lo, hi, dx(n,:), dt)
             else
                p0p => dataptr(p0_new_cart(n), i)
                call update_scal_3d_sphr(nstart, nstop, &
                                         sop(:,:,:,:), ng_so, snp(:,:,:,:), ng_sn, &
                                         sfpx(:,:,:,:), sfpy(:,:,:,:), sfpz(:,:,:,:), &
                                         ng_sf, fp(:,:,:,:), ng_f, &
                                         p0p(:,:,:,1), ng_p, lo, hi, dx(n,:), dt)
             end if
          end select
       end do

    end do

    ! restrict data and fill all ghost cells
    call ml_restrict_and_fill(nlevs,snew,mla%mba%rr,the_bc_level, &
                              icomp=nstart, &
                              bcomp=dm+nstart, &
                              nc=nstop-nstart+1, &
                              ng=snew(1)%ng)

    ! do the same for density if we updated the species
    if (nstart .eq. spec_comp .and. nstop .eq. (spec_comp+nspec-1)) then
       call ml_restrict_and_fill(nlevs,snew,mla%mba%rr,the_bc_level, &
                                 icomp=rho_comp, &
                                 bcomp=dm+rho_comp, &
                                 nc=1, &
                                 ng=snew(1)%ng)

    end if

    call destroy(bpt)

  end subroutine update_scal

  subroutine update_scal_1d(nstart,nstop,sold,ng_so,snew,ng_sn,sfluxx,ng_sf, &
                            force,ng_f,p0_new,lo,hi,dx,dt)

    use network,       only: nspec
    use eos_module,    only: eos, eos_input_rp
    use eos_type_module
    use probin_module, only: do_eos_h_above_cutoff, base_cutoff_density
    use variables,     only: spec_comp, rho_comp, rhoh_comp, temp_comp
    use pred_parameters
    use bl_constants_module

    integer           , intent(in   ) :: nstart, nstop, lo(:), hi(:)
    integer           , intent(in   ) :: ng_so, ng_sn, ng_sf, ng_f
    real (kind = dp_t), intent(in   ) ::   sold(lo(1)-ng_so:,:)
    real (kind = dp_t), intent(  out) ::   snew(lo(1)-ng_sn:,:)
    real (kind = dp_t), intent(in   ) :: sfluxx(lo(1)-ng_sf:,:)
    real (kind = dp_t), intent(in   ) ::  force(lo(1)-ng_f :,:)
    real (kind = dp_t), intent(in   ) :: p0_new(0:)
    real (kind = dp_t), intent(in   ) :: dt,dx(:)

    integer            :: i, comp, comp2
    real (kind = dp_t) :: divterm
    real (kind = dp_t) :: delta,frac,sumX
    logical            :: has_negative_species

    integer :: pt_index(MAX_SPACEDIM)

    type(eos_t) :: eos_state

    do comp = nstart, nstop

       do i=lo(1),hi(1)

          divterm = (sfluxx(i+1,comp) - sfluxx(i,comp))/dx(1)

          snew(i,comp) = sold(i,comp) + dt*(-divterm + force(i,comp))

       end do

    enddo

    if ( do_eos_h_above_cutoff .and. (nstart .eq. rhoh_comp) ) then

       do i = lo(1), hi(1)

          if (snew(i,rho_comp) .le. base_cutoff_density) then

             eos_state%rho = snew(i,rho_comp)
             eos_state%T   = sold(i,temp_comp)
             eos_state%p   = p0_new(i)
             eos_state%xn  =snew(i,spec_comp:spec_comp+nspec-1)/eos_state%rho

             pt_index(:) = (/i, -1, -1/)

             ! (rho,P) --> T,h
             call eos(eos_input_rp, eos_state, pt_index)

             snew(i,rhoh_comp) = snew(i,rho_comp) * eos_state%h

          end if

       enddo

    end if

    ! update density
    if (nstart .eq. spec_comp .and. nstop .eq. (spec_comp+nspec-1)) then

       snew(:,rho_comp) = sold(:,rho_comp)

       do i = lo(1), hi(1)

          has_negative_species = .false.

          ! define the update to rho as the sum of the updates to (rho X)_i
          do comp = nstart, nstop
             snew(i,rho_comp) = snew(i,rho_comp) + (snew(i,comp)-sold(i,comp))
             if (snew(i,comp) .lt. ZERO) has_negative_species = .true.
          enddo

          ! enforce a density floor
          if (snew(i,rho_comp) .lt. 0.5d0*base_cutoff_density) then
             do comp = nstart, nstop
                snew(i,comp) = snew(i,comp) * 0.5d0*base_cutoff_density/snew(i,rho_comp)
             end do
             snew(i,rho_comp) = 0.5d0*base_cutoff_density
          end if

          ! do not allow the species to leave here negative.
          if (has_negative_species) then
             do comp = nstart, nstop
                if (snew(i,comp) .lt. ZERO) then
                   delta = -snew(i,comp)
                   sumX = ZERO
                   do comp2 = nstart, nstop
                      if (comp2 .ne. comp .and. snew(i,comp2) .ge. ZERO) then
                         sumX = sumX + snew(i,comp2)
                      end if
                   enddo
                   do comp2 = nstart, nstop
                      if (comp2 .ne. comp .and. snew(i,comp2) .ge. ZERO) then
                         frac = snew(i,comp2) / sumX
                         snew(i,comp2) = snew(i,comp2) - frac * delta
                      end if
                   enddo
                   snew(i,comp) = ZERO
                end if
             end do
          end if

       enddo

    end if

  end subroutine update_scal_1d

  subroutine update_scal_2d(nstart,nstop,sold,ng_so,snew,ng_sn,sfluxx,sfluxy,ng_sf, &
                            force,ng_f,p0_new,lo,hi,dx,dt)

    use network,       only: nspec
    use eos_module,    only: eos, eos_input_rp
    use eos_type_module
    use probin_module, only: do_eos_h_above_cutoff, base_cutoff_density
    use variables,     only: spec_comp, rho_comp, rhoh_comp, temp_comp
    use pred_parameters
    use bl_constants_module

    integer           , intent(in   ) :: nstart, nstop, lo(:), hi(:)
    integer           , intent(in   ) :: ng_so, ng_sn, ng_sf, ng_f
    real (kind = dp_t), intent(in   ) ::   sold(lo(1)-ng_so:,lo(2)-ng_so:,:)
    real (kind = dp_t), intent(  out) ::   snew(lo(1)-ng_sn:,lo(2)-ng_sn:,:)
    real (kind = dp_t), intent(in   ) :: sfluxx(lo(1)-ng_sf:,lo(2)-ng_sf:,:)
    real (kind = dp_t), intent(in   ) :: sfluxy(lo(1)-ng_sf:,lo(2)-ng_sf:,:)
    real (kind = dp_t), intent(in   ) ::  force(lo(1)-ng_f :,lo(2)-ng_f :,:)
    real (kind = dp_t), intent(in   ) :: p0_new(0:)
    real (kind = dp_t), intent(in   ) :: dt,dx(:)

    integer            :: i, j, comp, comp2
    real (kind = dp_t) :: divterm
    real (kind = dp_t) :: delta,frac,sumX
    logical            :: has_negative_species

    integer :: pt_index(MAX_SPACEDIM)

    type(eos_t) :: eos_state

    do comp = nstart, nstop

       do j=lo(2),hi(2)
          do i=lo(1),hi(1)

             divterm = (sfluxx(i+1,j,comp) - sfluxx(i,j,comp))/dx(1) &
                     + (sfluxy(i,j+1,comp) - sfluxy(i,j,comp))/dx(2)

             snew(i,j,comp) = sold(i,j,comp) + dt*(-divterm + force(i,j,comp))

          end do
       end do

    enddo

    if ( do_eos_h_above_cutoff .and. (nstart .eq. rhoh_comp) ) then

       do j = lo(2), hi(2)
          do i = lo(1), hi(1)

             if (snew(i,j,rho_comp) .le. base_cutoff_density) then
                eos_state%rho = snew(i,j,rho_comp)
                eos_state%T   = sold(i,j,temp_comp)
                eos_state%p   = p0_new(j)
                eos_state%xn  = snew(i,j,spec_comp:spec_comp+nspec-1)/eos_state%rho

                pt_index(:) = (/i, j, -1/)

                ! (rho,P) --> T,h
                call eos(eos_input_rp, eos_state, pt_index)

                snew(i,j,rhoh_comp) = snew(i,j,rho_comp) * eos_state%h

             end if

          enddo
       enddo

    end if

    ! update density
    if (nstart .eq. spec_comp .and. nstop .eq. (spec_comp+nspec-1)) then

       snew(:,:,rho_comp) = sold(:,:,rho_comp)

       do j = lo(2), hi(2)
          do i = lo(1), hi(1)

             has_negative_species = .false.

             ! define the update to rho as the sum of the updates to (rho X)_i
             do comp = nstart, nstop
                snew(i,j,rho_comp) = snew(i,j,rho_comp) + (snew(i,j,comp)-sold(i,j,comp))
                if (snew(i,j,comp) .lt. ZERO) has_negative_species = .true.
             enddo

             ! enforce a density floor
             if (snew(i,j,rho_comp) .lt. 0.5d0*base_cutoff_density) then
                do comp = nstart, nstop
                   snew(i,j,comp) = snew(i,j,comp) * &
                        0.5d0*base_cutoff_density/snew(i,j,rho_comp)
                end do
                snew(i,j,rho_comp) = 0.5d0*base_cutoff_density
             end if

             ! do not allow the species to leave here negative.
             if (has_negative_species) then
                do comp = nstart, nstop
                   if (snew(i,j,comp) .lt. ZERO) then
                      delta = -snew(i,j,comp)
                      sumX = ZERO
                      do comp2 = nstart, nstop
                         if (comp2 .ne. comp .and. snew(i,j,comp2) .ge. ZERO) then
                            sumX = sumX + snew(i,j,comp2)
                         end if
                      enddo
                      do comp2 = nstart, nstop
                         if (comp2 .ne. comp .and. snew(i,j,comp2) .ge. ZERO) then
                            frac = snew(i,j,comp2) / sumX
                            snew(i,j,comp2) = snew(i,j,comp2) - frac * delta
                         end if
                      enddo
                      snew(i,j,comp) = ZERO
                   end if
                end do
             end if

          enddo
       enddo

    end if

  end subroutine update_scal_2d

  subroutine update_scal_3d_cart(nstart,nstop,sold,ng_so,snew,ng_sn,sfluxx,sfluxy,sfluxz, &
                                 ng_sf,force,ng_f,p0_new,lo,hi,dx,dt)

    use network,       only: nspec
    use eos_module,    only: eos, eos_input_rp
    use eos_type_module
    use probin_module, only: do_eos_h_above_cutoff, base_cutoff_density
    use variables,     only: spec_comp, rho_comp, rhoh_comp, temp_comp
    use pred_parameters
    use bl_constants_module

    integer           , intent(in   ) :: nstart, nstop, lo(:), hi(:)
    integer           , intent(in   ) :: ng_so, ng_sn, ng_sf, ng_f
    real (kind = dp_t), intent(in   ) ::   sold(lo(1)-ng_so:,lo(2)-ng_so:,lo(3)-ng_so:,:)
    real (kind = dp_t), intent(  out) ::   snew(lo(1)-ng_sn:,lo(2)-ng_sn:,lo(3)-ng_sn:,:)
    real (kind = dp_t), intent(in   ) :: sfluxx(lo(1)-ng_sf:,lo(2)-ng_sf:,lo(3)-ng_sf:,:)
    real (kind = dp_t), intent(in   ) :: sfluxy(lo(1)-ng_sf:,lo(2)-ng_sf:,lo(3)-ng_sf:,:)
    real (kind = dp_t), intent(in   ) :: sfluxz(lo(1)-ng_sf:,lo(2)-ng_sf:,lo(3)-ng_sf:,:)
    real (kind = dp_t), intent(in   ) ::  force(lo(1)-ng_f :,lo(2)-ng_f :,lo(3)-ng_f :,:)
    real (kind = dp_t), intent(in   ) :: p0_new(0:)
    real (kind = dp_t), intent(in   ) :: dt,dx(:)

    integer            :: i, j, k, comp, comp2
    real (kind = dp_t) :: divterm
    real (kind = dp_t) :: delta,frac,sumX
    logical            :: has_negative_species

    integer :: pt_index(MAX_SPACEDIM)

    type(eos_t) :: eos_state

    !$OMP PARALLEL PRIVATE(i,j,k,divterm,comp)
    do comp = nstart, nstop
       !$OMP DO
       do k = lo(3), hi(3)
          do j = lo(2), hi(2)
             do i = lo(1), hi(1)

                divterm = (sfluxx(i+1,j,k,comp) - sfluxx(i,j,k,comp))/dx(1) &
                        + (sfluxy(i,j+1,k,comp) - sfluxy(i,j,k,comp))/dx(2) &
                        + (sfluxz(i,j,k+1,comp) - sfluxz(i,j,k,comp))/dx(3)

                snew(i,j,k,comp) = sold(i,j,k,comp) + dt * (-divterm + force(i,j,k,comp))

             enddo
          enddo
       enddo
       !$OMP END DO NOWAIT
    end do
    !$OMP END PARALLEL

    if ( do_eos_h_above_cutoff .and. (nstart .eq. rhoh_comp) ) then

       !$OMP PARALLEL DO PRIVATE(i,j,k,eos_state,pt_index)
       do k = lo(3), hi(3)
          do j = lo(2), hi(2)
             do i = lo(1), hi(1)

                if (snew(i,j,k,rho_comp) .le. base_cutoff_density) then

                   eos_state%rho = snew(i,j,k,rho_comp)
                   eos_state%T   = sold(i,j,k,temp_comp)
                   eos_state%p   = p0_new(k)
                   eos_state%xn  = snew(i,j,k,spec_comp:spec_comp+nspec-1)/eos_state%rho

                   pt_index(:) = (/i, j, k/)

                   ! (rho,P) --> T,h
                   call eos(eos_input_rp, eos_state, pt_index)

                   snew(i,j,k,rhoh_comp) = snew(i,j,k,rho_comp) * eos_state%h

                end if

             enddo
          enddo
       enddo
       !$OMP END PARALLEL DO

    end if

    ! update density
    if (nstart .eq. spec_comp .and. nstop .eq. (spec_comp+nspec-1)) then

       snew(:,:,:,rho_comp) = sold(:,:,:,rho_comp)

       !$OMP PARALLEL DO PRIVATE(i,j,k,has_negative_species,comp,delta,sumX,comp2,frac)
       do k = lo(3), hi(3)
          do j = lo(2), hi(2)
             do i = lo(1), hi(1)

                has_negative_species = .false.

                ! define the update to rho as the sum of the updates to (rho X)_i
                do comp = nstart, nstop
                   snew(i,j,k,rho_comp) = snew(i,j,k,rho_comp) &
                        + (snew(i,j,k,comp)-sold(i,j,k,comp))
                   if (snew(i,j,k,comp) .lt. ZERO) has_negative_species = .true.
                enddo

                ! enforce a density floor
                if (snew(i,j,k,rho_comp) .lt. 0.5d0*base_cutoff_density) then
                   do comp = nstart, nstop
                      snew(i,j,k,comp) = snew(i,j,k,comp) * &
                           0.5d0*base_cutoff_density/snew(i,j,k,rho_comp)
                   end do
                   snew(i,j,k,rho_comp) = 0.5d0*base_cutoff_density
                end if

                ! do not allow the species to leave here negative.
                if (has_negative_species) then
                   do comp = nstart, nstop
                      if (snew(i,j,k,comp) .lt. ZERO) then
                         delta = -snew(i,j,k,comp)
                         sumX = ZERO
                         do comp2 = nstart, nstop
                            if (comp2 .ne. comp .and. snew(i,j,k,comp2) .ge. ZERO) then
                               sumX = sumX + snew(i,j,k,comp2)
                            end if
                         enddo
                         do comp2 = nstart, nstop
                            if (comp2 .ne. comp .and. snew(i,j,k,comp2) .ge. ZERO) then
                               frac = snew(i,j,k,comp2) / sumX
                               snew(i,j,k,comp2) = snew(i,j,k,comp2) - frac * delta
                            end if
                         enddo
                         snew(i,j,k,comp) = ZERO
                      end if
                   end do
                end if

             enddo
          enddo
       enddo
       !$OMP END PARALLEL DO

    end if

  end subroutine update_scal_3d_cart

  subroutine update_scal_3d_sphr(nstart,nstop,sold,ng_so,snew,ng_sn,sfluxx,sfluxy,sfluxz, &
                                 ng_sf,force,ng_f,p0_new_cart,ng_p,lo,hi,dx,dt)

    use network,       only: nspec
    use eos_module,    only: eos, eos_input_rp
    use eos_type_module
    use probin_module, only: do_eos_h_above_cutoff, base_cutoff_density
    use variables,     only: spec_comp, rho_comp, rhoh_comp, temp_comp
    use pred_parameters
    use bl_constants_module

    integer           , intent(in   ) :: nstart, nstop, lo(:), hi(:)
    integer           , intent(in   ) :: ng_so, ng_sn, ng_sf, ng_f, ng_p
    real (kind = dp_t), intent(in   ) ::    sold(lo(1)-ng_so:,lo(2)-ng_so:,lo(3)-ng_so:,:)
    real (kind = dp_t), intent(  out) ::    snew(lo(1)-ng_sn:,lo(2)-ng_sn:,lo(3)-ng_sn:,:)
    real (kind = dp_t), intent(in   ) ::  sfluxx(lo(1)-ng_sf:,lo(2)-ng_sf:,lo(3)-ng_sf:,:)
    real (kind = dp_t), intent(in   ) ::  sfluxy(lo(1)-ng_sf:,lo(2)-ng_sf:,lo(3)-ng_sf:,:)
    real (kind = dp_t), intent(in   ) ::  sfluxz(lo(1)-ng_sf:,lo(2)-ng_sf:,lo(3)-ng_sf:,:)
    real (kind = dp_t), intent(in   ) ::   force(lo(1)-ng_f :,lo(2)-ng_f :,lo(3)-ng_f :,:)
    real (kind = dp_t), intent(in   )::   p0_new_cart(lo(1)-ng_p :,lo(2)-ng_p :,lo(3)-ng_p :)
    real (kind = dp_t), intent(in   ) :: dt,dx(:)

    integer            :: i, j, k, comp, comp2
    real (kind = dp_t) :: divterm
    real (kind = dp_t) :: delta,frac,sumX
    logical            :: has_negative_species

    integer :: pt_index(MAX_SPACEDIM)
    type(eos_t) :: eos_state

    !
    ! is spherical
    !
    !$OMP PARALLEL PRIVATE(i,j,k,divterm,comp)
    do comp = nstart, nstop
       !$OMP DO
       do k = lo(3), hi(3)
          do j = lo(2), hi(2)
             do i = lo(1), hi(1)

                divterm = (sfluxx(i+1,j,k,comp) - sfluxx(i,j,k,comp))/dx(1) &
                        + (sfluxy(i,j+1,k,comp) - sfluxy(i,j,k,comp))/dx(2) &
                        + (sfluxz(i,j,k+1,comp) - sfluxz(i,j,k,comp))/dx(3)

                snew(i,j,k,comp) = sold(i,j,k,comp) + dt * (-divterm + force(i,j,k,comp))

             enddo
          enddo
       enddo
       !$OMP END DO NOWAIT
    end do
    !$OMP END PARALLEL

    if ( do_eos_h_above_cutoff .and. (nstart .eq. rhoh_comp) ) then

       !$OMP PARALLEL DO PRIVATE(i,j,k,eos_state,pt_index)
       do k = lo(3), hi(3)
          do j = lo(2), hi(2)
             do i = lo(1), hi(1)

                if (snew(i,j,k,rho_comp) .le. base_cutoff_density) then
                   eos_state%rho   = snew(i,j,k,rho_comp)
                   eos_state%T     = sold(i,j,k,temp_comp)
                   eos_state%p     = p0_new_cart(i,j,k)
                   eos_state%xn(:) = snew(i,j,k,spec_comp:spec_comp+nspec-1)/eos_state%rho

                   pt_index(:) = (/i, j, k/)

                   ! (rho,P) --> T,h
                   call eos(eos_input_rp, eos_state, pt_index)

                   snew(i,j,k,rhoh_comp) = snew(i,j,k,rho_comp) * eos_state%h

                end if

             enddo
          enddo
       enddo
       !$OMP END PARALLEL DO

    end if

    ! update density
    if (nstart .eq. spec_comp .and. nstop .eq. (spec_comp+nspec-1)) then

       snew(:,:,:,rho_comp) = sold(:,:,:,rho_comp)

       !$OMP PARALLEL DO PRIVATE(i,j,k,has_negative_species,comp,delta,sumX,comp2,frac)
       do k = lo(3), hi(3)
          do j = lo(2), hi(2)
             do i = lo(1), hi(1)

                has_negative_species = .false.

                ! define the update to rho as the sum of the updates to (rho X)_i
                do comp = nstart, nstop
                   snew(i,j,k,rho_comp) = snew(i,j,k,rho_comp) &
                        + (snew(i,j,k,comp)-sold(i,j,k,comp))
                   if (snew(i,j,k,comp) .lt. ZERO) has_negative_species = .true.
                enddo

                ! enforce a density floor
                if (snew(i,j,k,rho_comp) .lt. 0.5d0*base_cutoff_density) then
                   do comp = nstart, nstop
                      snew(i,j,k,comp) = snew(i,j,k,comp) * &
                           0.5d0*base_cutoff_density/snew(i,j,k,rho_comp)
                   end do
                   snew(i,j,k,rho_comp) = 0.5d0*base_cutoff_density
                end if

                ! do not allow the species to leave here negative.
                if (has_negative_species) then
                   do comp = nstart, nstop
                      if (snew(i,j,k,comp) .lt. ZERO) then
                         delta = -snew(i,j,k,comp)
                         sumX = ZERO
                         do comp2 = nstart, nstop
                            if (comp2 .ne. comp .and. snew(i,j,k,comp2) .ge. ZERO) then
                               sumX = sumX + snew(i,j,k,comp2)
                            end if
                         enddo
                         do comp2 = nstart, nstop
                            if (comp2 .ne. comp .and. snew(i,j,k,comp2) .ge. ZERO) then
                               frac = snew(i,j,k,comp2) / sumX
                               snew(i,j,k,comp2) = snew(i,j,k,comp2) - frac * delta
                            end if
                         enddo
                         snew(i,j,k,comp) = ZERO
                      end if
                   end do
                end if

             enddo
          enddo
       enddo
       !$OMP END PARALLEL DO
    end if

  end subroutine update_scal_3d_sphr

end module update_scal_module

max_stars_repo_name: sailoridy/MAESTRO

max_issues_repo_name: sailoridy/MAESTRO

max_forks_repo_name: sailoridy/MAESTRO


text:
      SUBROUTINE free_mem_funct3d
      USE vmec_main
      USE realspace
      USE vforces
      USE vacmod
      IMPLICIT NONE
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER :: istat1 = 0
C-----------------------------------------------

      IF (ALLOCATED(armn))
     1   DEALLOCATE (armn, azmn, brmn, bzmn, crmn, czmn, blmn, clmn,
     2   r1, ru, rv, z1, zu, zv, gcon, rcon, zcon, ru0, zu0,
     3   rcon0, zcon0, guu, guv, gvv, stat=istat1)
      IF (istat1 .ne. 0) STOP 'deallocation error#1 in funct3d'

      IF (ALLOCATED(brv))
     1   DEALLOCATE (brv, bphiv, bzv, bsqvac, stat=istat1)
      IF (istat1 .ne. 0) STOP 'deallocation error#2 in funct3d'

      IF (ALLOCATED(extra1))
     1   DEALLOCATE (extra1, extra2, extra3, extra4, stat=istat1)
      IF (istat1 .ne. 0) STOP 'deallocation error#3 in funct3d'

      END SUBROUTINE free_mem_funct3d

ext: f


text:
Subroutine FbarSpatial(AE, BE, MaxLRows, XLL, NDim, iDofT, NodElt, Sol0, Sol1, CommonPar, Param, JParam, DelT, DTm, Time)
!     ------------------------------------------------------------------
	use funcAux
	use loadingLib
	use finiteStrainLib
	use ptsGaussLib2

	implicit none
!~
!~ 	!   ===== SUBROUTINE ARGUMENTS  =======
	integer :: MaxLRows,Ndim,iDofT, NodElt! all integers
	Real*8 :: DelT, DTm,Time ! all reals
	Real*8 :: AE(MaxLRows,MaxLRows), BE(MaxLRows), XLL(*), Sol0(*), Sol1(*), Param(*), JParam(*), CommonPar(*)

!~ 	!   =====   END ARGUMENTS  =======

	Integer :: nG, i,j,k,l,m,n,p,q, A, B, ApRow, BpCol , App,Bqp, ipDim, ipDimE, ipDof,ipDofmod ! counters
	Integer :: constLaw
	integer , parameter ::  NdimE = 2! this code only makes sense in 2D and for triangles
	Real*8 :: PhiA,PhiB, dV , Det , bm(NdimE), MatPar(9)
	Real*8 :: dPhiA(NdimE), dPhiB(NdimE), dPhiB0(NdimE), detFs0, detFs
	Real*8 , dimension(NdimE,NdimE) :: GradUs , GradUs0, DGradUs, sigma, Fs, Fs0, Fbar
	Real*8 , dimension(NdimE,NdimE,NdimE,NdimE) :: Ds , tenQs
	real*8 , allocatable ::  SolU(:) , Xel(:) !! all have dimension NodG*NdimE
	integer :: pOrder , NGP, NodG, iFEMtype , iSimplex, iBubble
	Real*8 :: LoadPar(6) , he, energy
	integer :: iShiftU, iShiftDeltaU , LoadType , LoadTypeProg
	type(ptGaussClass) :: PtG, PtG0

	iShiftU = nint(CommonPar(1))
	iShiftDeltaU = nint(CommonPar(2))
	iFEMtype =  nint(commonPar(3))
	constLaw = nint(CommonPar(4))
	matPar(1:3) = commonPar(5:7)
	LoadType = nint(commonPar(8))
	LoadTypeProg = nint(commonPar(9))
	LoadPar(:) = CommonPar(10:15)

	call setFEMtype(iFEMtype,NodG,pOrder,NGP,iSimplex,iBubble)

	call getSliceAllocate(SolU,Sol1,1,NodG,iShiftU + 1 ,iShiftU + NdimE, iDofT)
	call getSliceAllocate(Xel,XLL,1,NodG,1 ,NdimE, Ndim)

	Xel = Xel + SolU

	call PtG%init(Xel,NodG,NdimE,NGP,pOrder,iBubble, iSimplex)

	nG = 1
	call PtG0%init(Xel,NodG,NdimE,1,pOrder,iBubble,iSimplex)
	call PtG0%calcGradU(GradUs0,SolU,nG)

	call setF(Fs0,LoadPar,Time,DelT,LoadTypeProg,LoadType)

	call computeFs(Fs0,detFs0,gradUs0,NdimE)

	Do nG = 1, NGP ! LoopGauss

		GradUs = 0.0d0

		call PtG%calcGradU(GradUs,SolU,nG)

		call setF(Fs,LoadPar,Time,DelT,LoadTypeProg,LoadType)

		call computeFs(Fs,detFs,gradUs,NdimE)

		Fbar = ((detFs0/detFs)**(1.0d0/2.0d0))*Fs !! for 2D

		call calcPpk(sigma,Fbar,NdimE,matPar,constLaw)
		call calcD(Ds,Fbar,NdimE,matPar,constLaw)
		call modifyToSpatial(Ds,sigma,Fbar,detFs0,NdimE)

		call build_tenQs(Ds,sigma,tenQs,NdimE)

		dV=ptG%dV(nG)

		Do A=1, NodG !LoopRow
			ApRow  = (A-1)*iDofT + iShiftDeltaU
			PhiA = PtG%Phi(A,nG)
			dPhiA= PtG%dPhi_G(:,A,nG)

			Do p = 1 , NdimE
				App = ApRow+p

				BE(App) = BE(App) - dot_product(sigma(p,:),dPhiA)*DV

				do B= 1 , NodG ! LoopCol ! not considering the simmetry
					BpCol  = (B-1)*iDofT + iShiftDeltaU
					PhiB  = PtG%Phi(B,nG)
					dPhiB = PtG%dPhi_G(:,B,nG)
					dPhiB0 = PtG0%dPhi_G(:,B,1)

					do q=1,NdimE
						Bqp = BpCol+q

						AE(App,Bqp) = AE(App,Bqp) + DGradPhiBdotGradPhiA(Ds,dPhiB,dPhiA,p,q,NdimE)*DV
						AE(App,Bqp) = AE(App,Bqp) + DGradPhiBdotGradPhiA(tenQs,dPhiB0 - dPhiB,dPhiA,p,q,NdimE)*DV
!~ 						write(0,*) AE(App,Bqp)
					end do ! loop Bq
				Enddo !LoopCol
			end do ! loop Ap
		Enddo !LoopRow
	EndDo !LoopGauss

 	if(iShiftU /= iShiftDeltaU) then
		Do A=1, NodG !LoopRow
			ApRow  = (A-1)*iDofT + iShiftU

			Do p = 1 , NdimE
				App = ApRow+p

				BE(App) = SolU((A-1)*NdimE + p)

				AE(App,App) = 1.0d0
				AE(App,App + iShiftDeltaU - iShiftU) = -1.0d0

			end do ! loop Ap
		end do !LoopRow
	end if

end subroutine

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Subroutine FbarSpatialS(Coupling,CommonPar,iDofT,Ndim,MaxLRows,iAdd)
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    use funcAux
    use ptsGaussLib2
    IMPLICIT NONE
    Integer :: MaxLRows,Ndim,iDofT, iAdd
    Integer   Coupling(MaxLRows,MaxLRows)
    Real*8 CommonPar(*)

    Integer :: A,B,q,p, ApRow, BpCol,iShiftU , iShiftDeltaU, Femtype, NodG
    integer , parameter :: NdimE = 2

    Coupling = 0

    iShiftU = nint(commonPar(1))
    iShiftDeltaU = nint(commonPar(2))
    femType = nint(commonPar(3))

    call setNodG(femtype, NodG)

    do A = 1,NodG
        ApRow = (A-1) * iDofT + iShiftDeltaU

		do B = 1,NodG
			BpCol = (B-1) * iDofT + iShiftDeltaU

			do p=1,nDimE
			do q=1,nDimE

				Coupling (ApRow + p, BpCol + q) = 1

			enddo
            enddo
		enddo
	enddo

 	if(iShiftU /= iShiftDeltaU) then
		do A = 1,NodG
			ApRow = (A-1) * iDofT + iShiftU
			do p=1,nDimE
				Coupling (ApRow + p, ApRow + p) = 1
				Coupling (ApRow + p, ApRow + p + iShiftDeltaU - iShiftU) = 1
			enddo
		enddo
	end if

	call numPrint(Coupling)

end Subroutine

subroutine computeFs(Fs,detFs,gradUs,NdimE)
	use funcAux
	implicit none

	real*8 , intent(inout) :: Fs(NdimE,NdimE)
	real*8, intent(out) :: detFs
	real*8, intent(in) :: gradUs(NdimE,NdimE)
	integer , intent(in) :: NdimE

	real*8 :: Fsinv(NdimE,NdimE), detFdummy

	Fs = Fs + deltaKron(1:NdimE,1:NdimE) !! Fs can be prescribed

	call matInv(FsInv,detFdummy,Fs)
	FsInv = Fsinv - gradUs

	call matInv(Fs,detFs,FsInv)
	detFs = 1.0d0/detFs

end subroutine

max_stars_repo_name: felipefr/Piola

max_issues_repo_name: felipefr/Piola

max_forks_repo_name: felipefr/Piola


text:
    !> Implements minimum restriction BC for continuum
    !!
    !!
	!! @param iShiftFluc = nint(CommonPar(1))
	!! @param iShiftLag = nint(CommonPar(2))
	!! @param pen  = CommonPar(3)
	!! @param eps  = commonPar(4)
	!!
    !! @author Rocha, Felipe Figueredo
!     ------------------------------------------------------------------
Subroutine minRestrictionRVE(AE, BE, MaxLRows, XLL, NDim, iDofT, NodELT, Sol0, &
						Sol1, CommonPar, Param, JParam, DelT, DTm, Time)
!     ------------------------------------------------------------------
	use funcAux

	IMPLICIT NONE

	!   ===== SUBROUTINE ARGUMENTS  =======
	integer :: MaxLRows,Ndim,iDofT,NodELT ! all integers
	Real*8 :: DelT, DTm,Time ! all reals
	! Reals Vectors and matrices
	Real*8  :: AE(MaxLRows,MaxLRows), BE(MaxLRows), XLL(*), Sol0(*), Sol1(*), &
		Param(*), JParam(*), CommonPar(*)

	!   =====   END ARGUMENTS  =======

	Integer :: i,j,k,l,m,n,p,q, A, B, ApRow, BpCol , App,Bqp! counters
	Integer :: NdofLag
	integer, parameter :: Ndim2D = 2, NodG = 2, NodLag = 3
	Real*8 :: eps, pen, lengthLine
	Real*8 :: normal0(Ndim2D) , vecU(Ndim2D)
	integer :: iShiftFluc, iShiftLag

	iShiftFluc = nint(CommonPar(1))
	iShiftLag = nint(CommonPar(2))
	pen = commonPar(3)
	eps = commonPar(4)

	NdofLag = Ndim2D*Ndim2D

	do i = 1 , Ndim2D
		vecU(i) = XLL(Ndim + i) - XLL(i)
	end do

	lengthLine = dsqrt(dot_product(vecU,vecU))
	normal0(1) = vecU(2)/lengthLine
	normal0(2) = -vecU(1)/lengthLine

	AE = 0.0d0
	BE = 0.0d0

	B = NodLag
	BpCol = (B-1)*iDofT + iShiftLag

	A = 1
	ApRow  = (A-1)*iDofT + iShiftFluc
	AE(BpCol + 1,ApRow + 1) = pen*0.5d0*lengthLine*normal0(1) ! l11
	AE(BpCol + 2,ApRow + 2) = pen*0.5d0*lengthLine*normal0(2) ! l12
	AE(BpCol + 3,ApRow + 1) = pen*0.5d0*lengthLine*normal0(1) ! l21
	AE(BpCol + 4,ApRow + 2) = pen*0.5d0*lengthLine*normal0(2) ! l22
	AE(BpCol + 1,ApRow + 1) = AE(ApRow + 1,BpCol + 1)
	AE(BpCol + 2,ApRow + 2) = AE(ApRow + 2,BpCol + 2)
	AE(BpCol + 3,ApRow + 1) = AE(ApRow + 1,BpCol + 3)
	AE(BpCol + 4,ApRow + 2) = AE(ApRow + 2,BpCol + 4)

	A = 2
	ApRow  = (A-1)*iDofT + iShiftFluc
	AE(BpCol + 1,ApRow + 1) = pen*0.5d0*lengthLine*normal0(1) ! l11
	AE(BpCol + 2,ApRow + 2) = pen*0.5d0*lengthLine*normal0(2) ! l12
	AE(BpCol + 3,ApRow + 1) = pen*0.5d0*lengthLine*normal0(1) ! l21
	AE(BpCol + 4,ApRow + 2) = pen*0.5d0*lengthLine*normal0(2) ! l22
	AE(BpCol + 1,ApRow + 1) = AE(ApRow + 1,BpCol + 1)
	AE(BpCol + 2,ApRow + 2) = AE(ApRow + 2,BpCol + 2)
	AE(BpCol + 3,ApRow + 1) = AE(ApRow + 1,BpCol + 3)
	AE(BpCol + 4,ApRow + 2) = AE(ApRow + 2,BpCol + 4)

	AE(BpCol + 1,BpCol + 1) = eps
	AE(BpCol + 2,BpCol + 2) = eps
	AE(BpCol + 3,BpCol + 3) = eps
	AE(BpCol + 4,BpCol + 4) = eps

	BE(BpCol + 1) = eps*Sol1(BpCol + 1)
	BE(BpCol + 2) = eps*Sol1(BpCol + 2)
	BE(BpCol + 3) = eps*Sol1(BpCol + 3)
	BE(BpCol + 4) = eps*Sol1(BpCol + 4)

end subroutine

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Subroutine minRestrictionRVES(Coupling,CommonPar,iDofT,Ndim,MaxLRows,iAdd)
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	use funcAux
	IMPLICIT NONE
	Integer :: MaxLRows,Ndim,iDofT, iAdd
	Integer   Coupling(MaxLRows,MaxLRows)
	Real*8 CommonPar(*)
	integer ::  p , q, App, ApRow, BpCol, Bqp, A,B, iShiftFluc, iShiftLag , NdofLag
	integer , parameter :: NodG = 2, NodLag = 3, Ndim2D = 2

	iShiftFluc = nint(CommonPar(1))
	iShiftLag = nint(CommonPar(2))

	NdofLag = Ndim2D*Ndim2D

	B = NodLag
	BpCol = (B-1)*iDofT + iShiftLag

	A = 1
	ApRow  = (A-1)*iDofT + iShiftFluc
	Coupling(BpCol + 1,ApRow + 1) = 1 ! l11
	Coupling(BpCol + 2,ApRow + 2) = 1 ! l12
	Coupling(BpCol + 3,ApRow + 1) = 1 ! l21
	Coupling(BpCol + 4,ApRow + 2) = 1 ! l22
	Coupling(BpCol + 1,ApRow + 1) = 1
	Coupling(BpCol + 2,ApRow + 2) = 1
	Coupling(BpCol + 3,ApRow + 1) = 1
	Coupling(BpCol + 4,ApRow + 2) = 1

	A = 2
	ApRow  = (A-1)*iDofT + iShiftFluc
	Coupling(BpCol + 1,ApRow + 1) = 1 ! l11
	Coupling(BpCol + 2,ApRow + 2) = 1 ! l12
	Coupling(BpCol + 3,ApRow + 1) = 1 ! l21
	Coupling(BpCol + 4,ApRow + 2) = 1 ! l22
	Coupling(BpCol + 1,ApRow + 1) = 1
	Coupling(BpCol + 2,ApRow + 2) = 1
	Coupling(BpCol + 3,ApRow + 1) = 1
	Coupling(BpCol + 4,ApRow + 2) = 1

	Coupling(BpCol + 1,BpCol + 1) = 1
	Coupling(BpCol + 2,BpCol + 2) = 1
	Coupling(BpCol + 3,BpCol + 3) = 1
	Coupling(BpCol + 4,BpCol + 4) = 1

end subroutine

max_stars_repo_name: felipefr/gpmaterials

max_issues_repo_name: felipefr/gpmaterials

max_forks_repo_name: felipefr/gpmaterials


text:
!===============================================================================
!> A simple task dispatcher, which handles a list of tasks that may be started
!> as background OpenMP tasks.  A task can be in any of 4 states:
!>
!> 1: ready for updating
!> 2: undergoing updating by some thread
!> 3: has been updated, but nbors not yet checked
!> 4: finished
!>
!> When starting, all tasks are normally in state 1, and the dispatcher starts
!> a number (omp%nthreads) of background tasks that update them, first changing
!> the state to 2.  When finished, the updating thread changes the state from 2
!> to 3. When the dispatcher finds a task in state 3, it runs a check on the
!> nbor tasks
!>
!> dispatcher%execute
!>   dispatcher%update
!>     link%check_task
!>       dispatcher%activate
!>         link%check_mpi
!>         link%update
!===============================================================================
MODULE dispatcher2_mod
  USE io_mod
  USE io_unit_mod
  USE omp_mod
  USE omp_timer_mod
  USE timer_mod
  USE trace_mod
  USE mpi_mesg_mod
  USE link_mod
  USE list_mod
  USE task_list_mod
  USE bits_mod
  USE task_mod
  USE experiment_mod
  implicit none
  private
  type, public :: dispatcher2_t
    integer   :: running
    type(task_list_t), pointer:: task_list => null()
  contains
    procedure :: execute
    procedure :: activate
    procedure :: update
  end type
  type(dispatcher2_t), public:: dispatcher2
CONTAINS

!===============================================================================
!> Initialize the list, by populating it from a task_list
!===============================================================================
SUBROUTINE execute (self, task_list, test)
  class(dispatcher2_t)    :: self
  type(task_list_t), pointer :: task_list
  logical:: test
  !.............................................................................
  class(link_t), pointer :: link
  real(8):: start
  !-----------------------------------------------------------------------------
  call trace_begin ('dispatcher2_t%execute')
  call task_list%startup
  self%running = 0
  self%task_list => task_list
  task_list%nq = task_list%n
  task_list%na = task_list%n
  link => task_list%head
  do while (associated(link))
    call link%task%set (bits%ready)
    link => link%next
  end do
  call tic (time=start)
  !$omp parallel
  !$omp master
  do while (self%task_list%na > 0 .and. wallclock() < io%job_seconds)
    call self%update (test)
  end do
  !$omp end master
  !$omp end parallel
  call trace_end
  call timer%print()
  call mpi_mesg%diagnostics(1)
  call toc ('wall time', timer%n_update, time=start)
END SUBROUTINE execute

!===============================================================================
!> Activate a task pointed to by a link, by starting a task%update as a
!> background OpenMP task.  Set the state of task to 1 when dormant, 2 when
!> ready, and 3 when active.
!===============================================================================
SUBROUTINE activate (self, link, test)
  class(dispatcher2_t)    :: self
  class(link_t), pointer :: link
  logical:: test
  !-----------------------------------------------------------------------------
  if (omp%nthreads==1) then
    link%task%state = 3
    call self%task_list%update (link, test)
    link%task%state = 1
  else
    !$omp task default(shared) firstprivate(link)
    !$omp atomic write
    link%task%state = 3
    call self%task_list%check_mpi ()
    call self%task_list%update (link, test)
    !$omp atomic write
    link%task%state = 1
    !$omp atomic
    self%task_list%nq = self%task_list%nq - 1
    !$omp end task
  end if
END SUBROUTINE activate

!===============================================================================
!> Update the list of running tasks.  For each task on the running list, if it
!> is no longer busy, check if any of its passive nbors has become executable.
!> Finally, check the task itself.
!===============================================================================
SUBROUTINE update (self, test)
  class(dispatcher2_t)    :: self
  logical                 :: test
  !.............................................................................
  class(link_t), pointer  :: link, nbor
  integer                 :: state
  logical                 :: ok
  integer, save           :: itimer=0
  !-----------------------------------------------------------------------------
  call trace_begin ('dispatcher2_t%update', itimer=itimer)
  link => self%task_list%head
  do while (associated(link))
    if (io%debug(2)) then
      print *, 'checking link', link%task%id, link%task%state
    end if
    !$omp atomic read
    state = link%task%state
    if (state==4) then
      link => link%next
      cycle
    end if
    if (link%task%time > io%end_time) then
      !$omp atomic
      self%task_list%na = self%task_list%na - 1
      !$omp atomic write
      link%task%state = 4
      cycle
    end if
    if (state<=1) then
      nbor => link%nbor
      do while (associated(nbor))
        !$omp atomic read
        state = nbor%link%task%state
        if (state==1) then
          if (io%debug(3)) then
            print '(i3,f12.6,a,i6,i3)', omp%thread, wallclock(), '       checking nbor', nbor%link%task%id, state
          end if
          call check_task (nbor%link, ok)
        end if
        nbor => nbor%next
      end do
      call check_task (link, ok)
    end if
    link => link%next
  end do
  call trace_end (itimer)
contains
  !-----------------------------------------------------------------------------
  ! If all nbors of a task are ahead in time activate the task.  Only the master
  ! thread is checking, so there should be no chance of starting multiple
  ! threads on the same update. The only state change performed by the many
  ! production threads is to go from state 2 to state 3, which is of no concern
  ! to the master thread, and then from state 3 to state 1, which will
  ! eventually be detected (typically much later) by the master thread. The
  ! delay before checking is not a disadvantage, since it actually increases
  ! the chances that the patch is ready for updating when it checked. At the end
  ! a change to state 4 indicates a finished task.
  !-----------------------------------------------------------------------------
  subroutine check_task (link, ok)
    class(link_t), pointer:: link, nbor
    class(task_t), pointer:: task
    logical:: ok
    integer:: state
    !---------------------------------------------------------------------------
    ok = .true.
    task => link%task
    nbor => link%nbor
    do while (ok.and.associated (nbor))
      ok = ok .and. nbor%task%is_ahead_of(task)
      nbor => nbor%next
    end do
    if (ok) then
      !$omp atomic read
      state = link%task%state
      if (state<=1) then
        if (state==1) then
          !$omp atomic
          self%task_list%nq = self%task_list%nq + 1
        end if
        !$omp atomic write
        link%task%state = 2
        call self%activate (link, test)
      else
         print '(i3,f12.6,a,i6,i3)', omp%thread, wallclock(), ' found unexpected state', link%task%id, link%task%state
      end if
    end if
  end subroutine check_task
END SUBROUTINE update

END MODULE dispatcher2_mod

max_stars_repo_name: applejwjcat/dispatch

max_issues_repo_name: applejwjcat/dispatch

max_forks_repo_name: applejwjcat/dispatch


!------------------------------------------------------------------------------!
!> Driver of the solution part of PB3D.
!------------------------------------------------------------------------------!
module driver_sol
#include <PB3D_macros.h>
#include <wrappers.h>
    use str_utilities
    use output_ops
    use messages
    use num_vars, only: dp, pi, max_str_ln, iu
    use grid_vars, only: grid_type
    use X_vars, only: X_2_type, modes_type, &
        &mds_X, mds_sol
    use vac_vars, only: vac_type
    use sol_vars, only: sol_type

    implicit none
    private
    public run_driver_sol

    ! global variables
#if ldebug
    logical :: debug_run_driver_sol = .false.                                   !< debug information for run_driver_sol \ldebug
    logical :: debug_interp_V = .false.                                         !< debug information for interp_v \ldebug
    integer :: ivs_stats(3)                                                     !< statistics of interp_v_spline (copy, linear, spline)
#endif

contains
    !> Main driver of PB3D solution part.
    !!
    !!  - sets up:
    !!      * \c grid_sol (only first Richardson level)
    !!      * \c sol
    !!  - writes to HDF5:
    !!      * \c grid_sol (only first Richardson level)
    !!      * \c sol
    !!  - deallocates:
    !!      * sol before setting up (but after guess)
    !!
    !! \return ierr
    integer function run_driver_sol(grid_eq,grid_X,grid_sol,X,vac,sol) &
        &result(ierr)

        use num_vars, only: EV_style, eq_style, rich_restart_lvl, rank, &
            &n_procs, X_grid_style, jump_to_sol
        use grid_vars, only: n_r_sol, n_alpha
        use PB3D_ops, only: reconstruct_PB3D_grid, reconstruct_PB3D_sol
        use SLEPC_ops, only: solve_EV_system_SLEPC
        use grid_ops, only: calc_norm_range, setup_grid_sol, &
            &print_output_grid, redistribute_output_grid
        use sol_ops, only: print_output_sol
        use rich_vars, only: rich_lvl
        use rich_ops, only: calc_rich_ex
        use vac_ops, only: calc_vac_res, print_output_vac
        use grid_utilities, only: trim_grid
        use X_ops, only: setup_modes, redistribute_output_X
#if ldebug
        use X_ops, only: print_debug_X_2
#endif

        character(*), parameter :: rout_name = 'run_driver_sol'

        ! input / output
        type(grid_type), intent(in) :: grid_eq                                  !< equilibrium grid
        type(grid_type), intent(in), target :: grid_X                           !< perturbation grid
        type(grid_type), intent(inout) :: grid_sol                              !< solution grid
        type(X_2_type), intent(in) :: X                                         !< integrated tensorial perturbation variables
        type(vac_type), intent(inout) :: vac                                    !< vacuum variables
        type(sol_type), intent(inout) :: sol                                    !< solution variables

        ! local variables
        type(grid_type) :: grid_sol_trim                                        ! trimmed solution grid
        type(grid_type) :: grid_X_rdst                                          ! redistributed perturbation grid
        type(grid_type) :: grid_X_sol                                           ! perturbation grid with solution normal part
        type(X_2_type) :: X_rdst                                                ! redistributed X
        type(X_2_type) :: X_sol                                                 ! interpolated X
        integer :: sol_limits(2)                                                ! min. and max. index of sol grid for this process
        integer :: rich_lvl_name                                                ! either the Richardson level or zero, to append to names
        real(dp), allocatable :: r_F_sol(:)                                     ! normal points in solution grid
        logical :: do_vac_ops                                                   ! whether specific calculations for vacuum are necessary
        character(len=max_str_ln) :: err_msg                                    ! error message

        ! initialize ierr
        ierr = 0

        ! set up solution grid if first level
        if (rich_lvl.eq.rich_restart_lvl) then
            ! Divide solution grid under group processes, calculating the limits
            ! and the normal coordinate.
            allocate(r_F_sol(n_r_sol))
            ierr = calc_norm_range('PB3D_sol',sol_limits=sol_limits,&
                &r_F_sol=r_F_sol)
            CHCKERR('')

            if (rich_lvl.eq.1) then
                call writo('Set up solution grid')
                call lvl_ud(1)

                call writo('Calculate the grid')
                call lvl_ud(1)
                ierr = setup_grid_sol(grid_eq,grid_X,grid_sol,r_F_sol,&
                    &sol_limits)
                CHCKERR('')
                call lvl_ud(-1)

                call writo('Write to output file')
                call lvl_ud(1)
                ierr = print_output_grid(grid_sol,'solution','sol',&
                    &remove_previous_arrs=(jump_to_sol.and.&
                    &(X_grid_style.eq.1 .or. X_grid_style.eq.3)))
                CHCKERR('')
                call lvl_ud(-1)

                ! set up modes
                ierr = setup_modes(mds_sol,grid_eq,grid_sol,plot_name='sol')
                CHCKERR('')

                call lvl_ud(-1)
            else
                ! restore solution grid and trim it
                ierr = reconstruct_PB3D_grid(grid_sol,'sol',&
                    &grid_limits=sol_limits)
                CHCKERR('')
                ierr = trim_grid(grid_sol,grid_sol_trim)
                CHCKERR('')

                ! set up modes
                ierr = setup_modes(mds_sol,grid_eq,grid_sol,plot_name='sol')
                CHCKERR('')

                ! reconstruct solution on trimmed grid
                ierr = reconstruct_PB3D_sol(mds_sol,grid_sol_trim,sol,'sol',&
                    &rich_lvl=rich_lvl-1)
                CHCKERR('')

                ! clean up
                call grid_sol_trim%dealloc()
            end if

            deallocate(r_F_sol)
        end if

        ! set up  whether Richardson level  has to be  appended to the  name and
        ! whether to do vacuum operations
        select case (eq_style)
            case (1)                                                            ! VMEC
                rich_lvl_name = rich_lvl                                        ! append richardson level
                do_vac_ops = .true.
            case (2)                                                            ! HELENA
                rich_lvl_name = 0                                               ! do not append name
                if (rich_lvl.eq.rich_restart_lvl) then
                    do_vac_ops = .true.
                else
                    do_vac_ops = .false.
                end if
        end select

        if (do_vac_ops) then
            ! calculate vacuum
            ierr = calc_vac_res(mds_sol,vac)
            CHCKERR('')

            call writo('Write to output file')
            call lvl_ud(1)
            CHCKERR('')
            if (rank.eq.n_procs-1) then
                ierr = print_output_vac(vac,'vac',rich_lvl=rich_lvl_name)
                CHCKERR('')
            end if
            call lvl_ud(-1)
        end if

        ! initialize interpolated X
        ierr = grid_X_sol%init([1,n_alpha,grid_sol%n(3)],&
            &[grid_sol%i_min,grid_sol%i_max],grid_sol%divided)
        CHCKERR('')
        grid_X_sol%r_F = grid_sol%r_F
        grid_X_sol%r_E = grid_sol%r_E
        grid_X_sol%loc_r_F = grid_sol%loc_r_F
        grid_X_sol%loc_r_E = grid_sol%loc_r_E

        select case (X_grid_style)
            case (1,3)                                                          ! equilibrium or enriched
                ! user output
                call writo('Redistribute and interpolate the perturbation &
                    &variables to solution grid')
                call lvl_ud(1)

                ! redistribute grid and X variables
                ierr = redistribute_output_grid(grid_X,grid_X_rdst)
                CHCKERR('')
                ierr = redistribute_output_X(mds_X,grid_X,grid_X_rdst,X,X_rdst)
                CHCKERR('')

                ! initialize
                call X_sol%init(mds_sol,grid_X_sol,is_field_averaged=.true.)

                ! interpolate
                ierr = interp_V(mds_X,grid_X_rdst,X_rdst,mds_sol,grid_X_sol,&
                    &X_sol)
                CHCKERR('')

                ! clean up
                call X_rdst%dealloc()
                call grid_X_rdst%dealloc()

                call lvl_ud(-1)
            case (2)                                                            ! solution
                ! user output
                call writo('Copy the perturbation variables to solution grid')
                call lvl_ud(1)

                ! copy
                call X%copy(mds_sol,grid_X,X_sol)

                call lvl_ud(-1)
        end select

#if ldebug
        ! write integrated  field-aligned tensorial  perturbation quantities
        ! to output and plot
        if (debug_run_driver_sol) then
            ierr = print_debug_X_2(mds_sol,grid_X_sol,X_sol)
            CHCKERR('')
        end if
#endif

        ! clean up
        call grid_X_sol%dealloc()

        ! solve the system
        call writo('Solving the system')
        call lvl_ud(1)
        select case (EV_style)
            case(1)                                                             ! SLEPC solver for EV problem
                ! solve the system
                ierr = solve_EV_system_SLEPC(mds_sol,grid_sol,X_sol,vac,sol)
                CHCKERR('')
            case default
                err_msg = 'No EV solver style associated with '//&
                    &trim(i2str(EV_style))
                ierr = 1
                CHCKERR(err_msg)
        end select
        call lvl_ud(-1)
        call writo('System solved')

        ! write solution variables to output
        ierr = print_output_sol(grid_sol,sol,'sol',rich_lvl=rich_lvl,&
            &remove_previous_arrs=&
            &(jump_to_sol.and.(X_grid_style.eq.1 .or. X_grid_style.eq.3)))
        CHCKERR('')

        ! calculate Richardson extrapolation factors if necessary
        call calc_rich_ex(sol%val)

        ! clean up
        call X_sol%dealloc()
    end function run_driver_sol

    !> \public  Interpolate  tensorial  perturbation  quantities  in  the  third
    !! dimension.
    !!
    !! The input grid should not be divided, though the output grid can be.
    !!
    !! The procedure  considers all  possible mode  number combinations.  For \c
    !! X_style 2  (fast), each  secondary mode  only lives  in a  certain normal
    !! range  of the  plasma. Therefore,  each secondary  mode pair  also has  a
    !! limited normal range, given by the overlap of the ranges of the members.
    !!
    !! The interpolated mode number combinations  have a normal range that might
    !! slightly differ from the input ranges, in which case extrapolation can be
    !! done if the method allows for it.
    !!
    !! \note If the input grid is too  coarse, but the interpolated grid is not,
    !! it is in  theory possible that there are mode  numbers and therefore mode
    !! number pairs  that exist  in the  input grid,  though they  do so  in the
    !! interpolated  grid.  In  this  case,  they  can  not  be  calculated  and
    !! are  currently  set  to  zero. However,  the  current  implementation  of
    !! setup_modes()  produces an  error if  less than  the full  range of  mode
    !! numbers is accurately present.
    integer function interp_V(mds_i,grid_i,X_i,mds_o,grid_o,X_o) result(ierr)
        use X_utilities, only: is_necessary_X, trim_modes
        use X_vars, only: n_mod_X
        use num_vars, only: rank
        use eq_vars, only: max_flux_F
        use num_utilities, only: c
        use sol_utilities, only: interp_V_spline
!#if ldebug
        !use num_vars, only: ex_plot_style
!#endif

        character(*), parameter :: rout_name = 'interp_V'

        ! input / output
        type(modes_type), intent(in), target :: mds_i                           !< general modes variables for input
        type(grid_type), intent(in) :: grid_i                                   !< grid at which \c X_i is tabulated
        type(X_2_type), intent(in), target :: X_i                               !< tensorial perturbation variable on input grid
        type(modes_type), intent(in), target :: mds_o                           !< general modes variables for output
        type(grid_type), intent(in) :: grid_o                                   !< grid at which \c X_o is interpolated
        type(X_2_type), intent(inout), target :: X_o                            !< interpolated tensorial perturbation variable

        ! local variables
        integer :: kd                                                           ! counter
        integer :: k, m                                                         ! counters for mode numbers
        integer :: c_loc_i(2)                                                   ! local input c for symmetric and asymmetric variables
        integer :: c_loc_o(2)                                                   ! local output c for symmetric and asymmetric variables
        integer :: n_mod_tot                                                    ! total amount of modes
        integer :: kdl_i(2), kdl_o(2)                                           ! limits on normal index for a mode combination
        integer :: id_lim_i(2), id_lim_o(2)                                     ! limits on total modes
        integer :: ivs_stat(6)                                                  ! local stats for interp_V_spline
        integer, pointer :: sec_i_loc(:,:), sec_o_loc(:,:)                      ! pointers to secondary mode variables
        real(dp), pointer :: r_i_loc(:), r_o_loc(:)                             ! local r_i and r_o for a mode combination
        complex(dp), pointer :: V_i(:,:,:), V_o(:,:,:)                          ! pointers to input and output PV_i and KV_i
        logical :: calc_this(2)                                                 ! whether this combination needs to be calculated
        logical :: extrap = .true.                                              ! whether extrapolation is used
        character(len=max_str_ln) :: err_msg                                    ! error message
#if ldebug
        integer :: km_id
        integer, allocatable :: norm_ext_i(:,:)                                 ! normal extent for input quantity mode combinations
        real(dp), allocatable :: r_loc_tot(:,:,:)                               ! r_i_loc and r_o_loc for all combinations
        character(len=max_str_ln) :: ivs_stats_names(3)                         ! names used in interp_V_spline statistics
        !complex(dp), allocatable :: V_plot(:,:)                                 ! for debug plotting of interpolated V
#endif

        ! initialize ierr
        ierr = 0

        ! test
        if (grid_i%n(2).ne.grid_o%n(2)) then
            ierr = 1
            err_msg = 'input and output grid are not compatible in the &
                &geodesic coordinate: '//trim(i2str(grid_i%n(2)))//' vs. '//&
                &trim(i2str(grid_o%n(2)))
            CHCKERR(err_msg)
        end if

        ! initialize to zero
        X_o%PV_0 = 0._dp
        X_o%PV_1 = 0._dp
        X_o%PV_2 = 0._dp
        X_o%KV_0 = 0._dp
        X_o%KV_1 = 0._dp
        X_o%KV_2 = 0._dp

        ! limit input modes to output modes
        ! (X grid should comprise sol grid)
        ierr = trim_modes(mds_i,mds_o,id_lim_i,id_lim_o)
        CHCKERR('')
        sec_i_loc => mds_i%sec(id_lim_i(1):id_lim_i(2),:)
        sec_o_loc => mds_o%sec(id_lim_o(1):id_lim_o(2),:)

        ! select all mode number combinations of interpolated grid
        n_mod_tot = size(sec_o_loc,1)
#if ldebug
        allocate(r_loc_tot(2,n_mod_tot**2,3))
        allocate(norm_ext_i(n_mod_tot,n_mod_tot))
        km_id = 0
        norm_ext_i = 0
#endif
        do m = 1,n_mod_tot
            do k = 1,n_mod_tot
#if ldebug
                ! test whether input and output mode numbers are consistent
                if (sec_i_loc(k,1).ne.sec_o_loc(k,1) .or. &
                    &sec_i_loc(m,1).ne.sec_o_loc(m,1)) then
                    ierr = 1
                    err_msg = 'For ('//trim(i2str(k))//','//trim(i2str(m))//&
                        &'), no consistency for modes in grid_i and grid_o'
                    CHCKERR(err_msg)
                end if
#endif

                ! set input and output grid normal limits for mode pair (k,m)
                kdl_i(1) = max(sec_i_loc(k,2),sec_i_loc(m,2))
                kdl_i(2) = min(sec_i_loc(k,3),sec_i_loc(m,3))
                kdl_o(1) = max(sec_o_loc(k,2),sec_o_loc(m,2))
                kdl_o(2) = min(sec_o_loc(k,3),sec_o_loc(m,3))

                ! limit to grid ranges
                kdl_i(1) = max(kdl_i(1),grid_i%i_min)
                kdl_i(2) = min(kdl_i(2),grid_i%i_max)
                kdl_o(1) = max(kdl_o(1),grid_o%i_min)
                kdl_o(2) = min(kdl_o(2),grid_o%i_max)

                ! limit output range to input range if not extrapolating
                if (.not.extrap) then
                    do kd = kdl_o(1),kdl_o(2)
                        if (grid_o%r_F(kd).ge.grid_i%r_F(kdl_i(1))) exit
                    end do
                    kdl_o(1) = kd

                    do kd = kdl_o(2),kdl_o(1),-1
                        if (grid_o%r_F(kd).le.grid_i%r_F(kdl_i(2))) exit
                    end do
                    kdl_o(2) = kd
                end if

                ! cycle if mode combination does not exist anywhere
                if (kdl_i(1).gt.kdl_i(2) .or. kdl_o(1).gt.kdl_o(2)) cycle

                ! convert limits to local
                kdl_i = kdl_i - grid_i%i_min + 1
                kdl_o = kdl_o - grid_o%i_min + 1

                ! get input ranges for this mode number
                r_i_loc => grid_i%loc_r_F(kdl_i(1):kdl_i(2))
                r_o_loc => grid_o%loc_r_F(kdl_o(1):kdl_o(2))

                ! check whether mode combination needs to be calculated
                calc_this(1) = is_necessary_X(.true.,&
                    &[sec_o_loc(k,1),sec_o_loc(m,1)])
                calc_this(2) = is_necessary_X(.false.,&
                    &[sec_o_loc(k,1),sec_o_loc(m,1)])

#if ldebug
                km_id = km_id + 1
                r_loc_tot(:,km_id,1) = [r_i_loc(1),r_i_loc(size(r_i_loc))]
                r_loc_tot(:,km_id,2) = [r_o_loc(1),r_o_loc(size(r_o_loc))]
                r_loc_tot(:,km_id,3) = km_id
                norm_ext_i(k,m) = size(r_i_loc)
                if (debug_interp_V) write(*,*) 'k, m', k, m, 'of', n_mod_tot, &
                        &'with calc_this = ', calc_this
#endif

                ! set up c_loc for input and output
                c_loc_i(1) = c([sec_i_loc(k,4),sec_i_loc(m,4)],.true.,n_mod_X)
                c_loc_i(2) = c([sec_i_loc(k,4),sec_i_loc(m,4)],.false.,n_mod_X)
                c_loc_o(1) = c([sec_o_loc(k,4),sec_o_loc(m,4)],.true.,n_mod_X)
                c_loc_o(2) = c([sec_o_loc(k,4),sec_o_loc(m,4)],.false.,n_mod_X)

                if (calc_this(1)) then
                    V_i => X_i%PV_0(:,:,kdl_i(1):kdl_i(2),c_loc_i(1))
                    V_o => X_o%PV_0(:,:,kdl_o(1):kdl_o(2),c_loc_o(1))
                    ierr = interp_V_spline(V_i,V_o,r_i_loc,r_o_loc,extrap,&
                        &ivs_stat(1))
                    CHCKERR('')

                    V_i => X_i%PV_2(:,:,kdl_i(1):kdl_i(2),c_loc_i(1))
                    V_o => X_o%PV_2(:,:,kdl_o(1):kdl_o(2),c_loc_o(1))
                    ierr = interp_V_spline(V_i,V_o,r_i_loc,r_o_loc,extrap,&
                        &ivs_stat(2))
                    CHCKERR('')

                    V_i => X_i%KV_0(:,:,kdl_i(1):kdl_i(2),c_loc_i(1))
                    V_o => X_o%KV_0(:,:,kdl_o(1):kdl_o(2),c_loc_o(1))
                    ierr = interp_V_spline(V_i,V_o,r_i_loc,r_o_loc,extrap,&
                        &ivs_stat(3))
                    CHCKERR('')

                    V_i => X_i%KV_2(:,:,kdl_i(1):kdl_i(2),c_loc_i(1))
                    V_o => X_o%KV_2(:,:,kdl_o(1):kdl_o(2),c_loc_o(1))
                    ierr = interp_V_spline(V_i,V_o,r_i_loc,r_o_loc,extrap,&
                        &ivs_stat(4))
                    CHCKERR('')

#if ldebug
                    if (debug_interp_V) then
                        !call writo('For [k,m] = ['//&
                            !&trim(i2str(sec_o_loc(k,1)))//','//&
                            !&trim(i2str(sec_o_loc(m,1)))//']:')
                        !call lvl_ud(1)
                        !call writo('input normal range: '//&
                            !&trim(i2str(kdl_i(1)))//'..'//trim(i2str(kdl_i(2))))
                        !call writo('output normal range: '//&
                            !&trim(i2str(kdl_o(1)))//'..'//trim(i2str(kdl_o(2))))
                        !allocate(V_plot(max(size(V_i,3),size(V_o,3)),4))
                        !V_plot(1:size(V_i,3),1) = V_i(1,1,:)
                        !V_plot(1:size(V_o,3),2) = V_o(1,1,:)
                        !V_plot(1:size(V_i,3),3) = r_i_loc/max_flux_F*2*pi
                        !V_plot(1:size(V_o,3),4) = r_o_loc/max_flux_F*2*pi
                        !V_plot(size(V_i,3)+1:size(V_plot,1),1) = &
                            !&V_plot(size(V_i,3),1)
                        !V_plot(size(V_o,3)+1:size(V_plot,1),2) = &
                            !&V_plot(size(V_o,3),2)
                        !V_plot(size(V_i,3)+1:size(V_plot,1),3) = &
                            !&V_plot(size(V_i,3),3)
                        !V_plot(size(V_o,3)+1:size(V_plot,1),4) = &
                            !&V_plot(size(V_o,3),4)
                        !call print_ex_2D(['V_i','V_o'],'RE_V',&
                            !&rp(V_plot(:,1:2)),x=rp(V_plot(:,3:4)))
                        !!call print_ex_2D(['V_i','V_o'],'IM_V',&
                            !!&ip(V_plot(:,1:2)),x=rp(V_plot(:,3:4)))
                        !deallocate(V_plot)
                        !if (ex_plot_style.eq.2) read(*,*)
                        !call lvl_ud(-1)
                    end if
#endif
                end if

                if (calc_this(2)) then
                    V_i => X_i%PV_1(:,:,kdl_i(1):kdl_i(2),c_loc_i(2))
                    V_o => X_o%PV_1(:,:,kdl_o(1):kdl_o(2),c_loc_o(2))
                    ierr = interp_V_spline(V_i,V_o,r_i_loc,r_o_loc,extrap,&
                        &ivs_stat(5))
                    CHCKERR('')

                    V_i => X_i%KV_1(:,:,kdl_i(1):kdl_i(2),c_loc_i(2))
                    V_o => X_o%KV_1(:,:,kdl_o(1):kdl_o(2),c_loc_o(2))
                    ierr = interp_V_spline(V_i,V_o,r_i_loc,r_o_loc,extrap,&
                        &ivs_stat(6))
                    CHCKERR('')
                end if
#if ldebug
                do kd = 1,size(ivs_stat)
                    !if (ivs_stat(kd).eq.1) call writo('Used direct copy for '//&
                        !&trim(i2str(kd))//' of mode combination ('//&
                        !&trim(i2str(k))//', '//trim(i2str(m))//') at ψ = '//&
                        !&trim(r2strt(r_loc_tot(1,km_id,1)/max_flux_F*2*pi)))
                    ivs_stats(ivs_stat(kd)) = ivs_stats(ivs_stat(kd)) + 1
                end do
#endif
            end do
        end do

#if ldebug
        call writo('interp_V_spline statistics:')
        call lvl_ud(1)
        ivs_stats_names(1) = 'direct copy'
        ivs_stats_names(2) = 'linear interpolation'
        ivs_stats_names(3) = 'spline interpolation'
        do kd = 1,3
            call writo(trim(ivs_stats_names(kd))//' was performed '//&
                &trim(i2str(ivs_stats(kd)))//' times: '//&
                &trim(r2strt(100._dp*ivs_stats(kd)/sum(ivs_stats)))//'%')
        end do
        if (ivs_stats(1).gt.0) then
            call writo('It is possible that more secondary modes and/or a &
                &lower max_njq_change are needed.',warning=.true.)
            call writo('If unphysical modes appear, consider changing these &
                &variables.',warning=.true.)
        end if
        call lvl_ud(-1)
        call print_ex_2D([''],'r_i_loc_'//trim(i2str(rank)),&
            &r_loc_tot(:,1:km_id,1)/max_flux_F*2*pi,x=r_loc_tot(:,1:km_id,3),&
            &draw=.false.)
        call draw_ex([''],'r_i_loc_'//trim(i2str(rank)),km_id,1,.false.,&
            &ex_plot_style=1)
        call print_ex_2D([''],'r_o_loc_'//trim(i2str(rank)),&
            &r_loc_tot(:,1:km_id,2)/max_flux_F*2*pi,x=r_loc_tot(:,1:km_id,3),&
            &draw=.false.)
        call draw_ex([''],'r_o_loc_'//trim(i2str(rank)),km_id,1,.false.,&
            &ex_plot_style=1)
        call plot_HDF5('normal extent','norm_ext_i_'//trim(i2str(rank)),1._dp*&
            &reshape(norm_ext_i*1._dp,[n_mod_tot,n_mod_tot,1]))
#endif

        ! clean up
        nullify(r_i_loc,r_o_loc)
        nullify(sec_i_loc,sec_o_loc)
        nullify(V_i,V_o)
    end function interp_V
end module driver_sol


text:
program modulo_01
real :: x1
integer :: i1
x1 = 3.3
print*, modulo(x1, 3.)
print*, modulo(x1, 2.)
print*, modulo(x1, 1.)
i1 = 4
print*, modulo(i1, 6)
print*, modulo(i1, 3)
print*, modulo(i1, 1)
end program

max_stars_repo_name: Thirumalai-Shaktivel/lfortran

max_issues_repo_name: Thirumalai-Shaktivel/lfortran

max_forks_repo_name: Thirumalai-Shaktivel/lfortran


text:
!  *********************************************************************
!  *                                                                   *
!  *                         subroutine genfld                         *
!  *                                                                   *
!  *********************************************************************
!  Single Precision Version 1.1
!  Written by Gordon A. Fenton, DalTech, Jan 31, 2001
!  Latest Update: Feb 22, 2005
!
!  PURPOSE  generates a soil property field
!
!  DESCRIPTION
!  This routine generates a soil property field. The field may be
!  deterministic, white noise, or a correlated random field.
!
!  ARGUMENTS
!
!   istat	unit number connected to a file to which the input data is to
!		echoed if echo is true. (input)
!
!    pfld	real array of size at least (nxe x nye) which, on output,
!		will contain the underlying soil property field (ie, the
!		normally distributed field, if the property is random).
!		(output)
!
!       p	real vector of length at least 7 containing the statistics
!		of the soil property field. In particular,
!		   p(1) = soil property mean
!		   p(2) = soil property standard deviation,
!		   p(3) = soil property distribution type;
!			  = 0.0 if property is deterministic (at mean value)
!			  = 1.0 if property is normally distributed
!			  = 2.0 if property is lognormally distributed (logn)
!			  = 3.0 if property is bounded
!			  = 4.0 if property is a function of phi
!		   p(4) = lower bound (bounded), or mean log-property (logn),
!			  or constant in functional relationship (f(phi))
!		   p(5) = upper bound (bounded), or sd of log-property (logn),
!			  or slope in functional relationship (f(phi))
!		   p(6) = m parameter (if bounded), or function type (f(phi))
!			  where in the last case,
!				p(6) = 0.0 for 1*phi
!				p(6) = 1.0 for sin(phi)
!				p(6) = 2.0 for tan(phi)
!		   p(7) = s parameter (if bounded)
!		If the soil property is bounded or functional, then p(1) and
!		p(2) are ignored and the parameters c(4) through c(7)
!		completely describe the distribution (or function). (input)
!
!     nxe	integer giving the number of elements describing the soil
!		mass in the x-direction (horizontally). (input)
!
!     nye	integer giving the number of elements describing the soil
!		mass in the y-direction (vertically). (input)
!
!      XL	real value containing the physical size of the soil mass in
!		the x-direction (horizontally). (input)
!
!      YL	real value containing the physical size of the soil mass in
!		the y-direction (vertically). (input)
!
!  varfnc	character string containing the name of the variance
!		function controlling the random fields. Possible variance
!		functions are
!		`dlavx2' - 2-D exponentially decaying (Markov) model
!			   requires X- and Y-direction scales of fluctuation
!		`dlafr2' - 2-D isotropic fractional Gaussian noise model
!			   requires (H,delta) as parameters. In this case,
!			   thx is H, and delta is the minimum element
!			   dimension.
!		`dlsep2' - 2-D separable (1D x 1D) Markov model
!			   requires X- and Y-direction scales of fluctuation
!		`dlsfr2' - 2-D separable fractional Gaussian noise model
!			   requires (H_x,H_y,delta) as parameters. In this
!			   case, thx is H_x, thy is H_y, and delta is the
!			   minimum element dimension.
!		`dlspx2' - 2-D separable Gaussian decaying model
!			   requires X- and Y-direction scales of fluctuation
!		(input)
!
!   kseed	integer giving the seed to be used to initialize the
!		pseudo-random number generator. Subsequent runs using
!		the same seed will result in the same sequence of random
!		numbers. (input)
!
!
!  REVISION HISTORY:
!  1.1	reordered numbering of f(phi) function types (0=1,1=sin,2=tan)
!	(this only changed the documentation above) (Feb 22/05)
!------------------------------------------------------------------------
      subroutine genfld(istat,pfld,p,nxe,nye,XL,YL,varfnc,kseed,liid)
      real pfld(nxe,*), p(*)
      real XL, YL
      integer istat, nxe, nye, kseed
      character*(*) varfnc
      logical liid
      real*8   dlavx2, dlsep2, dlspx2, dlafr2, dlsfr2
      external dlavx2, dlsep2, dlspx2, dlafr2, dlsfr2
      data zero/0.0/, one/1.0/
      data half/0.5/, thrpt5/3.5/

      if( p(3) .lt. half ) then			! property is deterministic
         do 10 j = 1, nye
         do 10 i = 1, nxe
            pfld(i,j) = zero			! we'll add mean back on later
  10     continue
      elseif( liid ) then			! white noise field
         do 20 j = 1, nye
         do 20 i = 1, nxe
            pfld(i,j) = gausv(one)
  20     continue
      elseif( p(3) .lt. thrpt5 ) then		! generate standard normal
         if( varfnc .eq. 'dlavx2' ) then
            call las2g(pfld,nxe,nye,XL,YL,dlavx2,kseed,0,istat)
         elseif( varfnc .eq. 'dlsep2' ) then
            call las2g(pfld,nxe,nye,XL,YL,dlsep2,kseed,0,istat)
         elseif( varfnc .eq. 'dlspx2' ) then
            call las2g(pfld,nxe,nye,XL,YL,dlspx2,kseed,0,istat)
         elseif( varfnc .eq. 'dlafr2' ) then
            call las2g(pfld,nxe,nye,XL,YL,dlafr2,kseed,0,istat)
         elseif( varfnc .eq. 'dlsfr2' ) then
            call las2g(pfld,nxe,nye,XL,YL,dlsfr2,kseed,0,istat)
         endif
      endif
!						all done
      return
      end

max_stars_repo_name: leemargetts/PoreFEM

max_issues_repo_name: leemargetts/PoreFEM

max_forks_repo_name: leemargetts/PoreFEM


! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
! Parameter Module File
!
! Generated by KPP-2.1 symbolic chemistry Kinetics PreProcessor
!       (http://www.cs.vt.edu/~asandu/Software/KPP)
! KPP is distributed under GPL, the general public licence
!       (http://www.gnu.org/copyleft/gpl.html)
! (C) 1995-1997, V. Damian & A. Sandu, CGRER, Univ. Iowa
! (C) 1997-2005, A. Sandu, Michigan Tech, Virginia Tech
!     With important contributions from:
!        M. Damian, Villanova University, USA
!        R. Sander, Max-Planck Institute for Chemistry, Mainz, Germany
!
! File                 : mozart_mosaic_4bin_vbs0_Parameters.f90
! Time                 : Mon Feb 27 11:08:41 2017
! Working directory    : /data/ksetigui/setigui/WRFV3_171015/chem/KPP/mechanisms/mozart_mosaic_4bin_vbs0
! Equation file        : mozart_mosaic_4bin_vbs0.kpp
! Output root filename : mozart_mosaic_4bin_vbs0
!
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

MODULE mozart_mosaic_4bin_vbs0_Parameters

  USE mozart_mosaic_4bin_vbs0_Precision
  PUBLIC
  SAVE

! NSPEC - Number of chemical species
  INTEGER, PARAMETER :: NSPEC = 91
! NVAR - Number of Variable species
  INTEGER, PARAMETER :: NVAR = 89
! NVARACT - Number of Active species
  INTEGER, PARAMETER :: NVARACT = 82
! NFIX - Number of Fixed species
  INTEGER, PARAMETER :: NFIX = 2
! NREACT - Number of reactions
  INTEGER, PARAMETER :: NREACT = 200
! NVARST - Starting of variables in conc. vect.
  INTEGER, PARAMETER :: NVARST = 1
! NFIXST - Starting of fixed in conc. vect.
  INTEGER, PARAMETER :: NFIXST = 90
! NONZERO - Number of nonzero entries in Jacobian
  INTEGER, PARAMETER :: NONZERO = 762
! LU_NONZERO - Number of nonzero entries in LU factoriz. of Jacobian
  INTEGER, PARAMETER :: LU_NONZERO = 858
! CNVAR - (NVAR+1) Number of elements in compressed row format
  INTEGER, PARAMETER :: CNVAR = 90
! NLOOKAT - Number of species to look at
  INTEGER, PARAMETER :: NLOOKAT = 0
! NMONITOR - Number of species to monitor
  INTEGER, PARAMETER :: NMONITOR = 0
! NMASS - Number of atoms to check mass balance
  INTEGER, PARAMETER :: NMASS = 1
! PI - Value of pi
  REAL(kind=dp), PARAMETER :: PI = 3.14159265358979

! Index declaration for variable species in C and VAR
!   VAR(ind_spc) = C(ind_spc)

  INTEGER, PARAMETER :: ind_SO4 = 1
  INTEGER, PARAMETER :: ind_NUME = 2
  INTEGER, PARAMETER :: ind_DEN = 3
  INTEGER, PARAMETER :: ind_BIOG1_c = 4
  INTEGER, PARAMETER :: ind_BIOG1_o = 5
  INTEGER, PARAMETER :: ind_SMPA = 6
  INTEGER, PARAMETER :: ind_VOCA = 7
  INTEGER, PARAMETER :: ind_SMPBB = 8
  INTEGER, PARAMETER :: ind_VOCBB = 9
  INTEGER, PARAMETER :: ind_NH3 = 10
  INTEGER, PARAMETER :: ind_C2H6 = 11
  INTEGER, PARAMETER :: ind_C3H8 = 12
  INTEGER, PARAMETER :: ind_BIGENE = 13
  INTEGER, PARAMETER :: ind_BIGALK = 14
  INTEGER, PARAMETER :: ind_N2O = 15
  INTEGER, PARAMETER :: ind_SO2 = 16
  INTEGER, PARAMETER :: ind_TOLUENE = 17
  INTEGER, PARAMETER :: ind_CRESOL = 18
  INTEGER, PARAMETER :: ind_H2O2 = 19
  INTEGER, PARAMETER :: ind_EO = 20
  INTEGER, PARAMETER :: ind_N2O5 = 21
  INTEGER, PARAMETER :: ind_XOH = 22
  INTEGER, PARAMETER :: ind_XOOH = 23
  INTEGER, PARAMETER :: ind_DMS = 24
  INTEGER, PARAMETER :: ind_C2H5OH = 25
  INTEGER, PARAMETER :: ind_MACROOH = 26
  INTEGER, PARAMETER :: ind_H2 = 27
  INTEGER, PARAMETER :: ind_C2H5OOH = 28
  INTEGER, PARAMETER :: ind_C3H7OOH = 29
  INTEGER, PARAMETER :: ind_ROOH = 30
  INTEGER, PARAMETER :: ind_ENEO2 = 31
  INTEGER, PARAMETER :: ind_MEKOOH = 32
  INTEGER, PARAMETER :: ind_CH3OOH = 33
  INTEGER, PARAMETER :: ind_HO2NO2 = 34
  INTEGER, PARAMETER :: ind_HYDRALD = 35
  INTEGER, PARAMETER :: ind_C2H4 = 36
  INTEGER, PARAMETER :: ind_ONIT = 37
  INTEGER, PARAMETER :: ind_TOLOOH = 38
  INTEGER, PARAMETER :: ind_EO2 = 39
  INTEGER, PARAMETER :: ind_CH3COOOH = 40
  INTEGER, PARAMETER :: ind_MEK = 41
  INTEGER, PARAMETER :: ind_CH3COOH = 42
  INTEGER, PARAMETER :: ind_GLYOXAL = 43
  INTEGER, PARAMETER :: ind_TERPOOH = 44
  INTEGER, PARAMETER :: ind_HNO3 = 45
  INTEGER, PARAMETER :: ind_POOH = 46
  INTEGER, PARAMETER :: ind_PAN = 47
  INTEGER, PARAMETER :: ind_CH4 = 48
  INTEGER, PARAMETER :: ind_ISOPOOH = 49
  INTEGER, PARAMETER :: ind_MPAN = 50
  INTEGER, PARAMETER :: ind_CH3OH = 51
  INTEGER, PARAMETER :: ind_TOLO2 = 52
  INTEGER, PARAMETER :: ind_BIGALD = 53
  INTEGER, PARAMETER :: ind_C10H16 = 54
  INTEGER, PARAMETER :: ind_O1D_CB4 = 55
  INTEGER, PARAMETER :: ind_ALKOOH = 56
  INTEGER, PARAMETER :: ind_MEKO2 = 57
  INTEGER, PARAMETER :: ind_CH3COCH3 = 58
  INTEGER, PARAMETER :: ind_PO2 = 59
  INTEGER, PARAMETER :: ind_O = 60
  INTEGER, PARAMETER :: ind_ISOPNO3 = 61
  INTEGER, PARAMETER :: ind_GLYALD = 62
  INTEGER, PARAMETER :: ind_ALKO2 = 63
  INTEGER, PARAMETER :: ind_C2H5O2 = 64
  INTEGER, PARAMETER :: ind_CO = 65
  INTEGER, PARAMETER :: ind_C3H7O2 = 66
  INTEGER, PARAMETER :: ind_ISOP = 67
  INTEGER, PARAMETER :: ind_ONITR = 68
  INTEGER, PARAMETER :: ind_C3H6 = 69
  INTEGER, PARAMETER :: ind_CH3CHO = 70
  INTEGER, PARAMETER :: ind_TERPO2 = 71
  INTEGER, PARAMETER :: ind_RO2 = 72
  INTEGER, PARAMETER :: ind_HYAC = 73
  INTEGER, PARAMETER :: ind_MACR = 74
  INTEGER, PARAMETER :: ind_CH3COCHO = 75
  INTEGER, PARAMETER :: ind_XO2 = 76
  INTEGER, PARAMETER :: ind_MVK = 77
  INTEGER, PARAMETER :: ind_MACRO2 = 78
  INTEGER, PARAMETER :: ind_CH2O = 79
  INTEGER, PARAMETER :: ind_ISOPO2 = 80
  INTEGER, PARAMETER :: ind_NO3 = 81
  INTEGER, PARAMETER :: ind_CH3O2 = 82
  INTEGER, PARAMETER :: ind_O3 = 83
  INTEGER, PARAMETER :: ind_OH = 84
  INTEGER, PARAMETER :: ind_MCO3 = 85
  INTEGER, PARAMETER :: ind_HO2 = 86
  INTEGER, PARAMETER :: ind_NO = 87
  INTEGER, PARAMETER :: ind_NO2 = 88
  INTEGER, PARAMETER :: ind_CH3CO3 = 89

! Index declaration for fixed species in C
!   C(ind_spc)

  INTEGER, PARAMETER :: ind_H2O = 90
  INTEGER, PARAMETER :: ind_M = 91

! Index declaration for fixed species in FIX
!    FIX(indf_spc) = C(ind_spc) = C(NVAR+indf_spc)

  INTEGER, PARAMETER :: indf_H2O = 1
  INTEGER, PARAMETER :: indf_M = 2

END MODULE mozart_mosaic_4bin_vbs0_Parameters


text:
MODULE neo_control
! Control parameters from inpt file
  USE neo_precision
  CHARACTER(120)                     :: in_file                             !SPH
  CHARACTER(120)                     :: out_file                            !SPH
  CHARACTER(120)                      :: cur_file
  INTEGER                            :: theta_n
  INTEGER                            :: phi_n
  INTEGER                            :: s_ind_in
  INTEGER                            :: write_progress
  INTEGER                            :: write_output_files
  INTEGER                            :: spline_test
  INTEGER                            :: max_m_mode, max_n_mode
  INTEGER                            :: lab_swi, inp_swi, ref_swi, eout_swi
  INTEGER                            :: no_fluxs
  INTEGER, DIMENSION(:), ALLOCATABLE :: fluxs_arr
END MODULE neo_control

max_stars_repo_name: joseluisvelasco/STELLOPT

max_issues_repo_name: joseluisvelasco/STELLOPT

max_forks_repo_name: joseluisvelasco/STELLOPT


text:
      SUBROUTINE SMCOMP (*,ZI,ZR,ZD)
C
C DRIVER PROGRAM FOR SYMMETRIC DECOMPOSITION.  SUBROUTINE SMCPH1 READS
C THE INPUT MATRIX AND STORES THE DATA EITHER IN MEMORY OR ON THE
C SPILL FILE.  SUBROUTINE SMCPH2 IS THEN CALLED TO PERFORM THE
C MATRIX DECOMPOSITION.
C
      REAL              ZR(4)
      INTEGER           ZI(4)
      INTEGER           MODULE(5), BEGN, END
      DOUBLE PRECISION  ZD(4)
      INCLUDE          'SMCOMX.COM'
C  mcb   - matrix control block for input matrix
C  lll   - matrix control block for lower triangular matrix
C  dbc   - dbc(1) = available scratch file, dbc(2-7) are not used
C  scr1, scr2, scr3 - three available scratch files
C  lcore - amount of open core available for use
C  ddr   - d.p. values of (real, imaginary) for scaled value of determinant
C  power - scale factor to apply to determinant, determinant=det * 10**power
C  mindd - d.p. value for minimum value of diagonal elements
C  chlsky - cholesky option when =1, i.e., form c matrix
C
      DATA  MODULE / 4HSMCO, 4HMP  , 3*4H    /
      DATA  BEGN   / 4HBEGN /
      DATA  END    / 4HEND  /
      IERROR = 0
      NCOL   = MCB(2)
      MODULE( 3 ) = BEGN
      STURM  = 0
      CALL CONMSG ( MODULE, 5, 0 )
      CALL SMCPH1 ( ZI, ZR, ZD )
      IF ( IERROR .EQ. 1 ) GO TO 701
      IF ( IERROR .NE. 0 ) GO TO 700
      CALL SMCPH2 ( ZI, ZR, ZD )
      IF ( IERROR .EQ. 1 ) GO TO 701
C
C print roots information if this is an eigenvalue problem, and keep
C two largest shift point data if several shift point movings are involved.
C
      IF ( SHFTPT .GT. 0. ) WRITE ( NOUT, 901 ) STURM, SHFTPT
901   FORMAT( 20X, I5, ' ROOTS BELOW ', 1P,E14.6 )
      IF ( STURM .NE. 0 ) GO TO 100
      IF ( KEEP  .LE. 0 ) GO TO 700
      STURM  = KEEP
      SHFTPT = PTSHFT
      GO TO 700
100   IF ( KEEP .GT. STURM ) GO TO 700
      JJ     = KEEP
      RS     = PTSHFT
      KEEP   = STURM
      PTSHFT = JJ
      SHFTPT = RS
700   MODULE( 3 ) = END
      CALL CONMSG ( MODULE, 5, 0 )
      IF ( IERROR .NE. 0 ) RETURN 1
      GO TO 777
701   CONTINUE
      MODULE( 3 ) = END
      CALL CONMSG ( MODULE, 5, 0 )
777   CONTINUE
      RETURN
      END


text:
module timedisc_probe_mod

use constants_mod

contains

pure subroutine timedisc_probe(quad,flux,state)

    use mesh_types_mod, only: quad_t
    use share_mod, only: rt => runtime

    use timedisc_const_mod

    type(quad_t), intent(in)    :: quad

# if PP_N_DIMS == 2
    real(dp), intent(in)        :: flux(N_NODES,N_NODES,N_VARS)
    real(dp), intent(out)       :: state(N_NODES,N_NODES,N_VARS)
# elif PP_N_DIMS == 3
    real(dp), intent(in)        :: flux(N_NODES,N_NODES,N_NODES,N_VARS)
    real(dp), intent(out)       :: state(N_NODES,N_NODES,N_NODES,N_VARS)
# endif

    associate(dt => rt%timestep, &
              uu => quad%pld%hydro%state, &
              r0 => quad%aux%timedisc%register)

    select case(rt%rkstage)

        case(1)
            state = uu + dt * flux

        case(2)
            state = 3.0_dp/4.0_dp * r0 + 1.0_dp/4.0_dp * uu + 1.0_dp/4.0_dp * dt * flux

        case(3)
            state = 1.0_dp/3.0_dp * r0 + 2.0_dp/3.0_dp * uu + 2.0_dp/3.0_dp * dt * flux

    end select
    end associate

end subroutine

end module

max_stars_repo_name: jmark/nemo

max_issues_repo_name: jmark/nemo

max_forks_repo_name: jmark/nemo


        !COMPILER-GENERATED INTERFACE MODULE: Fri Mar  6 15:34:21 2020
        ! This source file is for reference only and may not completely
        ! represent the generated interface used by the compiler.
        MODULE TSTIPO__genmod
          INTERFACE
            SUBROUTINE TSTIPO
            END SUBROUTINE TSTIPO
          END INTERFACE
        END MODULE TSTIPO__genmod


text:
program t
  call a
end program t

subroutine s
  print '(a)','s'
  return
  entry a
  print '(a)','a'
end subroutine s


text:
	SUBROUTINE IG_IFPT  ( rimnem, iret )
C************************************************************************
C* IG_IFPT								*
C*									*
C* This subroutine prints the interface values arrays to the decoder	*
C* log if the verbosity level is 3 or higher.				*
C*									*
C* IG_IFPT  ( RIMNEM, IRET )						*
C*                                                                      *
C* Input parameters:                                                    *
C*	RIMNEM (*)	CHAR*		Interface mnemonics for reals   *
C*									*
C* Output parameters:							*
C*	IRET		INTEGER		Return code 			*
C*					  0 = normal return 		*
C*									*
C**									*
C* Log:									*
C* C. Caruso Magee/NCEP 06/05                                           *
C************************************************************************
	INCLUDE		'GEMPRM.PRM'
	INCLUDE		'igcmn.cmn'
C*
	CHARACTER*(*)	rimnem (*)
C-----------------------------------------------------------------------
	iret = 0
        loglvl = 3
C
C*      Print all of the single-level real interface values
C*      that are not missing.
C
	CALL DC_IFPR  ( 3, rimnem, rivals, NRSLMN, ierfpr )
C
C*      Print all of the multi-level interface values.
C*      Print backscatter and related parameters.
C
        IF  ( rivals ( irnbks ) .gt. 0 )  THEN
C
            WRITE ( UNIT = logmsg, FMT = '( 6A10 )' )
     +          rimnem ( irnbks ),     rimnem ( irtobd (1) ),
     +          rimnem ( irqcb1 (1) ), rimnem ( irbks1 (1) ),
     +          rimnem ( irfos1 (1) ), rimnem ( irbks2 (1) )
            CALL DC_WLOG  ( loglvl, 'DC', 2, logmsg (1:60), ierwlg )
            WRITE ( UNIT = logmsg, FMT = '( 6A10 )' )
     +          rimnem ( irqcb2 (1) ), rimnem ( irelb1 (1) ),
     +          rimnem ( irhinb (1) ), rimnem ( irfos2 (1) ),
     +          rimnem ( irelb2 (1) ), rimnem ( irnvpb (1) )
            CALL DC_WLOG  ( loglvl, 'DC', 2, logmsg (1:60), ierwlg )
            WRITE ( UNIT = logmsg, FMT = '( 6F10.2 )' )
     +          rivals ( irnbks ),     rivals ( irtobd (1) ),
     +          rivals ( irqcb1 (1) ), rivals ( irbks1 (1) ),
     +          rivals ( irfos1 (1) ), rivals ( irbks2 (1) )
            CALL DC_WLOG  ( loglvl, 'DC', 2, logmsg (1:60), ierwlg )
            WRITE ( UNIT = logmsg, FMT = '( 6F10.2 )' )
     +          rivals ( irqcb2 (1) ), rivals ( irelb1 (1) ),
     +          rivals ( irhinb (1) ), rivals ( irfos2 (1) ),
     +          rivals ( irelb2 (1) ), rivals ( irnvpb (1) )
            CALL DC_WLOG  ( loglvl, 'DC', 2, logmsg (1:60), ierwlg )
C
            IF  ( rivals ( irnbks ) .gt. 1 )  THEN
                 DO jj = 2, rivals ( irnbks )
                     WRITE ( UNIT = logmsg,
     +                      FMT = '( 10X, 5F10.2 )' )
     +               rivals ( irtobd ( jj ) ), rivals ( irqcb1 ( jj ) ),
     +               rivals ( irbks1 ( jj ) ), rivals ( irfos1 ( jj ) ),
     +               rivals ( irbks2 ( jj ) )
                     CALL DC_WLOG  ( loglvl, 'DC', 2, logmsg (1:60),
     +                               ierwlg )
                     WRITE ( UNIT = logmsg,
     +                      FMT = '( 6F10.2 )' )
     +               rivals ( irqcb2 ( jj ) ), rivals ( irelb1 ( jj ) ),
     +               rivals ( irhinb ( jj ) ), rivals ( irfos2 ( jj ) ),
     +               rivals ( irelb2 ( jj ) ), rivals ( irnvpb ( jj ) )
                     CALL DC_WLOG  ( loglvl, 'DC', 2, logmsg (1:60),
     +                               ierwlg )
                 END DO
            END IF
C
        END IF
C
C*      Print mean frequency, brightness temperature, and assoc. parameters.
C
        IF  ( rivals ( irnmef ) .gt. 0 )  THEN
C
             WRITE ( UNIT = logmsg, FMT = '( 4A10 )' )
     +          rimnem ( irnmef ),
     +          rimnem ( irmefr (1) ),
     +          rimnem ( irqctm (1) ), rimnem ( irtmbr (1) )
             CALL DC_WLOG  ( loglvl, 'DC', 2, logmsg (1:40), ierwlg )
             WRITE ( UNIT = logmsg, FMT = '( 4F10.2 )' )
     +          rivals ( irnmef ),
     +          rivals ( irmefr (1) ),
     +          rivals ( irqctm (1) ), rivals ( irtmbr (1) )
            CALL DC_WLOG  ( loglvl, 'DC', 2, logmsg (1:40), ierwlg )
C
            IF  ( rivals ( irnmef ) .gt. 1 )  THEN
                 DO jj = 2, rivals ( irnmef )
                     WRITE ( UNIT = logmsg,
     +                      FMT = '( 10X, 3F10.2 )' )
     +               rivals ( irmefr ( jj ) ),
     +               rivals ( irqctm ( jj ) ), rivals ( irtmbr ( jj ) )
                    CALL DC_WLOG  ( loglvl, 'DC', 2, logmsg (1:40),
     +                              ierwlg )
                 END DO
            END IF
C
        END IF
C
C*      Print satellite derived wind calculation method and
C*      10 meter wind speed.
C
        IF  ( rivals ( irnsww ) .gt. 0 )  THEN
C
             WRITE ( UNIT = logmsg, FMT = '( 3A10 )' )
     +          rimnem ( irnsww ),
     +          rimnem ( irswcm (1) ), rimnem ( irws10 (1) )
             CALL DC_WLOG  ( loglvl, 'DC', 2, logmsg (1:30), ierwlg )
             WRITE ( UNIT = logmsg, FMT = '( 3F10.2 )' )
     +          rivals ( irnsww ),
     +          rivals ( irswcm (1) ), rivals ( irws10 (1) )
            CALL DC_WLOG  ( loglvl, 'DC', 2, logmsg (1:30), ierwlg )
C
            IF  ( rivals ( irnsww ) .gt. 1 )  THEN
                 DO jj = 2, rivals ( irnsww )
                     WRITE ( UNIT = logmsg,
     +                      FMT = '( 10X, 2F10.2 )' )
     +               rivals ( irswcm ( jj ) ), rivals ( irws10 ( jj ) )
                    CALL DC_WLOG  ( loglvl, 'DC', 2, logmsg (1:30),
     +                              ierwlg )
                 END DO
            END IF
C
        END IF
C
	RETURN
	END