text:
!
!
!	ToDo:
!
!	Introduce FFT
!	this decreases the DFT(discrete fourrier transform)
!	from O(n**2) to O(n log(n))
!
!
module wann_interp
	!	module for Wannier interpolation
	!
	!	the interpolation scheme from
	!			PRB 74, 195118 (2006)
	!	was used
	use constants,		only:		dp, fp_acc, i_dp, pi_dp
	use matrix_math,	only:		zheevd_wrapper, 		&
									zheevx_wrapper,			&
									matrix_comm,			&
									blas_matmul,			&
									is_equal_mat,			&
									is_herm_mat,			&
									is_skew_herm_mat
	use input_paras,	only:		debug_mode,	wf_centers,	&
									use_cart_velo,			&
									kubo_tol, a_latt
	use file_io,		only:		write_eig_binary,		&
									write_ham_binary,		&
									write_velo
	use mpi_community,	only:		mpi_root_id, mpi_id, mpi_nProcs, ierr
	use mpi
	use k_space,		only:		get_cart_kpt
	use omp_lib

	implicit none

	private
	public					::		get_wann_interp,		&
									get_kubo_curv

	contains

!public:
	subroutine get_wann_interp(		do_gauge_trafo, 									&
									H_real, r_real, 									&
									R_frac,	wf_centers,									&
									kpt_idx, kpt, 									&
									e_k, V_ka, A_ka, Om_kab								&
							)
		!
		!	interpolates the k-space:
		!			-	H_k	: 		hamiltonian
		!			-	H_ka:		k-space derivative of Ham
		!			-	A_ka:		Berry conncection (only if r_real given)
		!			-	Om_kab:		Berry curvature
		!the real space basis (H_real (real space Hamiltonian) and optionally r_real (real space postition operator)	)
		!
		!	see 	>>>	PRB 74, 195118 (2006)	<<<			for more details on Wannier interpolation
		!
		logical,						intent(in)				::	do_gauge_trafo
		complex(dp),					intent(in)				::	H_real(:,:,:)
		complex(dp),	allocatable, 	intent(inout)			::	r_real(:,:,:,:)
		integer,						intent(in)				::	kpt_idx
		real(dp),						intent(in)				::	R_frac(:,:), kpt(3)
		real(dp),		allocatable,	intent(in)				::	wf_centers(:,:)
		real(dp),						intent(out)				::	e_k(:)
		complex(dp),	allocatable,	intent(inout)			::	V_ka(:,:,:)
		complex(dp),	allocatable,	intent(inout)			::	A_ka(:,:,:), Om_kab(:,:,:,:)

		complex(dp),	allocatable								::	U_k(:,:), H_ka(:,:,:)
		!
		!
								allocate(	U_k(			size(H_real,1),		size(H_real,2)		)		)
		if(	allocated(V_ka)	)	allocate(	H_ka(	3	,	size(H_real,1),		size(H_real,2)	 	)		)
		!
		!
		!
		!ft onto k-space (W)-gauge
		call FT_R_to_k(H_real, r_real,  R_frac, wf_centers, kpt_idx, kpt, U_k,  H_ka, A_ka, Om_kab)
		!
		!
		!get energies (H)-gauge
		call zheevd_wrapper(U_k, e_k)
		!
		!
		!rotate back to (H)-gauge
		if (allocated( V_ka)) then
			if( do_gauge_trafo	)	call W_to_H_gaugeTRAFO(kpt_idx, e_k, U_k, H_ka, A_ka, Om_kab)
			call get_velo(e_k, H_ka, A_ka, 	V_ka)
			!	DEBUG
			if(debug_mode)	call check_H_gauge_herm(kpt_idx, kpt, A_ka, Om_kab, V_ka)
		end if

		!
		return
	end subroutine

	subroutine	get_kubo_curv(	en_k, V_ka, kubo_curv	)
		!
		!	use Kubo formula to calculate curvature of bands
		!
		!		curv_n = -2 sum_np	IM[		V^_   V^_	]	/	(E_np - E_n)**2
		!
		real(dp),						intent(in)			::	en_k(:)
		complex(dp),					intent(in)			::	V_ka(:,:,:)
		complex(dp),	allocatable,	intent(out)			::	kubo_curv(:,:,:)
		real(dp)											::	curv_nn(3,3),	dEE_nnp
		integer												::	n_wf, n, np,  j
		!
		n_wf		=	size(en_k,1)
		allocate(kubo_curv(3,3,n_wf))
		kubo_curv	=	cmplx(0.0_dp,0.0_dp,dp)
		!
		do n = 1, n_wf
			!
			curv_nn	=	0.0_dp
			do np = 1, n_wf
				if(np==n)					cycle
				!
				dEE_nnp				=	(	en_k(np) - en_k(n)	)**2
				if(	abs(dEE_nnp)< kubo_tol)	cycle
				!
				do j = 1, 3
					kubo_curv(:,j,n)	=	kubo_curv(:,j,n)	+	aimag(	V_ka(:,n,np) * V_ka(j,np,n)	)	/ dEE_nnp
				end do
			end do
		end do
		!
		kubo_curv	=	-2.0_dp * kubo_curv
		!
		return
	end subroutine

!private:
	subroutine FT_R_to_k(H_real, r_real, R_frac, atom_frac, kpt_idx, kpt, H_k,	H_ka, A_ka, Om_kab)
		!	interpolates real space Ham and position matrix to k-space,
		!	according to
		!		PRB 74, 195118 (2006)		EQ.(37)-(40)
		!
		!
		!	->	only the H_real, and H_k have to be allocated
		!	->	all other quantities are only calculated if allocated
		!
		complex(dp),					intent(in)				::	H_real(:,:,:)
		complex(dp),	allocatable, 	intent(inout)			::	r_real(:,:,:,:)
		real(dp),						intent(in)				::	R_frac(:,:), kpt(3)
		real(dp),		allocatable,	intent(in)				::	atom_frac(:,:)
		integer,						intent(in)				::	kpt_idx
		complex(dp),					intent(out)				::	H_k(:,:)
		complex(dp),	allocatable,	intent(inout)			::	H_ka(:,:,:), A_ka(:,:,:), Om_kab(:,:,:,:)
		real(dp),		allocatable								::	delta_at(:,:,:)
		real(dp)												::	dR(3), dR_cart(3),kpt_cart(3),ft_angle, two_pi
		complex(dp)												::	ft_phase
		logical													::	use_pos_op, do_en_grad
		integer    												::	sc, a, b, n_sc, n_wf, n, m
		!
		n_sc	=	size(R_frac,2)
		n_wf	=	size(H_real,1)
		two_pi	=	2.0_dp 	* pi_dp

		!jobs
		do_en_grad		= allocated(H_ka)
		use_pos_op		= allocated(A_ka) .and. allocated(r_real) .and. allocated(Om_kab)
		!
		!init
						H_k		= 	0.0_dp
		if(do_en_grad)	H_ka	= 	0.0_dp
		if(use_pos_op)	A_ka	= 	0.0_dp
		if(use_pos_op)	Om_kab	= 	0.0_dp
		!
		kpt_cart	=	get_cart_kpt(a_latt, kpt)

		!
		! I		INTERPOLATE HAM & VELO(s)
		!
		if( allocated(atom_frac)) then
			if(	kpt_idx <= mpi_nProcs	)then
				write(*,'(a,i7.7,a)')	"[#",mpi_id,";FT_R_to_k]: will use TIGHT-B FT CONVENTION"
			end if
			! ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!|
			!	CONSIDER ATOMIC POSITIONS IN FOURIER TRANSFORM															!|
			!	THIS IS THE 	"TIGHT BINDING" 	CONVENTION															!|
			!	---------																								!|
			allocate(	delta_at(size(atom_frac,1),size(atom_frac,2),size(atom_frac,2)))
			do n = 1, n_wf
				do m = 1, n_wf
					delta_at(:,n,m)	=	atom_frac(:,m) - atom_frac(:,n)
				end do
			end do
			!$OMP PARALLEL DEFAULT(none)								&
			!$OMP PRIVATE( 	dR, dR_cart,  ft_angle, ft_phase, m,n)		&
			!$OMP SHARED(  	use_cart_velo,kpt_cart, n_sc, n_wf,kpt_idx, R_frac, delta_at, a_latt, H_real, kpt, H_k, H_ka, do_en_grad, two_pi)
			!$OMP DO REDUCTION(+: H_k, H_ka)
			do sc = 1, n_sc																								!|
				do m = 1 , n_wf																							!|
					do n = 1, n_wf																						!|
						dR(:)		=	R_frac(:,sc) +	delta_at(:,n,m)													!|
						!	get ft angle																				!|
						if( use_cart_velo	) then																		!|
							dR_cart		=	intern_to_cart(a_latt, dR)													!|
							ft_angle	=	 			dot_product(	kpt_cart(:),	dR_cart(:)	)! cartesian coords	!|
						else																							!|
							ft_angle	=	two_pi * 	dot_product(		kpt(:) ,		dR(:)	)! internal coords	!|
						end if																							!|
						!																								!|
						!																								!|
						ft_phase	= 	cmplx(	cos(ft_angle), sin(ft_angle),	dp	)	!unit indepedent				!|
						!Hamilton operator																				!|
						H_k(n,m)	=			H_k(n,m)		+	ft_phase						* H_real(n,m,sc)	!|
						!																								!|
						!OPTIONAL energy gradients																		!|
						if( do_en_grad)		then																		!|
							!																							!|
							!	get dR_cart																				!|
							if( .not. use_cart_velo ) then																!|																					!|
								dR_cart(:)	=	two_pi  * dR(:)															!|
							end if																						!|
							!																							!|
							H_ka(:,n,m) 	=	H_ka(:,n,m)		+ ft_phase * cmplx(0.0_dp,dR_cart(:),dp)* H_real(n,m,sc)!|
						end if																							!|
					end do																								!|
				end do
			end do
			!$OMP END DO
			!$OMP END PARALLEL
			! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!|
			!
			!
		else
			if(	kpt_idx <= mpi_nProcs	)then
				write(*,'(a,i7.7,a)')	"[",mpi_id,"FT_R_to_k]: will use WANNIER FT CONVENTION"
			end if
			!
			! ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!|
			!	NEGLECT ATOMIC POSITIONS																				!|
			!	THIS IS THE 	"WANNIER"	CONVENTION																	!|
			!	---------
			!$OMP PARALLEL DEFAULT(none)						&
			!$OMP PRIVATE( dR_cart, ft_angle, ft_phase, a,m,n)	&
			!$OMP SHARED(use_cart_velo, n_sc, n_wf,  R_frac, a_latt, H_real, kpt, H_k, H_ka, do_en_grad, two_pi)
			!$OMP DO REDUCTION(+: H_k, H_ka)
			do sc = 1, n_sc																								!|
				ft_angle	=	two_pi * dot_product(	kpt(:),R_frac(:,sc)	)		! internal  coord					!|
				ft_phase	= 	cmplx(	cos(ft_angle), sin(ft_angle), dp)												!|
				!Hamilton operator																						!|
				H_k(:,:)	= 			H_k(:,:)		+	ft_phase								* H_real(:,:,sc)	!|
				!																										!|
				!OPTIONAL energy gradients																				!|
				if( do_en_grad)		then																				!|
					!	get dR_cart																						!|
					if( use_cart_velo ) then																			!|
						dR_cart(:)	= 	intern_to_cart( a_latt(:,:), R_frac(:,sc)	)									!|
					else																								!|
						dR_cart(:)	=	two_pi * R_frac(:,sc)															!|
					end if																								!|
					!																									!|
					do a = 1, 3																							!|
						H_ka(a,:,:) 	=		H_ka(a,:,:)		+	ft_phase * i_dp * 	dR_cart(a) * H_real(:,:,sc)		!|
					end do																								!|
				end if																									!|
			end do
			!$OMP END DO
			!$OMP END PARALLEL
			! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!|
		end if
		!
		!
		!
		!
		! II	INTERPOLATE CONN & CURV
		!
		if(use_pos_op) then
			if(	kpt_idx <= mpi_nProcs	)then
				write(*,'(a,i7.7,a)')	"[",mpi_id,"FT_R_to_k]: will interpolate conn & curv additonally"
			end if
			if( allocated(atom_frac)	) then
				! ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!|
				!	CONSIDER ATOMIC POSITIONS IN FOURIER TRANSFORM															!|
				!	THIS IS THE 	"TIGHT BINDING" 	CONVENTION															!|
				!	---------
				!$OMP PARALLEL DEFAULT(none) &
				!$OMP PRIVATE(m,n, dR, dR_cart, ft_angle, ft_phase, a, b) &
				!$OMP Shared(use_cart_velo, n_sc, n_wf, A_ka, Om_kab, a_latt, R_frac, delta_at, kpt,  r_real, two_pi)
				!$OMP DO REDUCTION(+: A_ka, Om_kab)
				do sc = 1, n_sc																								!|
					do m = 1 , n_wf																							!|
						do n = 1, n_wf																						!|
							dR(:)		=	R_frac(:,sc) +	delta_at(:,n,m)				! internal units					!|
							!																								!|
							ft_angle	=	2.0_dp * pi_dp * dot_product(kpt(:),dR(:))										!|
							ft_phase	= 	cmplx(	cos(ft_angle), sin(ft_angle)	,	dp	)								!|
							!																								!|
							!	get dR_cart																					!|
							if( use_cart_velo ) then																		!|
								dR_cart(:)	= 	intern_to_cart( a_latt(:,:), dR(:)	)										!|
							else																							!|
								dR_cart(:)	=	two_pi * dR(:)																!|
							end if																							!|
							!																								!|
							!																								!|
							do a = 1, 3																						!|
								!connection																					!|
								A_ka(a,n,m)			=	A_ka(a,n,m)		+	ft_phase				* r_real(a,n,m,sc)		!|
								!curvature																					!|
								do b = 1, 3																					!|
									Om_kab(a,b,n,m)	=	Om_kab(a,b,n,m)	&													!|
													 		+ 	ft_phase * i_dp * dR_cart(a) 	* r_real(b,n,m,sc)			!|
									Om_kab(a,b,n,m)	=	Om_kab(a,b,n,m)	& 													!|
															- 	ft_phase * i_dp * dR_cart(b) 	* r_real(a,n,m,sc)			!|
								end do																						!|
							end do																							!|
						end do																								!|
					end do																									!|
					!																										!|
					!																							!~~~~~~~~~~~!|
				end do
				!$OMP END DO
				!$OMP END PARALLEL
				! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!|
				!
				!
			else
				! ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!|
				!	NEGLECT ATOMIC POSITIONS																				!|
				!	THIS IS THE 	"WANNIER"	CONVENTION																	!|
				!	---------
				!$OMP PARALLEL DEFAULT(none) &
				!$OMP PRIVATE(dR_cart, ft_angle, ft_phase, a, b) &
				!$OMP Shared(use_cart_velo, n_sc, a_latt, R_frac, r_real, A_ka, Om_kab, kpt, two_pi)
				!$OMP DO REDUCTION(+: A_ka, Om_kab)
				do sc = 1, n_sc																								!|
					ft_angle	=	two_pi * dot_product(kpt(:),R_frac(:,sc))												!|
					ft_phase	= 	cmplx(	cos(ft_angle), sin(ft_angle),	dp	)											!|
					!	get dR_cart																							!|
					if( use_cart_velo ) then																				!|
						dR_cart(:)	= 	intern_to_cart( a_latt(:,:), R_frac(:,sc)	)										!|
					else																									!|
						dR_cart(:)	=	two_pi * R_frac(:,sc)																!|
					end if																									!|
					!																										!|
					do a = 1, 3																								!|
						!connection																							!|
						A_ka(a,:,:)			=	A_ka(a,:,:)		+	ft_phase		* r_real(a,:,:,sc)						!|
						!curvature																							!|
						do b = 1, 3																							!|
							Om_kab(a,b,:,:)	=	Om_kab(a,b,:,:)	&															!|
											 		+ 	ft_phase * i_dp * dR_cart(a) * r_real(b,:,:,sc)						!|
							Om_kab(a,b,:,:)	=	Om_kab(a,b,:,:)	& 															!|
													- 	ft_phase * i_dp * dR_cart(b) * r_real(a,:,:,sc)						!|
						end do																								!|
					end do																									!|
					!																										!|
					!																							!~~~~~~~~~~~!|
				end do
				!$OMP END DO
				!$OMP END PARALLEL
				! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~!|
			end if
		end if
		!
		if(debug_mode) then
			call debug_ft_phase(R_frac, delta_at, kpt, a_latt)
			call write_ham_binary(kpt_idx,	H_k)
			call check_W_gauge_herm(kpt, H_k, H_ka, A_ka, Om_kab)
		end if
		!
		return
	end subroutine

	pure function intern_to_cart(  a_latt, vec	) result(vec_cart)
		real(dp),	intent(in)		::		a_latt(3,3), vec(3)
		real(dp)					::		vec_cart(3)
		!integer						::		dim
		!
		!vec_cart	=	0.0_dp
		!do dim = 1, 3
		!	vec_cart(:)	=	vec_cart(:)	+	vec(dim)	*	a_latt(dim,:)
		!end do
		vec_cart(:)		=	matmul(		a_latt, vec(:)	)
		!
		return
	end function

	subroutine get_velo(en_k, H_ka, A_ka, V_ka)
		!	calc the (H)-gauge velocity matrix
		!
		!	PRB 74, 195118 (2006) EQ.(31)
		real(dp),						intent(in)		::		en_k(:)
		complex(dp),					intent(in)		::		H_ka(:,:,:)
		complex(dp),	allocatable,	intent(inout) 	::		A_ka(:,:,:)
		complex(dp),					intent(out)		::		V_ka(:,:,:)
		complex(dp)										::		cmplx_eDiff
		integer											::		m, n, n_wf
		!
		n_wf		=	size(en_k,1)
		!
		!
		if( allocated(A_ka)	) then
			!	if CONNECTION given
			!	-> apply correction to H_Ka
			!
			!$OMP PARALLEL DO DEFAULT(none) &
			!$OMP PRIVATE(n,cmplx_eDiff) &
			!$OMP SHARED(n_wf, en_k, H_Ka, A_ka, V_Ka)
			do m = 1, n_wf
				do n = 1, n_wf
					if(	n >	m )	then
						cmplx_eDiff	=	cmplx(	0.0_dp	,	en_k(m) - en_k(n),	dp)
						!
						V_ka(:,n,m)	= H_ka(:,n,m)	- cmplx_eDiff	* 	A_ka(:,n,m)
						V_ka(:,m,n)	= H_ka(:,m,n)	+ cmplx_eDiff	*	A_ka(:,m,n)
					end if
				end do
			end do
			!$OMP END PARALLEL DO
		else
			V_ka		=	H_ka
		end if
		!
		!
		return
	end subroutine

!
!
!
!
!
!
!
!	**************************************************************************************************************************************************
!	**************************************************************************************************************************************************
!	**************************************************************************************************************************************************
!
!						~~~~			GAUGE TRAFO ROUTINES				~~~~~~~~~~~~~~
!
!	**************************************************************************************************************************************************
!	**************************************************************************************************************************************************
!	**************************************************************************************************************************************************
!
	subroutine W_to_H_gaugeTRAFO(kpt_idx, e_k, U_k, H_ka, A_ka, Om_kab)
		!	see PRB 74, 195118 (2006) EQ. 21 - 31
		integer,						intent(in)		::	kpt_idx
		real(dp),						intent(in)		::	e_k(:)
		complex(dp),					intent(in)		::	U_k(:,:)
		complex(dp),	allocatable,	intent(inout)	::	H_ka(:,:,:)
		complex(dp),	allocatable, 	intent(inout)	::	A_ka(:,:,:), Om_kab(:,:,:,:)
		complex(dp),	allocatable						::	D_ka(:,:,:)
		!
		!
		!debug (W)-gauge
		if(debug_mode)	then
			call check_velo(U_k, H_ka)
			call write_eig_binary(kpt_idx,	U_k)
		end if
		!
		!
		!do 	(W) -> (Hbar)
		call rotate_gauge(U_k, H_ka, 	A_ka, Om_kab )
		if(debug_mode)	call check_Hbar_gauge_herm(H_ka, A_ka, Om_kab)
		!
		!
		!conn/curv	 (Hbar) -> (H)
		if( allocated(A_ka) )	then
			!	need gauge covar deriv (remember A_ka is not gauge covariant)
			allocate(		D_ka(	3,	size(H_ka,2),	size(H_ka,3))				)
			!
			!
			call get_gauge_covar_deriv(e_k, H_ka, D_ka)
			!
			call conn_gaugeTrafo(D_ka, A_ka)
			call curv_gaugeTrafo(D_ka, A_ka, Om_kab)
		end if
		!
		!
	end subroutine

	subroutine rotate_gauge(U_k, H_ka, A_ka, Om_kab)
		!	PRB 74, 195118 (2006)	EQ.(21)
		complex(dp),					intent(in)		::	U_k(:,:)
		complex(dp), 					intent(inout)	::	H_ka(:,:,:)
		complex(dp), 	allocatable,	intent(inout)	::	A_ka(:,:,:), Om_kab(:,:,:,:)
		complex(dp),	allocatable						::	U_dag(:,:), tmp(:,:), M_in(:,:)
		integer											::	a, b
		!
		allocate(	U_dag(		size(U_k,1), size(U_k,2)	))
		allocate(	tmp(		size(U_k,1), size(U_k,2)	))
		allocate(	M_in(		size(U_k,1), size(U_k,2)	))

		U_dag	=	conjg(	transpose(	U_k		))
		!
		!
		do a = 1, 3
			!
			!	VELOCITIES
										M_in			=	H_ka(a,:,:)
										!	BLAS
										!call	blas_matmul(	M_in, 	U_k,		tmp			)
										!call	blas_matmul(	U_dag, 	tmp,		M_in		)
										!	matmul
										tmp				=	blas_matmul(M_in, U_k)
										M_in			=	blas_matmul(U_dag,tmp)
										!
										H_ka(a,:,:)		=	M_in(:,:)
			!	CONNECTION
			if( allocated(A_ka))then
										M_in			=	A_ka(a,:,:)
										tmp(:,:)		=	blas_matmul(	U_dag		,	tmp		)
										A_ka(a,:,:)		=	blas_matmul(	tmp			, 	U_k		)
										!call	blas_matmul(	U_dag,	tmp		,	tmp		)
										!call	blas_matmul(	tmp,	U_k		,	A_ka(a,:,:))
			end if
			!	CURVATURE
			if( allocated(Om_kab)	)then
				do b = 1,3
										tmp				=	Om_kab(a,b,:,:)
										tmp				=	blas_matmul(	U_dag		,	tmp		)
										Om_kab(a,b,:,:)	=	blas_matmul(	tmp			,	U_k		)
				end do
			end if
			!
		end do
		!
		!
		return
	end subroutine

	function unit_rot(	M,	U, U_dag	)	result(M_rot)
		complex(dp),	intent(in)		::	M(:,:), U(:,:), U_dag(:,:)
		complex(dp),	allocatable		::	M_rot(:,:)
		complex(dp), save,	allocatable		::	tmp1(:,:), tmp2(:,:)

		return
	end function

	subroutine get_gauge_covar_deriv(e_k, H_ka,	D_ka )
		!	PRB 74, 195118 (2006)	EQ.(24)
		real(dp),			intent(in)		::	e_k(:)
		complex(dp),		intent(in)		::	H_ka(:,:,:)
		complex(dp),		intent(out)		::	D_ka(:,:,:)
		integer								::	m, n, a
		real(dp)							::	eDiff_mn, max_err
		!
		D_ka(:,:,:)	=	cmplx(0.0_dp, 0.0_dp, dp)
		!
		!$OMP PARALLEL DO DEFAULT(none) &
		!$OMP PRIVATE(n,eDiff_mn) &
		!$OMP SHARED(kubo_tol, e_k, H_Ka, D_ka)
		do m = 1, size(D_ka,3)
			do n = 1, size(D_ka,2)
				if(	n >	m )	then
					!
					!
					eDiff_mn	=	e_k(m)	- e_k(n)
					!
					if(abs(eDiff_mn) > 	kubo_tol	)	then
						D_ka(1:3,n,m)	=	H_ka(1:3,n,m) / 	eDiff_mn
						D_ka(1:3,m,n)	=	H_ka(1:3,m,n) /	( - eDiff_mn	)
					else
						write(*,'(a)',advance="no")	'[;get_gauge_covar_deriv]: '
						write(*,'(a,i6,a,i6)')		'WARNING degenerate bands detetected n=',n,' m=',m
					end if
					!
					!
				end if
			end do
		end do
		!$OMP END PARALLEL DO
		!
		if(debug_mode)	then
			do a = 1, 3
				!
				!
				if( .not. is_skew_herm_mat(D_ka(a,:,:), max_err)	)	then
					write(*,'(a,i1,a,f16.7)')	"[get_gauge_covar_deriv/DEBUG-MODE]:	WARNING D_(k,a=",a,&
																						") is not skew hermitian, max_err=", max_err
					if(.not. is_herm_mat(H_ka(a,:,:),max_err)	) then
						write(*,'(a,i1,a,f16.7)')	"[get_gauge_covar_deriv/DEBUG-MODE]:	WARNING H_(k,a=",a,&
																						") is not hermitian, max_err=", max_err
					end if
				end if
				!
				!
			end do
		end if
		!
		return
	end subroutine

	pure subroutine conn_gaugeTrafo(D_ka, A_ka)
		!	PRB 74, 195118 (2006)	EQ.(25)
		!
		!	Lapack
		!		https://software.intel.com/en-us/mkl-developer-reference-fortran-gemm#90EAA001-D4C8-4211-9EA0-B62F5ADE9CF0
		!		C :-
		complex(dp),		intent(in)		::	D_ka(:,:,:)
		complex(dp),		intent(inout)	::	A_ka(:,:,:)
		!
		!
		A_ka	=	A_ka	+ i_dp	*	D_ka
		!
		return
	end subroutine

	subroutine curv_gaugeTrafo(D_ka, A_ka, Om_kab)
		!	PRB 74, 195118 (2006)	EQ.(27)
		complex(dp),		intent(in)		::	D_ka(:,:,:), A_ka(:,:,:)
		complex(dp),		intent(inout)	::	Om_kab(:,:,:,:)
		complex(dp),	allocatable			::	mat_comm(:,:)
		integer								::	a, b
		!
		allocate(	mat_comm(	size(Om_kab,3),size(Om_kab,4)	)		)
		!
		do b = 1, 3
			do a = 1, 3
				Om_kab(a,b,:,:)	=	Om_kab(a,b,:,:)		-			matrix_comm(	D_ka(a,:,:), 	A_ka(b,:,:)		)
				!
				!
				Om_kab(a,b,:,:)	=	Om_kab(a,b,:,:)		+			matrix_comm(	D_ka(b,:,:), 	A_ka(a,:,:)		)
				!
				Om_kab(a,b,:,:)	=	Om_kab(a,b,:,:)		-	i_dp *	matrix_comm( D_ka(a,:,:), 	D_ka(b,:,:))
			end do
		end do
		!
		!
		return
	end subroutine

!
!
!
!
!
!
!
!	**************************************************************************************************************************************************
!	**************************************************************************************************************************************************
!	**************************************************************************************************************************************************
!
!						~~~~			DEBUGGING HELPER ROUTINES				~~~~~~~~~~~~~~
!
!	**************************************************************************************************************************************************
!	**************************************************************************************************************************************************
!	**************************************************************************************************************************************************
!

	subroutine debug_ft_phase(R_frac, delta_at, kpt, a_latt)
		real(dp),						intent(in)				::	R_frac(:,:), delta_at(:,:,:), kpt(3), a_latt(3,3)
		real(dp)												:: 	dR(3), dR_cart(3), kpt_cart(3), &
																	ft_agl, ft_agl_cart, two_pi
		complex(dp)												::	R_sum(3), R_sum_cart(3),			&
																	ft_phs, ft_phs_cart,&
																	ft_sum, ft_sum_cart, two_pi_i
		integer													::	m, n, sc, n_sc , n_wf
		!
		n_wf	=	size(delta_at,2)
		n_sc	=	size(R_frac,2)
		!
		kpt_cart(:)	=	get_cart_kpt(a_latt, kpt)
		!
		call MPI_BARRIER(MPI_COMM_WORLD, ierr)
		if(	mpi_id == mpi_root_id ) then
			write(*,*)	"[debug_ft_phase]: now check (n,m)=(",n,",",m,")	#sc=",n_sc	,"kpt=",kpt
			two_pi		=	2.0_dp * pi_dp
			two_pi_i	= 	two_pi * i_dp
			ft_sum 		=	0.0_dp
			ft_sum_cart	=	0.0_dp
			R_sum		=	0.0_dp
			R_sum_cart	=	0.0_dp
			!
			do m = 1 , n_wf
					do n = 1, n_wf
			do sc = 1, n_sc

						dR(:)		=	R_frac(:,sc)	+ delta_at(:,n,m)				! 	int
						dR_cart(:)	=	intern_to_cart( a_latt, dR )	!	cart
						!
						!	INTERNAL UNITS
						ft_agl		=	two_pi * dot_product(	kpt(:),	dR(:)	)
						ft_phs		=	cmplx( cos(ft_agl), sin(ft_agl), dp)
						ft_sum 		=	ft_sum 		+	ft_phs
						R_sum(:)	=	R_sum(:)	+	two_pi_i	* dR(:)	* ft_phs
						!	-------------------
						!
						!	CARTESIAN COORDS
						ft_agl_cart =	dot_product( kpt_cart(:), dR_cart(:) )
						ft_phs_cart = 	cmplx(	cos(ft_agl_cart), sin(ft_agl_cart), 	dp)
						ft_sum_cart =	ft_sum_cart + 	ft_phs_cart
						R_sum_cart(:)	=	R_sum_cart(:)	+	dR_cart(:)	* ft_phs_cart
						!	-------------------
						!
						end do
						!
			!if( abs(ft_sum - 1.0_dp) > 1e-4_dp) &
			write(*,*)	"[debug_ft_phase]:	WARNING      ft_sum=",ft_sum,'vs ',ft_sum_cart,"=ft_sum_cart"

			write(*,*)	"[debug_ft_phase]:	WARNING      R_sum=",R_sum,'.'
						write(*,*)	"[debug_ft_phase]:	WARNING      R_sum_cart=",R_sum_cart,'.'

					end do
				end do

		end if
		call MPI_BARRIER(MPI_COMM_WORLD, ierr)
		!
		return
	end subroutine

	logical function curv_is_herm( Om_kab, max_err)
		complex(dp), allocatable,	intent(in)		::	Om_kab(:,:,:,:)
		real(dp),			intent(out)		::	max_err
		integer								::	a, b
		!
		!
		if(allocated(Om_kab))then
			if( size(Om_kab,3) == size(Om_kab,4) )	then
				curv_is_herm	=	.true.
				do a = 1, 3
					do b = 1, 3
						if(.not. is_herm_mat(Om_kab(b,a,:,:),max_err)) 	then
							curv_is_herm = .false.
							!write(*,*)	"[conn_curv_is_herm/DEBUG-MODE]: Om_k,a=",b,",b=",a," is not hermitian"
						end if
					end do
				end do
			else
				curv_is_herm	= .false.
				stop " [curv_is_herm/DEBUG-MODE]: k-space matrices life on different basis sets"
			end if
		end if
		!
		return
	end function

	logical function velo_is_herm(V_ka, max_err)
		complex(dp),		intent(in)		::	V_ka(:,:,:)
		real(dp),			intent(out)		::	max_err
		integer								::	a
		!
		velo_is_herm	=	 ( size(V_ka,2) == size(V_ka,3) 	) .and. (size(V_ka,1)==3)
		!
		!
		if(velo_is_herm)	then
			do a = 1, 3
				if(.not. is_herm_mat(V_ka(a,:,:),max_err)) 	then
					velo_is_herm = .false.
				end if
			end do
		else
			velo_is_herm	= .false.
			stop " [velo_is_herm/DEBUG-MODE]: ERROR velo operator matrix not symmetric"
		end if
		!
		return
	end function

	subroutine check_velo(U_k, VW_ka)
		!	check if rotating forward and backwards wants is identity operation
		!
		complex(dp),		intent(in)		::	U_k(:,:), VW_ka(:,:,:)
		complex(dp),	allocatable			::	VW_new_ka(:,:,:), A_ka(:,:,:), Om_kab(:,:,:,:)
		!
		allocate(		VW_new_ka(size(VW_ka,1), size(VW_ka,2),size(VW_ka,3))	)
		!
		VW_new_ka	=	VW_ka
		!
		call rotate_gauge(U_k, VW_new_ka, A_ka, Om_kab)
		!
		!	now rotate back to wannier gauge
		call rotate_gauge(	conjg(transpose(U_k)),	VW_new_ka,	A_ka, Om_kab)
		!
		!
		if(			is_equal_mat(	1e-9_dp	, VW_ka(1,:,:), VW_new_ka(1,:,:))	&
			.and.	is_equal_mat(	1e-9_dp	, VW_ka(2,:,:), VW_new_ka(2,:,:))	&
			.and.	is_equal_mat(	1e-9_dp	, VW_ka(3,:,:), VW_new_ka(3,:,:))	&
			)	then
			write(*,*)	"[wann_interp/check_velo]:	SUCCESS gauge consistency seems fine "
		else
			write(*,*)	"[wann_interp/check_velo]:	WARNING gauge consistency not given "
		end if

		return
	end subroutine

	subroutine check_W_gauge_herm(kpt, H_k, H_ka, A_ka, Om_kab)
		real(dp),						intent(in)		::		kpt(3)
		complex(dp),					intent(in)		::		H_k(:,:)
		complex(dp),	allocatable, 	intent(in)		::		H_ka(:,:,:), A_ka(:,:,:), Om_kab(:,:,:,:)
		real(dp)										::		max_err
		character(len=44)								::		k_string
		character(len=51)								::		warn_msg
		character(len=26)								::		max_string
		character(len=32)								::		allo_lst
		logical 										::		ham_herm, velo_herm, conn_herm, curv_herm
		!
		warn_msg	=								'[check_w_gauge_herm/DEBUG-MODE]:	WARNING (W)-gauge '
		max_string	=								' IS NOT hermitian(max_err='
		write(k_string,'(a,f6.2,a,f6.2,a,f6.2,a)')	') at rel. kpt=( ',kpt(1),', ',kpt(2),', ',kpt(3),')'

		!
		!	CHECK HAMILTONIAN
		ham_herm	=	 is_herm_mat( H_k, max_err)
		if(.not. ham_herm) 				write(*,'(a,e16.7,a)')	warn_msg//"H_k"//max_string,	max_err ,k_string
		allo_lst	=	" ham, "
		!
		!	CHECK DERIVATIVE OF HAM
		if(allocated(H_ka)) then
			allo_lst	=	trim(allo_lst)	//	"velo, "
			velo_herm 	=	 velo_is_herm( H_ka, max_err)
			if(.not. velo_herm) 		write(*,'(a,e16.7,a)')	warn_msg//"H_ka"//max_string,	max_err ,k_string
		else
			velo_herm	=	.true.
		end if
		!
		!	CHECK CONNECTION
		if(allocated(A_ka)) then
			allo_lst	=	trim(allo_lst)	//	"conn, "
			conn_herm	=	velo_is_herm( A_ka,	max_err)
			if(.not. conn_herm)			write(*,'(a,e16.7,a)')	warn_msg//"A_ka"//max_string,	max_err ,k_string
		else
			conn_herm = .true.
		end if
		!
		!	CHECK CURVATURE
		if(allocated(Om_kab))	then
			allo_lst	=	trim(allo_lst)	//	"curv"
			curv_herm	= curv_is_herm( Om_kab,	max_err)
			if (.not. curv_herm) 		write(*,'(a,e16.7,a)')	warn_msg//"Om_kab"//max_string,	max_err ,k_string
		else
			curv_herm = .true.
		end if
		!
		!
		if( (ham_herm .and. velo_herm) .and. (conn_herm .and. curv_herm) ) then
			write(*,*)	"[check_w_gauge_herm/DEBUG-MODE]: SUCCESS (W)-gauge quantities ("//trim(allo_lst)//") are hermitian"
		end if
		!
		return
	end subroutine

	subroutine check_Hbar_gauge_herm(H_ka, A_ka, Om_kab)
		complex(dp),				intent(in)			::	H_ka(:,:,:)
		complex(dp), allocatable,	intent(in)			::	A_ka(:,:,:), Om_kab(:,:,:,:)
		real(dp)										::	max_err
		!
		if( .not. velo_is_herm(H_ka,max_err))	then
			write(*,'(a,f16.7)')	'[check_Hbar_gauge_herm/DEBUG-MODE]: WARNING	(Hbar)-gauge H_ka IS NOT hermitian, max_err=',max_err
		end if
		!
		if(allocated(A_ka)) then
			if( .not. velo_is_herm(A_ka,max_err))	then
				write(*,'(a,f16.7)')	'[check_Hbar_gauge_herm/DEBUG-MODE]: WARNING	(Hbar)-gauge A_ka IS NOT hermitian, max_err=',max_err
			end if
		end if
		!
		if(allocated(Om_kab))	then
			if( .not. curv_is_herm(Om_kab,max_err) )	then
				write(*,'(a,f16.7)')	'[check_Hbar_gauge_herm/DEBUG-MODE]: WARNING	(Hbar)-gauge Om_kab IS NOT hermitian, max_err=',max_err
			end if
		end if
		!
		!
		return
	end subroutine

	subroutine check_H_gauge_herm(kpt_idx, kpt, A_ka, Om_kab, V_ka)
		integer,						intent(in)			::	kpt_idx
		real(dp),						intent(in)			::	kpt(3)
		complex(dp),	allocatable,	intent(in)			::	A_ka(:,:,:), Om_kab(:,:,:,:)
		complex(dp),					intent(in)			::	V_ka(:,:,:)
		real(dp)											::	max_err
		character(len=31)									::	k_string
		character(len=24)									::	allo_lst
		logical												::	conn, curv, velo, is_herm
		!
		allo_lst	=	" "
		write(k_string,'(a,f6.2,a,f6.2,a,f6.2,a)')	'( ',kpt(1),', ',kpt(2),', ',kpt(3),') '
		is_herm	= .true.
		!
		!	VELOCITY
		allo_lst	=	trim(allo_lst) // "velo"
		velo		=	velo_is_herm( V_ka, max_err )
		is_herm 	=	is_herm .and. velo
		if(.not. velo) 		write(*,'(a,f16.7)')	"[check_H_gauge_herm/DEBUG-MODE]:	"								//	&
														"WARNING (H)-gauge V_KA IS NOT hermitian at rel. kpt= "			//	&
														k_string//"max_err=", max_err
		!
		!	CONNECTION
		if(allocated(A_ka)) then
			allo_lst	=	trim(allo_lst) // "conn, "
			conn	=	velo_is_herm(A_ka, max_err)
			is_herm =	conn
			if(.not. conn) 		write(*,'(a,f16.7)')	"[check_H_gauge_herm/DEBUG-MODE]:	"								//	&
															"WARNING (H)-gauge A_ka IS NOT hermitian at rel. kpt= "			//	&
															k_string//"max_err=", max_err
		else
			write(*,'(a,f16.7)')	"[check_H_gauge_herm/DEBUG-MODE]: NOTE	connection was not calculated"
		end if
		!
		!
		!	CURVATURE
		if(allocated(Om_kab)) then
			allo_lst	=	trim(allo_lst) // "curv, "
			curv	=	curv_is_herm( Om_kab, max_err)
			is_herm =	is_herm .and. curv
			if(.not. curv) 		write(*,'(a,f16.7)')	"[check_H_gauge_herm/DEBUG-MODE]:	"						 		//	&
															"WARNING (H)-gauge Om_kab IS NOT hermitian at rel. kpt= "		//	&
															k_string//"max_err=", max_err
		else
			write(*,'(a,f16.7)')	"[check_H_gauge_herm/DEBUG-MODE]: NOTE	curvature was not calculated"
		end if
		!
		!
		!	SUCCESS MESSAGE
		if(	is_herm ) 	write(*,'(a,i8)')	"[check_H_gauge_herm/DEBUG-MODE]:	"						//	&
															"SUCCESS (H)-gauge quantities ("//trim(allo_lst)//") are hermitian "	//  &
															" at  #kpt= ", kpt_idx
		!
		return
	end subroutine

end module wann_interp

max_stars_repo_name: Compiler-max/MEPinterp


text:
! RUN: %python %S/test_folding.py %s %flang_fc1
! Tests folding of SAME_TYPE_AS() and EXTENDS_TYPE_OF()
module m

  type :: t1
    real :: x
  end type
  type :: t2(k)
    integer, kind :: k
    real(kind=k) :: x
  end type
  type :: t3
    real :: x
  end type
  type, extends(t1) :: t4
    integer :: y
  end type

  type(t1) :: x1, y1
  type(t2(4)) :: x24, y24
  type(t2(8)) :: x28
  type(t3) :: x3
  type(t4) :: x4
  class(t1), allocatable :: a1
  class(t3), allocatable :: a3

  logical, parameter :: test_1 = same_type_as(x1, x1)
  logical, parameter :: test_2 = same_type_as(x1, y1)
  logical, parameter :: test_3 = same_type_as(x24, x24)
  logical, parameter :: test_4 = same_type_as(x24, y24)
  logical, parameter :: test_5 = same_type_as(x24, x28) ! ignores parameter
  logical, parameter :: test_6 = .not. same_type_as(x1, x3)
  logical, parameter :: test_7 = .not. same_type_as(a1, a3)

  logical, parameter :: test_11 = extends_type_of(x1, y1)
  logical, parameter :: test_12 = extends_type_of(x24, x24)
  logical, parameter :: test_13 = extends_type_of(x24, y24)
  logical, parameter :: test_14 = extends_type_of(x24, x28) ! ignores parameter
  logical, parameter :: test_15 = .not. extends_type_of(x1, x3)
  logical, parameter :: test_16 = .not. extends_type_of(a1, a3)
  logical, parameter :: test_17 = .not. extends_type_of(x1, x4)
  logical, parameter :: test_18 = extends_type_of(x4, x1)
end module

max_stars_repo_name: ornata/llvm-project

max_issues_repo_name: ornata/llvm-project

max_forks_repo_name: ornata/llvm-project


ext: f

max_stars_repo_name: manoelfranca/cilppp

max_issues_repo_name: manoelfranca/cilppp

max_forks_repo_name: manoelfranca/cilppp


MODULE dynldf_lap_blp
   !!======================================================================
   !!                   ***  MODULE  dynldf_lap_blp  ***
   !! Ocean dynamics:  lateral viscosity trend (laplacian and bilaplacian)
   !!======================================================================
   !! History : 3.7  ! 2014-01  (G. Madec, S. Masson)  Original code, re-entrant laplacian
   !!----------------------------------------------------------------------

   !!----------------------------------------------------------------------
   !!   dyn_ldf_lap   : update the momentum trend with the lateral viscosity using an iso-level   laplacian operator
   !!   dyn_ldf_blp   : update the momentum trend with the lateral viscosity using an iso-level bilaplacian operator
   !!----------------------------------------------------------------------
   USE oce            ! ocean dynamics and tracers
   USE dom_oce        ! ocean space and time domain
   USE ldfdyn         ! lateral diffusion: eddy viscosity coef.
   USE ldfslp         ! iso-neutral slopes
   USE zdf_oce        ! ocean vertical physics
   !
   USE in_out_manager ! I/O manager
   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)

   IMPLICIT NONE
   PRIVATE

   PUBLIC dyn_ldf_lap  ! called by dynldf.F90
   PUBLIC dyn_ldf_blp  ! called by dynldf.F90

   !! * Substitutions
   !!----------------------------------------------------------------------
   !!                   ***  vectopt_loop_substitute  ***
   !!----------------------------------------------------------------------
   !! ** purpose :   substitute the inner loop start/end indices with CPP macro
   !!                allow unrolling of do-loop (useful with vector processors)
   !!----------------------------------------------------------------------
   !!----------------------------------------------------------------------
   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
   !! $Id: vectopt_loop_substitute.h90 10068 2018-08-28 14:09:04Z nicolasmartin $
   !! Software governed by the CeCILL license (see ./LICENSE)
   !!----------------------------------------------------------------------
   !!----------------------------------------------------------------------
   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
   !! $Id: dynldf_lap_blp.F90 10068 2018-08-28 14:09:04Z nicolasmartin $
   !! Software governed by the CeCILL license (see ./LICENSE)
   !!----------------------------------------------------------------------
CONTAINS

   SUBROUTINE dyn_ldf_lap( kt, pub, pvb, pua, pva, kpass )
      !!----------------------------------------------------------------------
      !!                     ***  ROUTINE dyn_ldf_lap  ***
      !!
      !! ** Purpose :   Compute the before horizontal momentum diffusive
      !!      trend and add it to the general trend of momentum equation.
      !!
      !! ** Method  :   The Laplacian operator apply on horizontal velocity is
      !!      writen as :   grad_h( ahmt div_h(U )) - curl_h( ahmf curl_z(U) )
      !!
      !! ** Action : - pua, pva increased by the harmonic operator applied on pub, pvb.
      !!----------------------------------------------------------------------
      INTEGER                         , INTENT(in   ) ::   kt         ! ocean time-step index
      INTEGER                         , INTENT(in   ) ::   kpass      ! =1/2 first or second passage
      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pub, pvb   ! before velocity  [m/s]
      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pua, pva   ! velocity trend   [m/s2]
      !
      INTEGER  ::   ji, jj, jk   ! dummy loop indices
      REAL(wp) ::   zsign        ! local scalars
      REAL(wp) ::   zua, zva     ! local scalars
      REAL(wp), DIMENSION(jpi,jpj) ::   zcur, zdiv
      !!----------------------------------------------------------------------
      !
      IF( kt == nit000 .AND. lwp ) THEN
         WRITE(numout,*)
         WRITE(numout,*) 'dyn_ldf : iso-level harmonic (laplacian) operator, pass=', kpass
         WRITE(numout,*) '~~~~~~~ '
      ENDIF
      !
      IF( kpass == 1 ) THEN   ;   zsign =  1._wp      ! bilaplacian operator require a minus sign
      ELSE                    ;   zsign = -1._wp      !  (eddy viscosity coef. >0)
      ENDIF
      !
      !                                                ! ===============
      DO jk = 1, jpkm1                                 ! Horizontal slab
         !                                             ! ===============
         DO jj = 2, jpj
            DO ji = 2, jpi   ! vector opt.
               !                                      ! ahm * e3 * curl  (computed from 1 to jpim1/jpjm1)
!!gm open question here : e3f  at before or now ?    probably now...
!!gm note that ahmf has already been multiplied by fmask
               zcur(ji-1,jj-1) = ahmf(ji-1,jj-1,jk) * e3f_n(ji-1,jj-1,jk) * r1_e1e2f(ji-1,jj-1)       &
                  &     * (  e2v(ji  ,jj-1) * pvb(ji  ,jj-1,jk) - e2v(ji-1,jj-1) * pvb(ji-1,jj-1,jk)  &
                  &        - e1u(ji-1,jj  ) * pub(ji-1,jj  ,jk) + e1u(ji-1,jj-1) * pub(ji-1,jj-1,jk)  )
               !                                      ! ahm * div        (computed from 2 to jpi/jpj)
!!gm note that ahmt has already been multiplied by tmask
               zdiv(ji,jj)     = ahmt(ji,jj,jk) * r1_e1e2t(ji,jj) / e3t_b(ji,jj,jk)                                         &
                  &     * (  e2u(ji,jj)*e3u_b(ji,jj,jk) * pub(ji,jj,jk) - e2u(ji-1,jj)*e3u_b(ji-1,jj,jk) * pub(ji-1,jj,jk)  &
                  &        + e1v(ji,jj)*e3v_b(ji,jj,jk) * pvb(ji,jj,jk) - e1v(ji,jj-1)*e3v_b(ji,jj-1,jk) * pvb(ji,jj-1,jk)  )
            END DO
         END DO
         !
         DO jj = 2, jpjm1                             ! - curl( curl) + grad( div )
            DO ji = 2, jpim1   ! vector opt.
               pua(ji,jj,jk) = pua(ji,jj,jk) + zsign * (                                                 &
                  &              - ( zcur(ji  ,jj) - zcur(ji,jj-1) ) * r1_e2u(ji,jj) / e3u_n(ji,jj,jk)   &
                  &              + ( zdiv(ji+1,jj) - zdiv(ji,jj  ) ) * r1_e1u(ji,jj)                     )
                  !
               pva(ji,jj,jk) = pva(ji,jj,jk) + zsign * (                                                 &
                  &                ( zcur(ji,jj  ) - zcur(ji-1,jj) ) * r1_e1v(ji,jj) / e3v_n(ji,jj,jk)   &
                  &              + ( zdiv(ji,jj+1) - zdiv(ji  ,jj) ) * r1_e2v(ji,jj)                     )
            END DO
         END DO
         !                                             ! ===============
      END DO                                           !   End of slab
      !                                                ! ===============
      !
   END SUBROUTINE dyn_ldf_lap

   SUBROUTINE dyn_ldf_blp( kt, pub, pvb, pua, pva )
      !!----------------------------------------------------------------------
      !!                 ***  ROUTINE dyn_ldf_blp  ***
      !!
      !! ** Purpose :   Compute the before lateral momentum viscous trend
      !!              and add it to the general trend of momentum equation.
      !!
      !! ** Method  :   The lateral viscous trends is provided by a bilaplacian
      !!      operator applied to before field (forward in time).
      !!      It is computed by two successive calls to dyn_ldf_lap routine
      !!
      !! ** Action :   pta   updated with the before rotated bilaplacian diffusion
      !!----------------------------------------------------------------------
      INTEGER                         , INTENT(in   ) ::   kt         ! ocean time-step index
      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pub, pvb   ! before velocity fields
      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pua, pva   ! momentum trend
      !
      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zulap, zvlap   ! laplacian at u- and v-point
      !!----------------------------------------------------------------------
      !
      IF( kt == nit000 )  THEN
         IF(lwp) WRITE(numout,*)
         IF(lwp) WRITE(numout,*) 'dyn_ldf_blp : bilaplacian operator momentum '
         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
      ENDIF
      !
      zulap(:,:,:) = 0._wp
      zvlap(:,:,:) = 0._wp
      !
      CALL dyn_ldf_lap( kt, pub, pvb, zulap, zvlap, 1 )   ! rotated laplacian applied to ptb (output in zlap)
      !
      CALL lbc_lnk_multi( zulap, 'U', -1., zvlap, 'V', -1. )             ! Lateral boundary conditions
      !
      CALL dyn_ldf_lap( kt, zulap, zvlap, pua, pva, 2 )   ! rotated laplacian applied to zlap (output in pta)
      !
   END SUBROUTINE dyn_ldf_blp

   !!======================================================================
END MODULE dynldf_lap_blp


text:
address(The Memorial Union, 38.542292, 121.749440) is one of the main buildings on UC Davis Campus, and arguably the most publicly traveled building.  Often dubbed the MU, it contains a variety of interesting places and things. Students pay a yearly http://budget.ucdavis.edu/studentfees/documents/special/studentfeeoverview.pdf Memorial Union Fee of $85.50. Open 24 Hours during the week of finals as to provide students space to study.

 History

According to the http://vtour.ucdavis.edu/ UCD Vtour, it was dedicated in 1955 to honor 128 UC Davis students and alumni who lost their lives in Veterans military service.  Griffin Lounge houses the Golden Memory Book, which lists all of the names of former UC Davis students who died in the armed services.

At least part of the http://www.ucdavis.edu/history.html $1 million cost was raised by the thennew http://www.ucdavisalumni.com/mission.htm Cal Aggie Alumni Assocation (CAAA).  It was http://sunsite.berkeley.edu/uchistory/general_history/campuses/ucd/housing.html built on the site of West Hall, which was constructed in 1917 and housed 67 men until 1955.

When originally built the MU was only a http://ulibimage.ucdavis.edu/speccoll/east01/full/B9328.jpg single story building, housing a dining commons, kitchen, faculty club, student lounge, Alumni Association offices, a soda fountain and a student store.   But fund raising campaigns from 1961 to 1968 added the second, third and fourth stories along with the MU Games Area and the UC Davis Bookstore.  It was last renovated in 1996 to bring the seismic safety rating from poor to good, at a cost of $8.4 million.  The 4 story concrete lattice that is visible from the south side of the MU was added at that time to reinforce the building.

In spring of 2007, one of the largest flashmob flashmobs in US history occurred at the ATM machines.

 Areas in the MU

The MU hosts to several commercial operations, some owned and operated by UC Davis itself, other private businesses that lease space from UC Davis via the MU Auxillary Services department (homepage: http://wwwmuas.ucdavis.edu/).  Other space in the union is devoted to University administrative offices related to students, student government, and various meeting rooms.

 First Floor

  UC Davis Bookstore. We wouldnt be a campus without a place to gouge the hell out of us on textbook prices, right?
  Coffee House ASUCD Coffee House. The CoHo is the place to get some yums (cookies, cakes, etc). (I guess it http://wwwucdmag.ucdavis.edu/fall98/EndNotes.html used to be in East Hall!) http://www.coffeehouse.ucdavis.edu CoHo website.
  Aggie Card Office. This is the place to get an Aggie Card. The Registrar does not issue cards.
  Aggie Student Store. Affectionately known as the ASS. Sells candy and, more importantly, $1 Icees.
  Classical Notes/Campus Copies
  MU Information Center. A place for general information, maps, pamphlets and the Lost and Found
  Post Office. A small post office which recently moved from an even smaller location in the MU.
  STA Travel. A corporate travel agency that caters to students. Company home page: http://www.statravel.com/
  ASUCD Student Center
  The Corral  A store selling stamps, scantrons, blue books an other knickknacks.
  Griffin Lounge  A nice place to relax amidst the hustle and bustle.
  Memorial Union Early Polling Station A polling station open a week before most elections.
  ASUCD Bike Hub. This small store sells various bike accessories.
 Aggie ReStore. Here you can find a variety of things that have been donated or found throughout campus.
 Second Floor
  DeCarli Room
  Garrison Room
  Memorial Union Art Gallery
  Fielder Room
  Smith room
  MU II
  King Lounge
  Campus Unions Operations
  Student Affairs Event Services

 Third Floor
  ASUCD President Offices of the ASUCD President (ASUCD President, ASUCD Vice President, ASUCD Controller and Executive Office Staff)
  ASUCD Senate Office
  Angelina Malfitano Conference Room
  ASUCD Student Government Advisors Office
  SGAO ASUCD Student Government Administration Office
  ASUCD Student Services Office (Payroll, Accounts Payable & Receivable, ASUCD Business Manager)
  Mee Room
  Moss Room
  CalPIRG Office
  ASUCD Advocacy Units Office (City/County Affairs, Lobby Corps, University Affairs)

 Fourth Floor
  Center for Student Involvement
  Campus Events & Visitor Services
  Campus Unions

 Fifth Floor
  Sky Room

 Basement
Note: the MU Basement connects with Lower Freeborn Hall
  MU Games Area. Arcade, bowling and pool tables.

 Outdoor Meeting places
   The Memorial Union: Flagpole Flagpole
   The MU Patio
These are the default meeting places for many students, as they are more or less centrally located on campus and are well known to students by their first week at UCD.

 Other points of interest

  The UCD Bathrooms Cleanest Bathrooms in the building are on the second floor and up.  Avoid the first floor bathrooms if you can.
  There is a Bicycle Tire Pumps bike pump located on the west side of the Catering trailer near the bike racks.
  The Stanford Tree has been sighted here.
  The elevators and stairs can connect you to Lower Freeborn Hall Freeborn, and vice versa.
  The wind chimes on the north side are nice to listen to, particularly during the spring.

 Surveillance Cameras

On April 5, 2007, Campus Unions installed a number of surveillance cameras in the Memorial Union sending a strong message to everyone who enters the building, You are being watched.

 Former places in the MU
  University Haircutters. A very nice barbershop... on campus! Barbers enjoy talking about fly fishing. Now Closed as John has retired.
  Bookstore Special Events Room  unless this still exists?

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
       function dev_2(T)
        implicit none

        type(Tensor2), intent(in) :: T
        type(Tensor2) :: dev_2,Eye

        Eye = identity2(Eye)
        dev_2%ab = T%ab - tr(T)/3.d0*Eye%ab

       end function dev_2

       function dev_2s(T)
        implicit none

        type(Tensor2s), intent(in) :: T
        type(Tensor2s) :: dev_2s,Eye

        Eye = identity2(Eye)
        dev_2s%a6 = T%a6 - tr(T)/3.d0*Eye%a6

       end function dev_2s

max_stars_repo_name: kengwit/ttb

max_issues_repo_name: kengwit/ttb

max_forks_repo_name: kengwit/ttb


text:
!
!     CalculiX - A 3-dimensional finite element program
!              Copyright (C) 1998-2020 Guido Dhondt
!
!     This program is free software; you can redistribute it and/or
!     modify it under the terms of the GNU General Public License as
!     published by the Free Software Foundation(version 2);
!
!
!     This program is distributed in the hope that it will be useful,
!     but WITHOUT ANY WARRANTY; without even the implied warranty of
!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!     GNU General Public License for more details.
!
!     You should have received a copy of the GNU General Public License
!     along with this program; if not, write to the Free Software
!     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
!
      subroutine calcstressheatfluxfem(kon,lakon,ipkon,
     &  ielmat,ntmat_,vold,matname,mi,shcon,nshcon,
     &  turbulent,compressible,ipvar,var,sti,qfx,cocon,
     &  ncocon,ne,isti,iqfx)
!
!     calculating the stresses and heat flow in the integration points
!
      implicit none
!
      character*8 lakon(*)
      character*80 matname(*),amat
!
      integer kon(*),nelem,mi(*),ne,konl(20),ipkon(*),j,i1,i2,j1,ii,
     &  jj,
     &  indexe,isti,iqfx,kk,ielmat(mi(3),*),nshcon(*),ntmat_,nope,
     &  imat,ncocon(2,*),mint3d,k1,ipvar(*),index,turbulent,compressible
!
      real*8 shp(4,20),dvi,cond,cocon(0:6,ntmat_,*),div,
     &  shcon(0:3,ntmat_,*),voldl(0:mi(2),20),vold(0:mi(2),*),temp,
     &  tt(3,3),rho,t(3,3),vkl(3,3),xkin,unt,umt,
     &  var(*),sti(6,mi(1),*),dtem(3),qfx(3,mi(1),*)
!
      do nelem=1,ne
!
         if(ipkon(nelem).lt.0) cycle
         if(lakon(nelem)(1:1).ne.'F') cycle
         indexe=ipkon(nelem)
!
         imat=ielmat(1,nelem)
         amat=matname(imat)
!
         if(lakon(nelem)(4:4).eq.'4') then
            nope=4
            mint3d=1
         elseif(lakon(nelem)(4:4).eq.'6') then
            nope=6
            mint3d=2
         elseif(lakon(nelem)(4:5).eq.'8R') then
            nope=8
            mint3d=1
         elseif(lakon(nelem)(4:4).eq.'8') then
            nope=8
            mint3d=8
         elseif(lakon(nelem)(4:5).eq.'10') then
            nope=10
            mint3d=4
         elseif(lakon(nelem)(4:5).eq.'15') then
            nope=15
            mint3d=9
         elseif(lakon(nelem)(4:6).eq.'20R') then
            nope=20
            mint3d=8
         elseif(lakon(nelem)(4:4).eq.'2') then
            nope=20
            mint3d=27
         else
            nope=0
            mint3d=0
         endif
!
         do j=1,nope
            konl(j)=kon(indexe+j)
         enddo
!
!     storing the local temperature and velocity
!
         do i1=1,nope
            do i2=0,3
               voldl(i2,i1)=vold(i2,konl(i1))
            enddo
         enddo
!
!     computation of the matrix: loop over the Gauss points
!
         index=ipvar(nelem)
c                  write(*,*) 'index ',index
         do kk=1,mint3d
!
!     copying the shape functions and their derivatives from field var
!
            do jj=1,nope
               do ii=1,4
                  index=index+1
                  shp(ii,jj)=var(index)
c                  write(*,*) 'shp ',jj,ii,shp(ii,jj)
               enddo
            enddo
!
!     calculating of
!     the velocity gradient vkl
!     the temperature gradient dtem
!
            if(isti.gt.0) then
               do i1=1,3
                  do j1=1,3
                     vkl(i1,j1)=0.d0
                  enddo
               enddo
               do i1=1,nope
                  do j1=1,3
                     do k1=1,3
                        vkl(j1,k1)=vkl(j1,k1)+shp(k1,i1)*voldl(j1,i1)
                     enddo
                  enddo
               enddo
               if(compressible.eq.1) div=vkl(1,1)+vkl(2,2)+vkl(3,3)
            endif
c               write(*,*) 'calcstressheatflux ',kk,nelem,
c     &             ((vkl(i1,j1),j1=1,3),i1=1,3)
!
            if(iqfx.gt.0) then
               do i1=1,3
                  dtem(i1)=0.d0
               enddo
               do i1=1,nope
                  do j1=1,3
                     dtem(j1)=dtem(j1)+shp(j1,i1)*voldl(0,i1)
                  enddo
               enddo
            endif
!
!     storing shpv, vel and temp
!
            index=index+nope+2
            temp=var(index)
c               write(*,*) 'calcstressheatflux ',kk,nelem,
c     &             temp
!
!     determining the dissipative stress
!
            if(isti.gt.0) then
              call materialdata_dvifem(imat,ntmat_,temp,shcon,nshcon,
     &             dvi)
               do i1=1,3
                  do j1=i1,3
                     t(i1,j1)=vkl(i1,j1)+vkl(j1,i1)
                  enddo
                  if(compressible.eq.1) t(i1,i1)=t(i1,i1)-2.d0*div/3.d0
               enddo
!
!     calculating the stress
!
               if(turbulent.ne.0) then
!
                  rho=var(index+8)
                  xkin=var(index+10)
                  unt=var(index+12)
!
                  umt=unt*rho
!
!     calculating the turbulent stress
!
                  do i1=1,3
                     do j1=i1,3
                        tt(i1,j1)=umt*t(i1,j1)
                     enddo
                     tt(i1,i1)=tt(i1,i1)-2.d0*rho*xkin/3.d0
                  enddo
!
!     adding the viscous stress
!
                  do i1=1,3
                     do j1=i1,3
                        t(i1,j1)=dvi*t(i1,j1)+tt(i1,j1)
                     enddo
                  enddo
               else
!
                  do i1=1,3
                     do j1=i1,3
                        t(i1,j1)=dvi*t(i1,j1)
                     enddo
                  enddo
               endif
!
!     storing the total stress
!
               sti(1,kk,nelem)=t(1,1)
               sti(2,kk,nelem)=t(2,2)
               sti(3,kk,nelem)=t(3,3)
               sti(4,kk,nelem)=t(1,2)
               sti(5,kk,nelem)=t(1,3)
               sti(6,kk,nelem)=t(2,3)
c               write(*,*) 'calcstressheatflux ',kk,nelem,
c     &             (sti(i1,kk,nelem),i1=1,6)
            endif
!
!     storing the heat flow
!
            if(iqfx.gt.0) then
               call materialdata_cond(imat,ntmat_,temp,cocon,
     &               ncocon,cond)
               qfx(1,kk,nelem)=-cond*dtem(1)
               qfx(2,kk,nelem)=-cond*dtem(2)
               qfx(3,kk,nelem)=-cond*dtem(3)
            endif
!
            index=index+13
!
         enddo
      enddo
!
      return
      end

max_stars_repo_name: alleindrach/calculix-desktop

max_issues_repo_name: alleindrach/calculix-desktop

max_forks_repo_name: alleindrach/calculix-desktop


text:
c     ****************************************************************
c     *                                                              *
c     *                      subroutine gp_temps                     *
c     *                                                              *
c     *                       written by : rhd                       *
c     *                                                              *
c     *                   last modified : 8/17/2020 rhd              *
c     *                                                              *
c     *     compute increment of temperature at a gauss point for    *
c     *     all elements in the block.                               *
c     *     compute current temperature at end of step for           *
c     *     all elements of the block.                               *
c     *                                                              *
c     ****************************************************************
c
c
      subroutine gauss_pt_temps(
     &       dtemps_node_blk, gpn, etype, span, int_order,
     &       nnodel, gp_dtemps, temps_node_blk, gp_temps,
     &       temper_increment, temps_node_to_process,
     &       temperatures_init, temps_init_node_blk, gp_rtemps  )
      implicit none
      include 'param_def'
c
c              parameter declarations
c
      integer, intent(in) :: gpn, etype, span, int_order, nnodel
      double precision, intent(inout) :: gp_dtemps(*), gp_temps(*),
     &        gp_rtemps(*)
      double precision, intent(in) :: dtemps_node_blk(mxvl,*),
     &        temps_node_blk(mxvl,*), temps_init_node_blk(mxvl,*)
      logical, intent(in) :: temper_increment, temps_node_to_process,
     &                       temperatures_init
c
c              locally defined arrays-variables
c
      integer :: i, k, enode
      double precision :: sf(mxndel), xi, eta, zeta, weight
      double precision, parameter :: zero = 0.0d0
      logical, parameter :: local_debug = .false.
c
c              if there are incremental temperature changes
c              imposed in this load step, interpolate temperature
c              change at this gauss point for all elements
c              in the block from nodal values.
c
c              get the parametric coordinates for this
c              integration point. then get the nodal
c              shape functions evaluated at the point.
c
      if( temper_increment .or. temps_node_to_process .or.
     &     temperatures_init ) then
         call getgpts( etype, int_order, gpn, xi, eta, zeta, weight )
         call shapef( etype, xi, eta, zeta, sf(1) )
      end if
c
!DIR$ VECTOR ALIGNED
      do i = 1, span
         gp_dtemps(i) = zero
         gp_temps(i)  = zero
         gp_rtemps(i) = zero
      end do
c
      if( temper_increment ) then
         do enode = 1, nnodel
!DIR$ VECTOR ALIGNED
           do i = 1, span
             gp_dtemps(i) = gp_dtemps(i) +
     &                      sf(enode) * dtemps_node_blk(i,enode)
           end do
         end do
      end if
c
c              interpolate the current temperature at the end
c              of this load step for this gauss point for all elements
c              in the block from nodal values.
c
      if( temps_node_to_process ) then
        do enode = 1, nnodel
!DIR$ VECTOR ALIGNED
          do i = 1, span
            gp_temps(i) = gp_temps(i) +
     &                    sf(enode) * temps_node_blk(i,enode)
          end do
        end do
      end if
c
c              interpolate initial temperature
c              for this gauss point for all elements
c              in the block from nodal values.
c
      if( temperatures_init ) then
        do enode = 1, nnodel
!DIR$ VECTOR ALIGNED
          do i = 1, span
            gp_rtemps(i) = gp_rtemps(i) +
     &                     sf(enode) * temps_init_node_blk(i,enode)
          end do
        end do
      end if
c
      if ( .not. local_debug ) return
         write(*,*) '>> in  gauss_pt_temps:'
         write(*,*) '>> temps_init_node_block:'
         do i = 1, span
          write(*,*) '   element: ',i
          write(*,*) (k,temps_init_node_blk(i,k),k=1,nnodel)
         end do
         write(*,*) '>> gp_temps...'
         write(*,9000) (i,gp_temps(i),i=1,span)
         write(*,*) '>> gp_dtemps...'
         write(*,9000) (i,gp_dtemps(i),i=1,span)
         write(*,*) '>> gp_rtemps...'
         write(*,9000) (i,gp_rtemps(i),i=1,span)
      return
c
 9000 format(1x,i3,f15.6 )
      end
c     ****************************************************************
c     *                                                              *
c     *                      subroutine gp_temp_eps                  *
c     *                                                              *
c     *                       written by : rhd                       *
c     *                                                              *
c     *                   last modified : 8/25/2020 rhd              *
c     *                                                              *
c     *     compute incremental thermal strains at gauss point for   *
c     *     all elements of the block. subtract them from incr.      *
c     *     strains due to displacements.                            *
c     *                                                              *
c     ****************************************************************
c
c
      subroutine gp_temp_eps( span, deps, alpha_n1,  alpha_n,
     &      alpha_zero, temps_n1, dtemps, temps_initial, etype  )
      use main_data,  only : bar_types, link_types
      implicit none
      include 'param_def'
c
c              parameter declarations
c
      integer, intent(in)  :: span, etype
      double precision, intent(inout) :: deps(mxvl,*)
      double precision, intent(in) :: alpha_n1(mxvl,6),
     &                    alpha_n(mxvl,6), temps_n1(*),
     &                    temps_initial(*),
     &                    dtemps(*), alpha_zero
c
c              locals
c
      integer :: i, j
      double precision :: temps_n(mxvl), deps_theta(mxvl,6), tn1, tn
      double precision, parameter :: zero = 0.0d0
      logical :: is_bar_elem, is_link_elem
      logical, parameter :: local_debug = .false.
c
c              compute incremental thermal strain and subtract from
c              strain increment due to displacements.
c
      if( local_debug ) write(*,9100) (i,(deps(i,j),j=1,6),i=1,span)
      is_bar_elem  = bar_types(etype)
      is_link_elem = link_types(etype)
c
c              \Delta\veps_\theta =
c                  \alpha_{n+1) ( T_{n+1) - alpha_zero)
c                              -  \alpha_n (T_n - alpha_zero)
c
c              \alpha here can be anisotropic. alpha_zero adjustment
c              included but only applies to isotropic, temperature
c              dependent alpha values. default = 0
c
c              alpha_zero. same value for all elements, gpts
c              in block.
c
c              if the user sets a non-zero alpha_zero value
c              for temperature invariant alphas, the alpha_zero
c              terms for the incremental thermal strain cancel as
c              easily see below.
c
c              The user-specified initial temperature at t=0
c              passed in here but not used in current formulation.
c
!DIR$ VECTOR ALIGNED
      do i = 1, span
        temps_n(i) = temps_n1(i) - dtemps(i)
      end do
c
!DIR$ VECTOR ALIGNED
      do i = 1, span
       tn  = temps_n(i)  - alpha_zero
       tn1 = temps_n1(i) - alpha_zero
       deps_theta(i,1:6) = alpha_n1(i,1:6)*tn1 - alpha_n(i,1:6)*tn
      end do
!DIR$ VECTOR ALIGNED
      do i = 1, span
       deps(i,1:6) = deps(i,1:6) - deps_theta(i,1:6)
      end do
c
      if( is_bar_elem .or. is_link_elem ) deps(1:span,2:6) = zero
c
      if ( .not. local_debug ) return
      write(*,*) '>> thermal strain increments...'
      do i = 1, span
       write(*,9200) i, deps_theta(i,1:6)
      end do
      write(*,9000) (i,(deps(i,j),j=1,6),i=1,span)
      return
c
 9000 format(1x,'>> leaving  gp_temp_eps...',
     & /,'element strain increments after subtracting thermal',
     & ' increment...',/
     & (10x,i3,6e17.6) )
 9100 format(1x,'>> entering gp_temp_eps...',
     & /,'element strain increments before subtracting thermal',
     & ' increment...',/
     & (10x,i3,6e17.6) )
 9200 format(10x,i3,6e17.6)
c
      end

c     ****************************************************************
c     *                                                              *
c     *                      subroutine gp_temp_eps_n1               *
c     *                                                              *
c     *                       written by : rhd                       *
c     *                                                              *
c     *                   last modified : 8/27/2018 rhd              *
c     *                                                              *
c     *     compute total thermal strains at gauss point for         *
c     *     all elements of the block.                               *
c     *                                                              *
c     *   -- this routine not used anywhere in WARP3D. if used in    *
c     *      future it needs to be re-written to properly include    *
c     *      alpha_zero effects. Abaqus calls this term              *
c     *      *EXPANSION,ZERO                                         *
c     *                                                              *
c     ****************************************************************
c
c
      subroutine gp_temp_eps_n1( span, eps_theta_n1, alpha_n1,
     &                           gp_temps, gp_rtemps, etype  )
      use main_data,  only : bar_types, link_types
      implicit none
      include 'param_def'
c
c                      parameter declarations
c
      integer :: span, etype
      double precision :: eps_theta_n1(span,6), alpha_n1(mxvl,6),
     &                    gp_temps(*), gp_rtemps(*)
c
c                     locals
c
      integer :: i
      double precision :: dtn1
      double precision, parameter :: zero = 0.0d0
      logical :: is_bar_elem, is_link_elem
      logical, parameter :: local_debug = .false.
c
      if( local_debug ) write(*,9100)
c
c            compute (total) thermal strain at n+1
c
      is_bar_elem  = bar_types(etype)
      is_link_elem = link_types(etype)
c
c           instantaneous thermal strain is always just
c           current temperature * current  CTEs
c           adjust for reference temp.
c           \alpha here can be anisotropic
c
!DIR$ VECTOR ALIGNED
      do i = 1, span
         dtn1 = gp_temps(i) -  gp_rtemps(i)
         eps_theta_n1(i,1) =  alpha_n1(i,1)*dtn1
         eps_theta_n1(i,2) =  alpha_n1(i,2)*dtn1
         eps_theta_n1(i,3) =  alpha_n1(i,3)*dtn1
         eps_theta_n1(i,4) =  alpha_n1(i,4)*dtn1
         eps_theta_n1(i,5) =  alpha_n1(i,5)*dtn1
         eps_theta_n1(i,6) =  alpha_n1(i,6)*dtn1
      end do
c
      if( is_bar_elem .or. is_link_elem )
     &     eps_theta_n1(1:span,2:6) = zero
c
      if ( .not. local_debug ) return
      write(*,*) '>> thermal strains @ n+1...'
      do i = 1, span
       write(*,9200) i, eps_theta_n1(i,1:6)
      end do
      return
c
 9100 format(1x,'>> entering gp_temp_eps_n1...' )
 9200 format(10x,i3,6e17.6)
c
      end


text:
The UC Davis English Department is located in Voorhies Hall.  Their major program consists of  three different emphases: a General Emphasis, a Teaching Emphasis, and a Creative Writing emphasis.

If you have interest in pursuing a major in English, then check out the http://wwwenglish.ucdavis.edu/under_desc.htm#major Undergraduate Program requirements (or check the General Catalog).  For a full list of courses in the department, please consult the General Catalog.

Some people think that http://www.californiaaggie.com/media/paper981/news/2005/02/18/Features/Change.Your.Major1319638.shtml English is the best major at Davis. Some people think it will leave you unable to get a Employment job.  Regardless of your opinion, it is nice to have some Poetry pretty words in your life.

Yiyun Li, a member of the Department since 2008, was awarded a prestigious MacArthur Foundation fellowship in 2010.

 Advising
Staff Advisor: Lynda Jones
Office Hours: MF 812, 15
Office: Voorhies 177
Email: lmjones@ucdavis.edu

Peer Advisor: TBA
Office Hours: TBA
Office: Voorhies 112
Email: TBA

 Courses
Courses are divided between the English Departments http://wwwenglish.ucdavis.edu/coursched.htm literature and creative writing classes, for those interested in literature, and the http://writingprogram.ucdavis.edu/Expandedcoursedescription.htm composition courses that are taught as part of the University Writing Program which states on http://writingprogram.ucdavis.edu/ its official website that under the guidance of the Undergraduate Council of the U.C. Davis Academic Senate, the Writing Program has become an independent entity, separate from the English Department, but I believe most of the writing instructors are also part of the English department, though perhaps not all are.

English 180 (Childrens Literature)
    usually a very popular class with nonmajors, typical course reading lists will include books from fairy tales to Little Women to Harry Potter
    I took this class with Prof. Stenzel and would not recommend it. The course in general was unfocused, and Stenzel tries to fit too many texts into one quarter. Sometimes we would discuss a work for only half a class session or sometimes not at all. If you take this course, I would avoid Stenzel. Users/EyadDarras

English 101 (part of the Writing Program)
    I took it with Kerry Hanlon.  When I took the course it focused on creative writing assignments.  Our homework consisted of reading essays (taken from magazines like the New Yorker) and then answering questions on them.  Our assigned essays were to be modeled directly after one of the essays styles.  The essays we read were all gimmicky and focused on small stuff (but very well written and enjoyable to read)  the kind of things youd expect to read in a popular magazine.  Most unfortunately, the classroom I was in had laptops in it (which was surely seen as a positive, but I felt it was a distraction).  Additionally, I did not like the format my professor chose to take with the class.  Our __entire class discussion__ focused on answering the questions we were assigned for homework  questions we had already answered at home.  The discussions were not enlightening or meaningful, in my opinion.  There was also forced participation (10% of the grade) during class which I found to be childish and also, well, degrading to my course grade as I am very shy in class  especially when repeating material we had all already answered at home.

    I had very high hopes for this course.  As a Mathematics major I do not often get the opportunity to express myself creatively in essay form.  Unfortunately, I did not find the course to be stimulating.  I would recommend you take the competency test unless you can nab an awesome professor who will challenge you.  Users/PhilipNeustrom
    The composition courses are geared more towards making sure students grasp the mechanics of good writing, and while they do try to be as engaging and creative as possible, creative expression is not the main goal of the composition classes.  I would suggest the actual Creative Writing courses in Fiction, Poetry, and Drama from the English department. Users/IrenePark

    It also depends on the instructor. I took this course circa 2001 and it was fantastic.  Unfortunately, I am completely blanking on the instructors name right now, but she was not only engaging, but an effective teacher.  She also seemed to really care about us developing our writing style and finding our own voice.  We had a variety of writing assignments and she provided us with plenty of support and constructive feedback.  She didnt just give you a grade; she included plenty of specific comments.  The best part was that it was a tiny class.  I enjoyed it so much that I got hooked.  It was the impetus that I needed to double major and pursue an English degree.  I think this department is one of the best kept secrets at UC Davis.  I dont know whether people are aware of this wonderful program.  And, also keep in mind:  English professors often are the ones who write the best letters of recommendation.  This is something to consider if you are applying to a graduate or professional program.  Admissions committees want to know that their students write clearly and effectively.   Users/Curlygirl26

 Recommended Teachers
   http://wwwenglish.ucdavis.edu/faculty/clover/index.html Joshua Clover:  An Associate professor without a PhD, Clover writes for various magazines, mostly music, under a pseudonym while often teaching classes that involve poetry.  Hes the author of several poems published in the Best American Poetry collection, as well as a book entitled The Matrix, which analyzes the film as an allegory for turnofthemillennium American life and labor. Clover is an incredibly challenging professor, but its definitely worth it to take his classes if youre genuinely interested in expanding your knowledge in English and poetry.

   http://wwwenglish.ucdavis.edu/faculty/Freeman/Freeman.htm Elizabeth Freeman: Insanely intelligent and eloquent with her speech, Freeman is most notable for her inventive teaching methods, which really drill lessons into students heads.  A great distinguishing characteristic is that Freeman is a techieshe knows how to use her computer, unlike most English profs, and she incorporates her skill into teaching. She offers really clear and detailed lecture slides and even sets up a Jeopardy game for a final review session.

   http://wwwenglish.ucdavis.edu/faculty/hicks/hicks.htm Jack Hicks: A sophisticated and hilarious man, Hicks, too, has had a working history with publications including TV Guide and Rolling Stone.  He often teaches ENL 182  California Literature. Unlike Clover, Hicks is a pretty laidback man who realizes college is just something you have to put up with to get somewhere in life, so he isnt too hard on grading.

   http://wwwenglish.ucdavis.edu/faculty/vanleer/vanleer.htm David Van Leer: looking just a tad bit like the comedic writer Bruce Vilanch (and just as funny and gay), Van Leer has a wide array of talents and interests, ranging from early American literature (read: Puritan lit) to alt. sexualities in film. His film studies class (if its available) has one of the best selection of obscure, independent films ever; and only in his early American literature class will you realize that even in a truly Puritan environment, you really cant wipe out individuality and nonconventional sexualities. David (never Professor Van Leer) left UCDavis after fall quarter 2009.

 http://wwwenglish.ucdavis.edu/faculty/zender/default.htm Karl Zender (Retired): If youre into Faulkner, take his 177 course. This man knows more about Faulkner than I ever thought possible. He is sweet and personable and a fair grader and knows so much about Faulkner that even if you dont read a single page of the novels assigned (which you should!), you will come away feeling that you have learned something.

 Teachers to Avoid
  http://wwwenglish.ucdavis.edu/faculty/corin/default.htm Lucy Corin: For literature & ShortStory classes: lessons usually unorganized; grading generally arbitrary; very difficult to come out of the class being sure one has learned a thing. Next: ENL 100F (creative writing) is a class that you must petition to join by including a portfolio of your writing. This class is cherished by many English majors, begging the questions why a number of students (including wiki users) have dropped the class after a few weeks when taught by Lucy Corin. Rumor has it one girl left crying. It should be noted that Lucy Corin is a very smart and friendly person. She will likely be an incredible teacher in a few years, once shes gotten a grasp of her clearlydifficult job. For now, its a good idea to avoid. In the future, I recommend checking the various professor rating sites to see if things have improved.
     Rumor control: the girl left crying due to another student; Lucy was not the cause. Also, Ive taken 100F with Lucy twice and highly enjoyed it both quarters.
  http://wwwenglish.ucdavis.edu/faculty/schleiner/schleiner.htm Schleiner: If Shakespeare could be any worse than it is already, Schleiner is the one to make it happen. Led a pretty easy class, but incredibly boring; His text interpretations are also iffy, as are his acting skills.
     Rumor control: The above comment was obviously made by someone who doesnt like Shakespeare. That persons attitude toward the material may have more to do with their opinion of the course than the Professor.  Users/EyadDarras
  Waddington:  Now retired, but forever infamous for student hating. The man would get distracted and or irritated extremely easily by students, and would stop his lecture to glare at them silently with wrathful vengeance. If a student came in late, or a cell phone went off, or a student was laughing because he looks like an old rabbit (good thing he didnt see my doodling), he would throw a fit, shake all over and stare at them with desperate hatred. On second thought, I liked him. I just didnt like his grading scale.

 Senate Faculty

  http://www.english.ucdavis.edu/faculty/faculty.html Faculty Website

  Don Abbot
  http://english.ucdavis.edu/people/directory/gbloom Gina Bloom
  Nathan Brown
  Seeta Chaganti
  Joshua Clover
  Lucy Corin
  Joanne Feit Diehl
  Greg Dobbins
  Frances E. Dolan
  Margaret Ferguson
  http://english.ucdavis.edu/people/directory/ksfreder Kathleen Frederickson
  Lynn Freed Lynn R. Freed
  Users/BethFreeman Elizabeth Freeman
  http://english.ucdavis.edu/people/directory/dheard Danielle Heard
  Jack Hicks
  Pam Houston
  http://english.ucdavis.edu/people/directory/hlhsu Hsuan Hsu
  http://english.ucdavis.edu/people/directory/mjerng Mark Jerng
  Alessa Johns
  Richard A. Levin
  Desirée Martín
  Timothy Morton
  Marijane Osborn
  David Robertson (Retired)
  Catherine Robson
  Winfried Schleiner
  Scott C. Shershow
  Scott Simmon
  David Simpson
  Matthew Stratton
  Claire Waters
  Evan Watkins
  Joe Wenderoth
  Alan Williamson Alan B. Williamson
  Michael Ziser


text:
c
c***************************************************************************
c
      subroutine checkx(myid, i, j, k, omgrf,
     .   nkx1, nkx2, nky1, nky2, nphi1, nphi2,
     .   fdk, fek, ffk, u, v, w, xjx, xjy, xjz,
     .   nkdim1, nkdim2, mkdim1, mkdim2, lkdim1, lkdim2,
     .   nnodex, nnodey, nnodephi, nxmx, nymx, nphimx,
     .   xprime, yprime, phiprime)

      implicit none

      integer nkx1, nkx2, nky1, nky2, nphi1, nphi2,
     .    nkdim1, nkdim2, mkdim1, mkdim2, lkdim1, lkdim2,
     .    nnodex, nnodey, nnodephi, nxmx, nymx, nphimx,
     .    i, j, k, n, m, nphi, myid

      real xprime(nxmx), yprime(nymx), phiprime(nphimx)
      real eps0, omgrf

      complex xjx(nxmx, nymx, nphimx),
     .        xjy(nxmx, nymx, nphimx),
     .        xjz(nxmx, nymx, nphimx)

      complex fdk(nkdim1 : nkdim2, mkdim1 : mkdim2, lkdim1 : lkdim2),
     .        fek(nkdim1 : nkdim2, mkdim1 : mkdim2, lkdim1 : lkdim2),
     .        ffk(nkdim1 : nkdim2, mkdim1 : mkdim2, lkdim1 : lkdim2)

      complex term1, term2, term3

      complex xlhs, xrhs
      complex u(nkdim1 : nkdim2, mkdim1 : mkdim2, lkdim1 : lkdim2),
     .        v(nkdim1 : nkdim2, mkdim1 : mkdim2, lkdim1 : lkdim2),
     .        w(nkdim1 : nkdim2, mkdim1 : mkdim2, lkdim1 : lkdim2)
      complex zi

      zi = cmplx(0.0, 1.0)
      eps0 = 8.85e-12

      if (myid .eq. 0) then
         write(15, 910)
         write(15, 930)
         write(15, 920)
         write(15, 940)
         write(15, 920)

         write(6, 910)
         write(6, 930)
         write(6, 920)
         write(6, 940)
         write(6, 920)

         term1 = 0.0
         term2 = 0.0
         term3 = 0.0

         do nphi = nphi1, nphi2
            do n = nkx1, nkx2
               do m = nky1, nky2

                  term1 = term1 + fdk(n, m, nphi) * u(n, m, nphi)
                  term2 = term2 + fek(n, m, nphi) * v(n, m, nphi)
                  term3 = term3 + ffk(n, m, nphi) * w(n, m, nphi)

               end do
            end do
         end do

         xlhs = term1 + term2 + term3
         xrhs = -zi / omgrf / eps0 * xjx(i, j, k)
         write(15, 900) i, j, xprime(i), yprime(j), xlhs, xrhs
         write(6,  900) i, j, xprime(i), yprime(j), xlhs, xrhs
      endif

      return
  900 format(i5, i5, 1p9e12.3)
  940 format(1h , 4h   i, 5h    j, 4x, 6h x(i) , 6x, 6h y(j) ,
     1                    6x,  6hre lhs, 6x, 6him lhs, 6x,
     1                         6hre rhs, 6x, 6him rhs)
  910 format (1h1)
  920 format (1h0)
  930 format(3x, 14h u diagnostics)
      end

c
c***************************************************************************
c
      subroutine checky(myid, i, j, k, omgrf,
     .   nkx1, nkx2, nky1, nky2, nphi1, nphi2,
     .   fgk, fak, fpk, u, v, w, xjx, xjy, xjz,
     .   nkdim1, nkdim2, mkdim1, mkdim2, lkdim1, lkdim2,
     .   nnodex, nnodey, nnodephi, nxmx, nymx, nphimx,
     .   xprime, yprime, phiprime)

      implicit none

      integer nkx1, nkx2, nky1, nky2, nphi1, nphi2,
     .    nkdim1, nkdim2, mkdim1, mkdim2, lkdim1, lkdim2,
     .    nnodex, nnodey, nnodephi, nxmx, nymx, nphimx,
     .    i, j, k, n, m, nphi, myid

      real xprime(nxmx), yprime(nymx), phiprime(nphimx)
      real eps0, omgrf

      complex xjx(nxmx, nymx, nphimx),
     .        xjy(nxmx, nymx, nphimx),
     .        xjz(nxmx, nymx, nphimx)

      complex fgk(nkdim1 : nkdim2, mkdim1 : mkdim2, lkdim1 : lkdim2),
     .        fak(nkdim1 : nkdim2, mkdim1 : mkdim2, lkdim1 : lkdim2),
     .        fpk(nkdim1 : nkdim2, mkdim1 : mkdim2, lkdim1 : lkdim2)

      complex term1, term2, term3

      complex xlhs, xrhs
      complex u(nkdim1 : nkdim2, mkdim1 : mkdim2, lkdim1 : lkdim2),
     .        v(nkdim1 : nkdim2, mkdim1 : mkdim2, lkdim1 : lkdim2),
     .        w(nkdim1 : nkdim2, mkdim1 : mkdim2, lkdim1 : lkdim2)
      complex zi

      zi = cmplx(0.0, 1.0)
      eps0 = 8.85e-12

      if (myid .eq. 0) then
         write(15, 910)
         write(15, 930)
         write(15, 920)
         write(15, 940)
         write(15, 920)

         write(6, 910)
         write(6, 930)
         write(6, 920)
         write(6, 940)
         write(6, 920)

         term1 = 0.0
         term2 = 0.0
         term3 = 0.0

         do nphi = nphi1, nphi2
            do n = nkx1, nkx2
               do m = nky1, nky2

                  term1 = term1 + fgk(n, m, nphi) * u(n, m, nphi)
                  term2 = term2 + fak(n, m, nphi) * v(n, m, nphi)
                  term3 = term3 + fpk(n, m, nphi) * w(n, m, nphi)

               end do
            end do
         end do

         xlhs = term1 + term2 + term3
         xrhs = -zi / omgrf / eps0 * xjy(i, j, k)
         write(15, 900) i, j, xprime(i), yprime(j), xlhs, xrhs
         write(6,  900) i, j, xprime(i), yprime(j), xlhs, xrhs
      endif

      return
  900 format(i5, i5, 1p9e12.3)
  940 format(1h , 4h   i, 5h    j, 4x, 6h x(i) , 6x, 6h y(j) ,
     1                    6x,  6hre lhs, 6x, 6him lhs, 6x,
     1                         6hre rhs, 6x, 6him rhs)
  910 format (1h1)
  920 format (1h0)
  930 format(3x, 14h v diagnostics)
      end

c
c***************************************************************************
c
      subroutine checkz(myid, i, j, k, omgrf,
     .   nkx1, nkx2, nky1, nky2, nphi1, nphi2,
     .   frk, fqk, fsk, u, v, w, xjx, xjy, xjz,
     .   nkdim1, nkdim2, mkdim1, mkdim2, lkdim1, lkdim2,
     .   nnodex, nnodey, nnodephi, nxmx, nymx, nphimx,
     .   xprime, yprime, phiprime)

      implicit none

      integer nkx1, nkx2, nky1, nky2, nphi1, nphi2,
     .    nkdim1, nkdim2, mkdim1, mkdim2, lkdim1, lkdim2,
     .    nnodex, nnodey, nnodephi, nxmx, nymx, nphimx,
     .    i, j, k, n, m, nphi, myid

      real xprime(nxmx), yprime(nymx), phiprime(nphimx)
      real eps0, omgrf

      complex xjx(nxmx, nymx, nphimx),
     .        xjy(nxmx, nymx, nphimx),
     .        xjz(nxmx, nymx, nphimx)

      complex frk(nkdim1 : nkdim2, mkdim1 : mkdim2, lkdim1 : lkdim2),
     .        fqk(nkdim1 : nkdim2, mkdim1 : mkdim2, lkdim1 : lkdim2),
     .        fsk(nkdim1 : nkdim2, mkdim1 : mkdim2, lkdim1 : lkdim2)

      complex term1, term2, term3

      complex xlhs, xrhs
      complex u(nkdim1 : nkdim2, mkdim1 : mkdim2, lkdim1 : lkdim2),
     .        v(nkdim1 : nkdim2, mkdim1 : mkdim2, lkdim1 : lkdim2),
     .        w(nkdim1 : nkdim2, mkdim1 : mkdim2, lkdim1 : lkdim2)
      complex zi

      zi = cmplx(0.0, 1.0)
      eps0 = 8.85e-12

      if (myid .eq. 0) then
         write(15, 910)
         write(15, 930)
         write(15, 920)
         write(15, 940)
         write(15, 920)

         write(6, 910)
         write(6, 930)
         write(6, 920)
         write(6, 940)
         write(6, 920)

         term1 = 0.0
         term2 = 0.0
         term3 = 0.0

         do nphi = nphi1, nphi2
            do n = nkx1, nkx2
               do m = nky1, nky2

                  term1 = term1 + frk(n, m, nphi) * u(n, m, nphi)
                  term2 = term2 + fqk(n, m, nphi) * v(n, m, nphi)
                  term3 = term3 + fsk(n, m, nphi) * w(n, m, nphi)

               end do
            end do
         end do

         xlhs = term1 + term2 + term3
         xrhs = -zi / omgrf / eps0 * xjz(i, j, k)
         write(15, 900) i, j, xprime(i), yprime(j), xlhs, xrhs
         write(6,  900) i, j, xprime(i), yprime(j), xlhs, xrhs
      endif

      return
  900 format(i5, i5, 1p9e12.3)
  940 format(1h , 4h   i, 5h    j, 4x, 6h x(i) , 6x, 6h y(j) ,
     1                    6x,  6hre lhs, 6x, 6him lhs, 6x,
     1                         6hre rhs, 6x, 6him rhs)
  910 format (1h1)
  920 format (1h0)
  930 format(3x, 14h w diagnostics)
      end

c
c***************************************************************************
c


text:
! @@name:	declare_target.1f
! @@type:	F-free
! @@compilable:	yes
! @@linkable:	yes
! @@expect:	success
module module_fib
contains
   subroutine fib(N)
      integer :: N
      !$omp declare target
      !...
   end subroutine
end module
module params
integer :: THRESHOLD=1000000
end module
program my_fib
use params
use module_fib
   !$omp target if( N > THRESHOLD )
      call fib(N)
   !$omp end target
end program

max_stars_repo_name: kammerdienerb/flang

max_issues_repo_name: kammerdienerb/flang

max_forks_repo_name: kammerdienerb/flang


text:
       subroutine velect(iter,iconv,icorr,ispp,ifcore,
     1 nr,r,rab,zel,cdd,cdu,cdc,vod,vou,etot,y,yp,
     2 ypp,s1,s2,w)
c
c    velect generates the electronic output potential from
c    the electron charge density.  The ionic part is
c    added in dsolv1/dsolv2.
c
c  This version is a  revision of the exchange-correlation part
c  to introduce the interface ATOMXC which allows GGA calculations
c  with Perdew'96 and others xc functionals.
c  J.M. SOLER and L.C. BALBAS, january 1997.
c
c   version 5.62
c
       implicit double precision (a-h,o-z)
c
       character*1 ispp
       character*2 icorr
c
c  njtj  *** modification start  ***
       parameter (zero=0.D0,one=1.D0,pfive=.5D0,opf=1.5D0,pnn=.99D0)
Cray       parameter (zero=0.0,one=1.0,pfive=0.5,opf=1.5,pnn=0.99)
c  njtj  *** modification end  ***
c
      dimension r(nr),rab(nr),cdd(nr),cdu(nr),cdc(nr),
     1 vod(nr),vou(nr),etot(10),y(nr),yp(nr),ypp(nr),
     2 s1(nr),s2(nr),w(3*nr)

      parameter ( maxr = 1500 )
      dimension dens(maxr,2),vxc(maxr,2)

      pi=4*atan(one)
c
c------Machine dependent parameter-
c------Require exp(-2*expzer) to be within the range of the machine
c
c      expzer = 3.7D2
cApollo      expzer = 3.7D2
cSun      expzer = 3.7D2
cVax      expzer = 44.D0
Cray      expzer = 2.8E3
c
c      fit cd/r by splines
c
       y(1) = zero
       do 10 i=2,nr
         y(i) = (cdd(i)+cdu(i))/r(i)
 10    continue
       if (ifcore .eq. 2) then
         do 11 i=2,nr
           y(i) = y(i) + cdc(i)/r(i)
 11      continue
       endif
       isx = 0
       a1 = zero
       an = zero
       b1 = zero
       bn = zero
       nrm=nr
       call splift(r,y,yp,ypp,nrm,w,ierr,isx,a1,b1,an,bn)
       if(ierr.ne.1) then
         write(6,20000)ierr
         call ext(420+ierr)
       endif
20000  format(1x,'****** Error in splift ierr =',i2)
c
c      compute the integrals of cd/r and cd from
c      r(1)=0 to r(i)
c
       xlo = zero
       call spliq(r,y,yp,ypp,nrm,xlo,r,nrm,s2,ierr)
       if(ierr.ne.1) then
         write(6,20001)ierr
         call ext(440+ierr)
       endif
20001  format(1x,'****** Error in spliq ierr =',i2)
       do 20 i=1,nr
         ypp(i) = r(i)*ypp(i) + 2*yp(i)
         yp(i)  = r(i)*yp(i)  + y(i)
         y(i)   = r(i)*y(i)
 20    continue
       call spliq(r,y,yp,ypp,nrm,xlo,r,nrm,s1,ierr)
       if(ierr.ne.1) then
         write(6,20002)ierr
         call ext(460+ierr)
       endif
20002  format(1x,'****** Error in spliq ierr =',i2)
c
c      check normalization
c
       xnorm = zero
       if (ifcore .eq. 2 .and. iter .eq. 0 ) zel=s1(nr)
       if (zel .ne. zero) xnorm = zel/s1(nr)
       if (iter .gt. 3 .and. abs(zel-s1(nr)) .gt. 0.01) then
         if (zel .lt. s1(nr)+1.0 ) then
           write(6,24) iter,xnorm
 24    format(/,' warning *** charge density rescaled in',
     1 ' velect',/,' iteration number',i4,3x,
     2 'scaling factor =',f6.3,/)
         else
           xnorm=pnn*xnorm
           write(6,25) iter,xnorm
 25    format(/,' warning *** charge density partially rescaled in',
     1 ' velect',/,' iteration number',i4,3x,
     2 'scaling factor =',f6.3,/)
         endif
       endif
c
c      compute new hartree potential
c      renormalize the charge density
c
       do 30 i=2,nr
         vod(i) = 2 * xnorm*(s1(i)/r(i) + s2(nr) - s2(i))
         vou(i) = vod(i)
         cdd(i) = xnorm*cdd(i)
         cdu(i) = xnorm*cdu(i)
 30    continue
c
c      compute hartree contribution to total energy
c
       if (iconv .eq. 1) then
         ehart = zero
         ll = 4
         do 40 i=2,nr
           ehart = ehart+ll*(cdd(i)+cdu(i))*vod(i)*rab(i)
           ll = 6 - ll
 40      continue
         ehart = ehart / 6
       endif
c
c     Add exchange and correlation
c     This part is totally new. It is written to use
c     the package XC of J.M. Soler.
c     J.M. Soler and L.C. Balbas. January 1997
c
c     Compute dens(i,nspin) = density up, density down
c
      if (nr.gt.maxr) then
        write(6,*) 'velect: maxr too small. nr, maxr =', nr, maxr
        stop
      endif

      do 90 i=2,nr
        if (ispp .eq. 's') then
	  dens(i,1) = cdu(i)/(4.d0*pi*r(i)**2)
	  dens(i,2) = cdd(i)/(4.d0*pi*r(i)**2)
        else
          dens(i,1) = 0.5d0*(cdu(i) + cdd(i))/(4.d0*pi*r(i)**2)
          dens(i,2) = dens(i,1)
        endif
        if (ifcore .ge. 1) then
              dens(i,1) = dens(i,1) + 0.5d0 * cdc(i)/(4.d0*pi*r(i)**2)
              dens(i,2) = dens(i,2) + 0.5d0 * cdc(i)/(4.d0*pi*r(i)**2)
        endif
  90  continue

c
c     Extrapolate the density at r=0
c
      dens(1,1) = dens(2,1) - (dens(3,1)-dens(2,1))*r(2)/(r(3)-r(2))
      dens(1,2) = dens(2,2) - (dens(3,2)-dens(2,2))*r(2)/(r(3)-r(2))
c
c     Define 'irel' and 'nspin' for the interface ATOMXC
c
      if (ispp .eq. 'r') irel = 1
      if (ispp .ne. 'r') irel = 0
      nspin = 2
c
      r(1) = 0.0d0
           if (icorr .eq. 'ca') then
 	      call atomxc('LDA','ca',irel,nr,maxr,r,nspin,dens,
     .                     ex,ec,dx,dc,vxc)
	    elseif(icorr .eq. 'pw') then
              call atomxc('LDA','pw92',irel,nr,maxr,r,nspin,dens,
     .                    ex,ec,dx,dc,vxc)
            elseif(icorr .eq. 'pb') then
              call atomxc('GGA','pbe',irel,nr,maxr,r,nspin,dens,
     .                     ex,ec,dx,dc,vxc)
           endif

c
c   Add vxc to total potential and energies
c
      do 150 i=2,nr
        vou(i) = vou(i) + vxc(i,1)
        vod(i) = vod(i) + vxc(i,2)
  150 continue

clcb
      xccor=0.0d0
      ll = 4
      do 160 i=2,nr
        xccor = xccor + ll * rab(i)*
     .             (vxc(i,1)*cdu(i) + vxc(i,2)*cdd(i))
        ll = 6 - ll
  160 continue

      etot(4) = ehart
      etot(5) = xccor / 3
      etot(6) = xccor - 4*(ex + ec)
      etot(7) = ex + ec

clcb
c
c     Obtain total potential at r = 0
c
      vod(1) = vod(2) - (vod(3)-vod(2))*r(2)/(r(3)-r(2))
      vou(1) = vou(2) - (vou(3)-vou(2))*r(2)/(r(3)-r(2))

c     *** lcb-jms modification end ********************
      return
      end

max_stars_repo_name: qsnake/petot

max_issues_repo_name: qsnake/petot

max_forks_repo_name: qsnake/petot


text:
      SUBROUTINE DT_EVTOUT(Mode)

C***********************************************************************
C            MODE  = 1  plot content of complete DTEVT1 to out. unit   *
C                    3  plot entries of extended DTEVT1 (DTEVT2)       *
C                    4  plot entries of DTEVT1 and DTEVT2              *
C This version dated 11.12.94 is written by S. Roesler                 *
C***********************************************************************

      IMPLICIT NONE
      INTEGER irange , Mode
      SAVE

      INCLUDE 'inc/dtflka'

C event history

      INCLUDE 'inc/dtevt1'

      DIMENSION irange(NMXHKK)

      IF ( Mode.EQ.2 ) RETURN
      CALL DT_EVTPLO(irange,Mode)

      END SUBROUTINE

max_stars_repo_name: pzhristov/DPMJET

max_issues_repo_name: pzhristov/DPMJET

max_forks_repo_name: pzhristov/DPMJET


text:
c  99-08-20  add info to the output data file
c  99-06-28  increase number of AOs to 2000
c  99-04-23  add minor corrections in wording
C  99-02-23  add cntmos subroutine which counts mos and aos
C  99-02-10  read-in of filenames from standard i/o made formatted
C  98-07-02  skips part of the program if none of the vecs are to be removed
C                  same for aos; fixed vector file format (1P added)
C
C  Program for deleting a vector and/or ao from the $vec containing file
C   V.Kairys
C
C  Limits: 2000 Basis functions, 500 AOs or MOs can be deleted
C
      PROGRAM DELAOMO_RANGE
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      CHARACTER NEWFIL*255,OLDFIL*255,WORD*8
      DIMENSION NPLACE(500),NIFL(2000),EL(5),ARRAY(2000),ARNEW(2000),
     &NIFA(2000),NPLACA(500)
      WRITE(6,*)'Enter name of a file which contains a $VEC group'
      READ(5,1000)OLDFIL
      OPEN(7,FILE=OLDFIL,STATUS='OLD',FORM='FORMATTED')
      WRITE(6,*)'Enter destination file name'
      READ(5,1000)NEWFIL
      OPEN(8,FILE=NEWFIL,STATUS='NEW',FORM='FORMATTED')
c      WRITE(6,*)'How many vectors do you have in the vector file?'
c      READ(5,*)NVEC
c      WRITE(6,*)'What is the total number of basis functions?'
c      READ(5,*)NTOT
      call cntmos(7,nvec,ntot)
C   Calculates the number of lines for one vector
C   and also # of elements in the last line
      IF(MOD(NTOT,5).EQ.0)THEN
        NLIN=NTOT/5
        NELEMI=5
      ELSE
        NLIN=INT(NTOT/5)+1
        NELEMI=MOD(NTOT,5)
      END IF
C
      WRITE(6,*)'How many MOs do you want to delete? Enter 0 if none.'
      READ(5,*)NMOS
      WRITE(6,*)'How many AOs do you want to delete? Enter 0 if none.'
      READ(5,*)NAOS
      NLEFT=NTOT-NAOS
C   Calculates the number of lines in the new vector file
C   and also # of elements in the last line
      IF(MOD(NLEFT,5).EQ.0)THEN
        NLINL=NLEFT/5
        NELEM=5
      ELSE
        NLINL=INT(NLEFT/5)+1
        NELEM=MOD(NLEFT,5)
      END IF
C
C Initialize nplace and nifl
      DO 1 I=1,2000
1     NPLACE(I)=0
      DO 2 I=1,2000
2     NIFL(I)=0
      if(nmos.gt.0)then
      WRITE(6,*)'Write numbers of the MO-s which you'
      WRITE(6,*)' want to delete. Write all, separated by commas'
      WRITE(6,*)' For example, 1,5,15,27 <Enter>'
      READ(5,*)(NPLACE(I),I=1,NMOS)
      end if
C    if nifl=1 then the vector is to be removed
      DO 9 I=1,NMOS
9       NIFL(NPLACE(I))=1
C
C Initialize nplaca and nifa
      DO 21 I=1,500
21     NPLACA(I)=0
      DO 22 I=1,500
22     NIFA(I)=0
      if(naos.gt.0)then
      WRITE(6,*)'Enter the range of the AOs which you'
      WRITE(6,*)' want to delete. First enter the first member'
      WRITE(6,*)' and then the last, separated by a comma.'
      WRITE(6,*)' Example: 94,145 <Enter>'
      READ(5,*)NFIRST,NLAST
      naos1=nlast-nfirst+1
      if(naos1.ne.naos)then
        write(6,*)' Error: your range contains ',naos1,' AOs'
        goto 116
      end if
      do i=1,naos
         nplaca(i)=nfirst+i-1
      end do
c      READ(5,*)(NPLACA(I),I=1,NAOS)
      end if
C    if nifa=1 then the ao is to be removed
      DO 29 I=1,NAOS
29       NIFA(NPLACA(I))=1
C
3     FORMAT(I2,I3,1P,5E15.8)
C
C
   20  READ(7,9100,END=115)WORD
 9100  FORMAT(A8)
       IF(WORD.NE.' $VEC')GOTO 20
      WRITE(8,9200)OLDFIL,NAOS,NMOS
9200  FORMAT('THIS FILE IS CREATED FROM ',A15,' BY REMOVING ',I3,
     *       ' AO-S AND ',I3,' MO-S')
      WRITE(8,FMT='(5H $VEC)')
      IVC=0
      DO 8 I=1,NVEC
        IF(NIFL(I).EQ.1)THEN
         DO 11 J=1,NLIN
C  Read but no write!
11          READ(7,3)IN,IA,EL
        ELSE
         IVC=IVC+1
         DO 50 J=1,NLIN-1
50          READ(7,3)IN,IA,(ARRAY((J-1)*5+K),K=1,5)
c just one last line
         READ(7,3)IN,IA,(ARRAY((NLIN-1)*5+K),K=1,NELEMI)
c
c
c
c    redoing array into arnew
         IAR=0
         DO 51 K=1,NTOT
         IF(NIFA(K).EQ.1)GOTO 52
         IAR=IAR+1
         ARNEW(IAR)=ARRAY(K)
 52      CONTINUE
 51      CONTINUE
         if(iar.ne.nleft)then
           write(6,*)'ERROR: iar.neq.nleft'
         end if
c  writing arnew into a new vec file
         DO 53 IA=1,NLINL-1
53      WRITE(8,3)MOD(IVC,100),IA,(ARNEW((IA-1)*5+II),II=1,5)
c just one line more
        WRITE(8,3)MOD(IVC,100),NLINL,(ARNEW((NLINL-1)*5+II),II=1,NELEM)
        END IF
8     CONTINUE
      WRITE(8,FMT='(5H $END)')
      WRITE(6,*)'Finished...'
      GOTO 116
115   WRITE(6,*)' $VEC GROUP NOT FOUND'
      GOTO 116
116   CLOSE(7)
      CLOSE(8)
      STOP
1000  FORMAT(A255)
      END
C
      SUBROUTINE CNTMOS(IUNIT,NMOS,NAOS)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      CHARACTER CSTRNG*80,WORD*8
      DIMENSION DMY(5)
c
c
   20 READ(IUNIT,9100,END=115)WORD
      IF(WORD.NE.' $VEC')GOTO 20
c
c  scan records
c
      nmos=0
      numlv=0
      line=0
      iaprev=1
      nldone=0
      nfdone=0
      nleft=5
   30 READ (IUNIT,9068,END = 1400) CSTRNG
      line = line + 1
c  since len_trim isn't a standard fortran77 function, we have to do without it
      do  40 i=80,1,-1
          IF (CSTRNG(i:i) .NE. ' ') go to 50
 40   continue
 50   lentrm=i
c
ccc      write(6,*)'line',line,'lentrm',lentrm
      if(nfdone.eq.0 .and. lentrm.gt.10 .and. lentrm.lt.80) then
          nleft=(lentrm-5)/15
          nfdone=1
      end if
c
      if(lentrm.lt.10)then
c  end of vectors has been reached
        linem1=line-1
        if( mod(linem1,numlv) .ne. 0) then
           write(6,*)' ERROR: Vec File is irregular:'
           write(6,*)' the total number of lines cannot be divided by '
           write(6,*)' the number of lines per vector without remainder'
           go to 116
        end if
        nmos=linem1/numlv
        goto 114
      end if
c
      if(nldone.gt.0)goto 30
c   check out the first two integers in cstrng
      READ (UNIT=CSTRNG,FMT='(I2,I3,5E15.8)') IA,IB,DMY(1)
ccc     write(6,*)'iaprev',iaprev,'ia',ia
      if(iaprev.ne.ia)then
          numlv=line-1
          nldone=1
          naos=(numlv-1)*5+nleft
      end if
      iaprev=ia
      goto 30
cc
114   write(6,*)'THE PROGRAM FOUND ',NMOS,' MOS, ',naos,' AOS'
      GOTO 116
115   WRITE(6,*)' $VEC GROUP NOT FOUND'
      GOTO 116
1400  WRITE(6,*)' EOF FILE ENCOUNTERED'
 116  REWIND(IUNIT)
      RETURN
9068  FORMAT(A80)
9100  FORMAT(A8)
      END