text:
!##############################################################################
        subroutine boundu
!	Note 1: prescribed inflow always comes from East
!	Note 2: wall functions apply in ALL directions/if node inside vicous layer->no slip
!	Note 3: if LSM/LSMbase are selected, inflow/outflow BCs only applied BELOW surface/lid
!##############################################################################
        use vars
        use multidata
	  use mpi
	  use imb
        implicit none
        integer :: i,j,k,ib,ly,im,kk,jj
        integer :: is,ie,js,je,ks,ke,ktop
	  integer :: strlen,dummy
	  double precision :: Fwallu,dyy,dzz,dxx,up,ufric
	  character (LEN=29) :: filename,fileSEM
	  character (LEN=4) :: domain
	  character (LEN=5) :: name_end

        if (PERIODIC) call exchange_bc(1,pl_ex)
        do ly=0,pl_ex
        do ib=1,nbp
           dxx=dom(ib)%dx
           dyy=dom(ib)%dy
           dzz=dom(ib)%dz
           is=dom(ib)%isu; ie=dom(ib)%ieu
           js=dom(ib)%jsu; je=dom(ib)%jeu
           ks=dom(ib)%ksu; ke=dom(ib)%keu

! Boundary Conditions for u - i.e. shear at boundaries and/or diriclet conditions
!..............................................................................
!=== West ===> ..   4=wall  ..    1=Inflow	..	7=read inflow
!..............................................................................
        if (dom(ib)%iprev.lt.0) then

           if ((dom(ib)%bc_west.eq.4).or.(dom(ib)%bc_east.ge.61)) then
              do k=ks-1,ke+1; do j=js-1,je+1
                 dom(ib)%u(is-1-ly,j,k)= 0.0
              end do; end do

           else if (dom(ib)%bc_west.eq.1) then
              do k=ks-1,ke+1; do j=js-1,je+1
                if (L_LSM) then
			if (dom(ib)%phi(is,j,k) .ge. 0.d0) then
                    dom(ib)%u(is-1-ly,j,k)= ubulk
			else
                    dom(ib)%u(is-1-ly,j,k)= 0.d0
			end if
                elseif (L_LSMbase) then
			if (dom(ib)%zc(k) .le. length) then
                    dom(ib)%u(is-1-ly,j,k)= ubulk
			else
                    dom(ib)%u(is-1-ly,j,k)= 0.d0
			end if
		    else
                    dom(ib)%u(is-1-ly,j,k)= ubulk
		    end if
              end do; end do

          else if (dom(ib)%bc_west.eq.7) then					!Reading slices
        		write(name_end,'(I5)') ireadinlet
        		strlen=LEN(TRIM(ADJUSTL(name_end)))
       		name_end=REPEAT('0',(5-strlen))//
     &		TRIM(ADJUSTL(name_end))
			write(domain,'(I4)') dom_id(ib)
		  	strlen=LEN(TRIM(ADJUSTL(domain)))
		 	domain=REPEAT('0',(4-strlen))//
     &		TRIM(ADJUSTL(domain))
		    filename='inflow/Inlet_'//domain//'_'//name_end//'.dat'
		      open (unit=405, file=filename)
		      do k=ks-1,ke+1; do j=js-1,je+1
				read(405,*)dom(ib)%u(is-1-ly,j,k)
		      end do; end do
		      close (405)

          else if (dom(ib)%bc_west.eq.77) then 					!Mapping inflow
                write(domain,'(I4)') dom_id(ib)
                strlen=LEN(TRIM(ADJUSTL(domain)))
                domain=REPEAT('0',(4-strlen))//
     &                  TRIM(ADJUSTL(domain))
           filename='inflow/Mapping_'//domain//'.dat'
                open (unit=405, file=filename)
              do k=ks-1,ke+1; do j=js-1,je+1
                read(405,*)dom(ib)%u(is-1-ly,j,k)
              end do; end do
               close(405)

          else if (dom(ib)%bc_west.eq.8) then					!Reading SEM inlet
        		write(name_end,'(I5)') ireadinlet
        		strlen=LEN(TRIM(ADJUSTL(name_end)))
       		name_end=REPEAT('0',(5-strlen))//
     &			TRIM(ADJUSTL(name_end))
			write(domain,'(I4)') dom_id(ib)
        		strlen=LEN(TRIM(ADJUSTL(domain)))
       		domain=REPEAT('0',(4-strlen))//
     &			TRIM(ADJUSTL(domain))
		    fileSEM='inflow/Inlet_'//domain//'_'//name_end//'.dat'
            	open (unit=405, file=fileSEM)
			read(405,*)
			read(405,*)
              do k=ks-1,ke+1; do j=js-1,je+1
	   	 	read(405,*)up,dummy,dummy

	  IF (UPROF.eq.12) then 							!1/7th power law inlet condition
	   if (dom(ib)%yc(j).lt.((yen-yst)/2)) then
          dom(ib)%u(is-1-ly,j,k) = ubulk*(1.0d0+1.0d0/7.0d0)
     &	   *(DABS(2*dom(ib)%yc(j)/(yen-yst)))**(1./7.)
	   else
          dom(ib)%u(is-1-ly,j,k) =  ubulk*(1.0d0+1.0d0/7.0d0)
     &    *(DABS(2*((yen-yst)-dom(ib)%yc(j))/(yen-yst)))**(1.d0/7.d0)
	   endif
          dom(ib)%u(is-1-ly,j,k) = dom(ib)%u(is-1-ly,j,k)*(1.0d0+1./7.)
     &	   *(DABS(dom(ib)%zc(k)/(zen-zst)))**(1./7.) +up

	   ELSE IF (UPROF.eq.13) then 						!1/7th power law inlet condition. Only HORIZONTAL
	  if (dom(ib)%yc(j).lt.((yen-yst)/2)) then
          dom(ib)%u(is-1-ly,j,k) = ubulk*(1.0d0+1.0d0/7.0d0)
     &	   *(DABS(2*dom(ib)%yc(j)/(yen-yst)))**(1./7.)
	  else
          dom(ib)%u(is-1-ly,j,k) = ubulk*(1.0d0+1.0d0/7.0d0)
     &    *(DABS(2*((yen-yst)-dom(ib)%yc(j))/(yen-yst)))**(1.d0/7.d0)
	  endif
          dom(ib)%u(is-1-ly,j,k) = dom(ib)%u(is-1-ly,j,k)+up

	   ELSE IF (UPROF.eq.14) then 						!1/7th power law Only vertical
          dom(ib)%u(is-1-ly,j,k) = ubulk*(1.0d0+1./7.)
     &	   *(DABS(dom(ib)%zc(k)/(zen-zst)))**(1./7.)+up

	   ELSE IF (UPROF.eq.15) then 						!Logarithmic distribution
          dom(ib)%u(is-1-ly,j,k) = 0.0187d0*
     &  (1.d0/0.41d0*DLOG(DABS(dom(ib)%zc(k))*0.0187d0*10**6)+5.d0)+up

	  else
          dom(ib)%u(is-1-ly,j,k) = ubulk+up
	  endif
               end do ; end do
			close(405)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	   else if (dom(ib)%bc_west.eq.12) then 					!1/7th power law inlet condition
            do k=ks-1,ke+1; do j=js-1,je+1
	  if (dom(ib)%yc(j).lt.((yen-yst)/2)) then
          dom(ib)%u(is-1-ly,j,k) = ubulk*(1.0d0+1.0d0/7.0d0)
     &	   *(DABS(2*dom(ib)%yc(j)/(yen-yst)))**(1./7.)
	  else
          dom(ib)%u(is-1-ly,j,k) = ubulk*(1.0d0+1.0d0/7.0d0)
     &    *(DABS(2*((yen-yst)-dom(ib)%yc(j))/(yen-yst)))**(1.d0/7.d0)
	  endif
          dom(ib)%u(is-1-ly,j,k) = dom(ib)%u(is-1-ly,j,k)*(1.0d0+1./7.)
     &	   *(DABS(dom(ib)%zc(k)/(zen-zst)))**(1./7.)
	     enddo ; end do

	   else if (dom(ib)%bc_west.eq.13) then 					!1/7th power law inlet condition. Only horizontal
            do k=ks-1,ke+1; do j=js-1,je+1
	  if (dom(ib)%yc(j).lt.((yen-yst)/2)) then
          dom(ib)%u(is-1-ly,j,k) = ubulk*(1.0d0+1.0d0/7.0d0)
     &	   *(DABS(2*dom(ib)%yc(j)/(yen-yst)))**(1./7.)
	  else
          dom(ib)%u(is-1-ly,j,k) = ubulk*(1.0d0+1.0d0/7.0d0)
     &    *(DABS(2*((yen-yst)-dom(ib)%yc(j))/(yen-yst)))**(1.d0/7.d0)
	  endif
	     enddo ; end do

	   else if (dom(ib)%bc_west.eq.14) then 					!1/7th power law inlet condition. Only vertical
            do k=ks-1,ke+1; do j=js-1,je+1
          dom(ib)%u(is-1-ly,j,k) = ubulk*(1.0d0+1./7.)
     &	   *(DABS(dom(ib)%zc(k)/(zen-zst)))**(1./7.)
	     enddo ; end do

          else if (dom(ib)%bc_west.eq.15) THEN					!Logarithmic distribution
              do k=ks-1,ke+1; do j=js-1,je+1
                 dom(ib)%u(is-1-ly,j,k)=0.0187d0*
     &  (1.d0/0.41d0*DLOG(DABS(dom(ib)%zc(k))*0.0187d0*10**6)+5.d0)
              end do; end do

          else if (dom(ib)%bc_west.eq.17) THEN
              do k=ks-1,ke+1; do j=js-1,je+1
 			ufric=0.103594d0*ubulk+0.00568d0
                 dom(ib)%u(is-1-ly,j,k)=ufric*
     &  (1.d0/0.41d0*DLOG(DABS(dom(ib)%zc(k))*ufric*10**6))/LOG(10.d0)
              end do; end do

           end if
        end if
!...............................................................................
!=== East ===> ..  4=wall  ..   2=Outflow
!...............................................................................
        if (dom(ib)%inext.lt.0) then
           if ((dom(ib)%bc_east.eq.4).or.(dom(ib)%bc_east.ge.61)) then
              do k=ks-1,ke+1; do j=js-1,je+1
                 dom(ib)%u(ie+1+ly,j,k)= 0.0
              end do; end do

           else if (dom(ib)%bc_east.eq.2 .or.dom(ib)%bc_east.eq.21) then
              if (alfabc.eq.1)then
                do k=ks-1,ke+1; do j=js-1,je+1

                if (L_LSM) then

                    if (dom(ib)%bc_east.eq.2) then
                       dom(ib)%u(ie+1+ly,j,k)=dom(ib)%u(ie,j,k)
                    else
                       dom(ib)%u(ie+1+ly,j,k)=dom(ib)%uoo(ie+1+ly,j,k)
     & -dt*alfapr*(dom(ib)%uoo(ie+1+ly,j,k)-
     & dom(ib)%uoo(ie,j,k))/dom(ib)%dx
                    end if

		    else		!no LSM

                    if (dom(ib)%bc_east.eq.2) then
                       dom(ib)%u(ie+1+ly,j,k)=dom(ib)%u(ie,j,k)
                    else
                       dom(ib)%u(ie+1+ly,j,k)=dom(ib)%uoo(ie+1+ly,j,k)
     & -dt*alfapr*(dom(ib)%uoo(ie+1+ly,j,k)-
     & dom(ib)%uoo(ie,j,k))/dom(ib)%dx
                    end if

		    end if
                    end do; end do
		  endif
           end if
        end if
!...............................................................................
!=== South ===> ..   4=wall ..   3=Symmetry .. 6=Wall function
!...............................................................................
        if (dom(ib)%jprev.lt.0) then

           if (dom(ib)%bc_south.eq.4) then
              do k=ks-1,ke+1; do i=is-1,ie+1
                 dom(ib)%u(i,js-1-ly,k)= -dom(ib)%u(i,js+ly,k)
              end do; end do
           else if (dom(ib)%bc_south.eq.3) then
              do k=ks-1,ke+1; do i=is-1,ie+1
                 dom(ib)%u(i,js-1-ly,k)= dom(ib)%u(i,js+ly,k)
              end do; end do
	     else if (dom(ib)%bc_south.ge.61) then				!Wall functions
		if (ly.eq.0) then
	   		if (dom(ib)%bc_south.lt.63)
     &			call log_law(3,dom(ib)%bc_south)
			if (dom(ib)%bc_south.ge.63)
     &			call wall_function(3,dom(ib)%bc_south)
		      do k=ks-1,ke+1; do i=is-1,ie+1
				Fwallu = dom(ib)%tauws2(i,k)
     &					*dom(ib)%u(i,js,k)/dyy			!*Acell/Vcell
				Fwallu=SIGN(Fwallu,-dom(ib)%u(i,js,k))
				dom(ib)%u(i,js,k)=dom(ib)%u(i,js,k)
     &						+Fwallu*alfapr*dt
		           	dom(ib)%u(i,js-1,k)= dom(ib)%u(i,js,k)
		      end do; end do
		else
		      do k=ks-1,ke+1; do i=is-1,ie+1
		           	dom(ib)%u(i,js-1-ly,k)= dom(ib)%u(i,js,k)
		      end do; end do
		endif
	     endif
        end if
!.............................................................................
!=== North ===>  ..   4=wall ..   44=moving wall ..  3=Symmetry
!.............................................................................
        if (dom(ib)%jnext.lt.0) then
           if (dom(ib)%bc_north.eq.4) then
              do k=ks-1,ke+1; do i=is-1,ie+1
                 dom(ib)%u(i,je+1+ly,k)   = - dom(ib)%u(i,je-ly,k)
              end do; end do

           else if (dom(ib)%bc_north.eq.3) then
              do k=ks-1,ke+1; do i=is-1,ie+1
                 dom(ib)%u(i,je+1+ly,k)   =  dom(ib)%u(i,je-ly,k)
              end do; end do
	     else if (dom(ib)%bc_north.ge.61) then				!Wall functions
		if (ly.eq.0) then
	   		if (dom(ib)%bc_north.lt.63)
     &			call log_law(4,dom(ib)%bc_north)
			if (dom(ib)%bc_north.ge.63)
     &			call wall_function(4,dom(ib)%bc_north)
		      do k=ks-1,ke+1; do i=is-1,ie+1
				Fwallu = dom(ib)%tauwn2(i,k)
     &					*dom(ib)%u(i,je,k)/dyy			!*Acell/Vcell
				Fwallu=sign(Fwallu,-dom(ib)%u(i,je,k))
				dom(ib)%u(i,je,k)=dom(ib)%u(i,je,k)
     &						+Fwallu*alfapr*dt
		           	dom(ib)%u(i,je+1,k)= dom(ib)%u(i,je,k)
		      end do; end do
		else
		      do k=ks-1,ke+1; do i=is-1,ie+1
		           	dom(ib)%u(i,je+1+ly,k)= dom(ib)%u(i,je,k)
		      end do; end do
		endif
	     endif
        end if
!...............................................................................
!=== Bottom ===> ..   4=wall ..   3=Symmetry
!...............................................................................
        if (dom(ib)%kprev.lt.0) then
          if (dom(ib)%bc_bottom.eq.4) then
              do j=js-1,je+1; do i=is-1,ie+1
                 dom(ib)%u(i,j,ks-1-ly)= -dom(ib)%u(i,j,ks+ly)
              end do; end do

           else if (dom(ib)%bc_bottom.eq.1) then
              do j=js-1,je+1; do i=is-1,ie+1
                 dom(ib)%u(i,j,ks-1-ly)= 0.0
              end do; end do

           else if (dom(ib)%bc_bottom.eq.3) then
              do j=js-1,je+1; do i=is-1,ie+1
                 dom(ib)%u(i,j,ks-1-ly)= dom(ib)%u(i,j,ks+ly)
              end do; end do
	     else if (dom(ib)%bc_bottom.ge.61) then				!Wall functions Bruño2014
		if (ly.eq.0) then
	   		if (dom(ib)%bc_bottom.lt.63) then
     				call log_law(5,dom(ib)%bc_bottom)
			endif
			if (dom(ib)%bc_bottom.ge.63)
     &			call wall_function(5,dom(ib)%bc_bottom)
		      do j=js-1,je+1; do i=is-1,ie+1
				Fwallu = dom(ib)%tauwb2(i,j)
     &					*dom(ib)%u(i,j,ks)/dzz			!*Acell/Vcell??
				Fwallu = sign(Fwallu,-dom(ib)%u(i,j,ks))

				dom(ib)%u(i,j,ks)=dom(ib)%u(i,j,ks)
     &						+Fwallu*alfapr*dt

		           	dom(ib)%u(i,j,ks-1)= dom(ib)%u(i,j,ks)
		      end do; end do
		else
		      do j=js-1,je+1; do i=is-1,ie+1
		           	dom(ib)%u(i,j,ks-1-ly)= dom(ib)%u(i,j,ks)
		      end do; end do
		endif
           end if
        end if
!.............................................................................
!=== Top ===>  ..   4=wall ..     3=Symmetry
!.............................................................................
	  if (L_LSMbase) then
        if (dom(ib)%z(1).le.length .and. dom(ib)%z(ke).ge.length) then
	    k=1
          do while (dom(ib)%z(k).le.length)
	      k=k+1
	    end do
	    ktop=k-1
          do j=js-1,je+1; do i=is-1,ie+1
             dom(ib)%u(i,j,ktop) = dom(ib)%u(i,j,ktop-1) 	!THIS IS A SLIP CONDITION
             dom(ib)%u(i,j,ktop+1) = dom(ib)%u(i,j,ktop-2) 	!THIS IS A SLIP CONDITION
          enddo; enddo
	  !pablo:	03/19
	   	do k=ktop+2,ke	!Over the free-surface layer
		    do j=js-1,je+1; do i=is-1,ie+1
		       dom(ib)%u(i,j,k) = 0.0
		    enddo; enddo
		enddo
	  end if
        if (dom(ib)%z(1).ge.length) then
	    do k=1,ke
           do j=js-1,je+1; do i=is-1,ie+1
             dom(ib)%u(i,j,k) = 0.d0 !ENFORCE VELOCITIES TO ZERO IN THE AIR
          enddo; enddo; enddo
	  end if
	  end if

        if (dom(ib)%knext.lt.0) then
           if (dom(ib)%bc_top.eq.4) then
              do j=js-1,je+1; do i=is-1,ie+1
                 dom(ib)%u(i,j,ke+1+ly)   = - dom(ib)%u(i,j,ke-ly)
              end do; end do
           else if (dom(ib)%bc_top.eq.3) then
             do j=js-1,je+1; do i=is-1,ie+1
                 dom(ib)%u(i,j,ke+1+ly)= dom(ib)%u(i,j,ke-ly)
             end do; end do
	     else if (dom(ib)%bc_top.ge.61) then					!Wall functions Bruño2014
		if (ly.eq.0) then
	   		if (dom(ib)%bc_top.lt.63)
     &			call log_law(6,dom(ib)%bc_top)
			if (dom(ib)%bc_top.ge.63)
     &			call wall_function(6,dom(ib)%bc_top)
		      do j=js-1,je+1; do i=is-1,ie+1
				Fwallu = dom(ib)%tauwt2(i,j)
     &					*dom(ib)%u(i,j,ke)/dzz			!*Acell/Vcell
				Fwallu = sign(Fwallu,-dom(ib)%u(i,j,ke))
				dom(ib)%u(i,j,ke)=dom(ib)%u(i,j,ke)
     &						+Fwallu*alfapr*dt
                 		dom(ib)%u(i,j,ke+1) = dom(ib)%u(i,j,ke)
		      end do; end do
		else
		      do j=js-1,je+1; do i=is-1,ie+1
                 		dom(ib)%u(i,j,ke+1+ly) = dom(ib)%u(i,j,ke)
		      end do; end do
		endif
           end if

          if (trim(keyword).eq.'cavity') then
		  do j=js-1,je+1
		   do i=is-1,ie+1
!                 dom(ib)%u(i,j,ke-ly)  = 1.d0
                 dom(ib)%u(i,j,ke+1+ly)  = ubulk
              end do
		 	end do
	   	  endif
	  end if

        end do
        if(diff_sch.ne.3) call boundcoef(1)
        end do
        end subroutine
!#############################################################################
        subroutine boundv
!#############################################################################
        use vars
        use multidata
	  use mpi
	  use imb
        implicit none
        integer :: i,j,k,ib,ly,kk,jj
        integer :: is,ie,js,je,ks,ke,ktop
	  integer :: strlen
	  character (LEN=29):: filename,fileSEM
	  character (LEN=4) :: domain
	  character (LEN=5) :: name_end
	  double precision :: Fwallv,dxx,dzz,dummy,dyy

        if (PERIODIC) call exchange_bc(2,pl_ex)

        do ly=0,pl_ex
        do ib=1,nbp
	     dxx=dom(ib)%dx
	     dyy=dom(ib)%dy
	     dzz=dom(ib)%dz
           is=dom(ib)%isv; ie=dom(ib)%iev
           js=dom(ib)%jsv; je=dom(ib)%jev
           ks=dom(ib)%ksv; ke=dom(ib)%kev

! Boundary Conditions for V - i.e. shear at boundaries
!...............................................................................
!=== West ===>   4=wall   ..    1=Inflow
!...............................................................................
        if (dom(ib)%iprev.lt.0) then
           if (dom(ib)%bc_west.eq.4) then
              do k=ks-1,ke+1; do j=js-1,je+1
                 dom(ib)%v(is-1-ly,j,k)= -dom(ib)%v(is+ly,j,k)
              end do; end do

           else if (dom(ib)%bc_west.eq.1 .or. dom(ib)%bc_west.eq.12
     & .or. dom(ib)%bc_west.eq.14 .or. dom(ib)%bc_west.eq.15
     & .or. dom(ib)%bc_west.eq.17)then
              do k=ks-1,ke+1; do j=js-1,je+1
                 dom(ib)%v(is-1-ly,j,k)=0.0
              end do; end do

	     else if (dom(ib)%bc_west.ge.61) then					!Wall functions
		if (ly.eq.0) then
	   		if (dom(ib)%bc_west.lt.63)
     &			call log_law(1,dom(ib)%bc_west)
			if (dom(ib)%bc_west.ge.63)
     &			call wall_function(1,dom(ib)%bc_west)
		      do k=ks-1,ke+1; do j=js-1,je+1
				Fwallv = dom(ib)%tauww2(j,k)
     &					*dom(ib)%v(is,j,k)/dxx			!*Acell/Vcell
				Fwallv = sign(Fwallv,-dom(ib)%v(is,j,k))
				dom(ib)%v(is,j,k)=dom(ib)%v(is,j,k)
     &						+Fwallv*alfapr*dt
                 		dom(ib)%v(is-1,j,k)= dom(ib)%v(is,j,k)
		      end do; end do
		else
		      do k=ks-1,ke+1; do j=js-1,je+1
                 		dom(ib)%v(is-1-ly,j,k)= dom(ib)%v(is,j,k)
		      end do; end do
		endif

          else if (dom(ib)%bc_west.eq.7) then					!Rading slices
        		write(name_end,'(I5)') ireadinlet
        		strlen=LEN(TRIM(ADJUSTL(name_end)))
       		name_end=REPEAT('0',(5-strlen))//
     &		TRIM(ADJUSTL(name_end))
			write(domain,'(I4)') dom_id(ib)
		  	strlen=LEN(TRIM(ADJUSTL(domain)))
		 	domain=REPEAT('0',(4-strlen))//
     &		TRIM(ADJUSTL(domain))
		    filename='inflow/Inlet_'//domain//'_'//name_end//'.dat'
		      open (unit=405, file=filename)
		      do k=dom(ib)%ksu-1,dom(ib)%keu+1
			do j=dom(ib)%jsu-1,dom(ib)%jeu+1
				read(405,*)dummy
		      end do; end do
		      do k=ks-1,ke+1; do j=js-1,je+1
				read(405,*)dom(ib)%v(is-1-ly,j,k)
		      end do; end do
			close(405)

          else if (dom(ib)%bc_west.eq.77) then 					!Mapping inflow
                write(domain,'(I4)') dom_id(ib)
                strlen=LEN(TRIM(ADJUSTL(domain)))
                domain=REPEAT('0',(4-strlen))//
     &                  TRIM(ADJUSTL(domain))
           filename='inflow/Mapping_'//domain//'.dat'
                open (unit=405, file=filename)
		      do k=dom(ib)%ksu-1,dom(ib)%keu+1
			do j=dom(ib)%jsu-1,dom(ib)%jeu+1
				read(405,*)dummy
		      end do; end do
              do k=ks-1,ke+1; do j=js-1,je+1
			read(405,*)dom(ib)%v(is-1-ly,j,k)
              end do; end do
               close(405)

          else if (dom(ib)%bc_west.eq.8) then					!SEM
        		write(name_end,'(I5)') ireadinlet
        		strlen=LEN(TRIM(ADJUSTL(name_end)))
       		name_end=REPEAT('0',(5-strlen))//
     &			TRIM(ADJUSTL(name_end))
			write(domain,'(I4)') dom_id(ib)
        		strlen=LEN(TRIM(ADJUSTL(domain)))
       		domain=REPEAT('0',(4-strlen))//
     &			TRIM(ADJUSTL(domain))
		    fileSEM='inflow/Inlet_'//domain//'_'//name_end//'.dat'
            	open (unit=405, file=fileSEM)
			read(405,*)
			read(405,*)
              do k=ks-1,ke+1; do j=js-1,je+1
	   	   read(405,*)dummy,dom(ib)%v(is-1-ly,j,k),dummy
              end do; end do
			close(405)
           end if
        end if
!...............................................................................
!=== East ===>   4=wall   ..    2=Outflow
!...............................................................................
        if (dom(ib)%inext.lt.0) then
           if (dom(ib)%bc_east.eq.4) then
              do k=ks-1,ke+1; do j=js-1,je+1
                 dom(ib)%v(ie+1+ly,j,k)= -dom(ib)%v(ie-ly,j,k)
              end do; end do

           else if(dom(ib)%bc_east.eq.2 .or.dom(ib)%bc_east.eq.21) then
              if (alfabc.eq.1) then
                 do k=ks-1,ke+1; do j=js-1,je+1

                 if (L_LSM) then

			 if (dom(ib)%phi(ie,j,k) .ge. 0.0) then
                     if (dom(ib)%bc_east.eq.2) then
                       dom(ib)%v(ie+1+ly,j,k)= dom(ib)%v(ie,j,k)
                     else
                       dom(ib)%v(ie+1+ly,j,k)= dom(ib)%voo(ie+1+ly,j,k)
     & -dt*alfapr*(dom(ib)%voo(ie+1+ly,j,k)-
     & dom(ib)%voo(ie,j,k))/dom(ib)%dx
                     end if
			 else
			   dom(ib)%v(ie+1+ly,j,k)= dom(ib)%v(ie,j,k)
			 endif

		     else									!no LSM

                     if (dom(ib)%bc_east.eq.2) then
                       dom(ib)%v(ie+1+ly,j,k)= dom(ib)%v(ie,j,k)
                     else
                       dom(ib)%v(ie+1+ly,j,k)= dom(ib)%voo(ie+1+ly,j,k)
     & -dt*alfapr*(dom(ib)%voo(ie+1+ly,j,k)-
     & dom(ib)%voo(ie,j,k))/dom(ib)%dx
                     end if

		     endif
                 end do; end do
              end if
	     else if (dom(ib)%bc_east.ge.61) then					!Wall functions Bruño2014
		if (ly.eq.0) then
	   		if (dom(ib)%bc_east.lt.63)
     &			call log_law(2,dom(ib)%bc_east)
			if (dom(ib)%bc_east.ge.63)
     &			call wall_function(2,dom(ib)%bc_east)
		      do k=ks-1,ke+1; do j=js-1,je+1
				Fwallv = dom(ib)%tauwe2(j,k)
     &					*dom(ib)%v(ie,j,k)/dxx			!*Acell/Vcell
				Fwallv = sign(Fwallv,-dom(ib)%v(ie,j,k))
				dom(ib)%v(ie,j,k)=dom(ib)%v(ie,j,k)
     &						+Fwallv*alfapr*dt
                 		dom(ib)%v(ie+1,j,k)= dom(ib)%v(ie,j,k)
		      end do; end do
		else
		      do k=ks-1,ke+1; do j=js-1,je+1
                 		dom(ib)%v(ie+1+ly,j,k)= dom(ib)%v(ie,j,k)
		      end do; end do
		endif
           end if
        end if
!...............................................................................
!=== South ===>   4=wall   ..    3=Symmetry
!...............................................................................
        if (dom(ib)%jprev.lt.0) then
           if (dom(ib)%bc_south.eq.3) then
              do k=ks-1,ke+1; do i=is-1,ie+1
                 dom(ib)%v(i,js-1-ly,k)=0.0
              end do; end do

	else if ((dom(ib)%bc_south.eq.4).or.(dom(ib)%bc_south.ge.61)) then
              do k=ks-1,ke+1; do i=is-1,ie+1
                 dom(ib)%v(i,js-1-ly,k)=0.0
              end do; end do
           end if
        end if
!...............................................................................
!=== North ===>   4=wall   ..    3=Symmetry
!...............................................................................
        if (dom(ib)%jnext.lt.0) then
           if (dom(ib)%bc_north.eq.3) then
              do k=ks-1,ke+1; do i=is-1,ie+1
                 dom(ib)%v(i,je+1+ly,k)=0.0
              end do; end do

	else if ((dom(ib)%bc_north.eq.4).or.(dom(ib)%bc_north.ge.61)) then
              do k=ks-1,ke+1; do i=is-1,ie+1
                 dom(ib)%v(i,je+1+ly,k)=0.0
              end do; end do
           end if
        end if
!...............................................................................
!=== Bottom ===> ..   4=wall ..   3=Symmetry
!...............................................................................
        if (dom(ib)%kprev.lt.0) then
          if (dom(ib)%bc_bottom.eq.4) then
              do j=js-1,je+1; do i=is-1,ie+1
                 dom(ib)%v(i,j,ks-1-ly)= -dom(ib)%v(i,j,ks+ly)
              end do; end do

           else if (dom(ib)%bc_bottom.eq.1) then
              do j=js-1,je+1; do i=is-1,ie+1
                 dom(ib)%v(i,j,ks-1-ly)= 0.0
              end do; end do

           else if (dom(ib)%bc_bottom.eq.3) then
              do j=js-1,je+1; do i=is-1,ie+1
                 dom(ib)%v(i,j,ks-1-ly)= dom(ib)%v(i,j,ks+ly)
              end do; end do
	     else if (dom(ib)%bc_bottom.ge.61) then				!Wall functions Bruño2014
		if (ly.eq.0) then
	   		if (dom(ib)%bc_bottom.lt.63) then
     				call log_law(5,dom(ib)%bc_bottom)
			endif
			if (dom(ib)%bc_bottom.ge.63)
     &			call wall_function(5,dom(ib)%bc_bottom)
		      do j=js-1,je+1; do i=is-1,ie+1
				Fwallv = dom(ib)%tauwb2(i,j)
     &					*dom(ib)%v(i,j,ks)/dzz			!*Acell/Vcell
				Fwallv = sign(Fwallv,-dom(ib)%v(i,j,ks))
				dom(ib)%v(i,j,ks)=dom(ib)%v(i,j,ks)
     &						+Fwallv*alfapr*dt
                 		dom(ib)%v(i,j,ks-1)= dom(ib)%v(i,j,ks)
		      end do; end do
		else
		      do j=js-1,je+1; do i=is-1,ie+1
                 		dom(ib)%v(i,j,ks-1-ly)= dom(ib)%v(i,j,ks)
		      end do; end do
		endif
           end if
        end if
!...............................................................................
!=== Top ===>   4=wall   ..    3=Symmetry
!...............................................................................
	  if (L_LSMbase) then
        if (dom(ib)%z(1).le.length .and. dom(ib)%z(ke).ge.length) then
	    k=1
          do while (dom(ib)%z(k).le.length)
	      k=k+1
	    end do
	    ktop=k-1
          do j=js-1,je+1; do i=is-1,ie+1
            dom(ib)%v(i,j,ktop) = dom(ib)%v(i,j,ktop-1)
            dom(ib)%v(i,j,ktop+1) = dom(ib)%v(i,j,ktop-2)
	    enddo; enddo
	  end if
        if (dom(ib)%z(1).ge.length) then				! PABLO 03/19
	    do k=1,ke
          do j=js-1,je+1; do i=is-1,ie+1
            dom(ib)%v(i,j,k) = 0.0
	    enddo; enddo;enddo
	  end if
	  end if

        if (dom(ib)%knext.lt.0) then
           if (dom(ib)%bc_top.eq.4) then
              do j=js-1,je+1; do i=is-1,ie+1
                 dom(ib)%v(i,j,ke+1+ly)= -dom(ib)%v(i,j,ke-ly)
              end do; end do
           else if (dom(ib)%bc_top.eq.3) then
             do j=js-1,je+1; do i=is-1,ie+1
                 dom(ib)%v(i,j,ke+1+ly)= dom(ib)%v(i,j,ke-ly)
             end do; end do
	     else if (dom(ib)%bc_top.ge.61) then					!Wall functions Bruño2014
		if (ly.eq.0) then
	   		if (dom(ib)%bc_top.lt.63)
     &			call log_law(6,dom(ib)%bc_top)
			if (dom(ib)%bc_top.ge.63)
     &			call wall_function(6,dom(ib)%bc_top)

		      do j=js-1,je+1; do i=is-1,ie+1
				Fwallv = dom(ib)%tauwt2(i,j)
     &					*dom(ib)%v(i,j,ke)/dzz			!*Acell/Vcell
				Fwallv = sign(Fwallv,-dom(ib)%v(i,j,ke))
				dom(ib)%v(i,j,ke)=dom(ib)%v(i,j,ke)
     &						+Fwallv*alfapr*dt
                 		dom(ib)%v(i,j,ke+1)= dom(ib)%v(i,j,ke)
		      end do; end do
		else
		      do j=js-1,je+1; do i=is-1,ie+1
                 		dom(ib)%v(i,j,ke+1+ly)= dom(ib)%v(i,j,ke)
		      end do; end do
		endif
           end if !bc_top.eq.4

	     if (trim(keyword).eq.'cavity') then
		  do j=js-1,je+1 ; do i=is-1,ie+1
                dom(ib)%v(i,j,ke+1+ly) = dom(ib)%v(i,j,ke-ly)
          end do ; end do
	   	 endif
	  end if

        end do
        if(diff_sch.ne.3) call boundcoef(2)
        end do
        end subroutine
!#############################################################################
        subroutine boundw
!#############################################################################
        use vars
        use multidata
	  use mpi
	  use imb
        implicit none
        integer :: i,j,k,ib,ly,kk,jj,is,ie,js,je,ks,ke,ktop,strlen
	  character (LEN=29):: filename,fileSEM
	  character (LEN=4) :: domain
	  character (LEN=5) :: name_end
	  double precision :: Fwallw,dxx,dyy,dummy,dzz

        if (PERIODIC) call exchange_bc(3,pl_ex)

        do ly=0,pl_ex
        do ib=1,nbp
	     dxx=dom(ib)%dx
	     dyy=dom(ib)%dy
	     dzz=dom(ib)%dz
           is=dom(ib)%isw; ie=dom(ib)%iew
           js=dom(ib)%jsw; je=dom(ib)%jew
           ks=dom(ib)%ksw; ke=dom(ib)%kew

! Boundary Conditions for W - i.e. shear at boundaries
!...............................................................................
!=== West ===>   4=wall   ..    1=Inflow
!...............................................................................
        if (dom(ib)%iprev.lt.0) then
           if (dom(ib)%bc_west.eq.4) then
              do k=ks-1,ke+1; do j=js-1,je+1
                 dom(ib)%w(is-1-ly,j,k)= -dom(ib)%w(is+ly,j,k)
              end do; end do

           else if (dom(ib)%bc_west.eq.1 .or. dom(ib)%bc_west.eq.12
     &  .or. dom(ib)%bc_west.eq.14 .or. dom(ib)%bc_west.eq.15
     &  .or. dom(ib)%bc_west.eq.17)then
              do k=ks-1,ke+1; do j=js-1,je+1
                 dom(ib)%w(is-1-ly,j,k)= 0.d0
              end do; end do
	     else if (dom(ib)%bc_west.ge.61) then					!Wall functions Bruño2014
		if (ly.eq.0) then
	   		if (dom(ib)%bc_west.lt.63)
     &			call log_law(1,dom(ib)%bc_west)
			if (dom(ib)%bc_west.ge.63)
     &			call wall_function(1,dom(ib)%bc_west)

		      do k=ks-1,ke+1; do j=js-1,je+1
				Fwallw = dom(ib)%tauww2(j,k)
     &					*dom(ib)%w(is,j,k)/dxx			!*Acell/Vcell
				Fwallw = sign(Fwallw,-dom(ib)%w(is,j,k))
				dom(ib)%w(is,j,k)=dom(ib)%w(is,j,k)
     &						+Fwallw*alfapr*dt
                 		dom(ib)%w(is-1,j,k)= dom(ib)%w(is,j,k)
			enddo; enddo
		else
		      do k=ks-1,ke+1; do j=js-1,je+1
                 		dom(ib)%w(is-1-ly,j,k)= dom(ib)%w(is,j,k)
		      end do; end do
		endif
          else if (dom(ib)%bc_west.eq.7) then					!brunho2014 reading slices
        		write(name_end,'(I5)') ireadinlet
        		strlen=LEN(TRIM(ADJUSTL(name_end)))
       		name_end=REPEAT('0',(5-strlen))//
     &		TRIM(ADJUSTL(name_end))
			write(domain,'(I4)') dom_id(ib)
		  	strlen=LEN(TRIM(ADJUSTL(domain)))
		 	domain=REPEAT('0',(4-strlen))//
     &		TRIM(ADJUSTL(domain))
		    filename='inflow/Inlet_'//domain//'_'//name_end//'.dat'
		      open (unit=405, file=filename)
		      do k=dom(ib)%ksu-1,dom(ib)%keu+1
			do j=dom(ib)%jsu-1,dom(ib)%jeu+1
				read(405,*)dummy
		      end do; end do
		      do k=dom(ib)%ksv-1,dom(ib)%kev+1
			do j=dom(ib)%jsv-1,dom(ib)%jev+1
				read(405,*)dummy
		      end do; end do
		      do k=ks-1,ke+1; do j=js-1,je+1
				read(405,*)dom(ib)%w(is-1-ly,j,k)
			end do; end do
			close(405)

          else if (dom(ib)%bc_west.eq.77) then 					!Mapping inflow
                write(domain,'(I4)') dom_id(ib)
                strlen=LEN(TRIM(ADJUSTL(domain)))
                domain=REPEAT('0',(4-strlen))//
     &                  TRIM(ADJUSTL(domain))
           filename='inflow/Mapping_'//domain//'.dat'
                open (unit=405, file=filename)
		      do k=dom(ib)%ksu-1,dom(ib)%keu+1
			do j=dom(ib)%jsu-1,dom(ib)%jeu+1
				read(405,*)dummy
		      end do; end do
		      do k=dom(ib)%ksv-1,dom(ib)%kev+1
			do j=dom(ib)%jsv-1,dom(ib)%jev+1
				read(405,*)dummy
		      end do; end do
              do k=ks-1,ke+1; do j=js-1,je+1
			read(405,*)dom(ib)%w(is-1-ly,j,k)
              end do; end do
               close(405)

          else if (dom(ib)%bc_west.eq.8) then					!SEM
        		write(name_end,'(I5)') ireadinlet
        		strlen=LEN(TRIM(ADJUSTL(name_end)))
       		name_end=REPEAT('0',(5-strlen))//
     &			TRIM(ADJUSTL(name_end))
			write(domain,'(I4)') dom_id(ib)
        		strlen=LEN(TRIM(ADJUSTL(domain)))
       		domain=REPEAT('0',(4-strlen))//
     &			TRIM(ADJUSTL(domain))
		    fileSEM='inflow/Inlet_'//domain//'_'//name_end//'.dat'
            	open (unit=405, file=fileSEM)
			read(405,*)
			read(405,*)
              do k=ks-1,ke+1; do j=js-1,je+1
	         read(405,*)dummy,dummy,dom(ib)%w(is-1-ly,j,k)
              end do; end do
			close(405)

           end if
        end if
!...............................................................................
!=== East ===>   4=wall   ..    2=Outflow
!...............................................................................
        if (dom(ib)%inext.lt.0) then
           if (dom(ib)%bc_east.eq.4) then
              do k=ks-1,ke+1; do j=js-1,je+1
                 dom(ib)%w(ie+1+ly,j,k)= -dom(ib)%w(ie-ly,j,k)
              end do; end do

           elseif (dom(ib)%bc_east.eq.2 .or.dom(ib)%bc_east.eq.21) then
              if (alfabc.eq.1) then
                 do k=ks-1,ke+1; do j=js-1,je+1

                   if (L_LSM) then

			   if (dom(ib)%phi(ie,j,k) .ge. 0.0) then
                       if (dom(ib)%bc_east.eq.2) then
                       dom(ib)%w(ie+1+ly,j,k)= dom(ib)%w(ie,j,k)
                       else
                       dom(ib)%w(ie+1+ly,j,k)= dom(ib)%woo(ie+1+ly,j,k)
     & -dt*alfapr*(dom(ib)%woo(ie+1+ly,j,k)-
     & dom(ib)%woo(ie,j,k))/dom(ib)%dx
                       end if
			   else
                       dom(ib)%w(ie+1+ly,j,k)= dom(ib)%w(ie,j,k)
			   endif

			 else	!no LSM
                       if (dom(ib)%bc_east.eq.2) then
                       dom(ib)%w(ie+1+ly,j,k)= dom(ib)%w(ie,j,k)
                       else
                       dom(ib)%w(ie+1+ly,j,k)= dom(ib)%woo(ie+1+ly,j,k)
     & -dt*alfapr*(dom(ib)%woo(ie+1+ly,j,k)-
     & dom(ib)%woo(ie,j,k))/dom(ib)%dx
                       end if
			 endif

                 end do; end do
              end if
	     else if (dom(ib)%bc_east.ge.61) then					!Wall functions Bruño2014
		if (ly.eq.0) then
	   		if (dom(ib)%bc_east.lt.63)
     &			call log_law(2,dom(ib)%bc_east)
			if (dom(ib)%bc_east.ge.63)
     &			call wall_function(2,dom(ib)%bc_east)

		      do k=ks-1,ke+1; do j=js-1,je+1
				Fwallw = dom(ib)%tauwe2(j,k)
     &					*dom(ib)%w(ie,j,k)/dxx			!*Acell/Vcell
				Fwallw = sign(Fwallw,-dom(ib)%w(ie,j,k))
				dom(ib)%w(ie,j,k)=dom(ib)%w(ie,j,k)
     &						+Fwallw*alfapr*dt
                 		dom(ib)%w(ie+1,j,k)= dom(ib)%w(ie,j,k)
		      end do; end do
		else
		      do k=ks-1,ke+1; do j=js-1,je+1
                 		dom(ib)%w(ie+1+ly,j,k)= dom(ib)%w(ie,j,k)
		      end do; end do
		endif
           end if
        end if
!...............................................................................
!=== South ===> ..   4=wall ..   3=Symmetry
!...............................................................................
        if (dom(ib)%jprev.lt.0) then
           if (dom(ib)%bc_south.eq.4) then
              do k=ks-1,ke+1; do i=is-1,ie+1
                 dom(ib)%w(i,js-1-ly,k)= -dom(ib)%w(i,js+ly,k)
              end do; end do

           else if (dom(ib)%bc_south.eq.3) then
              do k=ks-1,ke+1; do i=is-1,ie+1
                 dom(ib)%w(i,js-1-ly,k)= dom(ib)%w(i,js+ly,k)
              end do; end do
	     else if (dom(ib)%bc_south.ge.61) then					!Wall functions Bruño2014
		if (ly.eq.0) then
	   		if (dom(ib)%bc_south.lt.63)
     &			call log_law(3,dom(ib)%bc_south)
			if (dom(ib)%bc_south.ge.63)
     &			call wall_function(3,dom(ib)%bc_south)
		      do k=ks-1,ke+1; do i=is-1,ie+1
				Fwallw = dom(ib)%tauws2(i,k)
     &					*dom(ib)%w(i,js,k)/dyy			!*Acell/Vcell
				Fwallw = sign(Fwallw,-dom(ib)%w(i,js,k))
				dom(ib)%w(i,js,k)=dom(ib)%w(i,js,k)
     &						+Fwallw*alfapr*dt
		           	dom(ib)%w(i,js-1,k)= dom(ib)%w(i,js,k)
		      end do; end do
		else
		      do k=ks-1,ke+1; do i=is-1,ie+1
		           	dom(ib)%w(i,js-1-ly,k)= dom(ib)%w(i,js,k)
		      end do; end do
		endif
           end if
        end if
!.............................................................................
!=== North ===>  ..   4=wall ..   44=moving wall ..  3=Symmetry
!.............................................................................
        if (dom(ib)%jnext.lt.0) then
           if (dom(ib)%bc_north.eq.4) then
              do k=ks-1,ke+1; do i=is-1,ie+1
                 dom(ib)%w(i,je+1+ly,k)   = - dom(ib)%w(i,je-ly,k)
              end do; end do
           else if (dom(ib)%bc_north.eq.3) then
              do k=ks-1,ke+1; do i=is-1,ie+1
                 dom(ib)%w(i,je+1+ly,k)   =  dom(ib)%w(i,je-ly,k)
              end do; end do
	     else if (dom(ib)%bc_north.ge.61) then					!Wall functions Bruño2014
		if (ly.eq.0) then
	   		if (dom(ib)%bc_north.lt.63)
     &			call log_law(4,dom(ib)%bc_north)
			if (dom(ib)%bc_north.ge.63)
     &			call wall_function(4,dom(ib)%bc_north)
		      do k=ks-1,ke+1; do i=is-1,ie+1
				Fwallw = dom(ib)%tauwn2(i,k)
     &					*dom(ib)%w(i,je,k)/dyy			!*Acell/Vcell
				Fwallw = sign(Fwallw,-dom(ib)%w(i,je,k))
				dom(ib)%w(i,je,k)=dom(ib)%w(i,je,k)
     &						+Fwallw*alfapr*dt
		           	dom(ib)%w(i,je+1,k)= dom(ib)%w(i,je,k)
		      end do; end do
		else
		      do k=ks-1,ke+1; do i=is-1,ie+1
				dom(ib)%w(i,je+1+ly,k)= dom(ib)%w(i,je,k)
		      end do; end do
		endif
           end if
        end if
!...............................................................................
!=== Bottom ===>   4=wall   ..    3=Symmetry
!...............................................................................
        if (dom(ib)%kprev.lt.0) then
           if (dom(ib)%bc_bottom.eq.3) then
              do j=js-1,je+1; do i=is-1,ie+1
                 dom(ib)%w(i,j,ks-1-ly)=0.0
              end do; end do

		else if ((dom(ib)%bc_bottom.eq.4).or.
     &		   (dom(ib)%bc_bottom.ge.61)) then
              do j=js-1,je+1; do i=is-1,ie+1
                 dom(ib)%w(i,j,ks-1-ly)=0.0
              end do; end do

           else if (dom(ib)%bc_bottom.eq.1) then
              do j=js-1,je+1; do i=is-1,ie+1
                 dom(ib)%w(i,j,ks-1-ly)=1.1
              end do; end do
           end if
        end if
!...............................................................................
!=== Top ===>   4=wall   ..    3=Symmetry
!...............................................................................
        if (L_LSMbase) then
        if (dom(ib)%z(1).le.length .and. dom(ib)%z(ke).ge.length) then
          k=1
          do while (dom(ib)%z(k).le.length)
            k=k+1
	    enddo
	    ktop=k-1
          do j=js-1,je+1; do i=is-1,ie+1
             dom(ib)%w(i,j,ktop) = dom(ib)%w(i,j,ktop-1)
	       dom(ib)%w(i,j,ktop+1) = dom(ib)%w(i,j,ktop-2)
!	       dom(ib)%w(i,j,ktop+1) = 0.0				!Pablo 03/19
	    enddo; enddo
	  end if
        if (dom(ib)%z(1).ge.length) then
	   do k=1,ke
          do j=js-1,je+1; do i=is-1,ie+1
	       dom(ib)%w(i,j,k) = 0.0
	    enddo; enddo;enddo
	  end if
	  end if
        if (dom(ib)%knext.lt.0) then
           if (dom(ib)%bc_top.eq.3) then
             do j=js-1,je+1; do i=is-1,ie+1
                 dom(ib)%w(i,j,ke+1+ly)=0.0
             end do; end do
           else if ((dom(ib)%bc_top.eq.4).or.
     &				(dom(ib)%bc_top.ge.61)) then
              do j=js-1,je+1; do i=is-1,ie+1
                 dom(ib)%w(i,j,ke+1+ly)=0.0
              end do; end do
           end if

	  	  if (trim(keyword).eq.'cavity') then
		  do j=js-1,je+1 ; do i=is-1,ie+1
                 dom(ib)%w(i,j,ke+1+ly)=0.0
          end do ; end do
	   	 endif
	  end if

        end do
        if(diff_sch.ne.3) call boundcoef(3)
        end do
        end subroutine
!##############################################################################
        subroutine boundcoef(op)
!##############################################################################
        use vars
        use multidata
        implicit none
        integer :: i,j,k,ib,op,is,ie,js,je,ks,ke
        double precision d,dxx,dyy,dzz
        double precision, pointer, dimension(:,:,:) :: fi

        do ib=1,nbp
           dxx=dom(ib)%dx*dom(ib)%dx
           dyy=dom(ib)%dy*dom(ib)%dy
           dzz=dom(ib)%dz*dom(ib)%dz

        select case (op)
           Case (1)
           is=dom(ib)%isu; ie=dom(ib)%ieu
           js=dom(ib)%jsu; je=dom(ib)%jeu
           ks=dom(ib)%ksu; ke=dom(ib)%keu
           fi => dom(ib)%u
           Case (2)
           is=dom(ib)%isv; ie=dom(ib)%iev
           js=dom(ib)%jsv; je=dom(ib)%jev
           ks=dom(ib)%ksv; ke=dom(ib)%kev
           fi => dom(ib)%v
           Case (3)
           is=dom(ib)%isw; ie=dom(ib)%iew
           js=dom(ib)%jsw; je=dom(ib)%jew
           ks=dom(ib)%ksw; ke=dom(ib)%kew
           fi => dom(ib)%w
        end select

! Boundary Conditions for u - i.e. shear at boundaries and/or dirichlet conditions
!..............................................................................
!=== West ===> ..   4=wall  ..    1=Inflow
!..............................................................................
        if (dom(ib)%iprev.lt.0) then
           do k=ks,ke; do j=js,je
           d=fac*dom(ib)%vis(is,j,k)/dxx
           dom(ib)%ap(is,j,k)=dom(ib)%ap(is,j,k)+d
           if(diff_sch.eq.1) then
              dom(ib)%su(is,j,k)=dom(ib)%su(is,j,k)+d*fi(is-1,j,k)
           else
              dom(ib)%su(is,j,k)=dom(ib)%su(is,j,k)+2.0*d*fi(is-1,j,k)
           end if
           end do; end do
        end if
!...............................................................................
!=== East ===> ..  4=wall  ..   2=Outflow
!...............................................................................
        if (dom(ib)%inext.lt.0) then
           do k=ks,ke; do j=js,je
           d=fac*dom(ib)%vis(ie,j,k)/dxx
           dom(ib)%ap(ie,j,k)=dom(ib)%ap(ie,j,k)+d
           if(diff_sch.eq.1) then
              dom(ib)%su(ie,j,k)=dom(ib)%su(ie,j,k)+d*fi(ie+1,j,k)
           else
              dom(ib)%su(ie,j,k)=dom(ib)%su(ie,j,k)+2.0*d*fi(ie+1,j,k)
           end if
           end do; end do
        end if
!...............................................................................
!=== South ===> ..   4=wall ..   3=Symmetry ..
!...............................................................................
        if (dom(ib)%jprev.lt.0) then
           do k=ks,ke; do i=is,ie
           d=fac*dom(ib)%vis(i,js,k)/dyy
           dom(ib)%ap(i,js,k)=dom(ib)%ap(i,js,k)+d
           if(diff_sch.eq.1) then
              dom(ib)%su(i,js,k)=dom(ib)%su(i,js,k)+d*fi(i,js-1,k)
           else
              dom(ib)%su(i,js,k)=dom(ib)%su(i,js,k)+2.0*d*fi(i,js-1,k)
           end if
           end do; end do
        end if
!.............................................................................
!=== North ===>  ..   4=wall ..   44=moving wall ..  3=Symmetry
!.............................................................................
        if (dom(ib)%jnext.lt.0) then
           do k=ks,ke; do i=is,ie
           d=fac*dom(ib)%vis(i,je,k)/dyy
           dom(ib)%ap(i,je,k)=dom(ib)%ap(i,je,k)+d
           if(diff_sch.eq.1) then
              dom(ib)%su(i,je,k)=dom(ib)%su(i,je,k)+d*fi(i,je+1,k)
           else
              dom(ib)%su(i,je,k)=dom(ib)%su(i,je,k)+2.0*d*fi(i,je+1,k)
           end if
           end do; end do
        end if
!...............................................................................
!=== Bottom ===> ..   4=wall ..   3=Symmetry
!...............................................................................
        if (dom(ib)%kprev.lt.0) then
           do j=js,je; do i=is,ie
           d=fac*dom(ib)%vis(i,j,ks)/dzz
           dom(ib)%ap(i,j,ks)= dom(ib)%ap(i,j,ks) + d
           if(diff_sch.eq.1) then
              dom(ib)%su(i,j,ks)=dom(ib)%su(i,j,ks)+d*fi(i,j,ks-1)
           else
              dom(ib)%su(i,j,ks)=dom(ib)%su(i,j,ks)+2.0*d*fi(i,j,ks-1)
           end if
           end do; end do
        end if
!.............................................................................
!=== Top ===>  ..   4=wall ..     3=Symmetry
!.............................................................................
        if (dom(ib)%knext.lt.0) then
           do j=js,je; do i=is,ie
           d=fac*dom(ib)%vis(i,j,ke)/dzz
           dom(ib)%ap(i,j,ke)=dom(ib)%ap(i,j,ke)+d
           if(diff_sch.eq.1) then
              dom(ib)%su(i,j,ke)=dom(ib)%su(i,j,ke)+d*fi(i,j,ke+1)
           else
              dom(ib)%su(i,j,ke)=dom(ib)%su(i,j,ke)+2.0*d*fi(i,j,ke+1)
           end if
           end do; end do
        end if

        end do

        end subroutine

ext: for


text:
program euler_094
   use EulerCommon
   integer(int64) :: L = 10**9
   integer(int64) :: a, c, s, p, m
   integer(int64) :: new(4)

   a = 1
   c = 1
   s = 0
   p = 0
   m = 1
   ! use generating function n(i) = 4*n(i-1) - n(i-2); cf. wiki on Heronian Triangles
   do while (p < L)
      ! if the next line doesn't work, get a better compiler
      new = [c, 4*c - a + 2*m, -m, s + p]
      a = new(1); c = new(2); m = new(3); s = new(4)
      p = 3*c - m
   end do

   print *,s

end program euler_094ents here

max_stars_repo_name: kylekanos/ProjectEuler

max_issues_repo_name: kylekanos/ProjectEuler

max_forks_repo_name: kylekanos/ProjectEuler


text:
! NOMENCLATURE
!------------------------------------------
!nj = NUMBER OF NODES
!nm = NUMBER OF MEMBERS
!npj = NUMBER OF LOADED JOINTS
!nmp = NUMBER OF MATERIAL PROPERTIES
!npm = NUMBER OF LOADED MEMBERS
!nxsp = NUMBER OF CROSS-SECTIONS
!neq = NUMBER OF EQUATIONS TO BE SOLVED
!npc = NUMBER OF LOADED CASES
!x,y = CO-ORDINATES
!idf = DEGREE OF FREEDOM
!zi = MOMENT OF INERTIA
!em = ELASTIC MODULOUS
!fj = APPLIED LOADS MATRIX
!xsp = CROSS-SECTIONAL PROPERTIES
!nxsp = NUMBER OF CROSS-SECTIONAL PROPERTIES
!mcp = MEMBER CONNECTIVITY
!acs = DISPLACEMENT VECTOR
!im = MOMENT OF INERTIA
!nip = NUMBER OF MOMENT OF INERTIAS
!xdst = POINT LOAD DISTANCE FROM LEFT
!PL = POINT LOAD MAGNITUDE
!UDL = UDL MAGNITUDE
!**************************************************************************************************

program FrameAnalysis

real,dimension(:,:),allocatable::a,fj,acs,ds,xdst,PL,UDL
real,dimension(:,:,:),allocatable::Aml
real,dimension(:),allocatable::x,y,em,xsp,q,im
integer,dimension(:),allocatable::idf
integer,dimension(:,:),allocatable:: mcp

!To compile f95 to executable: gfortran -o program/<executable_filename> program/FrameAnalysis.f95

!OPEN INPUT FILE AS UNIT 5
open(5,file="program/data_in.txt")

!OPEN OUTPUT FILE AS UNIT 8
open(8,file="program/data_out.txt")

!OPEN OUTPUT FILE AS UNIT 9 (JSON string)
open(9,file="program/data_out.json")

write(9,*)'{'

!**************************************************************************************************
! READ INPUT DATA
!**************************************************************************************************

!(1)INITIALIZE DATA
read(5,*)nj,nm,nmp,nxsp,nip,npc
neq=nj*3
allocate (acs(neq,npc),q(neq),a(neq,neq),ds(neq,npc),fj(neq,npc))
allocate (xdst(nm,npc),PL(nm,npc),UDL(nm,npc),Aml(6,nm,npc))
!(2)JOINT DATA
allocate (idf(neq),x(nj),y(nj))
do k=1,nj
  read(5,*)j,x(j),y(j),idf(3*j-2),idf(3*j-1),idf(3*j)
end do
!(3)MEMBER PROPERTIES
allocate (em(nmp))
do k=1,nmp
  read(5,*)j,em(j)
end do
!(4)CROSS-SECTIONAL PROPERTIES
allocate (xsp(nxsp))
do k=1,nxsp
  read(5,*)j,xsp(j)
end do
!(5)MOMENT OF INERTIAS
allocate(im(nip))
do k=1,nip
  read(5,*)j,im(j)
end do
!(6)MEMBER DATA
allocate (mcp(nm,5))
do k=1,nm
  read(5,*)j,mcp(j,1:5)
end do
!SET BANDWIDTH
nbw=0
nbt=0
do i=1,nm
  nbt=mcp(i,2)*3-(mcp(i,1)*3-2)+1
  if (nbt>nbw) then
    nbw=nbt
  end if
end do
!**************************************************************************************************
! PROGRAM INFORMATION
!**************************************************************************************************
write(8,*)
write(8,500)
write(8,505)'PLANE FRAME ANALYSIS PROGRAM'
write(8,505)'DEVELOPED BY: FREDERIC LAVOIE'
write(8,500)
write(8,*)
write(8,*)
!**************************************************************************************************
! PRINTOUT INPUT DATA
!**************************************************************************************************
write(8,500)
write(8,507)'INPUT DATA'
write(8,500)
write(8,*)
write(8,501)'NUMBER OF JOINTS =',nj
write(8,501)'NUMBER OF MEMBERS =',nm
! write(8,501)'NUMBER OF LOAD CASES =',npc
write(8,*)
write(8,*)
write(8,*)'JOINT DATA'
write(8,500)
write(8,503)
write(8,504)
do j=1,nj
  write(8,502)j,x(j),y(j),idf(3*j-2),idf(3*j-1),idf(3*j)
end do
write(8,*)
write(8,*)
write(8,*)'MEMBER PROPERTY DATA'
write(8,500)
write(8,*)
write(8,509)
write(8,510)
do j=1,nmp
  write(8,506)j,em(j)
end do
write(8,*)
write(8,*)
write(8,*)'CROSS-SECTIONAL PROPERTY DATA'
write(8,500)
write(8,*)
write(8,511)
do j=1,nxsp
  write(8,512)j,xsp(j)
end do
write(8,*)
write(8,*)
write(8,*)'MOMENT OF INERTIAS'
write(8,500)
write(8,*)
write(8,525)
write(8,524)
do j=1,nip
  write(8,526)j,im(j)
end do
write(8,*)
write(8,*)
write(8,*)'MEMBER DATA'
write(8,500)
write(8,*)
write(8,513)
write(8,514)
do j=1,nm
  write(8,508)j,mcp(j,1:5)
end do
write(8,*)
write(8,*)
write(8,*)'JOINT LOADS'
write(8,500)
write(8,*)

!**************************************************************************************************
! CALL SUBROUTINES
!**************************************************************************************************

call fstiff(x,y,em,im,xsp,mcp,idf,nj,nm,neq,nbw,a,nxsp,nmp,nip)

fj(:,:)=0
xdst(:,:)=0
PL(:,:)=0
UDL(:,:)=0
do i=1,npc
  read(5,*)npj,npm
  do k=1,npj
    read(5,*)j,fj(j*3-2,i),fj(j*3-1,i),fj(j*3,i)
  end do
  do k=1,npm
    read(5,*)j,xdst(j,i),PL(j,i),UDL(j,i)
  end do
  if (npc>1) then
    write(8,518)'LOAD CASE',i
    write(8,515)
    write(8,516)
    do j=1,nj
      write(8,517)j,fj(j*3-2,i),fj(j*3-1,i),fj(j*3,i)
    end do
    write(8,528)
    write(8,529)
    do j=1,nm
      write(8,527)j,xdst(j,i),PL(j,i),UDL(j,i)
    end do
    write(8,*)
    else
    write(8,515)
    write(8,516)
    do j=1,nj
      write(8,517)j,fj(j*3-2,1),fj(j*3-1,1),fj(j*3,1)
    end do
    write(8,528)
    write(8,529)
    do j=1,nm
      write(8,527)j,xdst(j,1),PL(j,1),UDL(j,1)
    end do
    write(8,*)
  end if
end do

call loads(x,y,nm,em,xsp,im,nmp,nip,nxsp,npc,UDL,xdst,PL,neq,fj,idf,mcp,acs,Aml)

call eqn(a,acs,nbw,neq,npc)

!**************************************************************************************************
! PRINTOUT PRIMARY UNKNOWNS
!**************************************************************************************************
write(8,*)
write(8,500)
write(8,519)'PRIMARY UNKNOWNS'
write(8,500)
write(8,*)
write(8,521)
write(8,522)
write(8,523)
write(9,*)' "primaryUnknowns":'
do i=1,npc
  write(9,*)'   ['
  if (npc>1) then
    write(8,*)
    write(8,518)'LOAD CASE',i

    do j=1,nj
      write(8,520)j,acs(j*3-2,i),acs(j*3-1,i),acs(j*3,i)
      if (j==nj) then
        write(9,600)'      [',acs(j*3-2,i),',',acs(j*3-1,i),',',acs(j*3,i),']'
      else
        write(9,600)'      [',acs(j*3-2,i),',',acs(j*3-1,i),',',acs(j*3,i),'],'
      end if
    end do

    else
    do j=1,nj
      write(8,520)j,acs(j*3-2,1),acs(j*3-1,1),acs(j*3,1)
      if (j==nj) then
        write(9,600)'      [',acs(j*3-2,1),',',acs(j*3-1,1),',',acs(j*3,1),']'
      else
        write(9,600)'      [',acs(j*3-2,1),',',acs(j*3-1,1),',',acs(j*3,1),'],'
      end if
    end do
  end if
  write(9,*)'   ],'
end do
write(9,*)

!**************************************************************************************************
! PRINTOUT SECONDARY UNKNOWNS
!**************************************************************************************************
write(8,*)
write(8,500)
write(8,519)'SECONDARY UNKNOWNS'
write(8,500)
write(8,*)
write(8,*)

call mforce(x,y,em,xsp,mcp,idf,nj,nm,neq,nmp,npc,acs,fj,nxsp,im,nip,Aml)

!**************************************************************************************************
! FORMAT STATEMENTS
!**************************************************************************************************
500 format(71('-'))
501 format(a25,i3)
502 format(i3,4x,f8.1,2x,f8.1,4x,b2,2x,b2,2x,b2)
503 format('JOINT',7x,'COORDINATES',9x,'D.O.F')
504 format('NUMBER',7x,'X',7x,'Y',8x,'X',3x,'Y',3x,'R')
505 format(17x,a35)
506 format(i3,2x,es10.1)
507 format(17x,a25)
508 format(4x,i3,7x,i3,5x,i3,8x,i3,11x,i3,11x,i3)
509 format(8x,'ELASTIC')
510 format(2x,'#',5x,'MODULUS')
511 format(2x,'#',6x,'AREA')
512 format(i3,2x,es10.1)
513 format(3x,'MEMBER',5x,'MEMBER  NODES',4x,' MATERIAL ',4x,'X-SECTIONAL',4x,' INERTIA')
514 format(3x,'NUMBER',5x,'j-END   k-END',4x,'PROPERTY #',4x,'PROPERTY # ',4x,'PROPERTY #')
515 format(2x,'JOINT',23x,'JOINT LOADS')
516 format(2x,'NUMBER',13x,'X',13x,'Y',10x,'MOMENT')
517 format(2x,i3,8x,f10.2,4x,f10.2,4x,f10.2)
518 format(a10,i3)
519 format(14x,a35)
520 format(4x,i3,5x,es14.3,4x,es14.3,4x,es14.3)
521 format(3x,'  NODE  ',21x,'JOINT')
522 format(3x,' NUMBER ',20x,'DISPLACEMENT')
523 format(20x,'X',17x,'Y',14x,'ROTATION')
524 format(2x,'#',5x,'INERTIA')
525 format(8x,'MOMENT OF')
526 format(i3,3x,es10.1)
527 format(2x,i3,8x,f10.2,4x,es10.2,5x,es10.2)
528 format(2x,'MEMBER',21x,'MEMBER LOADS')
529 format(2x,'NUMBER',10x,'X-DIST.',10x,'PL',10x,'UDL')
600 format(A,F25.4,A,F25.4,A,F25.4,A)

write(9,*)'}'
!**************************************************************************************************
!**************************************************************************************************
!**************************************************************************************************
!                            SUBROUTINES
!**************************************************************************************************
!**************************************************************************************************
!**************************************************************************************************
contains
!**************************************************************************************************
!**************************************************************************************************
subroutine fstiff(x,y,em,im,xsp,mcp,idf,nj,nm,neq,nbw,a,nxsp,nmp,nip)
real,dimension(nmp)::em
real,dimension(nip)::im
real,dimension(neq,neq)::a
real,dimension(nj)::x,y
real,dimension(nxsp)::xsp
real,dimension(6,6)::sm,r,rt,smrt,sms
integer,dimension(neq)::idf
integer,dimension(nm,5)::mcp
integer,dimension(6)::id
real::iim

sm(:,:)=0
r(:,:)=0
rt(:,:)=0
smrt(:,:)=0
sms(:,:)=0
a(:,:)=0

do i=1,nm
  xj=x(mcp(i,1))
  yj=y(mcp(i,1))
  xk=x(mcp(i,2))
  yk=y(mcp(i,2))
  eem=em(mcp(i,3))
  xxsp=xsp(mcp(i,4))
  iim=im(mcp(i,5))

  call fmstiff(xj,yj,xk,yk,eem,xxsp,iim,sm)
  call frotmat(xj,yj,xk,yk,r,rt)
  call multmat(6,6,6,rt,sm,smrt)
  call multmat(6,6,6,smrt,r,sms)

  id(1)=3*mcp(i,1)-2
  id(2)=3*mcp(i,1)-1
  id(3)=3*mcp(i,1)
  id(4)=3*mcp(i,2)-2
  id(5)=3*mcp(i,2)-1
  id(6)=3*mcp(i,2)

  !ASSEMBLE GLOBAL STIFFNESS MATRIX
  do nr=1,6
    do nc=1,6
  !CHECK IF STIFFNESS COEFFICIENT IS ABOVE DIAGONAL, ONLY
  !STORE VALUES ON OR ABOVE
    kn=id(nr)-id(nc)
    if(kn<0) goto 110
    ic=id(nc)
    ir=kn+1
    a(ir,ic)=a(ir,ic)+sms(nr,nc)
110 end do
  end do
end do

!APPLY JOINT RESTRAINT
do j=1,neq
  ic=j
  if(idf(ic)>0)then
    a(1,ic)=1.0
!PUT ZEROS INTO ROW AND COLUMN
    ir=ic
    do m=2,nbw
      a(m,ic)=0.0
      ir=ir-1
      if(ir>0) then
        a(m,ir)=0.0
      end if
    end do
  end if
end do

end subroutine fstiff
!**************************************************************************************************
!**************************************************************************************************
subroutine fmstiff(xa,ya,xb,yb,em,xsa,im,sm)
real,dimension(6,6)::sm
real::im
sm(:,:)=0.0
!COMPUTE LENGTH
dx=xb-xa
dy=yb-ya
d=(dx**2+dy**2)**0.5
!COMPUTE AE/L
!COMPUTE EI/L
aed=xsa*em/d
sm(1,1)=aed
sm(4,4)=aed
sm(1,4)=-aed
sm(4,1)=-aed
sm(2,2)=12*em*im/(d**3)
sm(5,5)=12*em*im/(d**3)
sm(2,5)=-12*em*im/(d**3)
sm(5,2)=-12*em*im/(d**3)
sm(3,2)=6*em*im/(d**2)
sm(2,3)=6*em*im/(d**2)
sm(6,2)=6*em*im/(d**2)
sm(2,6)=6*em*im/(d**2)
sm(5,3)=-6*em*im/(d**2)
sm(3,5)=-6*em*im/(d**2)
sm(6,5)=-6*em*im/(d**2)
sm(5,6)=-6*em*im/(d**2)
sm(3,3)=4*em*im/d
sm(6,6)=4*em*im/d
sm(6,3)=2*em*im/d
sm(3,6)=2*em*im/d

end subroutine fmstiff
!**************************************************************************************************
!**************************************************************************************************
subroutine frotmat(xa,ya,xb,yb,r,rt)
real,dimension(6,6)::r,rt
r(:,:)=0.0
rt(:,:)=0.0
!COMPUTE LENGTH
dx=xb-xa
dy=yb-ya
d=(dx**2+dy**2)**0.5
!FORM ROTATION MATRIX
r(1,1)=dx/d
r(1,2)=dy/d
r(2,1)=dy/d*(-1)
r(2,2)=dx/d
r(3,3)=1
r(4,4)=dx/d
r(4,5)=dy/d
r(5,4)=dy/d*(-1)
r(5,5)=dx/d
r(6,6)=1
!TRANSPOSE
do jj=1,6
  do ii=1,6
    rt(jj,ii)=r(ii,jj)
  end do
end do
end subroutine frotmat
!**************************************************************************************************
!**************************************************************************************************
subroutine multmat(m,n,np,a,b,c)
real,dimension(m,n)::a
real,dimension(n,np)::b
real,dimension(m,np)::c
c(:,:)=0.0
do mm=1,m
  do ji=1,np
    do kk=1,n
      c(mm,ji)=c(mm,ji)+a(mm,kk)*b(kk,ji)
    end do
  end do
end do
end subroutine multmat
!**************************************************************************************************
!**************************************************************************************************
subroutine loads(x,y,nm,em,xsp,im,nmp,nip,nxsp,npc,UDL,xdst,PL,neq,fj,idf,mcp,acs,Aml)

real::a,b,px,py,iim
real,dimension(neq,npc)::fj,acs,aams
integer,dimension(neq)::idf
integer,dimension(nm,5)::mcp
integer,dimension(6)::id
real,dimension(nj)::x,y
real,dimension(nmp)::em
real,dimension(nxsp)::xsp
real,dimension(nip)::im
real,dimension(nm,npc)::xdst,PL,UDL
real,dimension(6,1)::ams,am
real,dimension(6,nm,npc)::Aml
real,dimension(6,6)::r,rt

acs(:,:)=0
aams(:,:)=0
Aml(:,:,:)=0

do i=1,npc
  do j=1,nm

    am(:,:)=0
    ams(:,:)=0
    id(:)=0
    r(:,:)=0
    rt(:,:)=0

    xj=x(mcp(j,1))
    yj=y(mcp(j,1))
    xk=x(mcp(j,2))
    yk=y(mcp(j,2))
    eem=em(mcp(j,3))
    xxsp=xsp(mcp(j,4))
    iim=im(mcp(j,5))

    dx=xk-xj
    dy=yk-yj
    d=(dx**2+dy**2)**0.5

    id(1)=3*mcp(j,1)-2
    id(2)=3*mcp(j,1)-1
    id(3)=3*mcp(j,1)
    id(4)=3*mcp(j,2)-2
    id(5)=3*mcp(j,2)-1
    id(6)=3*mcp(j,2)

    px=PL(j,i)*dy/d
    py=PL(j,i)*dx/d
    a=xdst(j,i)
    b=d-a
    Aml(1,j,i)=Aml(1,j,i)-(px*b/d)
    Aml(2,j,i)=Aml(2,j,i)-((py*b**2)*(3*a+b)/d**3)
    Aml(3,j,i)=Aml(3,j,i)-((py*b**2*a)/d**2)
    Aml(4,j,i)=Aml(4,j,i)-(px*a/d)
    Aml(5,j,i)=Aml(5,j,i)-((py*a**2)*(a+3*b)/d**3)
    Aml(6,j,i)=Aml(6,j,i)+((py*a**2*b)/d**2)

    wx=UDL(j,i)*dy/d
    wy=UDL(j,i)*dx/d
    Aml(1,j,i)=Aml(1,j,i)-(wx*D/2)
    Aml(2,j,i)=Aml(2,j,i)-(wy*D/2)
    Aml(3,j,i)=Aml(3,j,i)-((wy*d**2)/12)
    Aml(4,j,i)=Aml(4,j,i)-(wx*d/2)
    Aml(5,j,i)=Aml(5,j,i)-(wy*d/2)
    Aml(6,j,i)=Aml(6,j,i)+((wy*d**2)/12)

    call frotmat(xj,yj,xk,yk,r,rt)

    do k=1,6
      am(k,1)=Aml(k,j,i)
    end do

    call multmat(6,6,1,rt,am,ams)

    do k=1,6
      aams(id(k),i)=aams(id(k),i)+ams(k,1)
    end do

  end do
end do

do k=1,npc
  do j=1,neq
    acs(j,k)=fj(j,k)-aams(j,k)
!CHECK RESTRAINT CONDITIONS
    if(idf(j)>0)then
      acs(j,k)=0.0
    end if
  end do
end do
end subroutine loads
!**************************************************************************************************
!**************************************************************************************************
subroutine eqn(a,b,m,n,x)
!   call eqn(a,acs,nbw,neq,npc)

integer:: k,l,nm,n,mr,j,jl,i,m1,max,ii,i1,kk,min,m,z,x
real:: pivot,c,sum
real,dimension(n,n):: a
real,dimension(n,x):: b

nm=n-1
do 200 k=1,nm
  pivot=a(1,k)
  mr=min(m,n-k+1)

  do 180 l=2,mr
      c=a(l,k)/pivot
!            -   0   +
     if (c) 150,180,150
150    i=0
       j=k+l-1

       do 160 jl=l,mr
         i=i+1
160      a(i,j)=a(i,j)-a(jl,k)*c
180      a(l,k)=c

200  continue
!z=1
do z=1,x
  !FORWARD SUBSTITUTION  (ld) * y = b TO STORE y IN b
  m1=m-1
  b(1,z)=b(1,z)/a(1,1)
  do i=2,n
    l=max(1,i-m1)
    sum=0.0
    ii=i+1
    i1=i-1
    do k=l,i1
      j=ii-k
      sum=sum+a(j,k)*b(k,z)*a(1,k)
    end do
    b(i,z)=(b(i,z)-sum)/a(1,i)
  end do
  !BACK SUBSTITUTION  l tr * x = y  TO STORE  x IN b.
  do kk=1,nm
    i=n-kk
    i1=min(n,i+m1)
    sum=0.0
    l=i+1
    do k=l,i1
      j=k-i+1
      sum=sum+a(j,i)*b(k,z)
    end do
    b(i,z)=b(i,z)-sum
  end do
end do
end subroutine eqn
!**************************************************************************************************
!**************************************************************************************************
subroutine mforce(x,y,em,xsp,mcp,idf,nj,nm,neq,nmp,npc,acs,fj,nxsp,im,nip,Aml)
real,dimension(nj)::x,y
integer,dimension(6)::id
real,dimension(neq,npc)::fj,acs,dms,am,ams
real,dimension(nmp)::em
real,dimension(nip)::im
real,dimension(nxsp)::xsp
real,dimension(6,nm,npc)::Aml
integer,dimension(neq)::idf
integer,dimension(nm,5)::mcp
real,dimension(6,6)::sm,r,rt,smr
real::iim
real,dimension(neq)::ar

write(8,105)
write(8,*)
write(9,*)' "secondaryUnknowns":'
do k=1,npc
  write(9,111)'    ['
  ar(:)=0
  do i=1,nm

    am(:,:)=0
    dms(:,:)=0
    ams(:,:)=0
    r(:,:)=0
    rt(:,:)=0
    sm(:,:)=0
    smr(:,:)=0

    xj=x(mcp(i,1))
    yj=y(mcp(i,1))
    xk=x(mcp(i,2))
    yk=y(mcp(i,2))
    eem=em(mcp(i,3))
    xxsp=xsp(mcp(i,4))
    iim=im(mcp(i,5))
    call fmstiff(xj,yj,xk,yk,eem,xxsp,iim,sm)
    call frotmat(xj,yj,xk,yk,r,rt)
    call multmat(6,6,6,sm,r,smr)
    id(1)=3*mcp(i,1)-2
    id(2)=3*mcp(i,1)-1
    id(3)=3*mcp(i,1)
    id(4)=3*mcp(i,2)-2
    id(5)=3*mcp(i,2)-1
    id(6)=3*mcp(i,2)
    do j=1,6
      dms(j,1)=dms(j,1)+acs(id(j),k)
    end do
    call multmat(6,6,1,smr,dms,am)
!   PRINT {Am}i
    do j=1,6
      am(j,1)=am(j,1)+Aml(j,i,k)
    end do
    write(8,104)i,am(1,1),am(2,1),am(3,1)
    write(8,108)am(4,1),am(5,1),am(6,1)
    write(9,110)'      [',am(1,1),',',am(2,1),',',am(3,1),',',am(4,1),',',am(5,1),',',am(6,1),'],'

!   SUM THE MEMBER FORCES INTO REACTIONS
    call multmat(6,6,1,rt,am,ams)
!   STORE FIXED END FORCES FOR SUPPORTED JOINTS
!   {Ar}
    do j=1,6
      if(idf(id(j))>0) then
        ar(id(j))=ar(id(j))+ams(j,1)
      end if
    end do
  end do
! INCLUDE JOINT LOADS THAT OCCURE AT SUPPORTS

  do j=1,neq
    if(idf(j)>0) then
      ar(j)=ar(j)-fj(j,k)
    end if
  end do
! PRINT OUT REACTIONS
  write(8,*)
  write(8,101)
  write(8,*)
  write(8,102)
  write(8,103)
  do j=1,nj
    if(j==nj) then
      if((idf(j*3-2)+idf(j*3-1)+idf(j*3))>0) then
        write(8,100)j,ar(j*3-2),ar(j*3-1),ar(j*3)
      end if
      write(9,112)'      [',j,',',ar(j*3-2),',',ar(j*3-1),',',ar(j*3),']'
    else
      if((idf(j*3-2)+idf(j*3-1)+idf(j*3))>0) then
        write(8,100)j,ar(j*3-2),ar(j*3-1),ar(j*3)
      end if
      write(9,112)'      [',j,',',ar(j*3-2),',',ar(j*3-1),',',ar(j*3),'],'
    end if
  end do
  write(8,*)
  write(8,106)
  write(8,*)
  if (npc>1)then
    if (k<npc)then
      write(8,107)
      write(8,*)
    end if
  end if
  write(9,*)'   ]'
end do

100 format(4x,i3,5x,es12.3,4x,es12.3,4x,es12.3)
101 format(3x,'SUPPORT  REACTIONS')
102 format(3x,'JOINT')
103 format(3x,'NUMBER',9X,'X',15x,'Y',13x,'MOMENT')
104 format(3x,'MEMBER',1x,i3,4x,es12.3,4x,es12.3,4x,es12.3)
105 format(3x,'MEMBER END FORCES')
106 format(71('-'))
107 format(3x,'MEMBER END FORCES')
108 format(17x,es12.3,4x,es12.3,4x,es12.3)
110 format(A,F25.4,A,F25.4,A,F25.4,A,F25.4,A,F25.4,A,F25.4,A)
111 format(A,i3,A)
112 format(A,i4,A,F25.4,A,F25.4,A,F25.4,A)
end subroutine mforce
!**************************************************************************************************
!**************************************************************************************************
end program

max_stars_repo_name: FredLavoie/structural-analysis

max_issues_repo_name: FredLavoie/structural-analysis

max_forks_repo_name: FredLavoie/structural-analysis


text:
c=======================================================================
c
c     subroutine SCHUR
c
c     Matrix Schur Decomposition
c
c     If A is a complex square matrix of dimention n, then there exists
c     a unitary matrix Q such that
c
c     Q'AQ = T = D + N
c
c     where ' is the conjugate transpose,
c     T is an upper triangular matrix,
c     D = Diag(l1,...,lp) (the li are eigenvalues of A),
c     and N is strictly upper triangular matrix.
c
c     Furthermore, Q can be chosen such that the eigenvalues li appear
c     in any order along the diagonal.
c
c-----------------------------------------------------------------------
      SUBROUTINE schur ( n, X, iwork, dwork, eigval, vpmat, info )
c-----------------------------------------------------------------------
c
c     INPUT
c            n      : matrix size                                integer
c            X      : input matrix (n*n)                          double
c
c     WORKSPACE
c            iwork  : vector (n)                                 integer
c            dwork  : vector ( n*(n+6) )                          double
c
c     OUTPUT
c            eigval : eigenvalues (n)                             double
c            vpmat  : eigenvectors (n*n)                          double
c            info   : diagnostic argument                        integer
c
c     CALL
c            YM     : copy a vectorized matrix in a vectorized matrix
c            DGEES  : Schur factorization (cf. LAPACK)
c
c-----------------------------------------------------------------------
c
      IMPLICIT NONE
c
c     arguments i/o
      INTEGER n, info
      DOUBLE PRECISION eigval(*), X(*), vpmat(n,*)
c
c     workspaces
      INTEGER iwork(*)
      DOUBLE PRECISION dwork(*)
c
c     local variables
      LOGICAL select
      INTEGER sdim, piw, pmat, pdw, pip, lwork, ilwork
      PARAMETER ( ilwork = 5 )
c
c     external subroutines
      EXTERNAL YM, DGEES
c
c-----------------------------------------------------------------------
c
c     initialisations
      info = 0
c
c     lwork  : size of workspace for DGEES
c              ( min 3*N , so ilwork = min 3 , more is best )
      lwork = n*ilwork
      sdim = 0
c
c     pointers for integer work space : iwork
c     ---------------------------------------
      piw = 1
c     piw   : pointer for DGEES  who needs (n)
c     pinext = piw + n
c     pinext : pointer for the next iwork array
c
c     Total size of iwork array ( n )
c
c     pointers for double precision work space  : work(nasmax,*)
c     ----------------------------------------------------------
      pmat = 1
c     pmat   : pointer for work  matrix so (n*n) more
      pdw = pmat + ( n*n )
c     pdw    : pointer for work matrix of DGEES so lwork=(n*5) more
      pip = pdw + lwork
c     pip    : pointer for imaginary part vector so (n) more
c     pdnext = pip + n
c     pdnext  : pointer of the next dwork array
c
c     Total size of work array = ( n*(n + 6) )
c
c-----------------------------------------------------------------------
c
c     saves matrix initial ( altered by DGEES )
      CALL YM ( n, n, X, dwork(pmat) )
c
c     DGEES (cf. LAPACK)
c     DGEES computes for an N-by-N real  nonsymmetric  matrix  A,
c     the eigenvalues, the real Schur form T, and, optionally, the
c     matrix of Schur vectors Z
      CALL DGEES( 'V', 'N', select, n, dwork(pmat), n,
     &            sdim, eigval, dwork(pip), vpmat,
     &            n, dwork(pdw), lwork, iwork(piw), info )
c
      RETURN
      END

ext: f

max_stars_repo_name: alpgodev/numx

max_issues_repo_name: alpgodev/numx

max_forks_repo_name: alpgodev/numx


text:
      SUBROUTINE HCOMPT
C
      IMPLICIT REAL*8(A-H,O-Z)
C
      COMMON/ADDMOM/HAWH(3),HAXWH(3),HELGM(3),HTIP(3)
C
      COMMON/CONSTS/PI,TWOPI,RAD
C
      COMMON/CSTVAL/TST,DUM9(40)
C
      COMMON/EXPAND/THRMPR(100),ITHRM(20)
C
      COMMON/HGIMBL/HGMB(3)
C
      COMMON/HOUTPT/IHCALC,IHREF,IHFLAG
C
      COMMON/HVCOMP/YB(3),YBD(3),CIN(3,3),CID,CIND(3)
C
      COMMON/HVCOUT/HBODY(3),HINERT(3),HMAG
C
      COMMON/IPOOL1/IDUM5,IDAMP,IDUM6(37)
C
      COMMON/MOMENT/IDUM1,IV,IDUM2(2),IW,IDUM3
C
      COMMON/OUTONE/OML,DUM5(9)
C
      COMMON/HWHEEL/HWM(3)
C
      COMMON/ISECBD/I2BDY
C
      COMMON/RPOOL1/DUM1(10),T,A(3,3),F(3,3),DUM2(10),OM(3),DUM3(81),
     1              PHID,PHI
C
      COMMON/RPOOL3/ZMS,YIZM(3,2)
C
      COMMON/RPOOL6/FDUM(3,3),CIY(3,3),CIZ(3,3),SAT(3,3),SZ1,SZ2,SZ3
C
      COMMON/RVISCS/VIV(3),DUM6(12),CIL
C
      COMMON/SBDOUT/B2DUM(34),H2B(3),Y2D(3)
C
      COMMON/XIN1/DUM7(3),ETTA,ZETTA,IDUM4
C
      DIMENSION HDUM(3),V1(3),V2(3),V3(3)
C
      OM1=OM(1)
      OM2=OM(2)
      OM3=OM(3)
      Y1=YB(1)
      Y2=YB(2)
      Y3=YB(3)
      HBODY(1)=(CIN(1,1)-(Y2*Y2+Y3*Y3)/ZMS)*OM1+(CIN(1,2)+Y1*Y2/ZMS)*OM2
     1        +(CIN(1,3)+Y1*Y3/ZMS)*OM3
      HBODY(2)=(CIN(2,1)+Y2*Y1/ZMS)*OM1+(CIN(2,2)-(Y3*Y3+Y1*Y1)/ZMS)*OM2
     1        +(CIN(2,3)+Y2*Y3/ZMS)*OM3
      HBODY(3)=(CIN(3,1)+Y3*Y1/ZMS)*OM1+(CIN(3,2)+Y3*Y2/ZMS)*OM2
     1        +(CIN(3,3)-(Y1*Y1+Y2*Y2)/ZMS)*OM3
C
      HBODY(1)=HBODY(1)+CIND(1)-(Y2*YBD(3)-Y3*YBD(2))/ZMS
      HBODY(2)=HBODY(2)+CIND(2)-(Y3*YBD(1)-Y1*YBD(3))/ZMS
      HBODY(3)=HBODY(3)+CIND(3)-(Y1*YBD(2)-Y2*YBD(1))/ZMS
C
      IF(IDAMP.EQ.0) GO TO 10
      CALL MATV(2,F,YB,V1)
      CALL MATV(2,F,YIZM(1,1),V2)
      HDUM(1)=PHID*(-CIZ(2,1)-(V2(2)-V1(2)/ZMS)*SZ1)
      HDUM(2)=PHID*(CIZ(3,3)+CIZ(1,1)+(V2(3)-V1(3)/ZMS)*SZ3+
     1(V2(1)-V1(1)/ZMS)*SZ1)
      HDUM(3)=PHID*(-CIZ(2,3)-(V2(2)-V1(2)/ZMS)*SZ3)
      CALL MATV(1,F,HDUM,V3)
      DO 5 I=1,3
      HBODY(I)=HBODY(I)+V3(I)
    5 CONTINUE
   10 CONTINUE
C
      IF(IV.EQ.0) GO TO 20
      WS1=OML*CIL
      DO 15 I=1,3
      HBODY(I)=HBODY(I)+VIV(I)*WS1
   15 CONTINUE
   20 CONTINUE
C
      IF(IW.EQ.0) GO TO 30
      DO 25 I=1,3
      HBODY(I)=HBODY(I)+HWM(I)
   25 CONTINUE
   30 CONTINUE
C
      IF(I2BDY.EQ.0) GO TO 32
C
      HBODY(1)=HBODY(1)+H2B(1)-(Y2*Y2D(3)-Y3*Y2D(2))/ZMS
      HBODY(2)=HBODY(2)+H2B(2)-(Y3*Y2D(1)-Y1*Y2D(3))/ZMS
      HBODY(3)=HBODY(3)+H2B(3)-(Y1*Y2D(2)-Y2*Y2D(1))/ZMS
C
   32 CONTINUE
C
      HBODY(1)=HBODY(1)+HGMB(1)+HAWH(1)+HAXWH(1)+HELGM(1)+HTIP(1)
      HBODY(2)=HBODY(2)+HGMB(2)+HAWH(2)+HAXWH(2)+HELGM(2)+HTIP(2)
      HBODY(3)=HBODY(3)+HGMB(3)+HAWH(3)+HAXWH(3)+HELGM(3)+HTIP(3)
C
      HMAG=0.0D0
      DO 35 I=1,3
      HINERT(I)=A(I,1)*HBODY(1)+A(I,2)*HBODY(2)+A(I,3)*HBODY(3)
      HMAG=HMAG+HINERT(I)*HINERT(I)
   35 CONTINUE
      HMAG=DSQRT(HMAG)
      IF(T.NE.TST) GO TO 50
      IF(IHREF.EQ.0) GO TO 50
      WS1=HINERT(1)/HMAG
      WS2=HINERT(2)/HMAG
      ZETTA=0.0D0
      IF(WS1.EQ.0.0D0.AND.WS2.EQ.0.0D0) GO TO 40
      ZETTA=DATAN2(WS2,WS1)/RAD
   40 CONTINUE
      WS1=HINERT(3)/HMAG
      ETTA=DARCOS(WS1)/RAD
   50 CONTINUE
      IF(IHCALC.EQ.0) IHFLAG=0
      IF(ITHRM(10).EQ.0) RETURN
      WRITE(9,9000) CIN
      WRITE(9,9000) CIND,YB,YBD
      WRITE(9,9000) OM,HBODY,HINERT
 9000 FORMAT(' HCOMPT',5X,1P9E13.5)
      RETURN
      END


text:
program sample
    use SeedClass
    implicit none

    type(Seed_) :: dry_seed,wet_seed

    call dry_seed%init(radius=2.0d0,width1=9.70d0,width2=8.20d0,width3=7.70d0)
    !call wet_seed%init(radius=2.0d0,width1=14.0d0,width2=10.0d0,width3=8.70d0)
    call dry_seed%createMesh(FileName="/home/haruka/test/seed_dry",withSTL=.true.,ElemType="Tetrahedra")
    !call wet_seed%createMesh(FileName="/home/haruka/test/seed_wet",withSTL=.true.,ElemType="Tetrahedra")

end program

max_stars_repo_name: kazulagi/plantFEM

max_issues_repo_name: kazulagi/plantFEM

max_forks_repo_name: kazulagi/plantFEM


text:
!--------------------------------------------------------------------------------
! Copyright (c) 2017 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
! This file is part of FLEUR and available as free software under the conditions
! of the MIT license as expressed in the LICENSE file in more detail.
!--------------------------------------------------------------------------------

MODULE m_init_wannier_defaults

CONTAINS

SUBROUTINE initWannierDefaults(wann)

   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   !!!
   !!!  This subroutine sets most of the attributes of the t_wann
   !!!  type to standard values.
   !!!
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

   USE m_types_setup

   IMPLICIT NONE

   TYPE(t_wann), INTENT(INOUT) :: wann

   wann%wan90version = 3 ! Set the standard to Wannier90-1.2

   wann%oc_num_orbs = 0
!        integer,allocatable :: oc_orbs(:)

   wann%l_perpmagatlres= .FALSE.
   wann%l_unformatted = .FALSE.
   wann%l_oc_f = .FALSE.
   wann%l_ndegen = .FALSE.
   wann%l_orbitalmom = .FALSE.
   wann%l_orbcomp = .FALSE.
   wann%l_orbcomprs = .FALSE.
   wann%l_denmat = .FALSE.
   wann%l_perturbrs = .FALSE.
   wann%l_perturb = .FALSE.
   wann%l_nedrho = .FALSE.
   wann%l_anglmomrs = .FALSE.
   wann%l_anglmom = .FALSE.
   wann%l_spindisp = .FALSE.
   wann%l_spindisprs = .FALSE.
   wann%l_socspicom = .FALSE.
   wann%l_socspicomrs = .FALSE.
   wann%l_offdiposoprs = .FALSE.
   wann%l_offdiposop = .FALSE.
   wann%l_torque = .FALSE.
   wann%l_torquers = .FALSE.
   wann%l_atomlist = .FALSE.

   wann%atomlist_num = 0 ! has to be initialize to atoms%nat or something smaller at some point
!        integer,allocatable :: atomlist(:)

   wann%l_berry = .FALSE.
   wann%l_perpmagrs = .FALSE.
   wann%l_perpmag = .FALSE.
   wann%l_perpmagat = .FALSE.
   wann%l_perpmagatrs = .FALSE.
   wann%l_socmatrs = .FALSE.
   wann%l_socmat = .FALSE.
   wann%l_socmatvecrs = .FALSE.
   wann%l_socmatvec = .FALSE.
   wann%l_soctomom = .FALSE.
   wann%l_kptsreduc2 = .FALSE.
   wann%l_nablapaulirs = .FALSE.
   wann%l_nablars = .FALSE.
   wann%l_surfcurr = .FALSE.
   wann%l_updown = .FALSE.
   wann%l_ahe = .FALSE.
   wann%l_she = .FALSE.
   wann%l_rmat = .FALSE.
   wann%l_nabla = .FALSE.
   wann%l_socodi = .FALSE.
   wann%l_pauli = .FALSE.
   wann%l_pauliat = .FALSE.
   wann%l_potmat = .FALSE.
   wann%l_projgen = .FALSE.
   wann%l_plot_symm = .FALSE.
   wann%l_socmmn0 = .FALSE.
   wann%l_bzsym = .FALSE.
   wann%l_hopping = .FALSE.
   wann%l_kptsreduc = .FALSE.
   wann%l_prepwan90 = .FALSE.
   wann%l_plot_umdat = .FALSE.
   wann%l_wann_plot = .FALSE.
   wann%l_bynumber = .FALSE.
   wann%l_stopopt = .FALSE.
   wann%l_matrixmmn = .FALSE.
   wann%l_matrixamn = .FALSE.
   wann%l_projmethod = .FALSE.
   wann%l_wannierize = .FALSE.
   wann%l_plotw90 = .FALSE.
   wann%l_byindex = .FALSE.
   wann%l_byenergy = .FALSE.
   wann%l_proj_plot = .FALSE.
   wann%l_bestproj = .FALSE.
   wann%l_ikptstart = .FALSE.
   wann%l_lapw = .FALSE.
   wann%l_plot_lapw = .FALSE.
   wann%l_fermi = .FALSE.
   wann%l_dipole = .FALSE.
   wann%l_dipole2 = .FALSE.
   wann%l_dipole3 = .FALSE.
   wann%l_mmn0 = .FALSE.
   wann%l_mmn0at = .FALSE.
   wann%l_manyfiles = .FALSE.
   wann%l_collectmanyfiles = .FALSE.
   wann%l_ldauwan = .FALSE.
   wann%l_lapw_kpts = .FALSE.
   wann%l_lapw_gfleur = .FALSE.
   wann%l_kpointgen = .FALSE.
   wann%l_w90kpointgen = .FALSE.
   wann%l_finishnocoplot = .FALSE.
   wann%l_finishgwf = .FALSE.
   wann%l_skipkov = .FALSE.
   wann%l_matrixuHu = .FALSE.
   wann%l_matrixuHu_dmi = .FALSE.

   wann%l_mmn0_unf_to_spn_unf= .FALSE.
   wann%l_mmn0_to_spn_unf= .FALSE.
   wann%l_mmn0_to_spn= .FALSE.
   wann%l_mmn0_to_spn2= .FALSE.
   wann%l_mmn0_unf_to_spn= .FALSE.
   wann%l_perpmag_unf_to_tor_unf= .FALSE.
   wann%l_perpmag_to_tor_unf= .FALSE.
   wann%l_perpmag_to_tor= .FALSE.
   wann%l_perpmag_unf_to_tor= .FALSE.
   wann%l_hsomtxvec_unf_to_lmpzsoc_unf= .FALSE.
   wann%l_hsomtxvec_to_lmpzsoc_unf= .FALSE.
   wann%l_hsomtxvec_to_lmpzsoc= .FALSE.
   wann%l_hsomtxvec_unf_to_lmpzsoc= .FALSE.
   wann%l_hsomtx_unf_to_hsoc_unf= .FALSE.
   wann%l_hsomtx_to_hsoc_unf= .FALSE.
   wann%l_hsomtx_to_hsoc= .FALSE.
   wann%l_hsomtx_unf_to_hsoc= .FALSE.

   wann%ikptstart = 1
   wann%band_min(1:2) = -1
   wann%band_max(1:2) = -1
   wann%gfthick = 0
   wann%gfcut = 0
   wann%unigrid(6) = 0
   wann%mhp(3) = 0

!---> gwf
   wann%l_ms = .FALSE.
   wann%l_sgwf = .FALSE.
   wann%l_socgwf = .FALSE.
   wann%l_gwf = .FALSE.
   wann%l_bs_comf = .FALSE.
   wann%l_exist = .FALSE.
   wann%l_opened = .FALSE.
   wann%l_cleverskip = .FALSE.
   wann%l_dim(3) = .FALSE.

   wann%scale_param = 0.0
   wann%aux_latt_const = 0.0
   wann%hdwf_t1 = 0.0
   wann%hdwf_t2 = 0.0
   wann%nparampts = 0
   wann%fn_eig = ''
   wann%param_file = ''

   wann%scale_param = 1.0
   wann%aux_latt_const = 8.0!5.5!5.45886450 !5.98136400 !8.0725882513951497 !5.4170 !1.0
   wann%param_file='qpts'
   wann%l_dim=.false.

END SUBROUTINE initWannierDefaults

END MODULE m_init_wannier_defaults

max_stars_repo_name: MRedies/FLEUR

max_issues_repo_name: MRedies/FLEUR

max_forks_repo_name: MRedies/FLEUR


text:
MODULE embm_data

  USE embm_lib
  USE embm_netcdf
  IMPLICIT NONE

CONTAINS

  ! This module reads netcdf restarts for EMBM
  SUBROUTINE inm_netcdf_embm
    USE genie_global, ONLY: write_status, gui_restart, istep_atm
    IMPLICIT NONE

    REAL, DIMENSION(maxi,maxj) :: temp_read, shum_read
    INTEGER :: iday, ifail, ncid
    LOGICAL :: lexist
    CHARACTER(LEN=200) :: fnamein
    REAL :: timestep
    INTEGER :: i, j
    REAL :: tmp_val(2), area

    timestep = 24.0 * 60.0 * 60.0 * yearlen / REAL(nyear * ndta)
    IF (gui_restart) THEN
       PRINT *, 'READING EMBM GUI RESTART FILE: gui_restart_embm.nc'
       fnamein = 'gui_restart_embm.nc'
    ELSE
       fnamein = TRIM(filenetin)
    END IF
    ifail = 0
    INQUIRE(FILE=TRIM(fnamein), EXIST=lexist)
    IF (.NOT. lexist) THEN
       PRINT *, ' Missing file ', TRIM(fnamein)
       ifail = 1
    END IF
    IF (ifail /= 0) THEN
       PRINT *, ' Correct error and try again '
       CALL write_status('ERRORED')
    END IF

    PRINT *, ' embm: Opening restart file for read: ', &
         & TRIM(filenetin)

    CALL open_file_nc(TRIM(fnamein), ncid)
    IF (gui_restart) &
         & CALL getscalari_data_nc(ncid, 'istep_atm', istep_atm, ifail)
    CALL getscalari_data_nc(ncid, 'ioffset', ioffset_rest, ifail)
    CALL getscalari_data_nc(ncid, 'iyear', iyear_rest, ifail)
    CALL getscalari_data_nc(ncid, 'imonth', imonth_rest, ifail)
    CALL getscalari_data_nc(ncid, 'iday', iday, ifail)
    CALL get2d_data_nc(ncid, 'air_temp', maxi, maxj, temp_read, ifail)
    CALL get2d_data_nc(ncid, 'humidity', maxi, maxj, shum_read, ifail)
    CALL close_file_nc(TRIM(filenetin), ncid)

    day_rest = iday
    ioffset_rest = MOD(ioffset_rest, NINT(yearlen))
    day_rest = day_rest + timestep / (24 * 60 * 60.0)
    ! This bit so that we don't get too far out in our count....
    ! Anchor to a day if we start drifting.
    ! Means timestep can never be less than 1/1000 of a day!!!!
    IF (ABS(iday - day_rest) <= 1.0E-3) day_rest = iday
    IF (day_rest >= 31) THEN
       day_rest = day_rest - 30
       imonth_rest = imonth_rest + 1
       IF (imonth_rest == 13) THEN
          imonth_rest = 1
          iyear_rest = iyear_rest + 1
       END IF
    END IF
    iday = NINT(day_rest)
    PRINT *, 'day in embm restart is now', day_rest
    PRINT *, 'iday in embm restart is now', iday

    tq(1,1:maxi,1:maxj) = temp_read(:,:)
    tq(2,1:maxi,1:maxj) = shum_read(:,:)

    WRITE (*,220) 'Avg T', 'Avg Q'
    tmp_val = 0
    area = 0.0
    DO j = 1, maxj
       DO i = 1, maxi
          area = area + ds(j)
          tmp_val = tmp_val + tq(:,i,j) * ds(j)
       END DO
    END DO
    WRITE (*,210) tmp_val(1) / area, (tmp_val(2) / area) * 1000.0
210 FORMAT(2f13.9)
220 FORMAT(2a13)

  END SUBROUTINE inm_netcdf_embm

  ! This module writes netcdf restarts for embm
  SUBROUTINE outm_netcdf_embm(istep)
    USE netcdf
    USE genie_global, ONLY: writing_gui_restarts
    IMPLICIT NONE
    INTEGER, INTENT(IN) :: istep

    REAL, DIMENSION(maxi,maxj) :: temp_write, shum_write
    REAL :: lons1(maxi), lats1(maxj)
    INTEGER :: i, j, ntempid, nshumid
    INTEGER :: nlon1id, nlongit1id, nlat1id, nlatit1id, nrecsid(1), ioffid
    INTEGER :: dim1pass(2)
    CHARACTER(LEN=200) :: fname

    ! For date and restarts...
    INTEGER :: iday, iyearid, imonthid, idayid, istepatmid
    CHARACTER(LEN=10) :: yearstring
    CHARACTER(LEN=2) :: monthstring, daystring
    CHARACTER(LEN=7) :: datestring

    ! For netcdf...
    INTEGER :: ncid
    REAL :: timestep

    REAL :: tmp_val(2), area

    ! output file format is yyyy_mm_dd
    timestep = yearlen / REAL(nyear * ndta)
    iday = NINT(day_rest)
    IF (MOD(istep, iwstp * ndta) == 0 .OR. writing_gui_restarts) THEN
       ! WRITE A RESTART.....
       lons1 = nclon1
       lats1 = nclat1
       temp_write(:,:) = REAL(tq(1,1:maxi,1:maxj))
       shum_write(:,:) = REAL(tq(2,1:maxi,1:maxj))

       WRITE (datestring,'(i7.7)') istep
       WRITE (yearstring,'(i10)') iyear_rest
       WRITE (monthstring,'(i2.2)') imonth_rest
       WRITE (daystring,'(i2.2)') iday

       IF (writing_gui_restarts) THEN
          fname = 'gui_restart_embm.nc'
       ELSE
          fname = TRIM(dirnetout) // '/embm_restart_' // &
               & TRIM(ADJUSTL(yearstring)) // '_' // &
               & monthstring // '_' // daystring // '.nc'
       END IF
       IF (debug_end) PRINT *, ' Opening netcdf restart file for write: ', TRIM(fname)
       CALL check_err(NF90_CREATE(TRIM(fname), NF90_CLOBBER, ncid))
       CALL check_err(NF90_DEF_DIM(ncid, 'nrecs', 1, nrecsid(1)))
       CALL check_err(NF90_DEF_DIM(ncid, 'longitude', maxi, nlon1id))
       CALL check_err(NF90_DEF_DIM(ncid, 'latitude', maxj, nlat1id))

       CALL check_err(NF90_DEF_VAR(ncid, 'longitude', &
            & NF90_REAL, (/ nlon1id /), nlongit1id))
       CALL check_err(NF90_DEF_VAR(ncid, 'latitude', &
            & NF90_REAL, (/ nlat1id /), nlatit1id))
       dim1pass = (/ nlon1id, nlat1id /)
       CALL check_err(NF90_DEF_VAR(ncid, 'ioffset', NF90_INT, nrecsid, ioffid))
       CALL check_err(NF90_DEF_VAR(ncid, 'iyear', NF90_INT, nrecsid, iyearid))
       CALL check_err(NF90_DEF_VAR(ncid, 'imonth', NF90_INT, nrecsid, imonthid))
       CALL check_err(NF90_DEF_VAR(ncid, 'iday', NF90_INT, nrecsid, idayid))
       IF (writing_gui_restarts) &
            & CALL check_err(NF90_DEF_VAR(ncid, 'istep_atm', NF90_INT, &
            &                             nrecsid, istepatmid))

       CALL check_err(NF90_DEF_VAR(ncid, 'air_temp', &
            & NF90_DOUBLE, dim1pass, ntempid))
       CALL check_err(NF90_DEF_VAR(ncid, 'humidity', &
            & NF90_DOUBLE, dim1pass, nshumid))
       CALL check_err(NF90_ENDDEF(ncid))

       CALL check_err(NF90_PUT_VAR(ncid, iyearid, iyear_rest))
       CALL check_err(NF90_PUT_VAR(ncid, imonthid, imonth_rest))
       CALL check_err(NF90_PUT_VAR(ncid, idayid, iday))
       CALL check_err(NF90_PUT_VAR(ncid, ioffid, ioffset_rest))
       IF (writing_gui_restarts) &
            & CALL check_err(NF90_PUT_VAR(ncid, istepatmid, istep))

       CALL check_err(NF90_PUT_VAR(ncid, nlongit1id, lons1))
       CALL check_err(NF90_PUT_VAR(ncid, nlatit1id, lats1))
       CALL check_err(NF90_PUT_VAR(ncid, ntempid, temp_write))
       CALL check_err(NF90_PUT_VAR(ncid, nshumid, shum_write))

       CALL check_err(NF90_CLOSE(ncid))

       IF (debug_end) WRITE (*,220) 'Avg T','Avg Q'
       tmp_val = 0.0
       area = 0.0
       DO j = 1, maxj
          DO i = 1, maxi
             area = area + ds(j)
             tmp_val = tmp_val + tq(:,i,j) * ds(j)
          END DO
       END DO
       IF (debug_end) WRITE (*,210) tmp_val(1) / area, (tmp_val(2) / area) * 1000.0

210    FORMAT(2f13.9)
220    FORMAT(2a13)
    END IF

    day_rest = day_rest + timestep
    ! This bit so that we don't get too far out in our count....
    ! Anchor to a day if we start drifting.
    ! Means timestep can never be less than 1/1000 of a day!!!!
    IF (ABS(iday - day_rest) <= 1e-3) THEN
       PRINT *, 'CORRECTING TIME-LAG! in outm_netcdf in genie-goldstein', &
            & iday, day_rest
       day_rest = iday
    END IF
    IF (day_rest >= 31) THEN
       day_rest = day_rest - 30
       imonth_rest = imonth_rest + 1
       IF (imonth_rest == 13) THEN
          imonth_rest = 1
          iyear_rest = iyear_rest + 1
       END IF
    END IF
    iday = NINT(day_rest)
  END SUBROUTINE outm_netcdf_embm

  subroutine outm_surf(unit, co2_in, albedo_ocn, usurf_ocn, latent_ocn, &
       & sensible_ocn, netsolar_ocn, netlong_ocn, evap_ocn, pptn_ocn, &
       & runoff_ocn, latent_atm, sensible_atm, netsolar_atm, netlong_atm, &
       & evap_atm, pptn_atm,dhght_sic, darea_sic)
    IMPLICIT NONE
    INTEGER, INTENT(IN) :: unit
    REAL, DIMENSION(maxi,maxj), INTENT(IN) :: &
         & co2_in, usurf_ocn, albedo_ocn, latent_ocn, sensible_ocn, &
         & netsolar_ocn, netlong_ocn, evap_ocn, pptn_ocn, runoff_ocn, &
         & latent_atm, sensible_atm, netsolar_atm, netlong_atm, &
         & evap_atm, pptn_atm, dhght_sic, darea_sic

    WRITE (unit,10) latent_ocn
    WRITE (unit,10) sensible_ocn
    WRITE (unit,10) netsolar_ocn
    WRITE (unit,10) netlong_ocn
    WRITE (unit,10) evap_ocn
    WRITE (unit,10) pptn_ocn
    WRITE (unit,10) runoff_ocn
    WRITE (unit,10) albedo_ocn
    WRITE (unit,10) usurf_ocn
    WRITE (unit,10) latent_atm
    WRITE (unit,10) sensible_atm
    WRITE (unit,10) netsolar_atm
    WRITE (unit,10) netlong_atm
    WRITE (unit,10) evap_atm
    WRITE (unit,10) pptn_atm
    WRITE (unit,10) co2_in
    WRITE (unit,10) dhght_sic
    WRITE (unit,10) darea_sic

10  FORMAT(e21.13)
  END SUBROUTINE outm_surf

END MODULE embm_data

max_stars_repo_name: derpycode/cgenie

max_issues_repo_name: derpycode/cgenie

max_forks_repo_name: derpycode/cgenie


text:
!C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!C        Chargement du système à partir d'un fichier
!C              (SYSTEM2.BIN) créé par matlab
!C------------------------------------------------------------
!C  return : ierr=0 si chargement ok.                 12.12.96
!C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

subroutine loadmat(A, IA, JA, b, xref, ierr)

    implicit none

    integer NN, NNELEM, IA(*), JA(*), ierr
    double precision  A(*), b(*), xref(*)
    integer :: ios, j

    open(UNIT = 1, FILE = 'system2.bin',                    &
         ACTION='read', STATUS = 'old', IOSTAT=ios)
    write(*,*) 'ios=', ios
    if (ios .ne. 0) then
        write(*,*) "file 'system2.bin' cannot be opened!"
        stop 1
    end if

    write(*,*)'Chargement du systeme : SYSTEM2.BIN'

    ! Dimensions du systטme :

    read(1,*) NN
    read(1,*) NNELEM
    write(*,*)'Dimension :',NN
    write(*,*)'NNZ       :',NNELEM

    !      if((NN.NE.NELEM)) then
    !         write(*,*)'Out of memory !'
    !         write(*,*)'Dim. max =',N
    !         write(*,*)'NNZ  max =',NELEM
    !         ierr=1
    !      else

    ! Chargement de la matrice :

    do j = 1,NNELEM
        read(1,*) A(j)
        read(1,*) JA(j)
    end do
    do j = 1,NN+1
        read(1,*) IA(j)
    end do

    ! Chargement du vecteur b :

    do j = 1,NN
        read(1,*) b(j)
    end do

    ! Chargement du vecteur xref :

    do j = 1,NN
        read(1,*) xref(j)
    end do

    ! Fermeture du fichier

    close (UNIT = 1)
    write(*,*)'Ok..'
    ierr=0

end subroutine loadmat

max_stars_repo_name: rboman/progs

max_issues_repo_name: rboman/progs

max_forks_repo_name: rboman/progs


text:
Chintan Desai is a candidate for the ASUCD Senate with LEAD L.E.A.D. running in Winter 2009 ASUCD election.

 Candidate Statement and Activities

Hello All!

My name is Chintan Desai and I am running for a seat on your ASUCD Senate this winter. As an Aggie, I have had the privilege of serving UC Davis in several different capacities since entering the school in the fall of 2006. I am running this quarter precisely because I am enthusiastic about this campus and am eager to help you all make the changes you want to improve your quality of life here at UC Davis!

If you see me oncampus, please, please come by and say Hi! I am eager to listen to the issues and ideas you believe will help make our campus a better place! LEAD is all about empowering students to create positive change on our campus and in our community, so dont be afraid to make your voice heard!

Sincerely,
Chintan Desai

Experience

  Resident Adviser  Ryerson 2 (The Deuce)  September 2007June 2008
  The California Aggie California Aggie  City News Reporter  January 2007June 2008
  Legislative Intern  Assembly District 33 (San Luis Obispo)  Office of Assemblyman Sam Blakeslee  June 2008September 2008
  University of London Study Abroad Participant  Immersion Program, Queen Mary University of London  September 2008December 2008

Senate Goals

Firstyear Outreach
  Establish the leadership councils currently in place in the housing areas as ASUCD advisory committees
  Work to provide senator office hours in DCs and other Student Housing areas
  Provide residents with ASUCD developed programs that RAs can utilize to educate freshmen on how to become involved
  Make it a requirement for ASUCD senators to regularly visit the dorms and give presentations

Campus Club Events
  Work with Campus Events to facilitate the process of putting on club events
  Make the eventplanning process more timeefficient, which would benefit smaller clubs with less resources
  Provide incentives for group leaders to become aware of the process of putting on events
  Establish an easy to follow online tutorial on how to put on events which would benefit any group

Athletic Events
  Make UC Davis football tickets available at Freeborn Hall, increasing the attendance for our Aggie sports teams!

Campus Entertainment
  With the ASUCD Coffee House closing for renovations next year, it is important to find new areas on campus for students to hangout between classes
  Establish the Gunrock Pub, located in the Silo Union building and currently open only during lunch hours, as such a place
  Work to make the Gunrock Pub a more studentfriendly location by expanding its hours and establishing programs where student clubs could use the facility for meetings

Parking
  Work with the TAPS advisory committee to move the free parking time from 10 PM to 9 PM on weekdays
  Continue to promote beneficial TAPS programs, such as The Comet, a payasyougo parking meter

Renters Advocacy
  Partner with Steven D. Lee, former City/County Affairs Director on ongoing projects related to renters advocacy:
  Expanding a renters education program, so that more students are knowledgeable about the complexities of renting an apartment in Davis
  Institutionalizing and regularly reviewing the Davis Model Lease, thus making it more applicable to current times which would benefit landlords and tenants
  Lobbying apartment complexes in Davis to renew leases later in the year (preferably after spring break), giving renters more flexibility in terms of planning for the next year

20090124 16:08:21 nbsp I love his Ostrowskian plan for a mock ASUCD Senate. Users/GregWebb

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!
!!!!   MIT License
!!!!
!!!!   ParaMonte: plain powerful parallel Monte Carlo library.
!!!!
!!!!   Copyright (C) 2012-present, The Computational Data Science Lab
!!!!
!!!!   This file is part of the ParaMonte library.
!!!!
!!!!   Permission is hereby granted, free of charge, to any person obtaining a
!!!!   copy of this software and associated documentation files (the "Software"),
!!!!   to deal in the Software without restriction, including without limitation
!!!!   the rights to use, copy, modify, merge, publish, distribute, sublicense,
!!!!   and/or sell copies of the Software, and to permit persons to whom the
!!!!   Software is furnished to do so, subject to the following conditions:
!!!!
!!!!   The above copyright notice and this permission notice shall be
!!!!   included in all copies or substantial portions of the Software.
!!!!
!!!!   THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
!!!!   EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
!!!!   MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
!!!!   IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,
!!!!   DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR
!!!!   OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE
!!!!   OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
!!!!
!!!!   ACKNOWLEDGMENT
!!!!
!!!!   ParaMonte is an honor-ware and its currency is acknowledgment and citations.
!!!!   As per the ParaMonte library license agreement terms, if you use any parts of
!!!!   this library for any purposes, kindly acknowledge the use of ParaMonte in your
!!!!   work (education/research/industry/development/...) by citing the ParaMonte
!!!!   library as described on this page:
!!!!
!!!!       https://github.com/cdslaborg/paramonte/blob/master/ACKNOWLEDGMENT.md
!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

module SpecMCMC_ProposalStartCovMat_mod

    use Constants_mod, only: RK
    implicit none

    character(*), parameter         :: MODULE_NAME = "@SpecMCMC_ProposalStartCovMat_mod"

    real(RK), allocatable           :: proposalStartCovMat(:,:) ! namelist input

    type                            :: ProposalStartCovMat_type
        logical                     :: isPresent
        real(RK), allocatable       :: Def(:,:)
        real(RK), allocatable       :: Val(:,:)
        real(RK)                    :: null
        character(:), allocatable   :: desc
    contains
        procedure, pass             :: set => setProposalStartCorMat, checkForSanity, nullifyNameListVar
    end type ProposalStartCovMat_type

    interface ProposalStartCovMat_type
        module procedure            :: constructProposalStartCovMat
    end interface ProposalStartCovMat_type

    private :: constructProposalStartCovMat, setProposalStartCorMat, checkForSanity, nullifyNameListVar

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

contains

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    function constructProposalStartCovMat(nd,methodName) result(self)
#if defined DLL_ENABLED && !defined CFI_ENABLED
        !DEC$ ATTRIBUTES DLLEXPORT :: constructProposalStartCovMat
#endif
        use Constants_mod, only: IK, NULL_RK
        use String_mod, only: num2str
        implicit none
        integer(IK), intent(in)         :: nd
        character(*), intent(in)        :: methodName
        type(ProposalStartCovMat_type)  :: self
        integer(IK)                     :: i
        self%isPresent = .false.
        allocate( self%Def(nd,nd) )
        self%Def = 0._RK
        do i = 1,nd
            self%Def(i,i) = 1._RK
        end do
        self%null   = NULL_RK
        self%desc   = &
        "proposalStartCovMat is a real-valued positive-definite matrix of size (ndim,ndim), where ndim is the dimension of the &
        &sampling space. It serves as the best-guess starting covariance matrix of the proposal distribution. &
        &To bring the sampling efficiency of " // methodName // " to within the desired requested range, the covariance matrix will &
        &be adaptively updated throughout the simulation, according to the user's requested schedule. If proposalStartCovMat &
        &is not provided by the user or it is completely missing from the input file, its value will be automatically computed &
        &via the input variables proposalStartCorMat and proposalStartStdVec (or via their default values, if not provided). &
        &The default value of proposalStartCovMat is an ndim-by-ndim Identity matrix."
    end function constructProposalStartCovMat

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    subroutine nullifyNameListVar(self,nd)
#if defined DLL_ENABLED && !defined CFI_ENABLED
        !DEC$ ATTRIBUTES DLLEXPORT :: nullifyNameListVar
#endif
        use Constants_mod, only: IK
        implicit none
        class(ProposalStartCovMat_type), intent(in) :: self
        integer(IK), intent(in)                     :: nd
        if (allocated(proposalStartCovMat)) deallocate(proposalStartCovMat)
        allocate(proposalStartCovMat(nd,nd))
        proposalStartCovMat = self%null
    end subroutine nullifyNameListVar

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    subroutine setProposalStartCorMat(self, proposalStartStdVec, proposalStartCorMat, proposalStartCovMat)
#if defined DLL_ENABLED && !defined CFI_ENABLED
        !DEC$ ATTRIBUTES DLLEXPORT :: setProposalStartCorMat
#endif
        use Statistics_mod, only: getCovMatFromCorMat
        use Constants_mod, only: RK, IK
        implicit none
        class(ProposalStartCovMat_type), intent(inout)  :: self
        real(RK), intent(in)                            :: proposalStartCorMat(:,:), proposalStartStdVec(:)
        real(RK), intent(in), optional                  :: proposalStartCovMat(:,:)
        integer(IK)                                     :: i, j, nd

        if (present(proposalStartCovMat)) then
            self%val = proposalStartCovMat
        else
            self%val = self%null
        end if

        self%isPresent = .false.
        nd = size(proposalStartCorMat(:,1))
        do i = 1, nd
            do j = 1, nd
                if (self%val(j,i)==self%null) then
                    self%val(j,i) = self%Def(j,i)
                else
                    self%isPresent = .true.
                end if
            end do
        end do

        if (self%isPresent) return

        self%val = getCovMatFromCorMat  ( nd = nd &
                                        , StdVec = proposalStartStdVec &
                                        , CorMat = proposalStartCorMat &
                                        )

    end subroutine setProposalStartCorMat

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    subroutine checkForSanity(self,Err,methodName,nd)
#if defined DLL_ENABLED && !defined CFI_ENABLED
        !DEC$ ATTRIBUTES DLLEXPORT :: checkForSanity
#endif
        use Constants_mod, only: IK, RK
        use Matrix_mod, only: isPosDef
        use String_mod, only: num2str
        use Err_mod, only: Err_type
        implicit none
        class(ProposalStartCovMat_type), intent(in) :: self
        integer(IK), intent(in)                     :: nd
        character(*), intent(in)                    :: methodName
        type(Err_type), intent(inout)               :: Err
        character(*), parameter                     :: PROCEDURE_NAME = "@checkForSanity()"
        if (.not.isPosDef(nd,self%val)) then
            Err%occurred = .true.
            Err%msg =   Err%msg // &
                        MODULE_NAME // PROCEDURE_NAME // ": Error occurred. &
                        &The input requested proposalStartCovMat for the proposal of " // methodName // &
                        " is not a positive-definite matrix.\n\n"
        end if
    end subroutine checkForSanity

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

end module SpecMCMC_ProposalStartCovMat_mod