diff --git a/Makefile b/Makefile index e57babf..a300718 100755 --- a/Makefile +++ b/Makefile @@ -1,3 +1,9 @@ +OSNF_DIR = osnf + +.PHONY: osnf cleanall +CLEANDIRS = $(OSNF_DIR) ./ + + DEBUG = -fbounds-check -g OPT =-O3 @@ -21,57 +27,48 @@ RANLIB = ranlib OBJ = o FFLAGS = $(OPT) $(DEBUG) -w -o FFLAGS2 = $(DEBUG) -w -O3 -o +VAR_TYPE = 1 # 0 single, 1 double main.exe : model_lib.a main.$(OBJ) variables.$(OBJ) initialisation.$(OBJ) \ mpi_module.$(OBJ) driver_code.$(OBJ) advection.$(OBJ) $(FOR2) $(FFLAGS2)main.exe main.$(OBJ) variables.$(OBJ) initialisation.$(OBJ) \ mpi_module.$(OBJ) driver_code.$(OBJ) advection.$(OBJ) -lm model_lib.a \ ${NETCDFLIB} -I ${NETCDFMOD} ${NETCDF_LIB} $(DEBUG) -model_lib.a : nrtype.$(OBJ) nr.$(OBJ) nrutil.$(OBJ) locate.$(OBJ) polint.$(OBJ) \ - rkqs.$(OBJ) rkck.$(OBJ) odeint.$(OBJ) zbrent.$(OBJ) \ - hygfx.$(OBJ) random.$(OBJ) - $(AR) rc model_lib.a nrutil.$(OBJ) locate.$(OBJ) polint.$(OBJ) \ - rkqs.$(OBJ) rkck.$(OBJ) odeint.$(OBJ) zbrent.$(OBJ) \ - hygfx.$(OBJ) random.$(OBJ) -locate.$(OBJ) : locate.f90 - $(FOR) locate.f90 $(FFLAGS)locate.$(OBJ) -polint.$(OBJ) : polint.f90 - $(FOR) polint.f90 $(FFLAGS)polint.$(OBJ) -nrtype.$(OBJ) : nrtype.f90 - $(FOR) nrtype.f90 $(FFLAGS)nrtype.$(OBJ) -nr.$(OBJ) : nr.f90 - $(FOR) nr.f90 $(FFLAGS)nr.$(OBJ) -nrutil.$(OBJ) : nrutil.f90 - $(FOR) nrutil.f90 $(FFLAGS)nrutil.$(OBJ) -rkqs.$(OBJ) : rkqs.f90 - $(FOR) rkqs.f90 $(FFLAGS)rkqs.$(OBJ) -rkck.$(OBJ) : rkck.f90 - $(FOR) rkck.f90 $(FFLAGS)rkck.$(OBJ) -odeint.$(OBJ) : odeint.f90 - $(FOR) odeint.f90 $(FFLAGS)odeint.$(OBJ) -zbrent.$(OBJ) : zbrent.f90 - $(FOR) zbrent.f90 $(FFLAGS2)zbrent.$(OBJ) -hygfx.$(OBJ) : hygfx.for - $(FOR) hygfx.for $(FFLAGS)hygfx.$(OBJ) -random.$(OBJ) : random.f90 - $(FOR) random.f90 $(FFLAGS)random.$(OBJ) -variables.$(OBJ) : variables.f90 - $(FOR) variables.f90 $(FFLAGS)variables.$(OBJ) -mpi_module.$(OBJ) : mpi_module.f90 - $(FOR) mpi_module.f90 $(FFLAGS)mpi_module.$(OBJ) -initialisation.$(OBJ) : initialisation.f90 random.$(OBJ) mpi_module.$(OBJ) +model_lib.a : osnf_code + $(AR) rc model_lib.a \ + $(OSNF_DIR)/numerics.$(OBJ) $(OSNF_DIR)/zeroin.$(OBJ) $(OSNF_DIR)/sfmin.$(OBJ) \ + $(OSNF_DIR)/fmin.$(OBJ) $(OSNF_DIR)/r1mach.$(OBJ) \ + $(OSNF_DIR)/d1mach.$(OBJ) $(OSNF_DIR)/dfsid1.$(OBJ) \ + $(OSNF_DIR)/poly_int.$(OBJ) $(OSNF_DIR)/find_pos.$(OBJ) \ + $(OSNF_DIR)/svode.$(OBJ) \ + $(OSNF_DIR)/slinpk.$(OBJ) $(OSNF_DIR)/vode.$(OBJ) \ + $(OSNF_DIR)/dlinpk.$(OBJ) $(OSNF_DIR)/vode_integrate.$(OBJ) \ + $(OSNF_DIR)/erfinv.$(OBJ) $(OSNF_DIR)/tridiagonal.$(OBJ) \ + $(OSNF_DIR)/hygfx.$(OBJ) $(OSNF_DIR)/random.$(OBJ) +variables.$(OBJ) : variables.f90 osnf_code + $(FOR) variables.f90 -I$(OSNF_DIR) $(FFLAGS)variables.$(OBJ) +mpi_module.$(OBJ) : mpi_module.f90 osnf_code + $(FOR) mpi_module.f90 -cpp -DVAR_TYPE=$(VAR_TYPE) -I$(OSNF_DIR) $(FFLAGS)mpi_module.$(OBJ) +initialisation.$(OBJ) : initialisation.f90 mpi_module.$(OBJ) osnf_code $(FOR) initialisation.f90 -I ${NETCDFMOD} \ - $(FFLAGS)initialisation.$(OBJ) -driver_code.$(OBJ) : driver_code.f90 advection.$(OBJ) mpi_module.$(OBJ) - $(FOR) driver_code.f90 -I ${NETCDFMOD} \ - $(FFLAGS)driver_code.$(OBJ) -advection.$(OBJ) : advection.f90 - $(FOR) advection.f90 $(FFLAGS)advection.$(OBJ) + -I$(OSNF_DIR) $(FFLAGS)initialisation.$(OBJ) +driver_code.$(OBJ) : driver_code.f90 advection.$(OBJ) mpi_module.$(OBJ) osnf_code + $(FOR) driver_code.f90 -I$(OSNF_DIR) -I ${NETCDFMOD} $(FFLAGS)driver_code.$(OBJ) +advection.$(OBJ) : advection.f90 osnf_code + $(FOR) advection.f90 -I$(OSNF_DIR) $(FFLAGS)advection.$(OBJ) main.$(OBJ) : main.f90 variables.$(OBJ) mpi_module.$(OBJ) initialisation.$(OBJ) \ driver_code.$(OBJ) advection.$(OBJ) $(FOR) main.f90 -I ${NETCDFMOD} $(FFLAGS)main.$(OBJ) -clean : +osnf_code: + $(MAKE) -C $(OSNF_DIR) + +clean: rm *.exe *.o *.mod *~ \ model_lib.a +cleanall: + for i in $(CLEANDIRS); do \ + $(MAKE) -C $$i clean; \ + done + diff --git a/advection.f90 b/advection.f90 index 74204d9..7594ae5 100644 --- a/advection.f90 +++ b/advection.f90 @@ -3,7 +3,7 @@ !>@brief !>advection routines module advection - use nrtype + use numerics_type private public :: lax_wendroff_sphere, lax_wendroff_ll, dissipation, smagorinsky contains @@ -26,25 +26,25 @@ module advection !>\f$ \frac{\partial \psi}{\partial t} + \frac{\partial u \psi}{\partial x} = 0 \f$ subroutine lax_wendroff_sphere(ip,jp,o_halo,dt,dx,dy,g,u,v,h,hs,re,theta,f_cor) - use nrtype + use numerics_type implicit none integer(i4b), intent(in) :: ip,jp,o_halo - real(sp), intent(in) :: dt, g, re - real(sp), intent(in), dimension(1-o_halo:ip+o_halo,1-o_halo:jp+o_halo) :: & + real(wp), intent(in) :: dt, g, re + real(wp), intent(in), dimension(1-o_halo:ip+o_halo,1-o_halo:jp+o_halo) :: & hs, f_cor, dx, dy - real(sp), dimension(1-o_halo:jp+o_halo), intent(in) :: theta - real(sp), intent(inout), dimension(1-o_halo:ip+o_halo,1-o_halo:jp+o_halo) :: & + real(wp), dimension(1-o_halo:jp+o_halo), intent(in) :: theta + real(wp), intent(inout), dimension(1-o_halo:ip+o_halo,1-o_halo:jp+o_halo) :: & h, u, v ! local variables: - real(sp), dimension(1-o_halo:ip+o_halo,1-o_halo:jp+o_halo) :: & + real(wp), dimension(1-o_halo:ip+o_halo,1-o_halo:jp+o_halo) :: & dy1, v1, h1, vh, uh, vh1, Ux, Uy, Vx, Vy, Vy2 - real(sp), dimension(1:ip,1:jp) :: & + real(wp), dimension(1:ip,1:jp) :: & uh_new, vh_new, h_new - real(sp), dimension(0:ip,1:jp) :: h_mid_xt, uh_mid_xt, Ux_mid_xt, Vx_mid_xt - real(sp), dimension(0:ip,1:jp) :: vh_mid_xt - real(sp), dimension(1:ip,0:jp) :: h_mid_yt, uh_mid_yt, vh_mid_yt, c_mid_yt, & + real(wp), dimension(0:ip,1:jp) :: h_mid_xt, uh_mid_xt, Ux_mid_xt, Vx_mid_xt + real(wp), dimension(0:ip,1:jp) :: vh_mid_xt + real(wp), dimension(1:ip,0:jp) :: h_mid_yt, uh_mid_yt, vh_mid_yt, c_mid_yt, & Uy_mid_yt, Vy_mid_yt, Vy_mid_yt2 integer(i4b) :: j, i @@ -55,7 +55,7 @@ subroutine lax_wendroff_sphere(ip,jp,o_halo,dt,dx,dy,g,u,v,h,hs,re,theta,f_cor) h1(:,j)=h(:,j)*cos(theta(j)) enddo do j=0,jp - c_mid_yt(:,j)=cos(0.5_sp*(theta(j+1)+theta(j))) + c_mid_yt(:,j)=cos(0.5_wp*(theta(j+1)+theta(j))) enddo uh=u *h @@ -63,27 +63,27 @@ subroutine lax_wendroff_sphere(ip,jp,o_halo,dt,dx,dy,g,u,v,h,hs,re,theta,f_cor) vh1=v1*h ! continuity equation (calculate mid-point values at 0.5*dt): - h_mid_xt = 0.5_sp*(h(1:ip+1,1:jp)+h(0:ip,1:jp)) & - -(0.5_sp*dt/(0.5_sp*(dx(1:ip+1,1:jp)+dx(0:ip,1:jp)))) & + h_mid_xt = 0.5_wp*(h(1:ip+1,1:jp)+h(0:ip,1:jp)) & + -(0.5_wp*dt/(0.5_wp*(dx(1:ip+1,1:jp)+dx(0:ip,1:jp)))) & *(uh(1:ip+1,1:jp)-uh(0:ip,1:jp)) - h_mid_yt = 0.5_sp*(h(1:ip,1:jp+1)+h(1:ip,0:jp)) & - -(0.5_sp*dt/(0.5_sp*(dy1(1:ip,1:jp+1)+dy1(1:ip,0:jp)))) & + h_mid_yt = 0.5_wp*(h(1:ip,1:jp+1)+h(1:ip,0:jp)) & + -(0.5_wp*dt/(0.5_wp*(dy1(1:ip,1:jp+1)+dy1(1:ip,0:jp)))) & *(vh1(1:ip,1:jp+1)-vh1(1:ip,0:jp)) ! v-phi, or u momentum equation (calculate mid-point values at 0.5*dt): - Ux = uh*u + g*h**2*(0.5_sp) + Ux = uh*u + g*h**2*(0.5_wp) Uy = uh*v1 - uh_mid_xt(0:ip,:) = 0.5_sp*(uh(1:ip+1,1:jp)+uh(0:ip,1:jp)) & - -(0.5_sp*dt/(0.5_sp*(dx(1:ip+1,1:jp)+dx(0:ip,1:jp))))* & + uh_mid_xt(0:ip,:) = 0.5_wp*(uh(1:ip+1,1:jp)+uh(0:ip,1:jp)) & + -(0.5_wp*dt/(0.5_wp*(dx(1:ip+1,1:jp)+dx(0:ip,1:jp))))* & (Ux(1:ip+1,1:jp)-Ux(0:ip,1:jp)) & - +0.125_sp*dt*(f_cor(1:ip+1,1:jp)+f_cor(0:ip,1:jp))* & + +0.125_wp*dt*(f_cor(1:ip+1,1:jp)+f_cor(0:ip,1:jp))* & (vh(1:ip+1,1:jp)+vh(0:ip,1:jp)) - uh_mid_yt = 0.5_sp*(uh(1:ip,1:jp+1)+uh(1:ip,0:jp)) & - -(0.5_sp*dt/(0.5_sp*(dy1(1:ip,1:jp+1)+dy1(1:ip,0:jp))))* & + uh_mid_yt = 0.5_wp*(uh(1:ip,1:jp+1)+uh(1:ip,0:jp)) & + -(0.5_wp*dt/(0.5_wp*(dy1(1:ip,1:jp+1)+dy1(1:ip,0:jp))))* & (Uy(1:ip,1:jp+1)-Uy(1:ip,0:jp)) & - +0.125_sp*dt*(f_cor(1:ip,1:jp+1)+f_cor(1:ip,0:jp))* & + +0.125_wp*dt*(f_cor(1:ip,1:jp+1)+f_cor(1:ip,0:jp))* & (vh(1:ip,1:jp+1)+vh(1:ip,0:jp)) @@ -91,15 +91,15 @@ subroutine lax_wendroff_sphere(ip,jp,o_halo,dt,dx,dy,g,u,v,h,hs,re,theta,f_cor) ! v-theta, or v momentum equation (calculate mid-point values at 0.5*dt): Vx = uh*v Vy = vh1*v - Vy2 = 0.5_sp*g*h**2 - vh_mid_xt(0:ip,1:jp) = 0.5_sp*(vh(1:ip+1,1:jp)+vh(0:ip,1:jp)) & - -(0.5_sp*dt/(0.5_sp*(dx(1:ip+1,1:jp)+dx(0:ip,1:jp))))*(Vx(1:ip+1,1:jp)-Vx(0:ip,1:jp)) & - -0.125_sp*dt*(f_cor(1:ip+1,1:jp)+f_cor(0:ip,1:jp))*(uh(1:ip+1,1:jp)+uh(0:ip,1:jp)) + Vy2 = 0.5_wp*g*h**2 + vh_mid_xt(0:ip,1:jp) = 0.5_wp*(vh(1:ip+1,1:jp)+vh(0:ip,1:jp)) & + -(0.5_wp*dt/(0.5_wp*(dx(1:ip+1,1:jp)+dx(0:ip,1:jp))))*(Vx(1:ip+1,1:jp)-Vx(0:ip,1:jp)) & + -0.125_wp*dt*(f_cor(1:ip+1,1:jp)+f_cor(0:ip,1:jp))*(uh(1:ip+1,1:jp)+uh(0:ip,1:jp)) - vh_mid_yt(1:ip,0:jp) = 0.5_sp*(vh(1:ip,1:jp+1)+vh(1:ip,0:jp)) & - -(0.5_sp*dt/(0.5_sp*(dy1(1:ip,1:jp+1)+dy1(1:ip,0:jp))))*(Vy(1:ip,1:jp+1)-Vy(1:ip,0:jp)) & - -(0.5_sp*dt/(dy(1:ip,0:jp)))*(Vy2(1:ip,1:jp+1)-Vy2(1:ip,0:jp)) & - -0.125_sp*dt*(f_cor(1:ip,1:jp+1)+f_cor(1:ip,0:jp))*(uh(1:ip,1:jp+1)+uh(1:ip,0:jp)) + vh_mid_yt(1:ip,0:jp) = 0.5_wp*(vh(1:ip,1:jp+1)+vh(1:ip,0:jp)) & + -(0.5_wp*dt/(0.5_wp*(dy1(1:ip,1:jp+1)+dy1(1:ip,0:jp))))*(Vy(1:ip,1:jp+1)-Vy(1:ip,0:jp)) & + -(0.5_wp*dt/(dy(1:ip,0:jp)))*(Vy2(1:ip,1:jp+1)-Vy2(1:ip,0:jp)) & + -0.125_wp*dt*(f_cor(1:ip,1:jp+1)+f_cor(1:ip,0:jp))*(uh(1:ip,1:jp+1)+uh(1:ip,0:jp)) ! calculate mid-point value of cos (theta) ! c_mid_yt=cos(0.5.*(THETA(:,2:end)+THETA(:,1:end-1))); @@ -108,42 +108,42 @@ subroutine lax_wendroff_sphere(ip,jp,o_halo,dt,dx,dy,g,u,v,h,hs,re,theta,f_cor) ! Now use the mid-point values to predict the values at the next timestep ! continuity: h_new = h(1:ip,1:jp) & - - (dt/(0.5_sp*(dx(1:ip,1:jp)+dx(0:ip-1,1:jp))))*(uh_mid_xt(1:ip,1:jp)-uh_mid_xt(0:ip-1,1:jp)) & - - (dt/(0.5_sp*(dy1(1:ip,1:jp)+dy1(1:ip,0:jp-1)))) * & + - (dt/(0.5_wp*(dx(1:ip,1:jp)+dx(0:ip-1,1:jp))))*(uh_mid_xt(1:ip,1:jp)-uh_mid_xt(0:ip-1,1:jp)) & + - (dt/(0.5_wp*(dy1(1:ip,1:jp)+dy1(1:ip,0:jp-1)))) * & (vh_mid_yt(1:ip,1:jp)*c_mid_yt(1:ip,1:jp)-vh_mid_yt(1:ip,0:jp-1)*c_mid_yt(1:ip,0:jp-1)) ! u-momentum equation: - Ux_mid_xt = uh_mid_xt*uh_mid_xt/h_mid_xt + 0.5_sp*g*h_mid_xt**2 + Ux_mid_xt = uh_mid_xt*uh_mid_xt/h_mid_xt + 0.5_wp*g*h_mid_xt**2 Uy_mid_yt = uh_mid_yt*vh_mid_yt/h_mid_yt*c_mid_yt uh_new = uh(1:ip,1:jp) & - - (dt/(0.5_sp*(dx(1:ip,1:jp)+dx(0:ip-1,1:jp))))* (Ux_mid_xt(1:ip,1:jp)-Ux_mid_xt(0:ip-1,1:jp)) & - - (dt/(0.5_sp*(dy1(1:ip,1:jp)+dy1(1:ip,0:jp-1))))*(Uy_mid_yt(1:ip,1:jp)-Uy_mid_yt(1:ip,0:jp-1)) + - (dt/(0.5_wp*(dx(1:ip,1:jp)+dx(0:ip-1,1:jp))))* (Ux_mid_xt(1:ip,1:jp)-Ux_mid_xt(0:ip-1,1:jp)) & + - (dt/(0.5_wp*(dy1(1:ip,1:jp)+dy1(1:ip,0:jp-1))))*(Uy_mid_yt(1:ip,1:jp)-Uy_mid_yt(1:ip,0:jp-1)) ! v-momentum equation: Vx_mid_xt = uh_mid_xt*vh_mid_xt/h_mid_xt Vy_mid_yt = vh_mid_yt*vh_mid_yt/h_mid_yt*c_mid_yt - Vy_mid_yt2 = 0.5_sp*g*h_mid_yt**2 + Vy_mid_yt2 = 0.5_wp*g*h_mid_yt**2 vh_new = vh(1:ip,1:jp) & - - (dt/(0.5_sp*(dx(1:ip,1:jp)+dx(0:ip-1,1:jp))))*(Vx_mid_xt(1:ip,1:jp)-Vx_mid_xt(0:ip-1,1:jp)) & - - (dt/(0.5_sp*(dy1(1:ip,1:jp)+dy1(1:ip,0:jp-1))))*(Vy_mid_yt(1:ip,1:jp)-Vy_mid_yt(1:ip,0:jp-1)) & + - (dt/(0.5_wp*(dx(1:ip,1:jp)+dx(0:ip-1,1:jp))))*(Vx_mid_xt(1:ip,1:jp)-Vx_mid_xt(0:ip-1,1:jp)) & + - (dt/(0.5_wp*(dy1(1:ip,1:jp)+dy1(1:ip,0:jp-1))))*(Vy_mid_yt(1:ip,1:jp)-Vy_mid_yt(1:ip,0:jp-1)) & - (dt/(dy(1:ip,0:jp-1) ))* & (Vy_mid_yt2(1:ip,1:jp)-Vy_mid_yt2(1:ip,0:jp-1)) ! add on Coriolis and contribution of orography to pressure gradient: - uh_new=uh_new +dt*.5_sp*(f_cor(1:ip,1:jp)*v(1:ip,1:jp) - & - g*(hs(2:ip+1,1:jp)-hs(0:ip-1,1:jp))/(2._sp*dx(1:ip,1:jp)))* & + uh_new=uh_new +dt*.5_wp*(f_cor(1:ip,1:jp)*v(1:ip,1:jp) - & + g*(hs(2:ip+1,1:jp)-hs(0:ip-1,1:jp))/(2._wp*dx(1:ip,1:jp)))* & (h(1:ip,1:jp)+h_new) - vh_new=vh_new -dt*.5_sp*(f_cor(1:ip,1:jp)*u(1:ip,1:jp) + & - g*(hs(1:ip,2:jp+1)-hs(1:ip,0:jp-1))/(2._sp*dy1(1:ip,1:jp)))* & + vh_new=vh_new -dt*.5_wp*(f_cor(1:ip,1:jp)*u(1:ip,1:jp) + & + g*(hs(1:ip,2:jp+1)-hs(1:ip,0:jp-1))/(2._wp*dy1(1:ip,1:jp)))* & (h(1:ip,1:jp)+h_new) -! uh_new=uh_new +dt*.5_sp*(f_cor(1:ip,1:jp)*(vh_mid_xt(1:ip,1:jp)+vh_mid_xt(0:ip-1,1:jp)) ) +! uh_new=uh_new +dt*.5_wp*(f_cor(1:ip,1:jp)*(vh_mid_xt(1:ip,1:jp)+vh_mid_xt(0:ip-1,1:jp)) ) ! -! vh_new=vh_new -dt*.5_sp*(f_cor(1:ip,1:jp)*(uh_mid_xt(1:ip,1:jp)+uh_mid_xt(0:ip-1,1:jp)) ) +! vh_new=vh_new -dt*.5_wp*(f_cor(1:ip,1:jp)*(uh_mid_xt(1:ip,1:jp)+uh_mid_xt(0:ip-1,1:jp)) ) ! re-calculate u and v. @@ -196,30 +196,30 @@ subroutine lax_wendroff_ll(ip,jp,o_halo,dt,g,u,v,h,hs,re,& theta,thetan,dtheta,dthetan, phi, phin, dphi, dphin, f_cor, & recqdq, recqdp, recqdp_s, recqdq_s, redq_s, redq, cq, cq_s) - use nrtype + use numerics_type implicit none integer(i4b), intent(in) :: ip,jp,o_halo - real(sp), intent(in) :: dt, g, re - real(sp), intent(in), dimension(1-o_halo:ip+o_halo,1-o_halo:jp+o_halo) :: & + real(wp), intent(in) :: dt, g, re + real(wp), intent(in), dimension(1-o_halo:ip+o_halo,1-o_halo:jp+o_halo) :: & hs, f_cor, & recqdq, recqdp, recqdp_s, recqdq_s, redq_s, redq, & cq, cq_s - real(sp), dimension(1-o_halo:jp+o_halo), intent(in) :: & + real(wp), dimension(1-o_halo:jp+o_halo), intent(in) :: & theta,thetan, dtheta, dthetan - real(sp), dimension(1-o_halo:ip+o_halo), intent(in) :: & + real(wp), dimension(1-o_halo:ip+o_halo), intent(in) :: & phi, phin, dphi, dphin - real(sp), intent(inout), dimension(1-o_halo:ip+o_halo,1-o_halo:jp+o_halo) :: & + real(wp), intent(inout), dimension(1-o_halo:ip+o_halo,1-o_halo:jp+o_halo) :: & h, u, v ! local variables: - real(sp), dimension(1-o_halo:ip+o_halo,1-o_halo:jp+o_halo) :: & + real(wp), dimension(1-o_halo:ip+o_halo,1-o_halo:jp+o_halo) :: & dy1, v1, h1, vh, uh, vh1, Ux, Uy, Vx, Vy, Vy2 - real(sp), dimension(1:ip,1:jp) :: & + real(wp), dimension(1:ip,1:jp) :: & uh_new, vh_new, h_new - real(sp), dimension(0:ip,1:jp) :: h_mid_xt, uh_mid_xt, Ux_mid_xt, Vx_mid_xt - real(sp), dimension(0:ip,1:jp) :: vh_mid_xt - real(sp), dimension(1:ip,0:jp) :: h_mid_yt, uh_mid_yt, vh_mid_yt, & + real(wp), dimension(0:ip,1:jp) :: h_mid_xt, uh_mid_xt, Ux_mid_xt, Vx_mid_xt + real(wp), dimension(0:ip,1:jp) :: vh_mid_xt + real(wp), dimension(1:ip,0:jp) :: h_mid_yt, uh_mid_yt, vh_mid_yt, & Uy_mid_yt, Vy_mid_yt, Vy_mid_yt2 integer(i4b) :: j, i @@ -232,27 +232,27 @@ subroutine lax_wendroff_ll(ip,jp,o_halo,dt,g,u,v,h,hs,re,& vh1=v1*h ! continuity equation (calculate mid-point values at 0.5*dt): - h_mid_xt = 0.5_sp*(h(1:ip+1,1:jp)+h(0:ip,1:jp)) & - -(0.5_sp*dt/(recqdp_s(0:ip,1:jp))) & + h_mid_xt = 0.5_wp*(h(1:ip+1,1:jp)+h(0:ip,1:jp)) & + -(0.5_wp*dt/(recqdp_s(0:ip,1:jp))) & *(uh(1:ip+1,1:jp)-uh(0:ip,1:jp)) - h_mid_yt = 0.5_sp*(h(1:ip,1:jp+1)+h(1:ip,0:jp)) & - -(0.5_sp*dt/(recqdq_s(1:ip,0:jp))) & + h_mid_yt = 0.5_wp*(h(1:ip,1:jp+1)+h(1:ip,0:jp)) & + -(0.5_wp*dt/(recqdq_s(1:ip,0:jp))) & *(vh1(1:ip,1:jp+1)-vh1(1:ip,0:jp)) ! v-phi, or u momentum equation (calculate mid-point values at 0.5*dt): - Ux = uh*u + g*h**2*(0.5_sp) + Ux = uh*u + g*h**2*(0.5_wp) Uy = uh*v1 - uh_mid_xt(0:ip,:) = 0.5_sp*(uh(1:ip+1,1:jp)+uh(0:ip,1:jp)) & - -(0.5_sp*dt/(recqdp_s(0:ip,1:jp)))* & + uh_mid_xt(0:ip,:) = 0.5_wp*(uh(1:ip+1,1:jp)+uh(0:ip,1:jp)) & + -(0.5_wp*dt/(recqdp_s(0:ip,1:jp)))* & (Ux(1:ip+1,1:jp)-Ux(0:ip,1:jp)) & - +0.125_sp*dt*(f_cor(1:ip+1,1:jp)+f_cor(0:ip,1:jp))* & + +0.125_wp*dt*(f_cor(1:ip+1,1:jp)+f_cor(0:ip,1:jp))* & (vh(1:ip+1,1:jp)+vh(0:ip,1:jp)) - uh_mid_yt = 0.5_sp*(uh(1:ip,1:jp+1)+uh(1:ip,0:jp)) & - -(0.5_sp*dt/(recqdq_s(1:ip,0:jp)))* & + uh_mid_yt = 0.5_wp*(uh(1:ip,1:jp+1)+uh(1:ip,0:jp)) & + -(0.5_wp*dt/(recqdq_s(1:ip,0:jp)))* & (Uy(1:ip,1:jp+1)-Uy(1:ip,0:jp)) & - +0.125_sp*dt*(f_cor(1:ip,1:jp+1)+f_cor(1:ip,0:jp))* & + +0.125_wp*dt*(f_cor(1:ip,1:jp+1)+f_cor(1:ip,0:jp))* & (vh(1:ip,1:jp+1)+vh(1:ip,0:jp)) @@ -260,15 +260,15 @@ subroutine lax_wendroff_ll(ip,jp,o_halo,dt,g,u,v,h,hs,re,& ! v-theta, or v momentum equation (calculate mid-point values at 0.5*dt): Vx = uh*v Vy = vh1*v - Vy2 = 0.5_sp*g*h**2 - vh_mid_xt(0:ip,1:jp) = 0.5_sp*(vh(1:ip+1,1:jp)+vh(0:ip,1:jp)) & - -(0.5_sp*dt/(recqdp_s(0:ip,1:jp)))*(Vx(1:ip+1,1:jp)-Vx(0:ip,1:jp)) & - -0.125_sp*dt*(f_cor(1:ip+1,1:jp)+f_cor(0:ip,1:jp))*(uh(1:ip+1,1:jp)+uh(0:ip,1:jp)) + Vy2 = 0.5_wp*g*h**2 + vh_mid_xt(0:ip,1:jp) = 0.5_wp*(vh(1:ip+1,1:jp)+vh(0:ip,1:jp)) & + -(0.5_wp*dt/(recqdp_s(0:ip,1:jp)))*(Vx(1:ip+1,1:jp)-Vx(0:ip,1:jp)) & + -0.125_wp*dt*(f_cor(1:ip+1,1:jp)+f_cor(0:ip,1:jp))*(uh(1:ip+1,1:jp)+uh(0:ip,1:jp)) - vh_mid_yt(1:ip,0:jp) = 0.5_sp*(vh(1:ip,1:jp+1)+vh(1:ip,0:jp)) & - -(0.5_sp*dt/(recqdq_s(1:ip,0:jp)))*(Vy(1:ip,1:jp+1)-Vy(1:ip,0:jp)) & - -(0.5_sp*dt/(redq_s(1:ip,0:jp)))*(Vy2(1:ip,1:jp+1)-Vy2(1:ip,0:jp)) & - -0.125_sp*dt*(f_cor(1:ip,1:jp+1)+f_cor(1:ip,0:jp))*(uh(1:ip,1:jp+1)+uh(1:ip,0:jp)) + vh_mid_yt(1:ip,0:jp) = 0.5_wp*(vh(1:ip,1:jp+1)+vh(1:ip,0:jp)) & + -(0.5_wp*dt/(recqdq_s(1:ip,0:jp)))*(Vy(1:ip,1:jp+1)-Vy(1:ip,0:jp)) & + -(0.5_wp*dt/(redq_s(1:ip,0:jp)))*(Vy2(1:ip,1:jp+1)-Vy2(1:ip,0:jp)) & + -0.125_wp*dt*(f_cor(1:ip,1:jp+1)+f_cor(1:ip,0:jp))*(uh(1:ip,1:jp+1)+uh(1:ip,0:jp)) ! calculate mid-point value of cos (theta) ! c_mid_yt=cos(0.5.*(THETA(:,2:end)+THETA(:,1:end-1))); @@ -283,7 +283,7 @@ subroutine lax_wendroff_ll(ip,jp,o_halo,dt,g,u,v,h,hs,re,& ! u-momentum equation: - Ux_mid_xt = uh_mid_xt*uh_mid_xt/h_mid_xt + 0.5_sp*g*h_mid_xt**2 + Ux_mid_xt = uh_mid_xt*uh_mid_xt/h_mid_xt + 0.5_wp*g*h_mid_xt**2 Uy_mid_yt = uh_mid_yt*vh_mid_yt/h_mid_yt*cq_s(1:ip,0:jp) uh_new = uh(1:ip,1:jp) & - (dt/(recqdp(0:ip-1,1:jp)))* (Ux_mid_xt(1:ip,1:jp)-Ux_mid_xt(0:ip-1,1:jp)) & @@ -293,7 +293,7 @@ subroutine lax_wendroff_ll(ip,jp,o_halo,dt,g,u,v,h,hs,re,& ! v-momentum equation: Vx_mid_xt = uh_mid_xt*vh_mid_xt/h_mid_xt Vy_mid_yt = vh_mid_yt*vh_mid_yt/h_mid_yt*cq_s(1:ip,0:jp) - Vy_mid_yt2 = 0.5_sp*g*h_mid_yt**2 + Vy_mid_yt2 = 0.5_wp*g*h_mid_yt**2 vh_new = vh(1:ip,1:jp) & - (dt/(recqdp(0:ip-1,1:jp)))*(Vx_mid_xt(1:ip,1:jp)-Vx_mid_xt(0:ip-1,1:jp)) & - (dt/(recqdq(1:ip,0:jp-1)))*(Vy_mid_yt(1:ip,1:jp)-Vy_mid_yt(1:ip,0:jp-1)) & @@ -302,11 +302,11 @@ subroutine lax_wendroff_ll(ip,jp,o_halo,dt,g,u,v,h,hs,re,& ! add on Coriolis and contribution of orography to pressure gradient: - uh_new=uh_new +dt*.5_sp*(f_cor(1:ip,1:jp)*v(1:ip,1:jp) - & + uh_new=uh_new +dt*.5_wp*(f_cor(1:ip,1:jp)*v(1:ip,1:jp) - & g*(hs(2:ip+1,1:jp)-hs(0:ip-1,1:jp))/(recqdp(1:ip,1:jp)+recqdp(0:ip-1,1:jp)))* & (h(1:ip,1:jp)+h_new) - vh_new=vh_new -dt*.5_sp*(f_cor(1:ip,1:jp)*u(1:ip,1:jp) + & + vh_new=vh_new -dt*.5_wp*(f_cor(1:ip,1:jp)*u(1:ip,1:jp) + & g*(hs(1:ip,2:jp+1)-hs(1:ip,0:jp-1))/(redq(1:ip,1:jp)+redq(1:ip,0:jp-1)))* & (h(1:ip,1:jp)+h_new) @@ -352,18 +352,18 @@ subroutine dissipation(ip,jp,o_halo,dt,f,delsq,re,& theta,thetan,dtheta,dthetan, phi, phin, dphi, dphin, & recq, cq_s, dp1, dq) - use nrtype + use numerics_type implicit none integer(i4b), intent(in) :: ip,jp,o_halo - real(sp), intent(in) :: dt, re - real(sp), dimension(1-o_halo:jp+o_halo), intent(in) :: & + real(wp), intent(in) :: dt, re + real(wp), dimension(1-o_halo:jp+o_halo), intent(in) :: & theta,thetan, dtheta, dthetan - real(sp), dimension(1-o_halo:ip+o_halo), intent(in) :: & + real(wp), dimension(1-o_halo:ip+o_halo), intent(in) :: & phi, phin, dphi, dphin - real(sp), intent(in), dimension(1-o_halo:ip+o_halo,1-o_halo:jp+o_halo) :: & + real(wp), intent(in), dimension(1-o_halo:ip+o_halo,1-o_halo:jp+o_halo) :: & f, & recq, cq_s, dp1, dq - real(sp), intent(inout), dimension(1:ip,1:jp) :: delsq + real(wp), intent(inout), dimension(1:ip,1:jp) :: delsq ! local variables: integer(i4b) :: j, i @@ -371,13 +371,13 @@ subroutine dissipation(ip,jp,o_halo,dt,f,delsq,re,& ! calculate del^2 using 2nd order difference ! (central difference of forward and backward): - delsq(1:ip,1:jp) =1._sp/(re*recq(1:ip,1:jp))* & + delsq(1:ip,1:jp) =1._wp/(re*recq(1:ip,1:jp))* & ( cq_s(1:ip,1:jp)*(f(1:ip,2:jp+1)-f(1:ip,1:jp))/dq(1:ip,1:jp) - & cq_s(1:ip,0:jp-1)*(f(1:ip,1:jp)-f(1:ip,0:jp-1))/dq(1:ip,0:jp-1) ) / & dq(1:ip,1:jp) delsq(1:ip,1:jp) = delsq(1:ip,1:jp) + & - 1._sp/(recq(1:ip,1:jp)**2._sp)* & + 1._wp/(recq(1:ip,1:jp)**2._wp)* & ( (f(2:ip+1,1:jp)-f(1:ip,1:jp))/dp1(1:ip,1:jp) - & (f(1:ip,1:jp)-f(0:ip-1,1:jp))/dp1(1:ip,0:jp-1) ) / & dp1(1:ip,1:jp) @@ -403,31 +403,31 @@ end subroutine dissipation subroutine smagorinsky(ip,jp,o_halo,cvis,u,v,vis,re,& recq, dp1, dq) - use nrtype + use numerics_type implicit none integer(i4b), intent(in) :: ip,jp,o_halo - real(sp), intent(in) :: cvis, re - real(sp), intent(in), dimension(1-o_halo:ip+o_halo,1-o_halo:jp+o_halo) :: & + real(wp), intent(in) :: cvis, re + real(wp), intent(in), dimension(1-o_halo:ip+o_halo,1-o_halo:jp+o_halo) :: & u,v, & recq, dp1, dq - real(sp), intent(inout), dimension(1:ip,1:jp) :: vis + real(wp), intent(inout), dimension(1:ip,1:jp) :: vis ! local variables: integer(i4b) :: j, i ! calculate viscosity using centred differences: - vis(1:ip,1:jp) = cvis**2._sp*re*recq(1:ip,1:jp)*dp1(1:ip,1:jp)*dq(1:ip,1:jp)* & + vis(1:ip,1:jp) = cvis**2._wp*re*recq(1:ip,1:jp)*dp1(1:ip,1:jp)*dq(1:ip,1:jp)* & sqrt( ( (u(2:ip+1,1:jp)-u(0:ip-1,1:jp))/ & - (recq(1:ip,1:jp)*(dp1(0:ip-1,1:jp)+dp1(1:ip,1:jp))) )**2._sp + & + (recq(1:ip,1:jp)*(dp1(0:ip-1,1:jp)+dp1(1:ip,1:jp))) )**2._wp + & ( (v(1:ip,2:jp+1)-v(1:ip,0:jp-1))/ & - (re*(dq(1:ip,1:jp)+dq(1:ip,0:jp-1))) )**2._sp + & - 0.5_sp*( & + (re*(dq(1:ip,1:jp)+dq(1:ip,0:jp-1))) )**2._wp + & + 0.5_wp*( & (u(1:ip,2:jp+1)-u(1:ip,0:jp-1))/ & (re*(dq(1:ip,1:jp)+dq(1:ip,0:jp-1)))+ & (v(2:ip+1,1:jp)-v(0:ip-1,1:jp))/ & (recq(1:ip,1:jp)*(dp1(1:ip,1:jp)+dp1(0:ip-1,1:jp))) & - )**2._sp ) + )**2._wp ) end subroutine smagorinsky diff --git a/driver_code.f90 b/driver_code.f90 index 95e524d..c48c06e 100644 --- a/driver_code.f90 +++ b/driver_code.f90 @@ -3,7 +3,7 @@ !>@brief !>drivers for the shallow water model module drivers - use nrtype + use numerics_type !use variables private public :: model_driver @@ -70,7 +70,7 @@ subroutine model_driver(ip,ipp, jp,jpp, ntim, f, & subgrid_model, viscous_dissipation, dissipate_h,vis, cvis, & vis_eq, lat_eq, & dims,id, world_process, rank, ring_comm) - use nrtype + use numerics_type use mpi_module use advection @@ -82,26 +82,26 @@ subroutine model_driver(ip,ipp, jp,jpp, ntim, f, & integer(i4b), intent(in) :: id, world_process, ring_comm, rank integer(i4b), dimension(2), intent(in) :: coords, dims character (len=*), intent(in) :: outputfile - real(sp), intent(in) :: f, re, g, rho, dt, output_interval - real(sp), dimension(1-o_halo:ipp+o_halo), intent(in) :: phi, phin, dphi, dphin - real(sp), dimension(1-o_halo:jpp+o_halo), intent(in) :: theta, thetan, u_nudge, & + real(wp), intent(in) :: f, re, g, rho, dt, output_interval + real(wp), dimension(1-o_halo:ipp+o_halo), intent(in) :: phi, phin, dphi, dphin + real(wp), dimension(1-o_halo:jpp+o_halo), intent(in) :: theta, thetan, u_nudge, & dtheta, dthetan - real(sp), dimension(1-o_halo:ipp+o_halo,1-o_halo:jpp+o_halo), & + real(wp), dimension(1-o_halo:ipp+o_halo,1-o_halo:jpp+o_halo), & intent(in) :: f_cor, hs, & recqdp, recqdp_s, recqdq_s, redq_s, redq, & recq, cq_s, cq, dp1, dq, recqdq - real(sp), dimension(1-o_halo:ipp+o_halo,1-o_halo:jpp+o_halo), & + real(wp), dimension(1-o_halo:ipp+o_halo,1-o_halo:jpp+o_halo), & intent(inout) :: h, u, v, height - real(sp), dimension(1-o_halo:ipp+o_halo,1-o_halo:jpp+o_halo), & + real(wp), dimension(1-o_halo:ipp+o_halo,1-o_halo:jpp+o_halo), & intent(in) :: dx, dy, x, y - real(sp), intent(in) :: vis, nudge_tau, cvis, lat_eq, vis_eq + real(wp), intent(in) :: vis, nudge_tau, cvis, lat_eq, vis_eq ! locals: integer(i4b) :: n, cur=1, j, error, rank2 - real(sp) :: time, time_last_output, output_time - real(sp), dimension(1-o_halo:ipp+o_halo,1-o_halo:jpp+o_halo) :: & + real(wp) :: time, time_last_output, output_time + real(wp), dimension(1-o_halo:ipp+o_halo,1-o_halo:jpp+o_halo) :: & u_old, v_old, h_old - real(sp), dimension(1:ipp,1:jpp) :: delsq, vort, visco + real(wp), dimension(1:ipp,1:jpp) :: delsq, vort, visco time_last_output=-output_interval @@ -116,11 +116,11 @@ subroutine model_driver(ip,ipp, jp,jpp, ntim, f, & !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! write netcdf variables ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - time=real(n-1,sp)*dt + time=real(n-1,wp)*dt if (time-time_last_output >= output_interval) then if (id==world_process) & print *,'output no ',cur,' at time (hrs) ', & - time/3600._sp,n,' steps of ',ntim + time/3600._wp,n,' steps of ',ntim !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! calculate diagnostics for output: vorticity, etc ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! @@ -168,21 +168,21 @@ subroutine model_driver(ip,ipp, jp,jpp, ntim, f, & ! mid-point rule: u(1:ipp,j)=u(1:ipp,j)+ & (u_nudge(j)- & - (0.5_sp*(u(1:ipp,j)+u_old(1:ipp,j)))/real(1,sp) ) & + (0.5_wp*(u(1:ipp,j)+u_old(1:ipp,j)))/real(1,wp) ) & /nudge_tau * dt ! v(1:ipp,j)=v(1:ipp,j)+& -! (0._sp- & -! sum(0.5_sp*(v(1:ipp,j)+v_old(1:ipp,j)))/real(ipp,sp) ) & +! (0._wp- & +! sum(0.5_wp*(v(1:ipp,j)+v_old(1:ipp,j)))/real(ipp,wp) ) & ! /nudge_tau *dt ! Derived by integrating du/dt=(u_nudge-u)/tau ! u(1:ipp,j)=u_nudge(j)- & ! (u_nudge(j)- & -! (u(1:ipp,j))/real(1,sp) ) * & +! (u(1:ipp,j))/real(1,wp) ) * & ! exp(-dt/nudge_tau ) -! v(1:ipp,j)=0._sp- & -! (0._sp- & -! (v(1:ipp,j))/real(1,sp) ) * & +! v(1:ipp,j)=0._wp- & +! (0._wp- & +! (v(1:ipp,j))/real(1,wp) ) * & ! exp(-dt/nudge_tau ) enddo @@ -211,7 +211,7 @@ subroutine model_driver(ip,ipp, jp,jpp, ntim, f, & !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! dissipate u ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - call dissipation(ipp,jpp,o_halo,dt,0.5_sp*(u_old+u), delsq,re,& + call dissipation(ipp,jpp,o_halo,dt,0.5_wp*(u_old+u), delsq,re,& theta,thetan,dtheta,dthetan, phi, phin, dphi, dphin, & recq, cq_s, dp1, dq) @@ -219,8 +219,8 @@ subroutine model_driver(ip,ipp, jp,jpp, ntim, f, & case (1) u(1:ipp,1:jpp)=u(1:ipp,1:jpp)+dt*delsq*vis case (2) - call smagorinsky(ipp,jpp,o_halo,cvis,0.5_sp*(u_old+u),& - 0.5_sp*(v_old+v),visco,re,recq, dp1, dq) + call smagorinsky(ipp,jpp,o_halo,cvis,0.5_wp*(u_old+u),& + 0.5_wp*(v_old+v),visco,re,recq, dp1, dq) u(1:ipp,1:jpp)=u(1:ipp,1:jpp)+dt*delsq*visco case default print *,'error subgrid ',subgrid_model @@ -236,7 +236,7 @@ subroutine model_driver(ip,ipp, jp,jpp, ntim, f, & !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! dissipate v ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - call dissipation(ipp,jpp,o_halo,dt,0.5_sp*(v_old+v), delsq,re,& + call dissipation(ipp,jpp,o_halo,dt,0.5_wp*(v_old+v), delsq,re,& theta,thetan,dtheta,dthetan, phi, phin, dphi, dphin, & recq, cq_s, dp1, dq) @@ -252,11 +252,11 @@ subroutine model_driver(ip,ipp, jp,jpp, ntim, f, & do j=1,jpp - if((theta(j) >-lat_eq*pi/180._sp) .and. & - (theta(j) < lat_eq*pi/180._sp)) then + if((theta(j) >-lat_eq*pi/180._wp) .and. & + (theta(j) < lat_eq*pi/180._wp)) then v(1:ipp,j)=v(1:ipp,j)+dt*delsq(1:ipp,j)*vis_eq* & - cos(theta(j)*90._sp/lat_eq) + cos(theta(j)*90._wp/lat_eq) endif enddo @@ -273,7 +273,7 @@ subroutine model_driver(ip,ipp, jp,jpp, ntim, f, & call exchange_halos(ring_comm, id, ipp, jpp, o_halo, h) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - call dissipation(ipp,jpp,o_halo,dt,0.5_sp*(h_old+h), delsq,re,& + call dissipation(ipp,jpp,o_halo,dt,0.5_wp*(h_old+h), delsq,re,& theta,thetan,dtheta,dthetan, phi, phin, dphi, dphin, & recq, cq_s, dp1, dq) @@ -300,14 +300,14 @@ subroutine model_driver(ip,ipp, jp,jpp, ntim, f, & !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! if(coords(2)==(dims(2)-1)) then -! v(:,jpp+1:jpp+o_halo)=0._sp +! v(:,jpp+1:jpp+o_halo)=0._wp ! do j=jpp+1,jpp+o_halo ! u(:,j)=u_nudge(j) ! enddo ! h(:,jpp+o_halo)=h(:,jpp) ! endif ! if(coords(2)==0) then -! v(:,1-o_halo)=0._sp +! v(:,1-o_halo)=0._wp ! do j=1-o_halo,0 ! u(:,j)=u_nudge(j) ! enddo @@ -363,18 +363,18 @@ subroutine diagnostics(ip,jp,o_halo,dt,u,v,vort,re,& theta,thetan,dtheta,dthetan, phi, phin, dphi, dphin, & recq, cq_s, cq, dp1, dq) - use nrtype + use numerics_type implicit none integer(i4b), intent(in) :: ip,jp,o_halo - real(sp), intent(in) :: dt, re - real(sp), dimension(1-o_halo:jp+o_halo), intent(in) :: & + real(wp), intent(in) :: dt, re + real(wp), dimension(1-o_halo:jp+o_halo), intent(in) :: & theta,thetan, dtheta, dthetan - real(sp), dimension(1-o_halo:ip+o_halo), intent(in) :: & + real(wp), dimension(1-o_halo:ip+o_halo), intent(in) :: & phi, phin, dphi, dphin - real(sp), intent(in), dimension(1-o_halo:ip+o_halo,1-o_halo:jp+o_halo) :: & + real(wp), intent(in), dimension(1-o_halo:ip+o_halo,1-o_halo:jp+o_halo) :: & u,v, & recq, cq_s, cq, dp1, dq - real(sp), intent(inout), dimension(1:ip,1:jp) :: vort + real(wp), intent(inout), dimension(1:ip,1:jp) :: vort ! local variables: integer(i4b) :: j, i @@ -382,7 +382,7 @@ subroutine diagnostics(ip,jp,o_halo,dt,u,v,vort,re,& ! calculate relative vorticity ! (central difference ): - vort(1:ip,1:jp) =-1._sp/(recq(1:ip,1:jp))* & + vort(1:ip,1:jp) =-1._wp/(recq(1:ip,1:jp))* & ( (u(1:ip,2:jp+1)*cq(1:ip,2:jp+1)- & u(1:ip,0:jp-1)*cq(1:ip,0:jp-1))/dq(1:ip,1:jp) - & (v(2:ip+1,1:jp)-v(0:ip-1,1:jp))/dp1(1:ip,1:jp) ) @@ -435,13 +435,13 @@ subroutine output(new_file,outputfile,n,ip,ipp,ipstart,jp,jpp,jpstart,o_halo, & logical, intent(inout) :: new_file character (len=*), intent(in) :: outputfile integer(i4b), intent(in) :: n, ip, ipp, ipstart, jp, jpp, jpstart, o_halo - real(sp), intent(in) :: time - real(sp), dimension(1-o_halo:ipp+o_halo), intent(in) :: phi - real(sp), dimension(1-o_halo:jpp+o_halo), intent(in) :: theta, u_nudge + real(wp), intent(in) :: time + real(wp), dimension(1-o_halo:ipp+o_halo), intent(in) :: phi + real(wp), dimension(1-o_halo:jpp+o_halo), intent(in) :: theta, u_nudge - real(sp), dimension(1-o_halo:ipp+o_halo,1-o_halo:jpp+o_halo), & + real(wp), dimension(1-o_halo:ipp+o_halo,1-o_halo:jpp+o_halo), & intent(in) :: f_cor, height, h, u, v - real(sp), dimension(1:ipp,1:jpp), intent(in) :: vort + real(wp), dimension(1:ipp,1:jpp), intent(in) :: vort integer(i4b), intent(in) :: id ,world_process, rank, ring_comm @@ -726,7 +726,7 @@ end subroutine output !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine check(status) use netcdf - use nrtype + use numerics_type integer(i4b), intent ( in) :: status if(status /= nf90_noerr) then diff --git a/hygfx.for b/hygfx.for deleted file mode 100755 index fb63507..0000000 --- a/hygfx.for +++ /dev/null @@ -1,307 +0,0 @@ - module hypergeo - private - public :: hygfx - contains - - SUBROUTINE HYGFX(A,B,C,X,HF) -C -C ==================================================== -C Purpose: Compute hypergeometric function F(a,b,c,x) -C Input : a --- Parameter -C b --- Parameter -C c --- Parameter, c <> 0,-1,-2,... -C x --- Argument ( x < 1 ) -C Output: HF --- F(a,b,c,x) -C Routines called: -C (1) GAMMA for computing gamma function -C (2) PSI for computing psi function -C ==================================================== -C - IMPLICIT DOUBLE PRECISION (A-H,O-Z) - LOGICAL L0,L1,L2,L3,L4,L5 - PI=3.141592653589793D0 - EL=.5772156649015329D0 - L0=C.EQ.INT(C).AND.C.LT.0.0 - L1=1.0D0-X.LT.1.0D-15.AND.C-A-B.LE.0.0 - L2=A.EQ.INT(A).AND.A.LT.0.0 - L3=B.EQ.INT(B).AND.B.LT.0.0 - L4=C-A.EQ.INT(C-A).AND.C-A.LE.0.0 - L5=C-B.EQ.INT(C-B).AND.C-B.LE.0.0 - IF (L0.OR.L1) THEN - WRITE(*,*)'The hypergeometric series is divergent' - RETURN - ENDIF - EPS=1.0D-15 - IF (X.GT.0.95) EPS=1.0D-8 - IF (X.EQ.0.0.OR.A.EQ.0.0.OR.B.EQ.0.0) THEN - HF=1.0D0 - RETURN - ELSE IF (1.0D0-X.EQ.EPS.AND.C-A-B.GT.0.0) THEN - CALL GAMMA(C,GC) - CALL GAMMA(C-A-B,GCAB) - CALL GAMMA(C-A,GCA) - CALL GAMMA(C-B,GCB) - HF=GC*GCAB/(GCA*GCB) - RETURN - ELSE IF (1.0D0+X.LE.EPS.AND.DABS(C-A+B-1.0).LE.EPS) THEN - G0=DSQRT(PI)*2.0D0**(-A) - CALL GAMMA(C,G1) - CALL GAMMA(1.0D0+A/2.0-B,G2) - CALL GAMMA(0.5D0+0.5*A,G3) - HF=G0*G1/(G2*G3) - RETURN - ELSE IF (L2.OR.L3) THEN - IF (L2) NM=INT(ABS(A)) - IF (L3) NM=INT(ABS(B)) - HF=1.0D0 - R=1.0D0 - DO 10 K=1,NM - R=R*(A+K-1.0D0)*(B+K-1.0D0)/(K*(C+K-1.0D0))*X -10 HF=HF+R - RETURN - ELSE IF (L4.OR.L5) THEN - IF (L4) NM=INT(ABS(C-A)) - IF (L5) NM=INT(ABS(C-B)) - HF=1.0D0 - R=1.0D0 - DO 15 K=1,NM - R=R*(C-A+K-1.0D0)*(C-B+K-1.0D0)/(K*(C+K-1.0D0))*X -15 HF=HF+R - HF=(1.0D0-X)**(C-A-B)*HF - RETURN - ENDIF - AA=A - BB=B - X1=X - IF (X.LT.0.0D0) THEN - X=X/(X-1.0D0) - IF (C.GT.A.AND.B.LT.A.AND.B.GT.0.0) THEN - A=BB - B=AA - ENDIF - B=C-B - ENDIF - IF (X.GE.0.75D0) THEN - GM=0.0D0 - IF (DABS(C-A-B-INT(C-A-B)).LT.1.0D-15) THEN - M=INT(C-A-B) - CALL GAMMA(A,GA) - CALL GAMMA(B,GB) - CALL GAMMA(C,GC) - CALL GAMMA(A+M,GAM) - CALL GAMMA(B+M,GBM) - CALL PSI(A,PA) - CALL PSI(B,PB) - IF (M.NE.0) GM=1.0D0 - DO 30 J=1,ABS(M)-1 -30 GM=GM*J - RM=1.0D0 - DO 35 J=1,ABS(M) -35 RM=RM*J - F0=1.0D0 - R0=1.0D0 - R1=1.0D0 - SP0=0.D0 - SP=0.0D0 - IF (M.GE.0) THEN - C0=GM*GC/(GAM*GBM) - C1=-GC*(X-1.0D0)**M/(GA*GB*RM) - DO 40 K=1,M-1 - R0=R0*(A+K-1.0D0)*(B+K-1.0)/(K*(K-M))*(1.0-X) -40 F0=F0+R0 - DO 45 K=1,M -45 SP0=SP0+1.0D0/(A+K-1.0)+1.0/(B+K-1.0)-1.0/K - F1=PA+PB+SP0+2.0D0*EL+DLOG(1.0D0-X) - DO 55 K=1,250 - SP=SP+(1.0D0-A)/(K*(A+K-1.0))+(1.0-B)/(K*(B+K-1.0)) - SM=0.0D0 - DO 50 J=1,M -50 SM=SM+(1.0D0-A)/((J+K)*(A+J+K-1.0))+1.0/ - & (B+J+K-1.0) - RP=PA+PB+2.0D0*EL+SP+SM+DLOG(1.0D0-X) - R1=R1*(A+M+K-1.0D0)*(B+M+K-1.0)/(K*(M+K))*(1.0-X) - F1=F1+R1*RP - IF (DABS(F1-HW).LT.DABS(F1)*EPS) GO TO 60 -55 HW=F1 -60 HF=F0*C0+F1*C1 - ELSE IF (M.LT.0) THEN - M=-M - C0=GM*GC/(GA*GB*(1.0D0-X)**M) - C1=-(-1)**M*GC/(GAM*GBM*RM) - DO 65 K=1,M-1 - R0=R0*(A-M+K-1.0D0)*(B-M+K-1.0)/(K*(K-M))*(1.0-X) -65 F0=F0+R0 - DO 70 K=1,M -70 SP0=SP0+1.0D0/K - F1=PA+PB-SP0+2.0D0*EL+DLOG(1.0D0-X) - DO 80 K=1,250 - SP=SP+(1.0D0-A)/(K*(A+K-1.0))+(1.0-B)/(K*(B+K-1.0)) - SM=0.0D0 - DO 75 J=1,M -75 SM=SM+1.0D0/(J+K) - RP=PA+PB+2.0D0*EL+SP-SM+DLOG(1.0D0-X) - R1=R1*(A+K-1.0D0)*(B+K-1.0)/(K*(M+K))*(1.0-X) - F1=F1+R1*RP - IF (DABS(F1-HW).LT.DABS(F1)*EPS) GO TO 85 -80 HW=F1 -85 HF=F0*C0+F1*C1 - ENDIF - ELSE - CALL GAMMA(A,GA) - CALL GAMMA(B,GB) - CALL GAMMA(C,GC) - CALL GAMMA(C-A,GCA) - CALL GAMMA(C-B,GCB) - CALL GAMMA(C-A-B,GCAB) - CALL GAMMA(A+B-C,GABC) - C0=GC*GCAB/(GCA*GCB) - C1=GC*GABC/(GA*GB)*(1.0D0-X)**(C-A-B) - HF=0.0D0 - R0=C0 - R1=C1 - DO 90 K=1,250 - R0=R0*(A+K-1.0D0)*(B+K-1.0)/(K*(A+B-C+K))*(1.0-X) - R1=R1*(C-A+K-1.0D0)*(C-B+K-1.0)/(K*(C-A-B+K)) - & *(1.0-X) - HF=HF+R0+R1 - IF (DABS(HF-HW).LT.DABS(HF)*EPS) GO TO 95 -90 HW=HF -95 HF=HF+C0+C1 - ENDIF - ELSE - A0=1.0D0 - IF (C.GT.A.AND.C.LT.2.0D0*A.AND. - & C.GT.B.AND.C.LT.2.0D0*B) THEN - A0=(1.0D0-X)**(C-A-B) - A=C-A - B=C-B - ENDIF - HF=1.0D0 - R=1.0D0 - DO 100 K=1,250 - R=R*(A+K-1.0D0)*(B+K-1.0D0)/(K*(C+K-1.0D0))*X - HF=HF+R - IF (DABS(HF-HW).LE.DABS(HF)*EPS) GO TO 105 -100 HW=HF -105 HF=A0*HF - ENDIF - IF (X1.LT.0.0D0) THEN - X=X1 - C0=1.0D0/(1.0D0-X)**AA - HF=C0*HF - ENDIF - A=AA - B=BB - IF (K.GT.120) WRITE(*,115) -115 FORMAT(1X,'Warning! You should check the accuracy') - RETURN - END - - - SUBROUTINE GAMMA(X,GA) -C -C ================================================== -C Purpose: Compute gamma function â(x) -C Input : x --- Argument of â(x) -C ( x is not equal to 0,-1,-2,úúú) -C Output: GA --- â(x) -C ================================================== -C - IMPLICIT DOUBLE PRECISION (A-H,O-Z) - DIMENSION G(26) - PI=3.141592653589793D0 - IF (X.EQ.INT(X)) THEN - IF (X.GT.0.0D0) THEN - GA=1.0D0 - M1=X-1 - DO 10 K=2,M1 -10 GA=GA*K - ELSE - GA=1.0D+300 - ENDIF - ELSE - IF (DABS(X).GT.1.0D0) THEN - Z=DABS(X) - M=INT(Z) - R=1.0D0 - DO 15 K=1,M -15 R=R*(Z-K) - Z=Z-M - ELSE - Z=X - ENDIF - DATA G/1.0D0,0.5772156649015329D0, - & -0.6558780715202538D0, -0.420026350340952D-1, - & 0.1665386113822915D0,-.421977345555443D-1, - & -.96219715278770D-2, .72189432466630D-2, - & -.11651675918591D-2, -.2152416741149D-3, - & .1280502823882D-3, -.201348547807D-4, - & -.12504934821D-5, .11330272320D-5, - & -.2056338417D-6, .61160950D-8, - & .50020075D-8, -.11812746D-8, - & .1043427D-9, .77823D-11, - & -.36968D-11, .51D-12, - & -.206D-13, -.54D-14, .14D-14, .1D-15/ - GR=G(26) - DO 20 K=25,1,-1 -20 GR=GR*Z+G(K) - GA=1.0D0/(GR*Z) - IF (DABS(X).GT.1.0D0) THEN - GA=GA*R - IF (X.LT.0.0D0) GA=-PI/(X*GA*DSIN(PI*X)) - ENDIF - ENDIF - RETURN - END - - - SUBROUTINE PSI(X,PS) -C -C ====================================== -C Purpose: Compute Psi function -C Input : x --- Argument of psi(x) -C Output: PS --- psi(x) -C ====================================== -C - IMPLICIT DOUBLE PRECISION (A-H,O-Z) - XA=DABS(X) - PI=3.141592653589793D0 - EL=.5772156649015329D0 - S=0.0D0 - IF (X.EQ.INT(X).AND.X.LE.0.0) THEN - PS=1.0D+300 - RETURN - ELSE IF (XA.EQ.INT(XA)) THEN - N=XA - DO 10 K=1 ,N-1 -10 S=S+1.0D0/K - PS=-EL+S - ELSE IF (XA+.5.EQ.INT(XA+.5)) THEN - N=XA-.5 - DO 20 K=1,N -20 S=S+1.0/(2.0D0*K-1.0D0) - PS=-EL+2.0D0*S-1.386294361119891D0 - ELSE - IF (XA.LT.10.0) THEN - N=10-INT(XA) - DO 30 K=0,N-1 -30 S=S+1.0D0/(XA+K) - XA=XA+N - ENDIF - X2=1.0D0/(XA*XA) - A1=-.8333333333333D-01 - A2=.83333333333333333D-02 - A3=-.39682539682539683D-02 - A4=.41666666666666667D-02 - A5=-.75757575757575758D-02 - A6=.21092796092796093D-01 - A7=-.83333333333333333D-01 - A8=.4432598039215686D0 - PS=DLOG(XA)-.5D0/XA+X2*(((((((A8*X2+A7)*X2+ - & A6)*X2+A5)*X2+A4)*X2+A3)*X2+A2)*X2+A1) - PS=PS-S - ENDIF - IF (X.LT.0.0) PS=PS-PI*DCOS(PI*X)/DSIN(PI*X)-1.0D0/X - RETURN - END - end module hypergeo diff --git a/initialisation.f90 b/initialisation.f90 index 8cb6352..4c2d292 100644 --- a/initialisation.f90 +++ b/initialisation.f90 @@ -3,7 +3,7 @@ !>@brief !>initialisation for the shallow water model module initialisation - use nrtype + use numerics_type private public :: allocate_and_set @@ -99,23 +99,23 @@ subroutine allocate_and_set(ipp,jpp,ntim, f, & rotation_period_hours, scale_height, slat, nlat, & slat_thresh, nlat_thresh, dims, id, comm2d) - use nrtype + use numerics_type use mpi use netcdf - use nr, only : locate, polint + use numerics, only : find_pos, poly_int use random, only : random_normal use mpi_module implicit none ! grid variables to be set integer(i4b), intent(inout) :: ipp, jpp, ntim, o_halo - real(sp), intent(inout) :: f, re, g, rho, dt - real(sp), intent(inout), allocatable, dimension(:,:) :: & + real(wp), intent(inout) :: f, re, g, rho, dt + real(wp), intent(inout), allocatable, dimension(:,:) :: & f_cor, h, hs, u, v, height, & dx, dy, x, y, & recqdp, recqdp_s, recqdq_s, redq_s, redq, & recq, cq_s, cq, dp1, dq, recqdq - real(sp), intent(inout), allocatable, dimension(:) :: & + real(wp), intent(inout), allocatable, dimension(:) :: & phi, theta, phin, thetan,u_nudge, & dphi, dtheta, dphin, dthetan @@ -123,12 +123,12 @@ subroutine allocate_and_set(ipp,jpp,ntim, f, & character (len=*), intent(in) :: inputfile logical, intent(in) :: add_random_height_noise, initially_geostrophic integer(i4b), intent(in) :: initial_winds, ip, jp - real(sp), intent(in) :: wind_factor, wind_shift, wind_reduce, runtime, & + real(wp), intent(in) :: wind_factor, wind_shift, wind_reduce, runtime, & dt_nm, grav, rho_nm, re_nm, & rotation_period_hours, scale_height, & slat_thresh, nlat_thresh, & u_jet, theta_jet, h_jet - real(sp), intent(inout) :: slat, nlat + real(wp), intent(inout) :: slat, nlat integer(i4b), dimension(2), intent(in) :: dims integer(i4b), dimension(2), intent(inout) :: coords integer(i4b), intent(inout) :: ipstart, jpstart @@ -136,11 +136,11 @@ subroutine allocate_and_set(ipp,jpp,ntim, f, & ! locals integer(i4b) :: iloc, error, AllocateStatus, ncid, varid1,varid2, dimid, nlats, & i, j - real(sp), dimension(:), allocatable :: latitude, wind - real(sp) :: var, dummy, delta_omega, slat_thresh2, nlat_thresh2 + real(wp), dimension(:), allocatable :: latitude, wind + real(wp) :: var, dummy, delta_omega, slat_thresh2, nlat_thresh2 ! for random number: - real(sp) :: r - real(sp), dimension(10,10) :: rs + real(wp) :: r + real(wp), dimension(10,10) :: rs integer(i4b) :: k, nbottom, ntop, tag1 integer(i4b), allocatable, dimension(:) :: seed @@ -158,7 +158,7 @@ subroutine allocate_and_set(ipp,jpp,ntim, f, & ! 2. scalar formulae: ntim=ceiling(runtime/dt) - f=2._sp*PI / (rotation_period_hours*3600._sp) + f=2._wp*PI / (rotation_period_hours*3600._wp) ! 3. arrays: @@ -172,13 +172,13 @@ subroutine allocate_and_set(ipp,jpp,ntim, f, & ! print *,'Coords of ',id,' are ',coords ! number of grid points in all but last: - ipp = floor(real(ip,sp)/real(dims(1),sp)) + ipp = floor(real(ip,wp)/real(dims(1),wp)) ipstart = ipp*(coords(1)) if(coords(1) == (dims(1)-1)) then ipp=ip-(dims(1)-1)*ipp ! number of grid points in last endif ! number of grid points in all but last: - jpp = floor(real(jp,sp)/real(dims(2),sp)) + jpp = floor(real(jp,wp)/real(dims(2),wp)) jpstart = jpp*(coords(2)) if(coords(2) == (dims(2)-1)) then jpp=jp-(dims(2)-1)*jpp ! number of grid points in last @@ -301,26 +301,26 @@ subroutine allocate_and_set(ipp,jpp,ntim, f, & ! interpolate to grid, do finite diffs, pass halos, etc - dphi=(2._sp*PI) / real(ip-1,sp) ! lon - dphin=(2._sp*PI) / real(ip-1,sp) ! lon + dphi=(2._wp*PI) / real(ip-1,wp) ! lon + dphin=(2._wp*PI) / real(ip-1,wp) ! lon ! set up longitude array: phi=dphi(1)*(/(i,i=ipstart+1-o_halo-1,ipstart+ipp+o_halo-1)/) - phin=phi+dphi/2._sp + phin=phi+dphi/2._wp nlat=min(nlat,nlat_thresh) slat=max(slat,slat_thresh) ! latitude array: dtheta=(nlat-slat) & - / real(jp-1,sp) * PI/180._sp ! lat + / real(jp-1,wp) * PI/180._wp ! lat ! deal with singularity at the poles: -! if((coords(2) == (dims(2)-1)) .and. (nlat .gt. 0._sp)) then -! dtheta(jpp) = 2._sp*(90._sp-nlat)*PI/180._sp +! if((coords(2) == (dims(2)-1)) .and. (nlat .gt. 0._wp)) then +! dtheta(jpp) = 2._wp*(90._wp-nlat)*PI/180._wp ! endif -! if((coords(2) == 0) .and. (slat .lt. 0._sp)) then -! dtheta(0) = 2._sp*(90._sp+slat)*PI/180._sp +! if((coords(2) == 0) .and. (slat .lt. 0._wp)) then +! dtheta(0) = 2._wp*(90._wp+slat)*PI/180._wp ! endif @@ -337,7 +337,7 @@ subroutine allocate_and_set(ipp,jpp,ntim, f, & endif if(coords(2) == 0) then - theta(0) = slat*PI/180._sp-dtheta(0) + theta(0) = slat*PI/180._wp-dtheta(0) endif do i=1,jpp+o_halo @@ -347,14 +347,14 @@ subroutine allocate_and_set(ipp,jpp,ntim, f, & if(ntop .ne. -1) then tag1=001 call MPI_Send(theta(jpp), 1, MPI_REAL8, ntop, & - tag1, MPI_COMM_WORLD, MPI_STATUS_IGNORE,error) + tag1, MPI_COMM_WORLD,error) endif -! theta=dtheta(10)*(/(i,i=jpstart+1-o_halo-1,jpstart+jpp+o_halo-1)/) + slat*PI/180._sp +! theta=dtheta(10)*(/(i,i=jpstart+1-o_halo-1,jpstart+jpp+o_halo-1)/) + slat*PI/180._wp - thetan=theta+dtheta/2._sp + thetan=theta+dtheta/2._wp do i=1-o_halo,jpp dthetan(i)=thetan(i+1)-thetan(i) enddo @@ -362,7 +362,7 @@ subroutine allocate_and_set(ipp,jpp,ntim, f, & ! if this is the top then set dtheta: if(ntop == -1 ) then - dthetan(jpp+1)=2._sp*(90._sp*PI/180._sp -dthetan(jpp)) + dthetan(jpp+1)=2._wp*(90._wp*PI/180._wp -dthetan(jpp)) endif if(ntop /= -1 ) then call MPI_Recv(dthetan(jpp+o_halo), 1, MPI_REAL8, ntop, & @@ -370,9 +370,9 @@ subroutine allocate_and_set(ipp,jpp,ntim, f, & endif if(nbottom /= -1 ) then call MPI_Send(dthetan(1), 1, MPI_REAL8, nbottom, & - tag1, MPI_COMM_WORLD, MPI_STATUS_IGNORE,error) + tag1, MPI_COMM_WORLD,error) endif -! dthetan=(nlat-slat) / real(jp-1,sp) * PI/180._sp ! lat +! dthetan=(nlat-slat) / real(jp-1,wp) * PI/180._wp ! lat !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! @@ -382,19 +382,19 @@ subroutine allocate_and_set(ipp,jpp,ntim, f, & !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! case (1) ! interpolate winds to u_nudge - delta_omega=2._sp*PI/(3600._sp)*(1._sp/10.656_sp-1._sp/rotation_period_hours) + delta_omega=2._wp*PI/(3600._wp)*(1._wp/10.656_wp-1._wp/rotation_period_hours) do i=1-o_halo,jpp+o_halo - iloc=locate(latitude(1:nlats),asin(sin(theta(i)))*180._sp/PI) + iloc=find_pos(latitude(1:nlats),asin(sin(theta(i)))*180._wp/PI) iloc=min(nlats-1,iloc) iloc=max(1,iloc) ! linear interp theta - call polint(latitude(iloc:iloc+1), wind(iloc:iloc+1), & - min(max( asin(sin(theta(i)))*180._sp/PI, & + call poly_int(latitude(iloc:iloc+1), wind(iloc:iloc+1), & + min(max( asin(sin(theta(i)))*180._wp/PI, & latitude(nlats)),latitude(1)), var,dummy) - if((theta(i)*180._sp/pi>90._sp) .or. & - (theta(i)*180._sp/pi<-90._sp) ) var=-var + if((theta(i)*180._wp/pi>90._wp) .or. & + (theta(i)*180._wp/pi<-90._wp) ) var=-var var=var+delta_omega*re*cos(theta(i)) @@ -407,8 +407,8 @@ subroutine allocate_and_set(ipp,jpp,ntim, f, & !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! case (2) do i=1-o_halo,jpp+o_halo - u_nudge(i)=u_jet* exp(-0.5_sp*(theta(i)-theta_jet*pi/180._sp)**2._sp & - / (h_jet*pi/180._sp)**2._sp ) + u_nudge(i)=u_jet* exp(-0.5_wp*(theta(i)-theta_jet*pi/180._wp)**2._wp & + / (h_jet*pi/180._wp)**2._wp ) enddo !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! @@ -419,7 +419,7 @@ subroutine allocate_and_set(ipp,jpp,ntim, f, & ! calculate Coriolis param: do i=1-o_halo,ipp+o_halo - f_cor(i,:)=2._sp*f*sin(theta) + f_cor(i,:)=2._wp*f*sin(theta) enddo ! if(coords(2)==0) f_cor(:,1-o_halo)=-f_cor(:,1-o_halo) ! if(coords(2)==dims(2)-1) f_cor(:,jpp+o_halo)=-f_cor(:,jpp+o_halo) @@ -441,13 +441,13 @@ subroutine allocate_and_set(ipp,jpp,ntim, f, & ! set u and v: - v(:,:)=0._sp + v(:,:)=0._wp do j=1-o_halo,jpp+o_halo u(:,j)=u_nudge(j) enddo ! set surface to zero. - hs(:,:)=0._sp + hs(:,:)=0._wp @@ -477,7 +477,7 @@ subroutine allocate_and_set(ipp,jpp,ntim, f, & do j=jpp,0,-1 height(:,j)=height(:,j+1)+ & - 0.25_sp*(f_cor(:,j+1)+f_cor(:,j))* & + 0.25_wp*(f_cor(:,j+1)+f_cor(:,j))* & (u(:,j+1)+u(:,j))* & re/g*(dtheta(j)) enddo @@ -488,7 +488,7 @@ subroutine allocate_and_set(ipp,jpp,ntim, f, & ipp+2*o_halo, & ! size of the data packet MPI_REAL8, & nbottom, & ! send to below (south) - tag1, MPI_COMM_WORLD, MPI_STATUS_IGNORE,error) + tag1, MPI_COMM_WORLD,error) endif !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! @@ -519,13 +519,13 @@ subroutine allocate_and_set(ipp,jpp,ntim, f, & if((i > ipstart) .and. (i <=ipstart+ipp) & .and. (j > jpstart) .and. (j <= jpstart+jpp) ) then - if ((theta(j-jpstart)*180._sp/PI) > 75._sp & - .and. (theta(j-jpstart)*180._sp/PI) < 80._sp) then + if ((theta(j-jpstart)*180._wp/PI) > 75._wp & + .and. (theta(j-jpstart)*180._wp/PI) < 80._wp) then height(i-ipstart,j-jpstart) = & height(i-ipstart,j-jpstart) + & - r*1000.e0_sp*0.6e5_sp/height(i-ipstart,j-jpstart) ! *& - !abs(f_cor(i-ipstart,j-jpstart))/3e-4_sp + r*1000.e0_wp*0.6e5_wp/height(i-ipstart,j-jpstart) ! *& + !abs(f_cor(i-ipstart,j-jpstart))/3e-4_wp endif endif @@ -544,13 +544,13 @@ subroutine allocate_and_set(ipp,jpp,ntim, f, & if((i > ipstart) .and. (i <=ipstart+ipp) & .and. (j > jpstart) .and. (j <= jpstart+jpp) ) then - if ((theta(j-jpstart)*180._sp/PI) > (theta_jet-h_jet*3._sp) & - .and. (theta(j-jpstart)*180._sp/PI) <(theta_jet+h_jet*3._sp)) then + if ((theta(j-jpstart)*180._wp/PI) > (theta_jet-h_jet*3._wp) & + .and. (theta(j-jpstart)*180._wp/PI) <(theta_jet+h_jet*3._wp)) then height(i-ipstart,j-jpstart) = & height(i-ipstart,j-jpstart) + & - r*1000.e0_sp*0.6e5_sp/height(i-ipstart,j-jpstart) ! *& - !abs(f_cor(i-ipstart,j-jpstart))/3e-4_sp + r*1000.e0_wp*0.6e5_wp/height(i-ipstart,j-jpstart) ! *& + !abs(f_cor(i-ipstart,j-jpstart))/3e-4_wp endif endif @@ -631,7 +631,7 @@ end subroutine allocate_and_set !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine check(status) use netcdf - use nrtype + use numerics_type integer(i4b), intent ( in) :: status if(status /= nf90_noerr) then diff --git a/locate.f90 b/locate.f90 deleted file mode 100755 index e7c1806..0000000 --- a/locate.f90 +++ /dev/null @@ -1,29 +0,0 @@ - FUNCTION locate(xx,x) - USE nrtype - IMPLICIT NONE - REAL(SP), DIMENSION(:), INTENT(IN) :: xx - REAL(SP), INTENT(IN) :: x - INTEGER(I4B) :: locate - INTEGER(I4B) :: n,jl,jm,ju - LOGICAL :: ascnd - n=size(xx) - ascnd = (xx(n) >= xx(1)) - jl=0 - ju=n+1 - do - if (ju-jl <= 1) exit - jm=(ju+jl)/2 - if (ascnd .eqv. (x >= xx(jm))) then - jl=jm - else - ju=jm - end if - end do - if (x == xx(1)) then - locate=1 - else if (x == xx(n)) then - locate=n-1 - else - locate=jl - end if - END FUNCTION locate diff --git a/main.f90 b/main.f90 index e961759..adf9c83 100644 --- a/main.f90 +++ b/main.f90 @@ -47,7 +47,6 @@ !>main programme reads in information, allocates arrays, then calls the model driver program main - use nrtype use variables use mpi use mpi_module diff --git a/mpi_module.f90 b/mpi_module.f90 index 70c389e..606de97 100644 --- a/mpi_module.f90 +++ b/mpi_module.f90 @@ -3,9 +3,16 @@ !>@brief !>mpi routines for shallow water model module mpi_module - use nrtype + use numerics_type use mpi + implicit none +#if VAR_TYPE==0 + integer(i4b), parameter :: MPIREAL=MPI_REAL4 +#endif +#if VAR_TYPE==1 + integer(i4b), parameter :: MPIREAL=MPI_REAL8 +#endif private public :: mpi_define, block_ring, exchange_halos @@ -43,7 +50,7 @@ end subroutine mpi_define subroutine exchange_halos(comm2d, id, ipp, jpp, o_halo, array) implicit none integer(i4b), intent(in) :: comm2d, id, ipp, jpp, o_halo - real(sp), intent(inout), & + real(wp), intent(inout), & dimension(1-o_halo:o_halo+ipp,1-o_halo:o_halo+jpp) :: array integer(i4b) :: error, nbrleft, nbrright, nbrbottom, nbrtop, tag1, & @@ -61,10 +68,10 @@ subroutine exchange_halos(comm2d, id, ipp, jpp, o_halo, array) if (nbrleft /= id) then tag1=010 ! send from left (specify destination): - call MPI_Issend(array(ipp+1-o_halo:ipp,1:jpp), jpp, MPI_REAL8, nbrright, & + call MPI_Issend(array(ipp+1-o_halo:ipp,1:jpp), jpp, MPIREAL, nbrright, & tag1, MPI_COMM_WORLD, request,error) ! receive from left (specify source): - call MPI_Recv(array(1-o_halo:0,1:jpp), jpp, MPI_REAL8, nbrleft, & + call MPI_Recv(array(1-o_halo:0,1:jpp), jpp, MPIREAL, nbrleft, & tag1, MPI_COMM_WORLD, status,error) call MPI_Wait(request, status, error) else @@ -76,10 +83,10 @@ subroutine exchange_halos(comm2d, id, ipp, jpp, o_halo, array) if (nbrright /= id) then tag1=010 ! send from right (specify destination): - call MPI_Issend(array(1:o_halo,1:jpp), jpp, MPI_REAL8, nbrleft, & + call MPI_Issend(array(1:o_halo,1:jpp), jpp, MPIREAL, nbrleft, & tag1, MPI_COMM_WORLD, request,error) ! receive from right (specify source): - call MPI_Recv(array(ipp+1:ipp+o_halo,1:jpp), jpp, MPI_REAL8, nbrright, & + call MPI_Recv(array(ipp+1:ipp+o_halo,1:jpp), jpp, MPIREAL, nbrright, & tag1, MPI_COMM_WORLD, status,error) call MPI_Wait(request, status, error) else @@ -90,10 +97,10 @@ subroutine exchange_halos(comm2d, id, ipp, jpp, o_halo, array) if ((nbrtop /= id)) then tag1=110 ! send from top (specify destination): - call MPI_Issend(array(1:ipp,jpp+1-o_halo:jpp), ipp, MPI_REAL8, nbrtop, & + call MPI_Issend(array(1:ipp,jpp+1-o_halo:jpp), ipp, MPIREAL, nbrtop, & tag1, MPI_COMM_WORLD, request,error) ! receive from top (specify source): - call MPI_Recv(array(1:ipp,1-o_halo:0), ipp, MPI_REAL8, nbrbottom, & + call MPI_Recv(array(1:ipp,1-o_halo:0), ipp, MPIREAL, nbrbottom, & tag1, MPI_COMM_WORLD, status,error) call MPI_Wait(request, status, error) else @@ -103,10 +110,10 @@ subroutine exchange_halos(comm2d, id, ipp, jpp, o_halo, array) if ((nbrbottom /= id)) then tag1=110 ! send from bottom (specify destination): - call MPI_Issend(array(1:ipp,1:o_halo), ipp, MPI_REAL8, nbrbottom, & + call MPI_Issend(array(1:ipp,1:o_halo), ipp, MPIREAL, nbrbottom, & tag1, MPI_COMM_WORLD, request,error) ! receive from bottom (specify source): - call MPI_Recv(array(1:ipp,jpp+1:jpp+o_halo), ipp, MPI_REAL8, nbrtop, & + call MPI_Recv(array(1:ipp,jpp+1:jpp+o_halo), ipp, MPIREAL, nbrtop, & tag1, MPI_COMM_WORLD, status,error) call MPI_Wait(request, status, error) else diff --git a/nr.f90 b/nr.f90 deleted file mode 100755 index 37dd921..0000000 --- a/nr.f90 +++ /dev/null @@ -1,3205 +0,0 @@ -MODULE nr - INTERFACE - SUBROUTINE airy(x,ai,bi,aip,bip) - USE nrtype - REAL(SP), INTENT(IN) :: x - REAL(SP), INTENT(OUT) :: ai,bi,aip,bip - END SUBROUTINE airy - END INTERFACE - INTERFACE - SUBROUTINE amebsa(p,y,pb,yb,ftol,func,iter,temptr) - USE nrtype - INTEGER(I4B), INTENT(INOUT) :: iter - REAL(SP), INTENT(INOUT) :: yb - REAL(SP), INTENT(IN) :: ftol,temptr - REAL(SP), DIMENSION(:), INTENT(INOUT) :: y,pb - REAL(SP), DIMENSION(:,:), INTENT(INOUT) :: p - INTERFACE - FUNCTION func(x) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(IN) :: x - REAL(SP) :: func - END FUNCTION func - END INTERFACE - END SUBROUTINE amebsa - END INTERFACE - INTERFACE - SUBROUTINE amoeba(p,y,ftol,func,iter) - USE nrtype - INTEGER(I4B), INTENT(OUT) :: iter - REAL(SP), INTENT(IN) :: ftol - REAL(SP), DIMENSION(:), INTENT(INOUT) :: y - REAL(SP), DIMENSION(:,:), INTENT(INOUT) :: p - INTERFACE - FUNCTION func(x) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(IN) :: x - REAL(SP) :: func - END FUNCTION func - END INTERFACE - END SUBROUTINE amoeba - END INTERFACE - INTERFACE - SUBROUTINE anneal(x,y,iorder) - USE nrtype - INTEGER(I4B), DIMENSION(:), INTENT(INOUT) :: iorder - REAL(SP), DIMENSION(:), INTENT(IN) :: x,y - END SUBROUTINE anneal - END INTERFACE - INTERFACE - SUBROUTINE asolve(b,x,itrnsp) - USE nrtype - REAL(DP), DIMENSION(:), INTENT(IN) :: b - REAL(DP), DIMENSION(:), INTENT(OUT) :: x - INTEGER(I4B), INTENT(IN) :: itrnsp - END SUBROUTINE asolve - END INTERFACE - INTERFACE - SUBROUTINE atimes(x,r,itrnsp) - USE nrtype - REAL(DP), DIMENSION(:), INTENT(IN) :: x - REAL(DP), DIMENSION(:), INTENT(OUT) :: r - INTEGER(I4B), INTENT(IN) :: itrnsp - END SUBROUTINE atimes - END INTERFACE - INTERFACE - SUBROUTINE avevar(data,ave,var) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(IN) :: data - REAL(SP), INTENT(OUT) :: ave,var - END SUBROUTINE avevar - END INTERFACE - INTERFACE - SUBROUTINE balanc(a) - USE nrtype - REAL(SP), DIMENSION(:,:), INTENT(INOUT) :: a - END SUBROUTINE balanc - END INTERFACE - INTERFACE - SUBROUTINE banbks(a,m1,m2,al,indx,b) - USE nrtype - INTEGER(I4B), INTENT(IN) :: m1,m2 - INTEGER(I4B), DIMENSION(:), INTENT(IN) :: indx - REAL(SP), DIMENSION(:,:), INTENT(IN) :: a,al - REAL(SP), DIMENSION(:), INTENT(INOUT) :: b - END SUBROUTINE banbks - END INTERFACE - INTERFACE - SUBROUTINE bandec(a,m1,m2,al,indx,d) - USE nrtype - INTEGER(I4B), INTENT(IN) :: m1,m2 - INTEGER(I4B), DIMENSION(:), INTENT(OUT) :: indx - REAL(SP), INTENT(OUT) :: d - REAL(SP), DIMENSION(:,:), INTENT(INOUT) :: a - REAL(SP), DIMENSION(:,:), INTENT(OUT) :: al - END SUBROUTINE bandec - END INTERFACE - INTERFACE - SUBROUTINE banmul(a,m1,m2,x,b) - USE nrtype - INTEGER(I4B), INTENT(IN) :: m1,m2 - REAL(SP), DIMENSION(:), INTENT(IN) :: x - REAL(SP), DIMENSION(:), INTENT(OUT) :: b - REAL(SP), DIMENSION(:,:), INTENT(IN) :: a - END SUBROUTINE banmul - END INTERFACE - INTERFACE - SUBROUTINE bcucof(y,y1,y2,y12,d1,d2,c) - USE nrtype - REAL(SP), INTENT(IN) :: d1,d2 - REAL(SP), DIMENSION(4), INTENT(IN) :: y,y1,y2,y12 - REAL(SP), DIMENSION(4,4), INTENT(OUT) :: c - END SUBROUTINE bcucof - END INTERFACE - INTERFACE - SUBROUTINE bcuint(y,y1,y2,y12,x1l,x1u,x2l,x2u,x1,x2,ansy,& - ansy1,ansy2) - USE nrtype - REAL(SP), DIMENSION(4), INTENT(IN) :: y,y1,y2,y12 - REAL(SP), INTENT(IN) :: x1l,x1u,x2l,x2u,x1,x2 - REAL(SP), INTENT(OUT) :: ansy,ansy1,ansy2 - END SUBROUTINE bcuint - END INTERFACE - INTERFACE beschb - SUBROUTINE beschb_s(x,gam1,gam2,gampl,gammi) - USE nrtype - REAL(DP), INTENT(IN) :: x - REAL(DP), INTENT(OUT) :: gam1,gam2,gampl,gammi - END SUBROUTINE beschb_s -!BL - SUBROUTINE beschb_v(x,gam1,gam2,gampl,gammi) - USE nrtype - REAL(DP), DIMENSION(:), INTENT(IN) :: x - REAL(DP), DIMENSION(:), INTENT(OUT) :: gam1,gam2,gampl,gammi - END SUBROUTINE beschb_v - END INTERFACE - INTERFACE bessi - FUNCTION bessi_s(n,x) - USE nrtype - INTEGER(I4B), INTENT(IN) :: n - REAL(SP), INTENT(IN) :: x - REAL(SP) :: bessi_s - END FUNCTION bessi_s -!BL - FUNCTION bessi_v(n,x) - USE nrtype - INTEGER(I4B), INTENT(IN) :: n - REAL(SP), DIMENSION(:), INTENT(IN) :: x - REAL(SP), DIMENSION(size(x)) :: bessi_v - END FUNCTION bessi_v - END INTERFACE - INTERFACE bessi0 - FUNCTION bessi0_s(x) - USE nrtype - REAL(SP), INTENT(IN) :: x - REAL(SP) :: bessi0_s - END FUNCTION bessi0_s -!BL - FUNCTION bessi0_v(x) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(IN) :: x - REAL(SP), DIMENSION(size(x)) :: bessi0_v - END FUNCTION bessi0_v - END INTERFACE - INTERFACE bessi1 - FUNCTION bessi1_s(x) - USE nrtype - REAL(SP), INTENT(IN) :: x - REAL(SP) :: bessi1_s - END FUNCTION bessi1_s -!BL - FUNCTION bessi1_v(x) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(IN) :: x - REAL(SP), DIMENSION(size(x)) :: bessi1_v - END FUNCTION bessi1_v - END INTERFACE - INTERFACE - SUBROUTINE bessik(x,xnu,ri,rk,rip,rkp) - USE nrtype - REAL(SP), INTENT(IN) :: x,xnu - REAL(SP), INTENT(OUT) :: ri,rk,rip,rkp - END SUBROUTINE bessik - END INTERFACE - INTERFACE bessj - FUNCTION bessj_s(n,x) - USE nrtype - INTEGER(I4B), INTENT(IN) :: n - REAL(SP), INTENT(IN) :: x - REAL(SP) :: bessj_s - END FUNCTION bessj_s -!BL - FUNCTION bessj_v(n,x) - USE nrtype - INTEGER(I4B), INTENT(IN) :: n - REAL(SP), DIMENSION(:), INTENT(IN) :: x - REAL(SP), DIMENSION(size(x)) :: bessj_v - END FUNCTION bessj_v - END INTERFACE - INTERFACE bessj0 - FUNCTION bessj0_s(x) - USE nrtype - REAL(SP), INTENT(IN) :: x - REAL(SP) :: bessj0_s - END FUNCTION bessj0_s -!BL - FUNCTION bessj0_v(x) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(IN) :: x - REAL(SP), DIMENSION(size(x)) :: bessj0_v - END FUNCTION bessj0_v - END INTERFACE - INTERFACE bessj1 - FUNCTION bessj1_s(x) - USE nrtype - REAL(SP), INTENT(IN) :: x - REAL(SP) :: bessj1_s - END FUNCTION bessj1_s -!BL - FUNCTION bessj1_v(x) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(IN) :: x - REAL(SP), DIMENSION(size(x)) :: bessj1_v - END FUNCTION bessj1_v - END INTERFACE - INTERFACE bessjy - SUBROUTINE bessjy_s(x,xnu,rj,ry,rjp,ryp) - USE nrtype - REAL(SP), INTENT(IN) :: x,xnu - REAL(SP), INTENT(OUT) :: rj,ry,rjp,ryp - END SUBROUTINE bessjy_s -!BL - SUBROUTINE bessjy_v(x,xnu,rj,ry,rjp,ryp) - USE nrtype - REAL(SP), INTENT(IN) :: xnu - REAL(SP), DIMENSION(:), INTENT(IN) :: x - REAL(SP), DIMENSION(:), INTENT(OUT) :: rj,rjp,ry,ryp - END SUBROUTINE bessjy_v - END INTERFACE - INTERFACE bessk - FUNCTION bessk_s(n,x) - USE nrtype - INTEGER(I4B), INTENT(IN) :: n - REAL(SP), INTENT(IN) :: x - REAL(SP) :: bessk_s - END FUNCTION bessk_s -!BL - FUNCTION bessk_v(n,x) - USE nrtype - INTEGER(I4B), INTENT(IN) :: n - REAL(SP), DIMENSION(:), INTENT(IN) :: x - REAL(SP), DIMENSION(size(x)) :: bessk_v - END FUNCTION bessk_v - END INTERFACE - INTERFACE bessk0 - FUNCTION bessk0_s(x) - USE nrtype - REAL(SP), INTENT(IN) :: x - REAL(SP) :: bessk0_s - END FUNCTION bessk0_s -!BL - FUNCTION bessk0_v(x) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(IN) :: x - REAL(SP), DIMENSION(size(x)) :: bessk0_v - END FUNCTION bessk0_v - END INTERFACE - INTERFACE bessk1 - FUNCTION bessk1_s(x) - USE nrtype - REAL(SP), INTENT(IN) :: x - REAL(SP) :: bessk1_s - END FUNCTION bessk1_s -!BL - FUNCTION bessk1_v(x) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(IN) :: x - REAL(SP), DIMENSION(size(x)) :: bessk1_v - END FUNCTION bessk1_v - END INTERFACE - INTERFACE bessy - FUNCTION bessy_s(n,x) - USE nrtype - INTEGER(I4B), INTENT(IN) :: n - REAL(SP), INTENT(IN) :: x - REAL(SP) :: bessy_s - END FUNCTION bessy_s -!BL - FUNCTION bessy_v(n,x) - USE nrtype - INTEGER(I4B), INTENT(IN) :: n - REAL(SP), DIMENSION(:), INTENT(IN) :: x - REAL(SP), DIMENSION(size(x)) :: bessy_v - END FUNCTION bessy_v - END INTERFACE - INTERFACE bessy0 - FUNCTION bessy0_s(x) - USE nrtype - REAL(SP), INTENT(IN) :: x - REAL(SP) :: bessy0_s - END FUNCTION bessy0_s -!BL - FUNCTION bessy0_v(x) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(IN) :: x - REAL(SP), DIMENSION(size(x)) :: bessy0_v - END FUNCTION bessy0_v - END INTERFACE - INTERFACE bessy1 - FUNCTION bessy1_s(x) - USE nrtype - REAL(SP), INTENT(IN) :: x - REAL(SP) :: bessy1_s - END FUNCTION bessy1_s -!BL - FUNCTION bessy1_v(x) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(IN) :: x - REAL(SP), DIMENSION(size(x)) :: bessy1_v - END FUNCTION bessy1_v - END INTERFACE - INTERFACE beta - FUNCTION beta_s(z,w) - USE nrtype - REAL(SP), INTENT(IN) :: z,w - REAL(SP) :: beta_s - END FUNCTION beta_s -!BL - FUNCTION beta_v(z,w) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(IN) :: z,w - REAL(SP), DIMENSION(size(z)) :: beta_v - END FUNCTION beta_v - END INTERFACE - INTERFACE betacf - FUNCTION betacf_s(a,b,x) - USE nrtype - REAL(SP), INTENT(IN) :: a,b,x - REAL(SP) :: betacf_s - END FUNCTION betacf_s -!BL - FUNCTION betacf_v(a,b,x) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(IN) :: a,b,x - REAL(SP), DIMENSION(size(x)) :: betacf_v - END FUNCTION betacf_v - END INTERFACE - INTERFACE betai - FUNCTION betai_s(a,b,x) - USE nrtype - REAL(SP), INTENT(IN) :: a,b,x - REAL(SP) :: betai_s - END FUNCTION betai_s -!BL - FUNCTION betai_v(a,b,x) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(IN) :: a,b,x - REAL(SP), DIMENSION(size(a)) :: betai_v - END FUNCTION betai_v - END INTERFACE - INTERFACE bico - FUNCTION bico_s(n,k) - USE nrtype - INTEGER(I4B), INTENT(IN) :: n,k - REAL(SP) :: bico_s - END FUNCTION bico_s -!BL - FUNCTION bico_v(n,k) - USE nrtype - INTEGER(I4B), DIMENSION(:), INTENT(IN) :: n,k - REAL(SP), DIMENSION(size(n)) :: bico_v - END FUNCTION bico_v - END INTERFACE - INTERFACE - FUNCTION bnldev(pp,n) - USE nrtype - REAL(SP), INTENT(IN) :: pp - INTEGER(I4B), INTENT(IN) :: n - REAL(SP) :: bnldev - END FUNCTION bnldev - END INTERFACE - INTERFACE - FUNCTION brent(ax,bx,cx,func,tol,xmin) - USE nrtype - REAL(SP), INTENT(IN) :: ax,bx,cx,tol - REAL(SP), INTENT(OUT) :: xmin - REAL(SP) :: brent - INTERFACE - FUNCTION func(x) - USE nrtype - REAL(SP), INTENT(IN) :: x - REAL(SP) :: func - END FUNCTION func - END INTERFACE - END FUNCTION brent - END INTERFACE - INTERFACE - SUBROUTINE broydn(x,check) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(INOUT) :: x - LOGICAL(LGT), INTENT(OUT) :: check - END SUBROUTINE broydn - END INTERFACE - INTERFACE - SUBROUTINE bsstep(y,dydx,x,htry,eps,yscal,hdid,hnext,derivs) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(INOUT) :: y - REAL(SP), DIMENSION(:), INTENT(IN) :: dydx,yscal - REAL(SP), INTENT(INOUT) :: x - REAL(SP), INTENT(IN) :: htry,eps - REAL(SP), INTENT(OUT) :: hdid,hnext - INTERFACE - SUBROUTINE derivs(x,y,dydx) - USE nrtype - REAL(SP), INTENT(IN) :: x - REAL(SP), DIMENSION(:), INTENT(IN) :: y - REAL(SP), DIMENSION(:), INTENT(OUT) :: dydx - END SUBROUTINE derivs - END INTERFACE - END SUBROUTINE bsstep - END INTERFACE - INTERFACE - SUBROUTINE caldat(julian,mm,id,iyyy) - USE nrtype - INTEGER(I4B), INTENT(IN) :: julian - INTEGER(I4B), INTENT(OUT) :: mm,id,iyyy - END SUBROUTINE caldat - END INTERFACE - INTERFACE - FUNCTION chder(a,b,c) - USE nrtype - REAL(SP), INTENT(IN) :: a,b - REAL(SP), DIMENSION(:), INTENT(IN) :: c - REAL(SP), DIMENSION(size(c)) :: chder - END FUNCTION chder - END INTERFACE - INTERFACE chebev - FUNCTION chebev_s(a,b,c,x) - USE nrtype - REAL(SP), INTENT(IN) :: a,b,x - REAL(SP), DIMENSION(:), INTENT(IN) :: c - REAL(SP) :: chebev_s - END FUNCTION chebev_s -!BL - FUNCTION chebev_v(a,b,c,x) - USE nrtype - REAL(SP), INTENT(IN) :: a,b - REAL(SP), DIMENSION(:), INTENT(IN) :: c,x - REAL(SP), DIMENSION(size(x)) :: chebev_v - END FUNCTION chebev_v - END INTERFACE - INTERFACE - FUNCTION chebft(a,b,n,func) - USE nrtype - REAL(SP), INTENT(IN) :: a,b - INTEGER(I4B), INTENT(IN) :: n - REAL(SP), DIMENSION(n) :: chebft - INTERFACE - FUNCTION func(x) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(IN) :: x - REAL(SP), DIMENSION(size(x)) :: func - END FUNCTION func - END INTERFACE - END FUNCTION chebft - END INTERFACE - INTERFACE - FUNCTION chebpc(c) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(IN) :: c - REAL(SP), DIMENSION(size(c)) :: chebpc - END FUNCTION chebpc - END INTERFACE - INTERFACE - FUNCTION chint(a,b,c) - USE nrtype - REAL(SP), INTENT(IN) :: a,b - REAL(SP), DIMENSION(:), INTENT(IN) :: c - REAL(SP), DIMENSION(size(c)) :: chint - END FUNCTION chint - END INTERFACE - INTERFACE - SUBROUTINE choldc(a,p) - USE nrtype - REAL(SP), DIMENSION(:,:), INTENT(INOUT) :: a - REAL(SP), DIMENSION(:), INTENT(OUT) :: p - END SUBROUTINE choldc - END INTERFACE - INTERFACE - SUBROUTINE cholsl(a,p,b,x) - USE nrtype - REAL(SP), DIMENSION(:,:), INTENT(IN) :: a - REAL(SP), DIMENSION(:), INTENT(IN) :: p,b - REAL(SP), DIMENSION(:), INTENT(INOUT) :: x - END SUBROUTINE cholsl - END INTERFACE - INTERFACE - SUBROUTINE chsone(bins,ebins,knstrn,df,chsq,prob) - USE nrtype - INTEGER(I4B), INTENT(IN) :: knstrn - REAL(SP), INTENT(OUT) :: df,chsq,prob - REAL(SP), DIMENSION(:), INTENT(IN) :: bins,ebins - END SUBROUTINE chsone - END INTERFACE - INTERFACE - SUBROUTINE chstwo(bins1,bins2,knstrn,df,chsq,prob) - USE nrtype - INTEGER(I4B), INTENT(IN) :: knstrn - REAL(SP), INTENT(OUT) :: df,chsq,prob - REAL(SP), DIMENSION(:), INTENT(IN) :: bins1,bins2 - END SUBROUTINE chstwo - END INTERFACE - INTERFACE - SUBROUTINE cisi(x,ci,si) - USE nrtype - REAL(SP), INTENT(IN) :: x - REAL(SP), INTENT(OUT) :: ci,si - END SUBROUTINE cisi - END INTERFACE - INTERFACE - SUBROUTINE cntab1(nn,chisq,df,prob,cramrv,ccc) - USE nrtype - INTEGER(I4B), DIMENSION(:,:), INTENT(IN) :: nn - REAL(SP), INTENT(OUT) :: chisq,df,prob,cramrv,ccc - END SUBROUTINE cntab1 - END INTERFACE - INTERFACE - SUBROUTINE cntab2(nn,h,hx,hy,hygx,hxgy,uygx,uxgy,uxy) - USE nrtype - INTEGER(I4B), DIMENSION(:,:), INTENT(IN) :: nn - REAL(SP), INTENT(OUT) :: h,hx,hy,hygx,hxgy,uygx,uxgy,uxy - END SUBROUTINE cntab2 - END INTERFACE - INTERFACE - FUNCTION convlv(data,respns,isign) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(IN) :: data - REAL(SP), DIMENSION(:), INTENT(IN) :: respns - INTEGER(I4B), INTENT(IN) :: isign - REAL(SP), DIMENSION(size(data)) :: convlv - END FUNCTION convlv - END INTERFACE - INTERFACE - FUNCTION correl(data1,data2) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(IN) :: data1,data2 - REAL(SP), DIMENSION(size(data1)) :: correl - END FUNCTION correl - END INTERFACE - INTERFACE - SUBROUTINE cosft1(y) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(INOUT) :: y - END SUBROUTINE cosft1 - END INTERFACE - INTERFACE - SUBROUTINE cosft2(y,isign) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(INOUT) :: y - INTEGER(I4B), INTENT(IN) :: isign - END SUBROUTINE cosft2 - END INTERFACE - INTERFACE - SUBROUTINE covsrt(covar,maska) - USE nrtype - REAL(SP), DIMENSION(:,:), INTENT(INOUT) :: covar - LOGICAL(LGT), DIMENSION(:), INTENT(IN) :: maska - END SUBROUTINE covsrt - END INTERFACE - INTERFACE - SUBROUTINE cyclic(a,b,c,alpha,beta,r,x) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(IN):: a,b,c,r - REAL(SP), INTENT(IN) :: alpha,beta - REAL(SP), DIMENSION(:), INTENT(OUT):: x - END SUBROUTINE cyclic - END INTERFACE - INTERFACE - SUBROUTINE daub4(a,isign) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(INOUT) :: a - INTEGER(I4B), INTENT(IN) :: isign - END SUBROUTINE daub4 - END INTERFACE - INTERFACE dawson - FUNCTION dawson_s(x) - USE nrtype - REAL(SP), INTENT(IN) :: x - REAL(SP) :: dawson_s - END FUNCTION dawson_s -!BL - FUNCTION dawson_v(x) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(IN) :: x - REAL(SP), DIMENSION(size(x)) :: dawson_v - END FUNCTION dawson_v - END INTERFACE - INTERFACE - FUNCTION dbrent(ax,bx,cx,func,dbrent_dfunc,tol,xmin) - USE nrtype - REAL(SP), INTENT(IN) :: ax,bx,cx,tol - REAL(SP), INTENT(OUT) :: xmin - REAL(SP) :: dbrent - INTERFACE - FUNCTION func(x) - USE nrtype - REAL(SP), INTENT(IN) :: x - REAL(SP) :: func - END FUNCTION func -!BL - FUNCTION dbrent_dfunc(x) - USE nrtype - REAL(SP), INTENT(IN) :: x - REAL(SP) :: dbrent_dfunc - END FUNCTION dbrent_dfunc - END INTERFACE - END FUNCTION dbrent - END INTERFACE - INTERFACE - SUBROUTINE ddpoly(c,x,pd) - USE nrtype - REAL(SP), INTENT(IN) :: x - REAL(SP), DIMENSION(:), INTENT(IN) :: c - REAL(SP), DIMENSION(:), INTENT(OUT) :: pd - END SUBROUTINE ddpoly - END INTERFACE - INTERFACE - FUNCTION decchk(string,ch) - USE nrtype - CHARACTER(1), DIMENSION(:), INTENT(IN) :: string - CHARACTER(1), INTENT(OUT) :: ch - LOGICAL(LGT) :: decchk - END FUNCTION decchk - END INTERFACE - INTERFACE - SUBROUTINE dfpmin(p,gtol,iter,fret,func,dfunc) - USE nrtype - INTEGER(I4B), INTENT(OUT) :: iter - REAL(SP), INTENT(IN) :: gtol - REAL(SP), INTENT(OUT) :: fret - REAL(SP), DIMENSION(:), INTENT(INOUT) :: p - INTERFACE - FUNCTION func(p) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(IN) :: p - REAL(SP) :: func - END FUNCTION func -!BL - FUNCTION dfunc(p) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(IN) :: p - REAL(SP), DIMENSION(size(p)) :: dfunc - END FUNCTION dfunc - END INTERFACE - END SUBROUTINE dfpmin - END INTERFACE - INTERFACE - FUNCTION dfridr(func,x,h,err) - USE nrtype - REAL(SP), INTENT(IN) :: x,h - REAL(SP), INTENT(OUT) :: err - REAL(SP) :: dfridr - INTERFACE - FUNCTION func(x) - USE nrtype - REAL(SP), INTENT(IN) :: x - REAL(SP) :: func - END FUNCTION func - END INTERFACE - END FUNCTION dfridr - END INTERFACE - INTERFACE - SUBROUTINE dftcor(w,delta,a,b,endpts,corre,corim,corfac) - USE nrtype - REAL(SP), INTENT(IN) :: w,delta,a,b - REAL(SP), INTENT(OUT) :: corre,corim,corfac - REAL(SP), DIMENSION(:), INTENT(IN) :: endpts - END SUBROUTINE dftcor - END INTERFACE - INTERFACE - SUBROUTINE dftint(func,a,b,w,cosint,sinint) - USE nrtype - REAL(SP), INTENT(IN) :: a,b,w - REAL(SP), INTENT(OUT) :: cosint,sinint - INTERFACE - FUNCTION func(x) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(IN) :: x - REAL(SP), DIMENSION(size(x)) :: func - END FUNCTION func - END INTERFACE - END SUBROUTINE dftint - END INTERFACE - INTERFACE - SUBROUTINE difeq(k,k1,k2,jsf,is1,isf,indexv,s,y) - USE nrtype - INTEGER(I4B), INTENT(IN) :: is1,isf,jsf,k,k1,k2 - INTEGER(I4B), DIMENSION(:), INTENT(IN) :: indexv - REAL(SP), DIMENSION(:,:), INTENT(OUT) :: s - REAL(SP), DIMENSION(:,:), INTENT(IN) :: y - END SUBROUTINE difeq - END INTERFACE - INTERFACE - FUNCTION eclass(lista,listb,n) - USE nrtype - INTEGER(I4B), DIMENSION(:), INTENT(IN) :: lista,listb - INTEGER(I4B), INTENT(IN) :: n - INTEGER(I4B), DIMENSION(n) :: eclass - END FUNCTION eclass - END INTERFACE - INTERFACE - FUNCTION eclazz(equiv,n) - USE nrtype - INTERFACE - FUNCTION equiv(i,j) - USE nrtype - LOGICAL(LGT) :: equiv - INTEGER(I4B), INTENT(IN) :: i,j - END FUNCTION equiv - END INTERFACE - INTEGER(I4B), INTENT(IN) :: n - INTEGER(I4B), DIMENSION(n) :: eclazz - END FUNCTION eclazz - END INTERFACE - INTERFACE - FUNCTION ei(x) - USE nrtype - REAL(SP), INTENT(IN) :: x - REAL(SP) :: ei - END FUNCTION ei - END INTERFACE - INTERFACE - SUBROUTINE eigsrt(d,v) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(INOUT) :: d - REAL(SP), DIMENSION(:,:), INTENT(INOUT) :: v - END SUBROUTINE eigsrt - END INTERFACE - INTERFACE elle - FUNCTION elle_s(phi,ak) - USE nrtype - REAL(SP), INTENT(IN) :: phi,ak - REAL(SP) :: elle_s - END FUNCTION elle_s -!BL - FUNCTION elle_v(phi,ak) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(IN) :: phi,ak - REAL(SP), DIMENSION(size(phi)) :: elle_v - END FUNCTION elle_v - END INTERFACE - INTERFACE ellf - FUNCTION ellf_s(phi,ak) - USE nrtype - REAL(SP), INTENT(IN) :: phi,ak - REAL(SP) :: ellf_s - END FUNCTION ellf_s -!BL - FUNCTION ellf_v(phi,ak) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(IN) :: phi,ak - REAL(SP), DIMENSION(size(phi)) :: ellf_v - END FUNCTION ellf_v - END INTERFACE - INTERFACE ellpi - FUNCTION ellpi_s(phi,en,ak) - USE nrtype - REAL(SP), INTENT(IN) :: phi,en,ak - REAL(SP) :: ellpi_s - END FUNCTION ellpi_s -!BL - FUNCTION ellpi_v(phi,en,ak) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(IN) :: phi,en,ak - REAL(SP), DIMENSION(size(phi)) :: ellpi_v - END FUNCTION ellpi_v - END INTERFACE - INTERFACE - SUBROUTINE elmhes(a) - USE nrtype - REAL(SP), DIMENSION(:,:), INTENT(INOUT) :: a - END SUBROUTINE elmhes - END INTERFACE - INTERFACE erf - FUNCTION erf_s(x) - USE nrtype - REAL(SP), INTENT(IN) :: x - REAL(SP) :: erf_s - END FUNCTION erf_s -!BL - FUNCTION erf_v(x) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(IN) :: x - REAL(SP), DIMENSION(size(x)) :: erf_v - END FUNCTION erf_v - END INTERFACE - INTERFACE erfc - FUNCTION erfc_s(x) - USE nrtype - REAL(SP), INTENT(IN) :: x - REAL(SP) :: erfc_s - END FUNCTION erfc_s -!BL - FUNCTION erfc_v(x) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(IN) :: x - REAL(SP), DIMENSION(size(x)) :: erfc_v - END FUNCTION erfc_v - END INTERFACE - INTERFACE erfcc - FUNCTION erfcc_s(x) - USE nrtype - REAL(SP), INTENT(IN) :: x - REAL(SP) :: erfcc_s - END FUNCTION erfcc_s -!BL - FUNCTION erfcc_v(x) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(IN) :: x - REAL(SP), DIMENSION(size(x)) :: erfcc_v - END FUNCTION erfcc_v - END INTERFACE - INTERFACE - SUBROUTINE eulsum(sum,term,jterm) - USE nrtype - REAL(SP), INTENT(INOUT) :: sum - REAL(SP), INTENT(IN) :: term - INTEGER(I4B), INTENT(IN) :: jterm - END SUBROUTINE eulsum - END INTERFACE - INTERFACE - FUNCTION evlmem(fdt,d,xms) - USE nrtype - REAL(SP), INTENT(IN) :: fdt,xms - REAL(SP), DIMENSION(:), INTENT(IN) :: d - REAL(SP) :: evlmem - END FUNCTION evlmem - END INTERFACE - INTERFACE expdev - SUBROUTINE expdev_s(harvest) - USE nrtype - REAL(SP), INTENT(OUT) :: harvest - END SUBROUTINE expdev_s -!BL - SUBROUTINE expdev_v(harvest) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(OUT) :: harvest - END SUBROUTINE expdev_v - END INTERFACE - INTERFACE - FUNCTION expint(n,x) - USE nrtype - INTEGER(I4B), INTENT(IN) :: n - REAL(SP), INTENT(IN) :: x - REAL(SP) :: expint - END FUNCTION expint - END INTERFACE - INTERFACE factln - FUNCTION factln_s(n) - USE nrtype - INTEGER(I4B), INTENT(IN) :: n - REAL(SP) :: factln_s - END FUNCTION factln_s -!BL - FUNCTION factln_v(n) - USE nrtype - INTEGER(I4B), DIMENSION(:), INTENT(IN) :: n - REAL(SP), DIMENSION(size(n)) :: factln_v - END FUNCTION factln_v - END INTERFACE - INTERFACE factrl - FUNCTION factrl_s(n) - USE nrtype - INTEGER(I4B), INTENT(IN) :: n - REAL(SP) :: factrl_s - END FUNCTION factrl_s -!BL - FUNCTION factrl_v(n) - USE nrtype - INTEGER(I4B), DIMENSION(:), INTENT(IN) :: n - REAL(SP), DIMENSION(size(n)) :: factrl_v - END FUNCTION factrl_v - END INTERFACE - INTERFACE - SUBROUTINE fasper(x,y,ofac,hifac,px,py,jmax,prob) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(IN) :: x,y - REAL(SP), INTENT(IN) :: ofac,hifac - INTEGER(I4B), INTENT(OUT) :: jmax - REAL(SP), INTENT(OUT) :: prob - REAL(SP), DIMENSION(:), POINTER :: px,py - END SUBROUTINE fasper - END INTERFACE - INTERFACE - SUBROUTINE fdjac(x,fvec,df) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(IN) :: fvec - REAL(SP), DIMENSION(:), INTENT(INOUT) :: x - REAL(SP), DIMENSION(:,:), INTENT(OUT) :: df - END SUBROUTINE fdjac - END INTERFACE - INTERFACE - SUBROUTINE fgauss(x,a,y,dyda) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(IN) :: x,a - REAL(SP), DIMENSION(:), INTENT(OUT) :: y - REAL(SP), DIMENSION(:,:), INTENT(OUT) :: dyda - END SUBROUTINE fgauss - END INTERFACE - INTERFACE - SUBROUTINE fit(x,y,a,b,siga,sigb,chi2,q,sig) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(IN) :: x,y - REAL(SP), INTENT(OUT) :: a,b,siga,sigb,chi2,q - REAL(SP), DIMENSION(:), OPTIONAL, INTENT(IN) :: sig - END SUBROUTINE fit - END INTERFACE - INTERFACE - SUBROUTINE fitexy(x,y,sigx,sigy,a,b,siga,sigb,chi2,q) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(IN) :: x,y,sigx,sigy - REAL(SP), INTENT(OUT) :: a,b,siga,sigb,chi2,q - END SUBROUTINE fitexy - END INTERFACE - INTERFACE - SUBROUTINE fixrts(d) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(INOUT) :: d - END SUBROUTINE fixrts - END INTERFACE - INTERFACE - FUNCTION fleg(x,n) - USE nrtype - REAL(SP), INTENT(IN) :: x - INTEGER(I4B), INTENT(IN) :: n - REAL(SP), DIMENSION(n) :: fleg - END FUNCTION fleg - END INTERFACE - INTERFACE - SUBROUTINE flmoon(n,nph,jd,frac) - USE nrtype - INTEGER(I4B), INTENT(IN) :: n,nph - INTEGER(I4B), INTENT(OUT) :: jd - REAL(SP), INTENT(OUT) :: frac - END SUBROUTINE flmoon - END INTERFACE - INTERFACE four1 -! SUBROUTINE four1_dp(data,isign) -! USE nrtype -! COMPLEX(DPC), DIMENSION(:), INTENT(INOUT) :: data -! INTEGER(I4B), INTENT(IN) :: isign -! END SUBROUTINE four1_dp -!BL - SUBROUTINE four1_sp(data,isign) - USE nrtype - COMPLEX(SPC), DIMENSION(:), INTENT(INOUT) :: data - INTEGER(I4B), INTENT(IN) :: isign - END SUBROUTINE four1_sp - END INTERFACE - INTERFACE - SUBROUTINE four1_alt(data,isign) - USE nrtype - COMPLEX(SPC), DIMENSION(:), INTENT(INOUT) :: data - INTEGER(I4B), INTENT(IN) :: isign - END SUBROUTINE four1_alt - END INTERFACE - INTERFACE - SUBROUTINE four1_gather(data,isign) - USE nrtype - COMPLEX(SPC), DIMENSION(:), INTENT(INOUT) :: data - INTEGER(I4B), INTENT(IN) :: isign - END SUBROUTINE four1_gather - END INTERFACE - INTERFACE - SUBROUTINE four2(data,isign) - USE nrtype - COMPLEX(SPC), DIMENSION(:,:), INTENT(INOUT) :: data - INTEGER(I4B),INTENT(IN) :: isign - END SUBROUTINE four2 - END INTERFACE - INTERFACE - SUBROUTINE four2_alt(data,isign) - USE nrtype - COMPLEX(SPC), DIMENSION(:,:), INTENT(INOUT) :: data - INTEGER(I4B), INTENT(IN) :: isign - END SUBROUTINE four2_alt - END INTERFACE - INTERFACE - SUBROUTINE four3(data,isign) - USE nrtype - COMPLEX(SPC), DIMENSION(:,:,:), INTENT(INOUT) :: data - INTEGER(I4B),INTENT(IN) :: isign - END SUBROUTINE four3 - END INTERFACE - INTERFACE - SUBROUTINE four3_alt(data,isign) - USE nrtype - COMPLEX(SPC), DIMENSION(:,:,:), INTENT(INOUT) :: data - INTEGER(I4B), INTENT(IN) :: isign - END SUBROUTINE four3_alt - END INTERFACE - INTERFACE - SUBROUTINE fourcol(data,isign) - USE nrtype - COMPLEX(SPC), DIMENSION(:,:), INTENT(INOUT) :: data - INTEGER(I4B), INTENT(IN) :: isign - END SUBROUTINE fourcol - END INTERFACE - INTERFACE - SUBROUTINE fourcol_3d(data,isign) - USE nrtype - COMPLEX(SPC), DIMENSION(:,:,:), INTENT(INOUT) :: data - INTEGER(I4B), INTENT(IN) :: isign - END SUBROUTINE fourcol_3d - END INTERFACE - INTERFACE - SUBROUTINE fourn_gather(data,nn,isign) - USE nrtype - COMPLEX(SPC), DIMENSION(:), INTENT(INOUT) :: data - INTEGER(I4B), DIMENSION(:), INTENT(IN) :: nn - INTEGER(I4B), INTENT(IN) :: isign - END SUBROUTINE fourn_gather - END INTERFACE - INTERFACE fourrow -! SUBROUTINE fourrow_dp(data,isign) -! USE nrtype -! COMPLEX(DPC), DIMENSION(:,:), INTENT(INOUT) :: data -! INTEGER(I4B), INTENT(IN) :: isign -! END SUBROUTINE fourrow_dp -!BL - SUBROUTINE fourrow_sp(data,isign) - USE nrtype - COMPLEX(SPC), DIMENSION(:,:), INTENT(INOUT) :: data - INTEGER(I4B), INTENT(IN) :: isign - END SUBROUTINE fourrow_sp - END INTERFACE - INTERFACE - SUBROUTINE fourrow_3d(data,isign) - USE nrtype - COMPLEX(SPC), DIMENSION(:,:,:), INTENT(INOUT) :: data - INTEGER(I4B), INTENT(IN) :: isign - END SUBROUTINE fourrow_3d - END INTERFACE - INTERFACE - FUNCTION fpoly(x,n) - USE nrtype - REAL(SP), INTENT(IN) :: x - INTEGER(I4B), INTENT(IN) :: n - REAL(SP), DIMENSION(n) :: fpoly - END FUNCTION fpoly - END INTERFACE - INTERFACE - SUBROUTINE fred2(a,b,t,f,w,g,ak) - USE nrtype - REAL(SP), INTENT(IN) :: a,b - REAL(SP), DIMENSION(:), INTENT(OUT) :: t,f,w - INTERFACE - FUNCTION g(t) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(IN) :: t - REAL(SP), DIMENSION(size(t)) :: g - END FUNCTION g -!BL - FUNCTION ak(t,s) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(IN) :: t,s - REAL(SP), DIMENSION(size(t),size(s)) :: ak - END FUNCTION ak - END INTERFACE - END SUBROUTINE fred2 - END INTERFACE - INTERFACE - FUNCTION fredin(x,a,b,t,f,w,g,ak) - USE nrtype - REAL(SP), INTENT(IN) :: a,b - REAL(SP), DIMENSION(:), INTENT(IN) :: x,t,f,w - REAL(SP), DIMENSION(size(x)) :: fredin - INTERFACE - FUNCTION g(t) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(IN) :: t - REAL(SP), DIMENSION(size(t)) :: g - END FUNCTION g -!BL - FUNCTION ak(t,s) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(IN) :: t,s - REAL(SP), DIMENSION(size(t),size(s)) :: ak - END FUNCTION ak - END INTERFACE - END FUNCTION fredin - END INTERFACE - INTERFACE - SUBROUTINE frenel(x,s,c) - USE nrtype - REAL(SP), INTENT(IN) :: x - REAL(SP), INTENT(OUT) :: s,c - END SUBROUTINE frenel - END INTERFACE - INTERFACE - SUBROUTINE frprmn(p,ftol,iter,fret) - USE nrtype - INTEGER(I4B), INTENT(OUT) :: iter - REAL(SP), INTENT(IN) :: ftol - REAL(SP), INTENT(OUT) :: fret - REAL(SP), DIMENSION(:), INTENT(INOUT) :: p - END SUBROUTINE frprmn - END INTERFACE - INTERFACE - SUBROUTINE ftest(data1,data2,f,prob) - USE nrtype - REAL(SP), INTENT(OUT) :: f,prob - REAL(SP), DIMENSION(:), INTENT(IN) :: data1,data2 - END SUBROUTINE ftest - END INTERFACE - INTERFACE - FUNCTION gamdev(ia) - USE nrtype - INTEGER(I4B), INTENT(IN) :: ia - REAL(SP) :: gamdev - END FUNCTION gamdev - END INTERFACE - INTERFACE gammln - FUNCTION gammln_s(xx) - USE nrtype - REAL(SP), INTENT(IN) :: xx - REAL(SP) :: gammln_s - END FUNCTION gammln_s -!BL - FUNCTION gammln_v(xx) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(IN) :: xx - REAL(SP), DIMENSION(size(xx)) :: gammln_v - END FUNCTION gammln_v - END INTERFACE - INTERFACE gammp - FUNCTION gammp_s(a,x) - USE nrtype - REAL(SP), INTENT(IN) :: a,x - REAL(SP) :: gammp_s - END FUNCTION gammp_s -!BL - FUNCTION gammp_v(a,x) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(IN) :: a,x - REAL(SP), DIMENSION(size(a)) :: gammp_v - END FUNCTION gammp_v - END INTERFACE - INTERFACE gammq - FUNCTION gammq_s(a,x) - USE nrtype - REAL(SP), INTENT(IN) :: a,x - REAL(SP) :: gammq_s - END FUNCTION gammq_s -!BL - FUNCTION gammq_v(a,x) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(IN) :: a,x - REAL(SP), DIMENSION(size(a)) :: gammq_v - END FUNCTION gammq_v - END INTERFACE - INTERFACE gasdev - SUBROUTINE gasdev_s(harvest) - USE nrtype - REAL(SP), INTENT(OUT) :: harvest - END SUBROUTINE gasdev_s -!BL - SUBROUTINE gasdev_v(harvest) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(OUT) :: harvest - END SUBROUTINE gasdev_v - END INTERFACE - INTERFACE - SUBROUTINE gaucof(a,b,amu0,x,w) - USE nrtype - REAL(SP), INTENT(IN) :: amu0 - REAL(SP), DIMENSION(:), INTENT(INOUT) :: a,b - REAL(SP), DIMENSION(:), INTENT(OUT) :: x,w - END SUBROUTINE gaucof - END INTERFACE - INTERFACE - SUBROUTINE gauher(x,w) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(OUT) :: x,w - END SUBROUTINE gauher - END INTERFACE - INTERFACE - SUBROUTINE gaujac(x,w,alf,bet) - USE nrtype - REAL(SP), INTENT(IN) :: alf,bet - REAL(SP), DIMENSION(:), INTENT(OUT) :: x,w - END SUBROUTINE gaujac - END INTERFACE - INTERFACE - SUBROUTINE gaulag(x,w,alf) - USE nrtype - REAL(SP), INTENT(IN) :: alf - REAL(SP), DIMENSION(:), INTENT(OUT) :: x,w - END SUBROUTINE gaulag - END INTERFACE - INTERFACE - SUBROUTINE gauleg(x1,x2,x,w) - USE nrtype - REAL(SP), INTENT(IN) :: x1,x2 - REAL(SP), DIMENSION(:), INTENT(OUT) :: x,w - END SUBROUTINE gauleg - END INTERFACE - INTERFACE - SUBROUTINE gaussj(a,b) - USE nrtype - REAL(SP), DIMENSION(:,:), INTENT(INOUT) :: a,b - END SUBROUTINE gaussj - END INTERFACE - INTERFACE gcf - FUNCTION gcf_s(a,x,gln) - USE nrtype - REAL(SP), INTENT(IN) :: a,x - REAL(SP), OPTIONAL, INTENT(OUT) :: gln - REAL(SP) :: gcf_s - END FUNCTION gcf_s -!BL - FUNCTION gcf_v(a,x,gln) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(IN) :: a,x - REAL(SP), DIMENSION(:), OPTIONAL, INTENT(OUT) :: gln - REAL(SP), DIMENSION(size(a)) :: gcf_v - END FUNCTION gcf_v - END INTERFACE - INTERFACE - FUNCTION golden(ax,bx,cx,func,tol,xmin) - USE nrtype - REAL(SP), INTENT(IN) :: ax,bx,cx,tol - REAL(SP), INTENT(OUT) :: xmin - REAL(SP) :: golden - INTERFACE - FUNCTION func(x) - USE nrtype - REAL(SP), INTENT(IN) :: x - REAL(SP) :: func - END FUNCTION func - END INTERFACE - END FUNCTION golden - END INTERFACE - INTERFACE gser - FUNCTION gser_s(a,x,gln) - USE nrtype - REAL(SP), INTENT(IN) :: a,x - REAL(SP), OPTIONAL, INTENT(OUT) :: gln - REAL(SP) :: gser_s - END FUNCTION gser_s -!BL - FUNCTION gser_v(a,x,gln) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(IN) :: a,x - REAL(SP), DIMENSION(:), OPTIONAL, INTENT(OUT) :: gln - REAL(SP), DIMENSION(size(a)) :: gser_v - END FUNCTION gser_v - END INTERFACE - INTERFACE - SUBROUTINE hqr(a,wr,wi) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(OUT) :: wr,wi - REAL(SP), DIMENSION(:,:), INTENT(INOUT) :: a - END SUBROUTINE hqr - END INTERFACE - INTERFACE - SUBROUTINE hunt(xx,x,jlo) - USE nrtype - INTEGER(I4B), INTENT(INOUT) :: jlo - REAL(SP), INTENT(IN) :: x - REAL(SP), DIMENSION(:), INTENT(IN) :: xx - END SUBROUTINE hunt - END INTERFACE - INTERFACE - SUBROUTINE hypdrv(s,ry,rdyds) - USE nrtype - REAL(SP), INTENT(IN) :: s - REAL(SP), DIMENSION(:), INTENT(IN) :: ry - REAL(SP), DIMENSION(:), INTENT(OUT) :: rdyds - END SUBROUTINE hypdrv - END INTERFACE - INTERFACE - FUNCTION hypgeo(a,b,c,z) - USE nrtype - COMPLEX(SPC), INTENT(IN) :: a,b,c,z - COMPLEX(SPC) :: hypgeo - END FUNCTION hypgeo - END INTERFACE - INTERFACE - SUBROUTINE hypser(a,b,c,z,series,deriv) - USE nrtype - COMPLEX(SPC), INTENT(IN) :: a,b,c,z - COMPLEX(SPC), INTENT(OUT) :: series,deriv - END SUBROUTINE hypser - END INTERFACE - INTERFACE - FUNCTION icrc(crc,buf,jinit,jrev) - USE nrtype - CHARACTER(1), DIMENSION(:), INTENT(IN) :: buf - INTEGER(I2B), INTENT(IN) :: crc,jinit - INTEGER(I4B), INTENT(IN) :: jrev - INTEGER(I2B) :: icrc - END FUNCTION icrc - END INTERFACE - INTERFACE - FUNCTION igray(n,is) - USE nrtype - INTEGER(I4B), INTENT(IN) :: n,is - INTEGER(I4B) :: igray - END FUNCTION igray - END INTERFACE - INTERFACE - RECURSIVE SUBROUTINE index_bypack(arr,index,partial) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(IN) :: arr - INTEGER(I4B), DIMENSION(:), INTENT(INOUT) :: index - INTEGER, OPTIONAL, INTENT(IN) :: partial - END SUBROUTINE index_bypack - END INTERFACE - INTERFACE indexx - SUBROUTINE indexx_sp(arr,index) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(IN) :: arr - INTEGER(I4B), DIMENSION(:), INTENT(OUT) :: index - END SUBROUTINE indexx_sp - SUBROUTINE indexx_i4b(iarr,index) - USE nrtype - INTEGER(I4B), DIMENSION(:), INTENT(IN) :: iarr - INTEGER(I4B), DIMENSION(:), INTENT(OUT) :: index - END SUBROUTINE indexx_i4b - END INTERFACE - INTERFACE - FUNCTION interp(uc) - USE nrtype - REAL(DP), DIMENSION(:,:), INTENT(IN) :: uc - REAL(DP), DIMENSION(2*size(uc,1)-1,2*size(uc,1)-1) :: interp - END FUNCTION interp - END INTERFACE - INTERFACE - FUNCTION rank(indx) - USE nrtype - INTEGER(I4B), DIMENSION(:), INTENT(IN) :: indx - INTEGER(I4B), DIMENSION(size(indx)) :: rank - END FUNCTION rank - END INTERFACE - INTERFACE - FUNCTION irbit1(iseed) - USE nrtype - INTEGER(I4B), INTENT(INOUT) :: iseed - INTEGER(I4B) :: irbit1 - END FUNCTION irbit1 - END INTERFACE - INTERFACE - FUNCTION irbit2(iseed) - USE nrtype - INTEGER(I4B), INTENT(INOUT) :: iseed - INTEGER(I4B) :: irbit2 - END FUNCTION irbit2 - END INTERFACE - INTERFACE - SUBROUTINE jacobi(a,d,v,nrot) - USE nrtype - INTEGER(I4B), INTENT(OUT) :: nrot - REAL(SP), DIMENSION(:), INTENT(OUT) :: d - REAL(SP), DIMENSION(:,:), INTENT(INOUT) :: a - REAL(SP), DIMENSION(:,:), INTENT(OUT) :: v - END SUBROUTINE jacobi - END INTERFACE - INTERFACE - SUBROUTINE jacobn(x,y,dfdx,dfdy) - USE nrtype - REAL(SP), INTENT(IN) :: x - REAL(SP), DIMENSION(:), INTENT(IN) :: y - REAL(SP), DIMENSION(:), INTENT(OUT) :: dfdx - REAL(SP), DIMENSION(:,:), INTENT(OUT) :: dfdy - END SUBROUTINE jacobn - END INTERFACE - INTERFACE - FUNCTION julday(mm,id,iyyy) - USE nrtype - INTEGER(I4B), INTENT(IN) :: mm,id,iyyy - INTEGER(I4B) :: julday - END FUNCTION julday - END INTERFACE - INTERFACE - SUBROUTINE kendl1(data1,data2,tau,z,prob) - USE nrtype - REAL(SP), INTENT(OUT) :: tau,z,prob - REAL(SP), DIMENSION(:), INTENT(IN) :: data1,data2 - END SUBROUTINE kendl1 - END INTERFACE - INTERFACE - SUBROUTINE kendl2(tab,tau,z,prob) - USE nrtype - REAL(SP), DIMENSION(:,:), INTENT(IN) :: tab - REAL(SP), INTENT(OUT) :: tau,z,prob - END SUBROUTINE kendl2 - END INTERFACE - INTERFACE - FUNCTION kermom(y,m) - USE nrtype - REAL(DP), INTENT(IN) :: y - INTEGER(I4B), INTENT(IN) :: m - REAL(DP), DIMENSION(m) :: kermom - END FUNCTION kermom - END INTERFACE - INTERFACE - SUBROUTINE ks2d1s(x1,y1,quadvl,d1,prob) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(IN) :: x1,y1 - REAL(SP), INTENT(OUT) :: d1,prob - INTERFACE - SUBROUTINE quadvl(x,y,fa,fb,fc,fd) - USE nrtype - REAL(SP), INTENT(IN) :: x,y - REAL(SP), INTENT(OUT) :: fa,fb,fc,fd - END SUBROUTINE quadvl - END INTERFACE - END SUBROUTINE ks2d1s - END INTERFACE - INTERFACE - SUBROUTINE ks2d2s(x1,y1,x2,y2,d,prob) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(IN) :: x1,y1,x2,y2 - REAL(SP), INTENT(OUT) :: d,prob - END SUBROUTINE ks2d2s - END INTERFACE - INTERFACE - SUBROUTINE ksone(data,func,d,prob) - USE nrtype - REAL(SP), INTENT(OUT) :: d,prob - REAL(SP), DIMENSION(:), INTENT(INOUT) :: data - INTERFACE - FUNCTION func(x) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(IN) :: x - REAL(SP), DIMENSION(size(x)) :: func - END FUNCTION func - END INTERFACE - END SUBROUTINE ksone - END INTERFACE - INTERFACE - SUBROUTINE kstwo(data1,data2,d,prob) - USE nrtype - REAL(SP), INTENT(OUT) :: d,prob - REAL(SP), DIMENSION(:), INTENT(IN) :: data1,data2 - END SUBROUTINE kstwo - END INTERFACE - INTERFACE - SUBROUTINE laguer(a,x,its) - USE nrtype - INTEGER(I4B), INTENT(OUT) :: its - COMPLEX(SPC), INTENT(INOUT) :: x - COMPLEX(SPC), DIMENSION(:), INTENT(IN) :: a - END SUBROUTINE laguer - END INTERFACE - INTERFACE - SUBROUTINE lfit(x,y,sig,a,maska,covar,chisq,funcs) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(IN) :: x,y,sig - REAL(SP), DIMENSION(:), INTENT(INOUT) :: a - LOGICAL(LGT), DIMENSION(:), INTENT(IN) :: maska - REAL(SP), DIMENSION(:,:), INTENT(INOUT) :: covar - REAL(SP), INTENT(OUT) :: chisq - INTERFACE - SUBROUTINE funcs(x,arr) - USE nrtype - REAL(SP),INTENT(IN) :: x - REAL(SP), DIMENSION(:), INTENT(OUT) :: arr - END SUBROUTINE funcs - END INTERFACE - END SUBROUTINE lfit - END INTERFACE - INTERFACE - SUBROUTINE linbcg(b,x,itol,tol,itmax,iter,err) - USE nrtype - REAL(DP), DIMENSION(:), INTENT(IN) :: b - REAL(DP), DIMENSION(:), INTENT(INOUT) :: x - INTEGER(I4B), INTENT(IN) :: itol,itmax - REAL(DP), INTENT(IN) :: tol - INTEGER(I4B), INTENT(OUT) :: iter - REAL(DP), INTENT(OUT) :: err - END SUBROUTINE linbcg - END INTERFACE - INTERFACE - SUBROUTINE linmin(p,xi,fret) - USE nrtype - REAL(SP), INTENT(OUT) :: fret - REAL(SP), DIMENSION(:), TARGET, INTENT(INOUT) :: p,xi - END SUBROUTINE linmin - END INTERFACE - INTERFACE - SUBROUTINE lnsrch(xold,fold,g,p,x,f,stpmax,check,func) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(IN) :: xold,g - REAL(SP), DIMENSION(:), INTENT(INOUT) :: p - REAL(SP), INTENT(IN) :: fold,stpmax - REAL(SP), DIMENSION(:), INTENT(OUT) :: x - REAL(SP), INTENT(OUT) :: f - LOGICAL(LGT), INTENT(OUT) :: check - INTERFACE - FUNCTION func(x) - USE nrtype - REAL(SP) :: func - REAL(SP), DIMENSION(:), INTENT(IN) :: x - END FUNCTION func - END INTERFACE - END SUBROUTINE lnsrch - END INTERFACE - INTERFACE - FUNCTION locate(xx,x) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(IN) :: xx - REAL(SP), INTENT(IN) :: x - INTEGER(I4B) :: locate - END FUNCTION locate - END INTERFACE - INTERFACE - FUNCTION lop(u) - USE nrtype - REAL(DP), DIMENSION(:,:), INTENT(IN) :: u - REAL(DP), DIMENSION(size(u,1),size(u,1)) :: lop - END FUNCTION lop - END INTERFACE - INTERFACE - SUBROUTINE lubksb(a,indx,b) - USE nrtype - REAL(SP), DIMENSION(:,:), INTENT(IN) :: a - INTEGER(I4B), DIMENSION(:), INTENT(IN) :: indx - REAL(SP), DIMENSION(:), INTENT(INOUT) :: b - END SUBROUTINE lubksb - END INTERFACE - INTERFACE - SUBROUTINE ludcmp(a,indx,d) - USE nrtype - REAL(SP), DIMENSION(:,:), INTENT(INOUT) :: a - INTEGER(I4B), DIMENSION(:), INTENT(OUT) :: indx - REAL(SP), INTENT(OUT) :: d - END SUBROUTINE ludcmp - END INTERFACE - INTERFACE - SUBROUTINE machar(ibeta,it,irnd,ngrd,machep,negep,iexp,minexp,& - maxexp,eps,epsneg,xmin,xmax) - USE nrtype - INTEGER(I4B), INTENT(OUT) :: ibeta,iexp,irnd,it,machep,maxexp,& - minexp,negep,ngrd - REAL(SP), INTENT(OUT) :: eps,epsneg,xmax,xmin - END SUBROUTINE machar - END INTERFACE - INTERFACE - SUBROUTINE medfit(x,y,a,b,abdev) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(IN) :: x,y - REAL(SP), INTENT(OUT) :: a,b,abdev - END SUBROUTINE medfit - END INTERFACE - INTERFACE - SUBROUTINE memcof(data,xms,d) - USE nrtype - REAL(SP), INTENT(OUT) :: xms - REAL(SP), DIMENSION(:), INTENT(IN) :: data - REAL(SP), DIMENSION(:), INTENT(OUT) :: d - END SUBROUTINE memcof - END INTERFACE - INTERFACE - SUBROUTINE mgfas(u,maxcyc) - USE nrtype - REAL(DP), DIMENSION(:,:), INTENT(INOUT) :: u - INTEGER(I4B), INTENT(IN) :: maxcyc - END SUBROUTINE mgfas - END INTERFACE - INTERFACE - SUBROUTINE mglin(u,ncycle) - USE nrtype - REAL(DP), DIMENSION(:,:), INTENT(INOUT) :: u - INTEGER(I4B), INTENT(IN) :: ncycle - END SUBROUTINE mglin - END INTERFACE - INTERFACE - SUBROUTINE midexp(funk,aa,bb,s,n) - USE nrtype - REAL(SP), INTENT(IN) :: aa,bb - REAL(SP), INTENT(INOUT) :: s - INTEGER(I4B), INTENT(IN) :: n - INTERFACE - FUNCTION funk(x) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(IN) :: x - REAL(SP), DIMENSION(size(x)) :: funk - END FUNCTION funk - END INTERFACE - END SUBROUTINE midexp - END INTERFACE - INTERFACE - SUBROUTINE midinf(funk,aa,bb,s,n) - USE nrtype - REAL(SP), INTENT(IN) :: aa,bb - REAL(SP), INTENT(INOUT) :: s - INTEGER(I4B), INTENT(IN) :: n - INTERFACE - FUNCTION funk(x) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(IN) :: x - REAL(SP), DIMENSION(size(x)) :: funk - END FUNCTION funk - END INTERFACE - END SUBROUTINE midinf - END INTERFACE - INTERFACE - SUBROUTINE midpnt(func,a,b,s,n) - USE nrtype - REAL(SP), INTENT(IN) :: a,b - REAL(SP), INTENT(INOUT) :: s - INTEGER(I4B), INTENT(IN) :: n - INTERFACE - FUNCTION func(x) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(IN) :: x - REAL(SP), DIMENSION(size(x)) :: func - END FUNCTION func - END INTERFACE - END SUBROUTINE midpnt - END INTERFACE - INTERFACE - SUBROUTINE midsql(funk,aa,bb,s,n) - USE nrtype - REAL(SP), INTENT(IN) :: aa,bb - REAL(SP), INTENT(INOUT) :: s - INTEGER(I4B), INTENT(IN) :: n - INTERFACE - FUNCTION funk(x) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(IN) :: x - REAL(SP), DIMENSION(size(x)) :: funk - END FUNCTION funk - END INTERFACE - END SUBROUTINE midsql - END INTERFACE - INTERFACE - SUBROUTINE midsqu(funk,aa,bb,s,n) - USE nrtype - REAL(SP), INTENT(IN) :: aa,bb - REAL(SP), INTENT(INOUT) :: s - INTEGER(I4B), INTENT(IN) :: n - INTERFACE - FUNCTION funk(x) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(IN) :: x - REAL(SP), DIMENSION(size(x)) :: funk - END FUNCTION funk - END INTERFACE - END SUBROUTINE midsqu - END INTERFACE - INTERFACE - RECURSIVE SUBROUTINE miser(func,regn,ndim,npts,dith,ave,var) - USE nrtype - INTERFACE - FUNCTION func(x) - USE nrtype - REAL(SP) :: func - REAL(SP), DIMENSION(:), INTENT(IN) :: x - END FUNCTION func - END INTERFACE - REAL(SP), DIMENSION(:), INTENT(IN) :: regn - INTEGER(I4B), INTENT(IN) :: ndim,npts - REAL(SP), INTENT(IN) :: dith - REAL(SP), INTENT(OUT) :: ave,var - END SUBROUTINE miser - END INTERFACE - INTERFACE - SUBROUTINE mmid(y,dydx,xs,htot,nstep,yout,derivs) - USE nrtype - INTEGER(I4B), INTENT(IN) :: nstep - REAL(SP), INTENT(IN) :: xs,htot - REAL(SP), DIMENSION(:), INTENT(IN) :: y,dydx - REAL(SP), DIMENSION(:), INTENT(OUT) :: yout - INTERFACE - SUBROUTINE derivs(x,y,dydx) - USE nrtype - REAL(SP), INTENT(IN) :: x - REAL(SP), DIMENSION(:), INTENT(IN) :: y - REAL(SP), DIMENSION(:), INTENT(OUT) :: dydx - END SUBROUTINE derivs - END INTERFACE - END SUBROUTINE mmid - END INTERFACE - INTERFACE - SUBROUTINE mnbrak(ax,bx,cx,fa,fb,fc,func) - USE nrtype - REAL(SP), INTENT(INOUT) :: ax,bx - REAL(SP), INTENT(OUT) :: cx,fa,fb,fc - INTERFACE - FUNCTION func(x) - USE nrtype - REAL(SP), INTENT(IN) :: x - REAL(SP) :: func - END FUNCTION func - END INTERFACE - END SUBROUTINE mnbrak - END INTERFACE - INTERFACE - SUBROUTINE mnewt(ntrial,x,tolx,tolf,usrfun) - USE nrtype - INTEGER(I4B), INTENT(IN) :: ntrial - REAL(SP), INTENT(IN) :: tolx,tolf - REAL(SP), DIMENSION(:), INTENT(INOUT) :: x - INTERFACE - SUBROUTINE usrfun(x,fvec,fjac) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(IN) :: x - REAL(SP), DIMENSION(:), INTENT(OUT) :: fvec - REAL(SP), DIMENSION(:,:), INTENT(OUT) :: fjac - END SUBROUTINE usrfun - END INTERFACE - END SUBROUTINE mnewt - END INTERFACE - INTERFACE - SUBROUTINE moment(data,ave,adev,sdev,var,skew,curt) - USE nrtype - REAL(SP), INTENT(OUT) :: ave,adev,sdev,var,skew,curt - REAL(SP), DIMENSION(:), INTENT(IN) :: data - END SUBROUTINE moment - END INTERFACE - INTERFACE - SUBROUTINE mp2dfr(a,s,n,m) - USE nrtype - INTEGER(I4B), INTENT(IN) :: n - INTEGER(I4B), INTENT(OUT) :: m - CHARACTER(1), DIMENSION(:), INTENT(INOUT) :: a - CHARACTER(1), DIMENSION(:), INTENT(OUT) :: s - END SUBROUTINE mp2dfr - END INTERFACE - INTERFACE - SUBROUTINE mpdiv(q,r,u,v,n,m) - USE nrtype - CHARACTER(1), DIMENSION(:), INTENT(OUT) :: q,r - CHARACTER(1), DIMENSION(:), INTENT(IN) :: u,v - INTEGER(I4B), INTENT(IN) :: n,m - END SUBROUTINE mpdiv - END INTERFACE - INTERFACE - SUBROUTINE mpinv(u,v,n,m) - USE nrtype - CHARACTER(1), DIMENSION(:), INTENT(OUT) :: u - CHARACTER(1), DIMENSION(:), INTENT(IN) :: v - INTEGER(I4B), INTENT(IN) :: n,m - END SUBROUTINE mpinv - END INTERFACE - INTERFACE - SUBROUTINE mpmul(w,u,v,n,m) - USE nrtype - CHARACTER(1), DIMENSION(:), INTENT(IN) :: u,v - CHARACTER(1), DIMENSION(:), INTENT(OUT) :: w - INTEGER(I4B), INTENT(IN) :: n,m - END SUBROUTINE mpmul - END INTERFACE - INTERFACE - SUBROUTINE mppi(n) - USE nrtype - INTEGER(I4B), INTENT(IN) :: n - END SUBROUTINE mppi - END INTERFACE - INTERFACE - SUBROUTINE mprove(a,alud,indx,b,x) - USE nrtype - REAL(SP), DIMENSION(:,:), INTENT(IN) :: a,alud - INTEGER(I4B), DIMENSION(:), INTENT(IN) :: indx - REAL(SP), DIMENSION(:), INTENT(IN) :: b - REAL(SP), DIMENSION(:), INTENT(INOUT) :: x - END SUBROUTINE mprove - END INTERFACE - INTERFACE - SUBROUTINE mpsqrt(w,u,v,n,m) - USE nrtype - CHARACTER(1), DIMENSION(:), INTENT(OUT) :: w,u - CHARACTER(1), DIMENSION(:), INTENT(IN) :: v - INTEGER(I4B), INTENT(IN) :: n,m - END SUBROUTINE mpsqrt - END INTERFACE - INTERFACE - SUBROUTINE mrqcof(x,y,sig,a,maska,alpha,beta,chisq,funcs) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(IN) :: x,y,a,sig - REAL(SP), DIMENSION(:), INTENT(OUT) :: beta - REAL(SP), DIMENSION(:,:), INTENT(OUT) :: alpha - REAL(SP), INTENT(OUT) :: chisq - LOGICAL(LGT), DIMENSION(:), INTENT(IN) :: maska - INTERFACE - SUBROUTINE funcs(x,a,yfit,dyda) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(IN) :: x,a - REAL(SP), DIMENSION(:), INTENT(OUT) :: yfit - REAL(SP), DIMENSION(:,:), INTENT(OUT) :: dyda - END SUBROUTINE funcs - END INTERFACE - END SUBROUTINE mrqcof - END INTERFACE - INTERFACE - SUBROUTINE mrqmin(x,y,sig,a,maska,covar,alpha,chisq,funcs,alamda) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(IN) :: x,y,sig - REAL(SP), DIMENSION(:), INTENT(INOUT) :: a - REAL(SP), DIMENSION(:,:), INTENT(OUT) :: covar,alpha - REAL(SP), INTENT(OUT) :: chisq - REAL(SP), INTENT(INOUT) :: alamda - LOGICAL(LGT), DIMENSION(:), INTENT(IN) :: maska - INTERFACE - SUBROUTINE funcs(x,a,yfit,dyda) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(IN) :: x,a - REAL(SP), DIMENSION(:), INTENT(OUT) :: yfit - REAL(SP), DIMENSION(:,:), INTENT(OUT) :: dyda - END SUBROUTINE funcs - END INTERFACE - END SUBROUTINE mrqmin - END INTERFACE - INTERFACE - SUBROUTINE newt(x,check) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(INOUT) :: x - LOGICAL(LGT), INTENT(OUT) :: check - END SUBROUTINE newt - END INTERFACE - INTERFACE - SUBROUTINE odeint(ystart,x1,x2,eps,h1,hmin,derivs,rkqs) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(INOUT) :: ystart - REAL(SP), INTENT(IN) :: x1,x2,eps,h1,hmin - INTERFACE - SUBROUTINE derivs(x,y,dydx) - USE nrtype - REAL(SP), INTENT(IN) :: x - REAL(SP), DIMENSION(:), INTENT(IN) :: y - REAL(SP), DIMENSION(:), INTENT(OUT) :: dydx - END SUBROUTINE derivs -!BL - SUBROUTINE rkqs(y,dydx,x,htry,eps,yscal,hdid,hnext,derivs) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(INOUT) :: y - REAL(SP), DIMENSION(:), INTENT(IN) :: dydx,yscal - REAL(SP), INTENT(INOUT) :: x - REAL(SP), INTENT(IN) :: htry,eps - REAL(SP), INTENT(OUT) :: hdid,hnext - INTERFACE - SUBROUTINE derivs(x,y,dydx) - USE nrtype - REAL(SP), INTENT(IN) :: x - REAL(SP), DIMENSION(:), INTENT(IN) :: y - REAL(SP), DIMENSION(:), INTENT(OUT) :: dydx - END SUBROUTINE derivs - END INTERFACE - END SUBROUTINE rkqs - END INTERFACE - END SUBROUTINE odeint - END INTERFACE - INTERFACE - SUBROUTINE orthog(anu,alpha,beta,a,b) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(IN) :: anu,alpha,beta - REAL(SP), DIMENSION(:), INTENT(OUT) :: a,b - END SUBROUTINE orthog - END INTERFACE - INTERFACE - SUBROUTINE pade(cof,resid) - USE nrtype - REAL(DP), DIMENSION(:), INTENT(INOUT) :: cof - REAL(SP), INTENT(OUT) :: resid - END SUBROUTINE pade - END INTERFACE - INTERFACE - FUNCTION pccheb(d) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(IN) :: d - REAL(SP), DIMENSION(size(d)) :: pccheb - END FUNCTION pccheb - END INTERFACE - INTERFACE - SUBROUTINE pcshft(a,b,d) - USE nrtype - REAL(SP), INTENT(IN) :: a,b - REAL(SP), DIMENSION(:), INTENT(INOUT) :: d - END SUBROUTINE pcshft - END INTERFACE - INTERFACE - SUBROUTINE pearsn(x,y,r,prob,z) - USE nrtype - REAL(SP), INTENT(OUT) :: r,prob,z - REAL(SP), DIMENSION(:), INTENT(IN) :: x,y - END SUBROUTINE pearsn - END INTERFACE - INTERFACE - SUBROUTINE period(x,y,ofac,hifac,px,py,jmax,prob) - USE nrtype - INTEGER(I4B), INTENT(OUT) :: jmax - REAL(SP), INTENT(IN) :: ofac,hifac - REAL(SP), INTENT(OUT) :: prob - REAL(SP), DIMENSION(:), INTENT(IN) :: x,y - REAL(SP), DIMENSION(:), POINTER :: px,py - END SUBROUTINE period - END INTERFACE - INTERFACE plgndr - FUNCTION plgndr_s(l,m,x) - USE nrtype - INTEGER(I4B), INTENT(IN) :: l,m - REAL(SP), INTENT(IN) :: x - REAL(SP) :: plgndr_s - END FUNCTION plgndr_s -!BL - FUNCTION plgndr_v(l,m,x) - USE nrtype - INTEGER(I4B), INTENT(IN) :: l,m - REAL(SP), DIMENSION(:), INTENT(IN) :: x - REAL(SP), DIMENSION(size(x)) :: plgndr_v - END FUNCTION plgndr_v - END INTERFACE - INTERFACE - FUNCTION poidev(xm) - USE nrtype - REAL(SP), INTENT(IN) :: xm - REAL(SP) :: poidev - END FUNCTION poidev - END INTERFACE - INTERFACE - FUNCTION polcoe(x,y) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(IN) :: x,y - REAL(SP), DIMENSION(size(x)) :: polcoe - END FUNCTION polcoe - END INTERFACE - INTERFACE - FUNCTION polcof(xa,ya) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(IN) :: xa,ya - REAL(SP), DIMENSION(size(xa)) :: polcof - END FUNCTION polcof - END INTERFACE - INTERFACE - SUBROUTINE poldiv(u,v,q,r) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(IN) :: u,v - REAL(SP), DIMENSION(:), INTENT(OUT) :: q,r - END SUBROUTINE poldiv - END INTERFACE - INTERFACE - SUBROUTINE polin2(x1a,x2a,ya,x1,x2,y,dy) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(IN) :: x1a,x2a - REAL(SP), DIMENSION(:,:), INTENT(IN) :: ya - REAL(SP), INTENT(IN) :: x1,x2 - REAL(SP), INTENT(OUT) :: y,dy - END SUBROUTINE polin2 - END INTERFACE - INTERFACE - SUBROUTINE polint(xa,ya,x,y,dy) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(IN) :: xa,ya - REAL(SP), INTENT(IN) :: x - REAL(SP), INTENT(OUT) :: y,dy - END SUBROUTINE polint - END INTERFACE - INTERFACE - SUBROUTINE powell(p,xi,ftol,iter,fret) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(INOUT) :: p - REAL(SP), DIMENSION(:,:), INTENT(INOUT) :: xi - INTEGER(I4B), INTENT(OUT) :: iter - REAL(SP), INTENT(IN) :: ftol - REAL(SP), INTENT(OUT) :: fret - END SUBROUTINE powell - END INTERFACE - INTERFACE - FUNCTION predic(data,d,nfut) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(IN) :: data,d - INTEGER(I4B), INTENT(IN) :: nfut - REAL(SP), DIMENSION(nfut) :: predic - END FUNCTION predic - END INTERFACE - INTERFACE - FUNCTION probks(alam) - USE nrtype - REAL(SP), INTENT(IN) :: alam - REAL(SP) :: probks - END FUNCTION probks - END INTERFACE - INTERFACE psdes - SUBROUTINE psdes_s(lword,rword) - USE nrtype - INTEGER(I4B), INTENT(INOUT) :: lword,rword - END SUBROUTINE psdes_s -!BL - SUBROUTINE psdes_v(lword,rword) - USE nrtype - INTEGER(I4B), DIMENSION(:), INTENT(INOUT) :: lword,rword - END SUBROUTINE psdes_v - END INTERFACE - INTERFACE - SUBROUTINE pwt(a,isign) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(INOUT) :: a - INTEGER(I4B), INTENT(IN) :: isign - END SUBROUTINE pwt - END INTERFACE - INTERFACE - SUBROUTINE pwtset(n) - USE nrtype - INTEGER(I4B), INTENT(IN) :: n - END SUBROUTINE pwtset - END INTERFACE - INTERFACE pythag -! FUNCTION pythag_dp(a,b) -! USE nrtype -! REAL(DP), INTENT(IN) :: a,b -! REAL(DP) :: pythag_dp -! END FUNCTION pythag_dp -!BL - FUNCTION pythag_sp(a,b) - USE nrtype - REAL(SP), INTENT(IN) :: a,b - REAL(SP) :: pythag_sp - END FUNCTION pythag_sp - END INTERFACE - INTERFACE - SUBROUTINE pzextr(iest,xest,yest,yz,dy) - USE nrtype - INTEGER(I4B), INTENT(IN) :: iest - REAL(SP), INTENT(IN) :: xest - REAL(SP), DIMENSION(:), INTENT(IN) :: yest - REAL(SP), DIMENSION(:), INTENT(OUT) :: yz,dy - END SUBROUTINE pzextr - END INTERFACE - INTERFACE - SUBROUTINE qrdcmp(a,c,d,sing) - USE nrtype - REAL(SP), DIMENSION(:,:), INTENT(INOUT) :: a - REAL(SP), DIMENSION(:), INTENT(OUT) :: c,d - LOGICAL(LGT), INTENT(OUT) :: sing - END SUBROUTINE qrdcmp - END INTERFACE - INTERFACE - FUNCTION qromb(func,a,b) - USE nrtype - REAL(SP), INTENT(IN) :: a,b - REAL(SP) :: qromb - INTERFACE - FUNCTION func(x) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(IN) :: x - REAL(SP), DIMENSION(size(x)) :: func - END FUNCTION func - END INTERFACE - END FUNCTION qromb - END INTERFACE - INTERFACE - FUNCTION qromo(func,a,b,choose) - USE nrtype - REAL(SP), INTENT(IN) :: a,b - REAL(SP) :: qromo - INTERFACE - FUNCTION func(x) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(IN) :: x - REAL(SP), DIMENSION(size(x)) :: func - END FUNCTION func - END INTERFACE - INTERFACE - SUBROUTINE choose(funk,aa,bb,s,n) - USE nrtype - REAL(SP), INTENT(IN) :: aa,bb - REAL(SP), INTENT(INOUT) :: s - INTEGER(I4B), INTENT(IN) :: n - INTERFACE - FUNCTION funk(x) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(IN) :: x - REAL(SP), DIMENSION(size(x)) :: funk - END FUNCTION funk - END INTERFACE - END SUBROUTINE choose - END INTERFACE - END FUNCTION qromo - END INTERFACE - INTERFACE - SUBROUTINE qroot(p,b,c,eps) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(IN) :: p - REAL(SP), INTENT(INOUT) :: b,c - REAL(SP), INTENT(IN) :: eps - END SUBROUTINE qroot - END INTERFACE - INTERFACE - SUBROUTINE qrsolv(a,c,d,b) - USE nrtype - REAL(SP), DIMENSION(:,:), INTENT(IN) :: a - REAL(SP), DIMENSION(:), INTENT(IN) :: c,d - REAL(SP), DIMENSION(:), INTENT(INOUT) :: b - END SUBROUTINE qrsolv - END INTERFACE - INTERFACE - SUBROUTINE qrupdt(r,qt,u,v) - USE nrtype - REAL(SP), DIMENSION(:,:), INTENT(INOUT) :: r,qt - REAL(SP), DIMENSION(:), INTENT(INOUT) :: u - REAL(SP), DIMENSION(:), INTENT(IN) :: v - END SUBROUTINE qrupdt - END INTERFACE - INTERFACE - FUNCTION qsimp(func,a,b) - USE nrtype - REAL(SP), INTENT(IN) :: a,b - REAL(SP) :: qsimp - INTERFACE - FUNCTION func(x) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(IN) :: x - REAL(SP), DIMENSION(size(x)) :: func - END FUNCTION func - END INTERFACE - END FUNCTION qsimp - END INTERFACE - INTERFACE - FUNCTION qtrap(func,a,b) - USE nrtype - REAL(SP), INTENT(IN) :: a,b - REAL(SP) :: qtrap - INTERFACE - FUNCTION func(x) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(IN) :: x - REAL(SP), DIMENSION(size(x)) :: func - END FUNCTION func - END INTERFACE - END FUNCTION qtrap - END INTERFACE - INTERFACE - SUBROUTINE quadct(x,y,xx,yy,fa,fb,fc,fd) - USE nrtype - REAL(SP), INTENT(IN) :: x,y - REAL(SP), DIMENSION(:), INTENT(IN) :: xx,yy - REAL(SP), INTENT(OUT) :: fa,fb,fc,fd - END SUBROUTINE quadct - END INTERFACE - INTERFACE - SUBROUTINE quadmx(a) - USE nrtype - REAL(SP), DIMENSION(:,:), INTENT(OUT) :: a - END SUBROUTINE quadmx - END INTERFACE - INTERFACE - SUBROUTINE quadvl(x,y,fa,fb,fc,fd) - USE nrtype - REAL(SP), INTENT(IN) :: x,y - REAL(SP), INTENT(OUT) :: fa,fb,fc,fd - END SUBROUTINE quadvl - END INTERFACE - INTERFACE - FUNCTION ran(idum) - INTEGER(selected_int_kind(9)), INTENT(INOUT) :: idum - REAL :: ran - END FUNCTION ran - END INTERFACE - INTERFACE ran0 - SUBROUTINE ran0_s(harvest) - USE nrtype - REAL(SP), INTENT(OUT) :: harvest - END SUBROUTINE ran0_s -!BL - SUBROUTINE ran0_v(harvest) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(OUT) :: harvest - END SUBROUTINE ran0_v - END INTERFACE - INTERFACE ran1 - SUBROUTINE ran1_s(harvest) - USE nrtype - REAL(SP), INTENT(OUT) :: harvest - END SUBROUTINE ran1_s -!BL - SUBROUTINE ran1_v(harvest) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(OUT) :: harvest - END SUBROUTINE ran1_v - END INTERFACE - INTERFACE ran2 - SUBROUTINE ran2_s(harvest) - USE nrtype - REAL(SP), INTENT(OUT) :: harvest - END SUBROUTINE ran2_s -!BL - SUBROUTINE ran2_v(harvest) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(OUT) :: harvest - END SUBROUTINE ran2_v - END INTERFACE - INTERFACE ran3 - SUBROUTINE ran3_s(harvest) - USE nrtype - REAL(SP), INTENT(OUT) :: harvest - END SUBROUTINE ran3_s -!BL - SUBROUTINE ran3_v(harvest) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(OUT) :: harvest - END SUBROUTINE ran3_v - END INTERFACE - INTERFACE - SUBROUTINE ratint(xa,ya,x,y,dy) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(IN) :: xa,ya - REAL(SP), INTENT(IN) :: x - REAL(SP), INTENT(OUT) :: y,dy - END SUBROUTINE ratint - END INTERFACE - INTERFACE - SUBROUTINE ratlsq(func,a,b,mm,kk,cof,dev) - USE nrtype - REAL(DP), INTENT(IN) :: a,b - INTEGER(I4B), INTENT(IN) :: mm,kk - REAL(DP), DIMENSION(:), INTENT(OUT) :: cof - REAL(DP), INTENT(OUT) :: dev - INTERFACE - FUNCTION func(x) - USE nrtype - REAL(DP), DIMENSION(:), INTENT(IN) :: x - REAL(DP), DIMENSION(size(x)) :: func - END FUNCTION func - END INTERFACE - END SUBROUTINE ratlsq - END INTERFACE - INTERFACE ratval - FUNCTION ratval_s(x,cof,mm,kk) - USE nrtype - REAL(DP), INTENT(IN) :: x - INTEGER(I4B), INTENT(IN) :: mm,kk - REAL(DP), DIMENSION(mm+kk+1), INTENT(IN) :: cof - REAL(DP) :: ratval_s - END FUNCTION ratval_s -!BL - FUNCTION ratval_v(x,cof,mm,kk) - USE nrtype - REAL(DP), DIMENSION(:), INTENT(IN) :: x - INTEGER(I4B), INTENT(IN) :: mm,kk - REAL(DP), DIMENSION(mm+kk+1), INTENT(IN) :: cof - REAL(DP), DIMENSION(size(x)) :: ratval_v - END FUNCTION ratval_v - END INTERFACE - INTERFACE rc - FUNCTION rc_s(x,y) - USE nrtype - REAL(SP), INTENT(IN) :: x,y - REAL(SP) :: rc_s - END FUNCTION rc_s -!BL - FUNCTION rc_v(x,y) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(IN) :: x,y - REAL(SP), DIMENSION(size(x)) :: rc_v - END FUNCTION rc_v - END INTERFACE - INTERFACE rd - FUNCTION rd_s(x,y,z) - USE nrtype - REAL(SP), INTENT(IN) :: x,y,z - REAL(SP) :: rd_s - END FUNCTION rd_s -!BL - FUNCTION rd_v(x,y,z) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(IN) :: x,y,z - REAL(SP), DIMENSION(size(x)) :: rd_v - END FUNCTION rd_v - END INTERFACE - INTERFACE realft -! SUBROUTINE realft_dp(data,isign,zdata) -! USE nrtype -! REAL(DP), DIMENSION(:), INTENT(INOUT) :: data -! INTEGER(I4B), INTENT(IN) :: isign -! COMPLEX(DPC), DIMENSION(:), OPTIONAL, TARGET :: zdata -! END SUBROUTINE realft_dp -!BL - SUBROUTINE realft_sp(data,isign,zdata) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(INOUT) :: data - INTEGER(I4B), INTENT(IN) :: isign - COMPLEX(SPC), DIMENSION(:), OPTIONAL, TARGET :: zdata - END SUBROUTINE realft_sp - END INTERFACE - INTERFACE - RECURSIVE FUNCTION recur1(a,b) RESULT(u) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(IN) :: a,b - REAL(SP), DIMENSION(size(a)) :: u - END FUNCTION recur1 - END INTERFACE - INTERFACE - FUNCTION recur2(a,b,c) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(IN) :: a,b,c - REAL(SP), DIMENSION(size(a)) :: recur2 - END FUNCTION recur2 - END INTERFACE - INTERFACE - SUBROUTINE relax(u,rhs) - USE nrtype - REAL(DP), DIMENSION(:,:), INTENT(INOUT) :: u - REAL(DP), DIMENSION(:,:), INTENT(IN) :: rhs - END SUBROUTINE relax - END INTERFACE - INTERFACE - SUBROUTINE relax2(u,rhs) - USE nrtype - REAL(DP), DIMENSION(:,:), INTENT(INOUT) :: u - REAL(DP), DIMENSION(:,:), INTENT(IN) :: rhs - END SUBROUTINE relax2 - END INTERFACE - INTERFACE - FUNCTION resid(u,rhs) - USE nrtype - REAL(DP), DIMENSION(:,:), INTENT(IN) :: u,rhs - REAL(DP), DIMENSION(size(u,1),size(u,1)) :: resid - END FUNCTION resid - END INTERFACE - INTERFACE rf - FUNCTION rf_s(x,y,z) - USE nrtype - REAL(SP), INTENT(IN) :: x,y,z - REAL(SP) :: rf_s - END FUNCTION rf_s -!BL - FUNCTION rf_v(x,y,z) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(IN) :: x,y,z - REAL(SP), DIMENSION(size(x)) :: rf_v - END FUNCTION rf_v - END INTERFACE - INTERFACE rj - FUNCTION rj_s(x,y,z,p) - USE nrtype - REAL(SP), INTENT(IN) :: x,y,z,p - REAL(SP) :: rj_s - END FUNCTION rj_s -!BL - FUNCTION rj_v(x,y,z,p) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(IN) :: x,y,z,p - REAL(SP), DIMENSION(size(x)) :: rj_v - END FUNCTION rj_v - END INTERFACE - INTERFACE - SUBROUTINE rk4(y,dydx,x,h,yout,derivs) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(IN) :: y,dydx - REAL(SP), INTENT(IN) :: x,h - REAL(SP), DIMENSION(:), INTENT(OUT) :: yout - INTERFACE - SUBROUTINE derivs(x,y,dydx) - USE nrtype - REAL(SP), INTENT(IN) :: x - REAL(SP), DIMENSION(:), INTENT(IN) :: y - REAL(SP), DIMENSION(:), INTENT(OUT) :: dydx - END SUBROUTINE derivs - END INTERFACE - END SUBROUTINE rk4 - END INTERFACE - INTERFACE - SUBROUTINE rkck(y,dydx,x,h,yout,yerr,derivs) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(IN) :: y,dydx - REAL(SP), INTENT(IN) :: x,h - REAL(SP), DIMENSION(:), INTENT(OUT) :: yout,yerr - INTERFACE - SUBROUTINE derivs(x,y,dydx) - USE nrtype - REAL(SP), INTENT(IN) :: x - REAL(SP), DIMENSION(:), INTENT(IN) :: y - REAL(SP), DIMENSION(:), INTENT(OUT) :: dydx - END SUBROUTINE derivs - END INTERFACE - END SUBROUTINE rkck - END INTERFACE - INTERFACE - SUBROUTINE rkdumb(vstart,x1,x2,nstep,derivs) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(IN) :: vstart - REAL(SP), INTENT(IN) :: x1,x2 - INTEGER(I4B), INTENT(IN) :: nstep - INTERFACE - SUBROUTINE derivs(x,y,dydx) - USE nrtype - REAL(SP), INTENT(IN) :: x - REAL(SP), DIMENSION(:), INTENT(IN) :: y - REAL(SP), DIMENSION(:), INTENT(OUT) :: dydx - END SUBROUTINE derivs - END INTERFACE - END SUBROUTINE rkdumb - END INTERFACE - INTERFACE - SUBROUTINE rkqs(y,dydx,x,htry,eps,yscal,hdid,hnext,derivs) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(INOUT) :: y - REAL(SP), DIMENSION(:), INTENT(IN) :: dydx,yscal - REAL(SP), INTENT(INOUT) :: x - REAL(SP), INTENT(IN) :: htry,eps - REAL(SP), INTENT(OUT) :: hdid,hnext - INTERFACE - SUBROUTINE derivs(x,y,dydx) - USE nrtype - REAL(SP), INTENT(IN) :: x - REAL(SP), DIMENSION(:), INTENT(IN) :: y - REAL(SP), DIMENSION(:), INTENT(OUT) :: dydx - END SUBROUTINE derivs - END INTERFACE - END SUBROUTINE rkqs - END INTERFACE - INTERFACE - SUBROUTINE rlft2(data,spec,speq,isign) - USE nrtype - REAL(SP), DIMENSION(:,:), INTENT(INOUT) :: data - COMPLEX(SPC), DIMENSION(:,:), INTENT(OUT) :: spec - COMPLEX(SPC), DIMENSION(:), INTENT(OUT) :: speq - INTEGER(I4B), INTENT(IN) :: isign - END SUBROUTINE rlft2 - END INTERFACE - INTERFACE - SUBROUTINE rlft3(data,spec,speq,isign) - USE nrtype - REAL(SP), DIMENSION(:,:,:), INTENT(INOUT) :: data - COMPLEX(SPC), DIMENSION(:,:,:), INTENT(OUT) :: spec - COMPLEX(SPC), DIMENSION(:,:), INTENT(OUT) :: speq - INTEGER(I4B), INTENT(IN) :: isign - END SUBROUTINE rlft3 - END INTERFACE - INTERFACE - SUBROUTINE rotate(r,qt,i,a,b) - USE nrtype - REAL(SP), DIMENSION(:,:), TARGET, INTENT(INOUT) :: r,qt - INTEGER(I4B), INTENT(IN) :: i - REAL(SP), INTENT(IN) :: a,b - END SUBROUTINE rotate - END INTERFACE - INTERFACE - SUBROUTINE rsolv(a,d,b) - USE nrtype - REAL(SP), DIMENSION(:,:), INTENT(IN) :: a - REAL(SP), DIMENSION(:), INTENT(IN) :: d - REAL(SP), DIMENSION(:), INTENT(INOUT) :: b - END SUBROUTINE rsolv - END INTERFACE - INTERFACE - FUNCTION rstrct(uf) - USE nrtype - REAL(DP), DIMENSION(:,:), INTENT(IN) :: uf - REAL(DP), DIMENSION((size(uf,1)+1)/2,(size(uf,1)+1)/2) :: rstrct - END FUNCTION rstrct - END INTERFACE - INTERFACE - FUNCTION rtbis(func,x1,x2,xacc) - USE nrtype - REAL(SP), INTENT(IN) :: x1,x2,xacc - REAL(SP) :: rtbis - INTERFACE - FUNCTION func(x) - USE nrtype - REAL(SP), INTENT(IN) :: x - REAL(SP) :: func - END FUNCTION func - END INTERFACE - END FUNCTION rtbis - END INTERFACE - INTERFACE - FUNCTION rtflsp(func,x1,x2,xacc) - USE nrtype - REAL(SP), INTENT(IN) :: x1,x2,xacc - REAL(SP) :: rtflsp - INTERFACE - FUNCTION func(x) - USE nrtype - REAL(SP), INTENT(IN) :: x - REAL(SP) :: func - END FUNCTION func - END INTERFACE - END FUNCTION rtflsp - END INTERFACE - INTERFACE - FUNCTION rtnewt(funcd,x1,x2,xacc) - USE nrtype - REAL(SP), INTENT(IN) :: x1,x2,xacc - REAL(SP) :: rtnewt - INTERFACE - SUBROUTINE funcd(x,fval,fderiv) - USE nrtype - REAL(SP), INTENT(IN) :: x - REAL(SP), INTENT(OUT) :: fval,fderiv - END SUBROUTINE funcd - END INTERFACE - END FUNCTION rtnewt - END INTERFACE - INTERFACE - FUNCTION rtsafe(funcd,x1,x2,xacc) - USE nrtype - REAL(SP), INTENT(IN) :: x1,x2,xacc - REAL(SP) :: rtsafe - INTERFACE - SUBROUTINE funcd(x,fval,fderiv) - USE nrtype - REAL(SP), INTENT(IN) :: x - REAL(SP), INTENT(OUT) :: fval,fderiv - END SUBROUTINE funcd - END INTERFACE - END FUNCTION rtsafe - END INTERFACE - INTERFACE - FUNCTION rtsec(func,x1,x2,xacc) - USE nrtype - REAL(SP), INTENT(IN) :: x1,x2,xacc - REAL(SP) :: rtsec - INTERFACE - FUNCTION func(x) - USE nrtype - REAL(SP), INTENT(IN) :: x - REAL(SP) :: func - END FUNCTION func - END INTERFACE - END FUNCTION rtsec - END INTERFACE - INTERFACE - SUBROUTINE rzextr(iest,xest,yest,yz,dy) - USE nrtype - INTEGER(I4B), INTENT(IN) :: iest - REAL(SP), INTENT(IN) :: xest - REAL(SP), DIMENSION(:), INTENT(IN) :: yest - REAL(SP), DIMENSION(:), INTENT(OUT) :: yz,dy - END SUBROUTINE rzextr - END INTERFACE - INTERFACE - FUNCTION savgol(nl,nrr,ld,m) - USE nrtype - INTEGER(I4B), INTENT(IN) :: nl,nrr,ld,m - REAL(SP), DIMENSION(nl+nrr+1) :: savgol - END FUNCTION savgol - END INTERFACE - INTERFACE - SUBROUTINE scrsho(func) - USE nrtype - INTERFACE - FUNCTION func(x) - USE nrtype - REAL(SP), INTENT(IN) :: x - REAL(SP) :: func - END FUNCTION func - END INTERFACE - END SUBROUTINE scrsho - END INTERFACE - INTERFACE - FUNCTION select(k,arr) - USE nrtype - INTEGER(I4B), INTENT(IN) :: k - REAL(SP), DIMENSION(:), INTENT(INOUT) :: arr - REAL(SP) :: select - END FUNCTION select - END INTERFACE - INTERFACE - FUNCTION select_bypack(k,arr) - USE nrtype - INTEGER(I4B), INTENT(IN) :: k - REAL(SP), DIMENSION(:), INTENT(INOUT) :: arr - REAL(SP) :: select_bypack - END FUNCTION select_bypack - END INTERFACE - INTERFACE - SUBROUTINE select_heap(arr,heap) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(IN) :: arr - REAL(SP), DIMENSION(:), INTENT(OUT) :: heap - END SUBROUTINE select_heap - END INTERFACE - INTERFACE - FUNCTION select_inplace(k,arr) - USE nrtype - INTEGER(I4B), INTENT(IN) :: k - REAL(SP), DIMENSION(:), INTENT(IN) :: arr - REAL(SP) :: select_inplace - END FUNCTION select_inplace - END INTERFACE - INTERFACE - SUBROUTINE simplx(a,m1,m2,m3,icase,izrov,iposv) - USE nrtype - REAL(SP), DIMENSION(:,:), INTENT(INOUT) :: a - INTEGER(I4B), INTENT(IN) :: m1,m2,m3 - INTEGER(I4B), INTENT(OUT) :: icase - INTEGER(I4B), DIMENSION(:), INTENT(OUT) :: izrov,iposv - END SUBROUTINE simplx - END INTERFACE - INTERFACE - SUBROUTINE simpr(y,dydx,dfdx,dfdy,xs,htot,nstep,yout,derivs) - USE nrtype - REAL(SP), INTENT(IN) :: xs,htot - REAL(SP), DIMENSION(:), INTENT(IN) :: y,dydx,dfdx - REAL(SP), DIMENSION(:,:), INTENT(IN) :: dfdy - INTEGER(I4B), INTENT(IN) :: nstep - REAL(SP), DIMENSION(:), INTENT(OUT) :: yout - INTERFACE - SUBROUTINE derivs(x,y,dydx) - USE nrtype - REAL(SP), INTENT(IN) :: x - REAL(SP), DIMENSION(:), INTENT(IN) :: y - REAL(SP), DIMENSION(:), INTENT(OUT) :: dydx - END SUBROUTINE derivs - END INTERFACE - END SUBROUTINE simpr - END INTERFACE - INTERFACE - SUBROUTINE sinft(y) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(INOUT) :: y - END SUBROUTINE sinft - END INTERFACE - INTERFACE - SUBROUTINE slvsm2(u,rhs) - USE nrtype - REAL(DP), DIMENSION(3,3), INTENT(OUT) :: u - REAL(DP), DIMENSION(3,3), INTENT(IN) :: rhs - END SUBROUTINE slvsm2 - END INTERFACE - INTERFACE - SUBROUTINE slvsml(u,rhs) - USE nrtype - REAL(DP), DIMENSION(3,3), INTENT(OUT) :: u - REAL(DP), DIMENSION(3,3), INTENT(IN) :: rhs - END SUBROUTINE slvsml - END INTERFACE - INTERFACE - SUBROUTINE sncndn(uu,emmc,sn,cn,dn) - USE nrtype - REAL(SP), INTENT(IN) :: uu,emmc - REAL(SP), INTENT(OUT) :: sn,cn,dn - END SUBROUTINE sncndn - END INTERFACE - INTERFACE - FUNCTION snrm(sx,itol) - USE nrtype - REAL(DP), DIMENSION(:), INTENT(IN) :: sx - INTEGER(I4B), INTENT(IN) :: itol - REAL(DP) :: snrm - END FUNCTION snrm - END INTERFACE - INTERFACE - SUBROUTINE sobseq(x,init) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(OUT) :: x - INTEGER(I4B), OPTIONAL, INTENT(IN) :: init - END SUBROUTINE sobseq - END INTERFACE - INTERFACE - SUBROUTINE solvde(itmax,conv,slowc,scalv,indexv,nb,y) - USE nrtype - INTEGER(I4B), INTENT(IN) :: itmax,nb - REAL(SP), INTENT(IN) :: conv,slowc - REAL(SP), DIMENSION(:), INTENT(IN) :: scalv - INTEGER(I4B), DIMENSION(:), INTENT(IN) :: indexv - REAL(SP), DIMENSION(:,:), INTENT(INOUT) :: y - END SUBROUTINE solvde - END INTERFACE - INTERFACE - SUBROUTINE sor(a,b,c,d,e,f,u,rjac) - USE nrtype - REAL(DP), DIMENSION(:,:), INTENT(IN) :: a,b,c,d,e,f - REAL(DP), DIMENSION(:,:), INTENT(INOUT) :: u - REAL(DP), INTENT(IN) :: rjac - END SUBROUTINE sor - END INTERFACE - INTERFACE - SUBROUTINE sort(arr) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(INOUT) :: arr - END SUBROUTINE sort - END INTERFACE - INTERFACE - SUBROUTINE sort2(arr,slave) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(INOUT) :: arr,slave - END SUBROUTINE sort2 - END INTERFACE - INTERFACE - SUBROUTINE sort3(arr,slave1,slave2) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(INOUT) :: arr,slave1,slave2 - END SUBROUTINE sort3 - END INTERFACE - INTERFACE - SUBROUTINE sort_bypack(arr) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(INOUT) :: arr - END SUBROUTINE sort_bypack - END INTERFACE - INTERFACE - SUBROUTINE sort_byreshape(arr) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(INOUT) :: arr - END SUBROUTINE sort_byreshape - END INTERFACE - INTERFACE - SUBROUTINE sort_heap(arr) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(INOUT) :: arr - END SUBROUTINE sort_heap - END INTERFACE - INTERFACE - SUBROUTINE sort_pick(arr) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(INOUT) :: arr - END SUBROUTINE sort_pick - END INTERFACE - INTERFACE - SUBROUTINE sort_radix(arr) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(INOUT) :: arr - END SUBROUTINE sort_radix - END INTERFACE - INTERFACE - SUBROUTINE sort_shell(arr) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(INOUT) :: arr - END SUBROUTINE sort_shell - END INTERFACE - INTERFACE - SUBROUTINE spctrm(p,k,ovrlap,unit,n_window) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(OUT) :: p - INTEGER(I4B), INTENT(IN) :: k - LOGICAL(LGT), INTENT(IN) :: ovrlap - INTEGER(I4B), OPTIONAL, INTENT(IN) :: n_window,unit - END SUBROUTINE spctrm - END INTERFACE - INTERFACE - SUBROUTINE spear(data1,data2,d,zd,probd,rs,probrs) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(IN) :: data1,data2 - REAL(SP), INTENT(OUT) :: d,zd,probd,rs,probrs - END SUBROUTINE spear - END INTERFACE - INTERFACE sphbes - SUBROUTINE sphbes_s(n,x,sj,sy,sjp,syp) - USE nrtype - INTEGER(I4B), INTENT(IN) :: n - REAL(SP), INTENT(IN) :: x - REAL(SP), INTENT(OUT) :: sj,sy,sjp,syp - END SUBROUTINE sphbes_s -!BL - SUBROUTINE sphbes_v(n,x,sj,sy,sjp,syp) - USE nrtype - INTEGER(I4B), INTENT(IN) :: n - REAL(SP), DIMENSION(:), INTENT(IN) :: x - REAL(SP), DIMENSION(:), INTENT(OUT) :: sj,sy,sjp,syp - END SUBROUTINE sphbes_v - END INTERFACE - INTERFACE - SUBROUTINE splie2(x1a,x2a,ya,y2a) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(IN) :: x1a,x2a - REAL(SP), DIMENSION(:,:), INTENT(IN) :: ya - REAL(SP), DIMENSION(:,:), INTENT(OUT) :: y2a - END SUBROUTINE splie2 - END INTERFACE - INTERFACE - FUNCTION splin2(x1a,x2a,ya,y2a,x1,x2) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(IN) :: x1a,x2a - REAL(SP), DIMENSION(:,:), INTENT(IN) :: ya,y2a - REAL(SP), INTENT(IN) :: x1,x2 - REAL(SP) :: splin2 - END FUNCTION splin2 - END INTERFACE - INTERFACE - SUBROUTINE spline(x,y,yp1,ypn,y2) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(IN) :: x,y - REAL(SP), INTENT(IN) :: yp1,ypn - REAL(SP), DIMENSION(:), INTENT(OUT) :: y2 - END SUBROUTINE spline - END INTERFACE - INTERFACE - FUNCTION splint(xa,ya,y2a,x) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(IN) :: xa,ya,y2a - REAL(SP), INTENT(IN) :: x - REAL(SP) :: splint - END FUNCTION splint - END INTERFACE - INTERFACE sprsax - SUBROUTINE sprsax_dp(sa,x,b) - USE nrtype - TYPE(sprs2_dp), INTENT(IN) :: sa - REAL(DP), DIMENSION (:), INTENT(IN) :: x - REAL(DP), DIMENSION (:), INTENT(OUT) :: b - END SUBROUTINE sprsax_dp -!BL - SUBROUTINE sprsax_sp(sa,x,b) - USE nrtype - TYPE(sprs2_sp), INTENT(IN) :: sa - REAL(SP), DIMENSION (:), INTENT(IN) :: x - REAL(SP), DIMENSION (:), INTENT(OUT) :: b - END SUBROUTINE sprsax_sp - END INTERFACE - INTERFACE sprsdiag - SUBROUTINE sprsdiag_dp(sa,b) - USE nrtype - TYPE(sprs2_dp), INTENT(IN) :: sa - REAL(DP), DIMENSION(:), INTENT(OUT) :: b - END SUBROUTINE sprsdiag_dp -!BL - SUBROUTINE sprsdiag_sp(sa,b) - USE nrtype - TYPE(sprs2_sp), INTENT(IN) :: sa - REAL(SP), DIMENSION(:), INTENT(OUT) :: b - END SUBROUTINE sprsdiag_sp - END INTERFACE - INTERFACE sprsin - SUBROUTINE sprsin_sp(a,thresh,sa) - USE nrtype - REAL(SP), DIMENSION(:,:), INTENT(IN) :: a - REAL(SP), INTENT(IN) :: thresh - TYPE(sprs2_sp), INTENT(OUT) :: sa - END SUBROUTINE sprsin_sp -!BL - SUBROUTINE sprsin_dp(a,thresh,sa) - USE nrtype - REAL(DP), DIMENSION(:,:), INTENT(IN) :: a - REAL(DP), INTENT(IN) :: thresh - TYPE(sprs2_dp), INTENT(OUT) :: sa - END SUBROUTINE sprsin_dp - END INTERFACE - INTERFACE - SUBROUTINE sprstp(sa) - USE nrtype - TYPE(sprs2_sp), INTENT(INOUT) :: sa - END SUBROUTINE sprstp - END INTERFACE - INTERFACE sprstx - SUBROUTINE sprstx_dp(sa,x,b) - USE nrtype - TYPE(sprs2_dp), INTENT(IN) :: sa - REAL(DP), DIMENSION (:), INTENT(IN) :: x - REAL(DP), DIMENSION (:), INTENT(OUT) :: b - END SUBROUTINE sprstx_dp -!BL - SUBROUTINE sprstx_sp(sa,x,b) - USE nrtype - TYPE(sprs2_sp), INTENT(IN) :: sa - REAL(SP), DIMENSION (:), INTENT(IN) :: x - REAL(SP), DIMENSION (:), INTENT(OUT) :: b - END SUBROUTINE sprstx_sp - END INTERFACE - INTERFACE - SUBROUTINE stifbs(y,dydx,x,htry,eps,yscal,hdid,hnext,derivs) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(INOUT) :: y - REAL(SP), DIMENSION(:), INTENT(IN) :: dydx,yscal - REAL(SP), INTENT(IN) :: htry,eps - REAL(SP), INTENT(INOUT) :: x - REAL(SP), INTENT(OUT) :: hdid,hnext - INTERFACE - SUBROUTINE derivs(x,y,dydx) - USE nrtype - REAL(SP), INTENT(IN) :: x - REAL(SP), DIMENSION(:), INTENT(IN) :: y - REAL(SP), DIMENSION(:), INTENT(OUT) :: dydx - END SUBROUTINE derivs - END INTERFACE - END SUBROUTINE stifbs - END INTERFACE - INTERFACE - SUBROUTINE stiff(y,dydx,x,htry,eps,yscal,hdid,hnext,derivs) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(INOUT) :: y - REAL(SP), DIMENSION(:), INTENT(IN) :: dydx,yscal - REAL(SP), INTENT(INOUT) :: x - REAL(SP), INTENT(IN) :: htry,eps - REAL(SP), INTENT(OUT) :: hdid,hnext - INTERFACE - SUBROUTINE derivs(x,y,dydx) - USE nrtype - REAL(SP), INTENT(IN) :: x - REAL(SP), DIMENSION(:), INTENT(IN) :: y - REAL(SP), DIMENSION(:), INTENT(OUT) :: dydx - END SUBROUTINE derivs - END INTERFACE - END SUBROUTINE stiff - END INTERFACE - INTERFACE - SUBROUTINE stoerm(y,d2y,xs,htot,nstep,yout,derivs) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(IN) :: y,d2y - REAL(SP), INTENT(IN) :: xs,htot - INTEGER(I4B), INTENT(IN) :: nstep - REAL(SP), DIMENSION(:), INTENT(OUT) :: yout - INTERFACE - SUBROUTINE derivs(x,y,dydx) - USE nrtype - REAL(SP), INTENT(IN) :: x - REAL(SP), DIMENSION(:), INTENT(IN) :: y - REAL(SP), DIMENSION(:), INTENT(OUT) :: dydx - END SUBROUTINE derivs - END INTERFACE - END SUBROUTINE stoerm - END INTERFACE - INTERFACE svbksb -! SUBROUTINE svbksb_dp(u,w,v,b,x) -! USE nrtype -! REAL(DP), DIMENSION(:,:), INTENT(IN) :: u,v -! REAL(DP), DIMENSION(:), INTENT(IN) :: w,b -! REAL(DP), DIMENSION(:), INTENT(OUT) :: x -! END SUBROUTINE svbksb_dp -!BL - SUBROUTINE svbksb_sp(u,w,v,b,x) - USE nrtype - REAL(SP), DIMENSION(:,:), INTENT(IN) :: u,v - REAL(SP), DIMENSION(:), INTENT(IN) :: w,b - REAL(SP), DIMENSION(:), INTENT(OUT) :: x - END SUBROUTINE svbksb_sp - END INTERFACE - INTERFACE svdcmp -! SUBROUTINE svdcmp_dp(a,w,v) -! USE nrtype -! REAL(DP), DIMENSION(:,:), INTENT(INOUT) :: a -! REAL(DP), DIMENSION(:), INTENT(OUT) :: w -! REAL(DP), DIMENSION(:,:), INTENT(OUT) :: v -! END SUBROUTINE svdcmp_dp -!BL - SUBROUTINE svdcmp_sp(a,w,v) - USE nrtype - REAL(SP), DIMENSION(:,:), INTENT(INOUT) :: a - REAL(SP), DIMENSION(:), INTENT(OUT) :: w - REAL(SP), DIMENSION(:,:), INTENT(OUT) :: v - END SUBROUTINE svdcmp_sp - END INTERFACE - INTERFACE - SUBROUTINE svdfit(x,y,sig,a,v,w,chisq,funcs) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(IN) :: x,y,sig - REAL(SP), DIMENSION(:), INTENT(OUT) :: a,w - REAL(SP), DIMENSION(:,:), INTENT(OUT) :: v - REAL(SP), INTENT(OUT) :: chisq - INTERFACE - FUNCTION funcs(x,n) - USE nrtype - REAL(SP), INTENT(IN) :: x - INTEGER(I4B), INTENT(IN) :: n - REAL(SP), DIMENSION(n) :: funcs - END FUNCTION funcs - END INTERFACE - END SUBROUTINE svdfit - END INTERFACE - INTERFACE - SUBROUTINE svdvar(v,w,cvm) - USE nrtype - REAL(SP), DIMENSION(:,:), INTENT(IN) :: v - REAL(SP), DIMENSION(:), INTENT(IN) :: w - REAL(SP), DIMENSION(:,:), INTENT(OUT) :: cvm - END SUBROUTINE svdvar - END INTERFACE - INTERFACE - FUNCTION toeplz(r,y) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(IN) :: r,y - REAL(SP), DIMENSION(size(y)) :: toeplz - END FUNCTION toeplz - END INTERFACE - INTERFACE - SUBROUTINE tptest(data1,data2,t,prob) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(IN) :: data1,data2 - REAL(SP), INTENT(OUT) :: t,prob - END SUBROUTINE tptest - END INTERFACE - INTERFACE - SUBROUTINE tqli(d,e,z) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(INOUT) :: d,e - REAL(SP), DIMENSION(:,:), OPTIONAL, INTENT(INOUT) :: z - END SUBROUTINE tqli - END INTERFACE - INTERFACE - SUBROUTINE trapzd(func,a,b,s,n) - USE nrtype - REAL(SP), INTENT(IN) :: a,b - REAL(SP), INTENT(INOUT) :: s - INTEGER(I4B), INTENT(IN) :: n - INTERFACE - FUNCTION func(x) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(IN) :: x - REAL(SP), DIMENSION(size(x)) :: func - END FUNCTION func - END INTERFACE - END SUBROUTINE trapzd - END INTERFACE - INTERFACE - SUBROUTINE tred2(a,d,e,novectors) - USE nrtype - REAL(SP), DIMENSION(:,:), INTENT(INOUT) :: a - REAL(SP), DIMENSION(:), INTENT(OUT) :: d,e - LOGICAL(LGT), OPTIONAL, INTENT(IN) :: novectors - END SUBROUTINE tred2 - END INTERFACE -! On a purely serial machine, for greater efficiency, remove -! the generic name tridag from the following interface, -! and put it on the next one after that. - INTERFACE tridag - RECURSIVE SUBROUTINE tridag_par(a,b,c,r,u) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(IN) :: a,b,c,r - REAL(SP), DIMENSION(:), INTENT(OUT) :: u - END SUBROUTINE tridag_par - END INTERFACE - INTERFACE - SUBROUTINE tridag_ser(a,b,c,r,u) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(IN) :: a,b,c,r - REAL(SP), DIMENSION(:), INTENT(OUT) :: u - END SUBROUTINE tridag_ser - END INTERFACE - INTERFACE - SUBROUTINE ttest(data1,data2,t,prob) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(IN) :: data1,data2 - REAL(SP), INTENT(OUT) :: t,prob - END SUBROUTINE ttest - END INTERFACE - INTERFACE - SUBROUTINE tutest(data1,data2,t,prob) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(IN) :: data1,data2 - REAL(SP), INTENT(OUT) :: t,prob - END SUBROUTINE tutest - END INTERFACE - INTERFACE - SUBROUTINE twofft(data1,data2,fft1,fft2) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(IN) :: data1,data2 - COMPLEX(SPC), DIMENSION(:), INTENT(OUT) :: fft1,fft2 - END SUBROUTINE twofft - END INTERFACE - INTERFACE - FUNCTION vander(x,q) - USE nrtype - REAL(DP), DIMENSION(:), INTENT(IN) :: x,q - REAL(DP), DIMENSION(size(x)) :: vander - END FUNCTION vander - END INTERFACE - INTERFACE - SUBROUTINE vegas(region,func,init,ncall,itmx,nprn,tgral,sd,chi2a) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(IN) :: region - INTEGER(I4B), INTENT(IN) :: init,ncall,itmx,nprn - REAL(SP), INTENT(OUT) :: tgral,sd,chi2a - INTERFACE - FUNCTION func(pt,wgt) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(IN) :: pt - REAL(SP), INTENT(IN) :: wgt - REAL(SP) :: func - END FUNCTION func - END INTERFACE - END SUBROUTINE vegas - END INTERFACE - INTERFACE - SUBROUTINE voltra(t0,h,t,f,g,ak) - USE nrtype - REAL(SP), INTENT(IN) :: t0,h - REAL(SP), DIMENSION(:), INTENT(OUT) :: t - REAL(SP), DIMENSION(:,:), INTENT(OUT) :: f - INTERFACE - FUNCTION g(t) - USE nrtype - REAL(SP), INTENT(IN) :: t - REAL(SP), DIMENSION(:), POINTER :: g - END FUNCTION g -!BL - FUNCTION ak(t,s) - USE nrtype - REAL(SP), INTENT(IN) :: t,s - REAL(SP), DIMENSION(:,:), POINTER :: ak - END FUNCTION ak - END INTERFACE - END SUBROUTINE voltra - END INTERFACE - INTERFACE - SUBROUTINE wt1(a,isign,wtstep) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(INOUT) :: a - INTEGER(I4B), INTENT(IN) :: isign - INTERFACE - SUBROUTINE wtstep(a,isign) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(INOUT) :: a - INTEGER(I4B), INTENT(IN) :: isign - END SUBROUTINE wtstep - END INTERFACE - END SUBROUTINE wt1 - END INTERFACE - INTERFACE - SUBROUTINE wtn(a,nn,isign,wtstep) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(INOUT) :: a - INTEGER(I4B), DIMENSION(:), INTENT(IN) :: nn - INTEGER(I4B), INTENT(IN) :: isign - INTERFACE - SUBROUTINE wtstep(a,isign) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(INOUT) :: a - INTEGER(I4B), INTENT(IN) :: isign - END SUBROUTINE wtstep - END INTERFACE - END SUBROUTINE wtn - END INTERFACE - INTERFACE - FUNCTION wwghts(n,h,kermom) - USE nrtype - INTEGER(I4B), INTENT(IN) :: n - REAL(SP), INTENT(IN) :: h - REAL(SP), DIMENSION(n) :: wwghts - INTERFACE - FUNCTION kermom(y,m) - USE nrtype - REAL(DP), INTENT(IN) :: y - INTEGER(I4B), INTENT(IN) :: m - REAL(DP), DIMENSION(m) :: kermom - END FUNCTION kermom - END INTERFACE - END FUNCTION wwghts - END INTERFACE - INTERFACE - SUBROUTINE zbrac(func,x1,x2,succes) - USE nrtype - REAL(SP), INTENT(INOUT) :: x1,x2 - LOGICAL(LGT), INTENT(OUT) :: succes - INTERFACE - FUNCTION func(x) - USE nrtype - REAL(SP), INTENT(IN) :: x - REAL(SP) :: func - END FUNCTION func - END INTERFACE - END SUBROUTINE zbrac - END INTERFACE - INTERFACE - SUBROUTINE zbrak(func,x1,x2,n,xb1,xb2,nb) - USE nrtype - INTEGER(I4B), INTENT(IN) :: n - INTEGER(I4B), INTENT(OUT) :: nb - REAL(SP), INTENT(IN) :: x1,x2 - REAL(SP), DIMENSION(:), POINTER :: xb1,xb2 - INTERFACE - FUNCTION func(x) - USE nrtype - REAL(SP), INTENT(IN) :: x - REAL(SP) :: func - END FUNCTION func - END INTERFACE - END SUBROUTINE zbrak - END INTERFACE - INTERFACE - FUNCTION zbrent(func,x1,x2,tol) - USE nrtype - REAL(SP), INTENT(IN) :: x1,x2,tol - REAL(SP) :: zbrent - INTERFACE - FUNCTION func(x) - USE nrtype - REAL(SP), INTENT(IN) :: x - REAL(SP) :: func - END FUNCTION func - END INTERFACE - END FUNCTION zbrent - END INTERFACE - INTERFACE - SUBROUTINE zrhqr(a,rtr,rti) - USE nrtype - REAL(SP), DIMENSION(:), INTENT(IN) :: a - REAL(SP), DIMENSION(:), INTENT(OUT) :: rtr,rti - END SUBROUTINE zrhqr - END INTERFACE - INTERFACE - FUNCTION zriddr(func,x1,x2,xacc) - USE nrtype - REAL(SP), INTENT(IN) :: x1,x2,xacc - REAL(SP) :: zriddr - INTERFACE - FUNCTION func(x) - USE nrtype - REAL(SP), INTENT(IN) :: x - REAL(SP) :: func - END FUNCTION func - END INTERFACE - END FUNCTION zriddr - END INTERFACE - INTERFACE - SUBROUTINE zroots(a,roots,polish) - USE nrtype - COMPLEX(SPC), DIMENSION(:), INTENT(IN) :: a - COMPLEX(SPC), DIMENSION(:), INTENT(OUT) :: roots - LOGICAL(LGT), INTENT(IN) :: polish - END SUBROUTINE zroots - END INTERFACE -END MODULE nr diff --git a/nrtype.f90 b/nrtype.f90 deleted file mode 100755 index bda9374..0000000 --- a/nrtype.f90 +++ /dev/null @@ -1,30 +0,0 @@ -MODULE nrtype - INTEGER, PARAMETER :: I4B = SELECTED_INT_KIND(9) - INTEGER, PARAMETER :: I2B = SELECTED_INT_KIND(4) - INTEGER, PARAMETER :: I1B = SELECTED_INT_KIND(2) - INTEGER, PARAMETER :: SP = KIND(1.0d0) ! TO MAKE SINGLE, REMOVE THE `D0'' - INTEGER, PARAMETER :: DP = KIND(1.0D0) - INTEGER, PARAMETER :: SPC = KIND((1.0,1.0)) - INTEGER, PARAMETER :: DPC = KIND((1.0D0,1.0D0)) - INTEGER, PARAMETER :: LGT = KIND(.true.) - REAL(SP), PARAMETER :: PI=3.141592653589793238462643383279502884197_sp - REAL(SP), PARAMETER :: PIO2=1.57079632679489661923132169163975144209858_sp - REAL(SP), PARAMETER :: TWOPI=6.283185307179586476925286766559005768394_sp - REAL(SP), PARAMETER :: SQRT2=1.41421356237309504880168872420969807856967_sp - REAL(SP), PARAMETER :: EULER=0.5772156649015328606065120900824024310422_sp - REAL(DP), PARAMETER :: PI_D=3.141592653589793238462643383279502884197_dp - REAL(DP), PARAMETER :: PIO2_D=1.57079632679489661923132169163975144209858_dp - REAL(DP), PARAMETER :: TWOPI_D=6.283185307179586476925286766559005768394_dp - TYPE sprs2_sp - INTEGER(I4B) :: n,len - REAL(SP), DIMENSION(:), POINTER :: val - INTEGER(I4B), DIMENSION(:), POINTER :: irow - INTEGER(I4B), DIMENSION(:), POINTER :: jcol - END TYPE sprs2_sp - TYPE sprs2_dp - INTEGER(I4B) :: n,len - REAL(DP), DIMENSION(:), POINTER :: val - INTEGER(I4B), DIMENSION(:), POINTER :: irow - INTEGER(I4B), DIMENSION(:), POINTER :: jcol - END TYPE sprs2_dp -END MODULE nrtype diff --git a/nrutil.f90 b/nrutil.f90 deleted file mode 100755 index b9dd9ae..0000000 --- a/nrutil.f90 +++ /dev/null @@ -1,1164 +0,0 @@ -MODULE nrutil - USE nrtype - IMPLICIT NONE - INTEGER(I4B), PARAMETER :: NPAR_ARTH=16,NPAR2_ARTH=8 - INTEGER(I4B), PARAMETER :: NPAR_GEOP=4,NPAR2_GEOP=2 - INTEGER(I4B), PARAMETER :: NPAR_CUMSUM=16 - INTEGER(I4B), PARAMETER :: NPAR_CUMPROD=8 - INTEGER(I4B), PARAMETER :: NPAR_POLY=8 - INTEGER(I4B), PARAMETER :: NPAR_POLYTERM=8 - INTERFACE array_copy -! MODULE PROCEDURE array_copy_r, array_copy_d, array_copy_i - MODULE PROCEDURE array_copy_r, array_copy_i - END INTERFACE - INTERFACE swap - MODULE PROCEDURE swap_i,swap_r,swap_rv,swap_c, & - swap_cv,swap_cm,swap_z,swap_zv,swap_zm, & - masked_swap_rs,masked_swap_rv,masked_swap_rm - END INTERFACE - INTERFACE reallocate - MODULE PROCEDURE reallocate_rv,reallocate_rm,& - reallocate_iv,reallocate_im,reallocate_hv - END INTERFACE - INTERFACE imaxloc - MODULE PROCEDURE imaxloc_r,imaxloc_i - END INTERFACE - INTERFACE assert - MODULE PROCEDURE assert1,assert2,assert3,assert4,assert_v - END INTERFACE - INTERFACE assert_eq - MODULE PROCEDURE assert_eq2,assert_eq3,assert_eq4,assert_eqn - END INTERFACE - INTERFACE arth -! MODULE PROCEDURE arth_r, arth_d, arth_i - MODULE PROCEDURE arth_r, arth_i - END INTERFACE - INTERFACE geop -! MODULE PROCEDURE geop_r, geop_d, geop_i, geop_c, geop_dv - MODULE PROCEDURE geop_r, geop_i, geop_c, geop_dv - END INTERFACE - INTERFACE cumsum - MODULE PROCEDURE cumsum_r,cumsum_i - END INTERFACE - INTERFACE poly -! MODULE PROCEDURE poly_rr,poly_rrv,poly_dd,poly_ddv,& -! poly_rc,poly_cc,poly_msk_rrv,poly_msk_ddv - MODULE PROCEDURE poly_rr,poly_rrv,& - poly_rc,poly_cc,poly_msk_rrv - END INTERFACE - INTERFACE poly_term - MODULE PROCEDURE poly_term_rr,poly_term_cc - END INTERFACE - INTERFACE outerprod -! MODULE PROCEDURE outerprod_r,outerprod_d - MODULE PROCEDURE outerprod_r - END INTERFACE - INTERFACE outerdiff -! MODULE PROCEDURE outerdiff_r,outerdiff_d,outerdiff_i - MODULE PROCEDURE outerdiff_r,outerdiff_i - END INTERFACE - INTERFACE scatter_add -! MODULE PROCEDURE scatter_add_r,scatter_add_d - MODULE PROCEDURE scatter_add_r - END INTERFACE - INTERFACE scatter_max -! MODULE PROCEDURE scatter_max_r,scatter_max_d - MODULE PROCEDURE scatter_max_r - END INTERFACE - INTERFACE diagadd - MODULE PROCEDURE diagadd_rv,diagadd_r - END INTERFACE - INTERFACE diagmult - MODULE PROCEDURE diagmult_rv,diagmult_r - END INTERFACE - INTERFACE get_diag -! MODULE PROCEDURE get_diag_rv, get_diag_dv - MODULE PROCEDURE get_diag_rv - END INTERFACE - INTERFACE put_diag - MODULE PROCEDURE put_diag_rv, put_diag_r - END INTERFACE -CONTAINS -!BL - SUBROUTINE array_copy_r(src,dest,n_copied,n_not_copied) - REAL(SP), DIMENSION(:), INTENT(IN) :: src - REAL(SP), DIMENSION(:), INTENT(OUT) :: dest - INTEGER(I4B), INTENT(OUT) :: n_copied, n_not_copied - n_copied=min(size(src),size(dest)) - n_not_copied=size(src)-n_copied - dest(1:n_copied)=src(1:n_copied) - END SUBROUTINE array_copy_r -!BL - SUBROUTINE array_copy_d(src,dest,n_copied,n_not_copied) - REAL(DP), DIMENSION(:), INTENT(IN) :: src - REAL(DP), DIMENSION(:), INTENT(OUT) :: dest - INTEGER(I4B), INTENT(OUT) :: n_copied, n_not_copied - n_copied=min(size(src),size(dest)) - n_not_copied=size(src)-n_copied - dest(1:n_copied)=src(1:n_copied) - END SUBROUTINE array_copy_d -!BL - SUBROUTINE array_copy_i(src,dest,n_copied,n_not_copied) - INTEGER(I4B), DIMENSION(:), INTENT(IN) :: src - INTEGER(I4B), DIMENSION(:), INTENT(OUT) :: dest - INTEGER(I4B), INTENT(OUT) :: n_copied, n_not_copied - n_copied=min(size(src),size(dest)) - n_not_copied=size(src)-n_copied - dest(1:n_copied)=src(1:n_copied) - END SUBROUTINE array_copy_i -!BL -!BL - SUBROUTINE swap_i(a,b) - INTEGER(I4B), INTENT(INOUT) :: a,b - INTEGER(I4B) :: dum - dum=a - a=b - b=dum - END SUBROUTINE swap_i -!BL - SUBROUTINE swap_r(a,b) - REAL(SP), INTENT(INOUT) :: a,b - REAL(SP) :: dum - dum=a - a=b - b=dum - END SUBROUTINE swap_r -!BL - SUBROUTINE swap_rv(a,b) - REAL(SP), DIMENSION(:), INTENT(INOUT) :: a,b - REAL(SP), DIMENSION(SIZE(a)) :: dum - dum=a - a=b - b=dum - END SUBROUTINE swap_rv -!BL - SUBROUTINE swap_c(a,b) - COMPLEX(SPC), INTENT(INOUT) :: a,b - COMPLEX(SPC) :: dum - dum=a - a=b - b=dum - END SUBROUTINE swap_c -!BL - SUBROUTINE swap_cv(a,b) - COMPLEX(SPC), DIMENSION(:), INTENT(INOUT) :: a,b - COMPLEX(SPC), DIMENSION(SIZE(a)) :: dum - dum=a - a=b - b=dum - END SUBROUTINE swap_cv -!BL - SUBROUTINE swap_cm(a,b) - COMPLEX(SPC), DIMENSION(:,:), INTENT(INOUT) :: a,b - COMPLEX(SPC), DIMENSION(size(a,1),size(a,2)) :: dum - dum=a - a=b - b=dum - END SUBROUTINE swap_cm -!BL - SUBROUTINE swap_z(a,b) - COMPLEX(DPC), INTENT(INOUT) :: a,b - COMPLEX(DPC) :: dum - dum=a - a=b - b=dum - END SUBROUTINE swap_z -!BL - SUBROUTINE swap_zv(a,b) - COMPLEX(DPC), DIMENSION(:), INTENT(INOUT) :: a,b - COMPLEX(DPC), DIMENSION(SIZE(a)) :: dum - dum=a - a=b - b=dum - END SUBROUTINE swap_zv -!BL - SUBROUTINE swap_zm(a,b) - COMPLEX(DPC), DIMENSION(:,:), INTENT(INOUT) :: a,b - COMPLEX(DPC), DIMENSION(size(a,1),size(a,2)) :: dum - dum=a - a=b - b=dum - END SUBROUTINE swap_zm -!BL - SUBROUTINE masked_swap_rs(a,b,mask) - REAL(SP), INTENT(INOUT) :: a,b - LOGICAL(LGT), INTENT(IN) :: mask - REAL(SP) :: swp - if (mask) then - swp=a - a=b - b=swp - end if - END SUBROUTINE masked_swap_rs -!BL - SUBROUTINE masked_swap_rv(a,b,mask) - REAL(SP), DIMENSION(:), INTENT(INOUT) :: a,b - LOGICAL(LGT), DIMENSION(:), INTENT(IN) :: mask - REAL(SP), DIMENSION(size(a)) :: swp - where (mask) - swp=a - a=b - b=swp - end where - END SUBROUTINE masked_swap_rv -!BL - SUBROUTINE masked_swap_rm(a,b,mask) - REAL(SP), DIMENSION(:,:), INTENT(INOUT) :: a,b - LOGICAL(LGT), DIMENSION(:,:), INTENT(IN) :: mask - REAL(SP), DIMENSION(size(a,1),size(a,2)) :: swp - where (mask) - swp=a - a=b - b=swp - end where - END SUBROUTINE masked_swap_rm -!BL -!BL - FUNCTION reallocate_rv(p,n) - REAL(SP), DIMENSION(:), POINTER :: p, reallocate_rv - INTEGER(I4B), INTENT(IN) :: n - INTEGER(I4B) :: nold,ierr - allocate(reallocate_rv(n),stat=ierr) - if (ierr /= 0) call & - nrerror('reallocate_rv: problem in attempt to allocate memory') - if (.not. associated(p)) RETURN - nold=size(p) - reallocate_rv(1:min(nold,n))=p(1:min(nold,n)) - deallocate(p) - END FUNCTION reallocate_rv -!BL - FUNCTION reallocate_iv(p,n) - INTEGER(I4B), DIMENSION(:), POINTER :: p, reallocate_iv - INTEGER(I4B), INTENT(IN) :: n - INTEGER(I4B) :: nold,ierr - allocate(reallocate_iv(n),stat=ierr) - if (ierr /= 0) call & - nrerror('reallocate_iv: problem in attempt to allocate memory') - if (.not. associated(p)) RETURN - nold=size(p) - reallocate_iv(1:min(nold,n))=p(1:min(nold,n)) - deallocate(p) - END FUNCTION reallocate_iv -!BL - FUNCTION reallocate_hv(p,n) - CHARACTER(1), DIMENSION(:), POINTER :: p, reallocate_hv - INTEGER(I4B), INTENT(IN) :: n - INTEGER(I4B) :: nold,ierr - allocate(reallocate_hv(n),stat=ierr) - if (ierr /= 0) call & - nrerror('reallocate_hv: problem in attempt to allocate memory') - if (.not. associated(p)) RETURN - nold=size(p) - reallocate_hv(1:min(nold,n))=p(1:min(nold,n)) - deallocate(p) - END FUNCTION reallocate_hv -!BL - FUNCTION reallocate_rm(p,n,m) - REAL(SP), DIMENSION(:,:), POINTER :: p, reallocate_rm - INTEGER(I4B), INTENT(IN) :: n,m - INTEGER(I4B) :: nold,mold,ierr - allocate(reallocate_rm(n,m),stat=ierr) - if (ierr /= 0) call & - nrerror('reallocate_rm: problem in attempt to allocate memory') - if (.not. associated(p)) RETURN - nold=size(p,1) - mold=size(p,2) - reallocate_rm(1:min(nold,n),1:min(mold,m))=& - p(1:min(nold,n),1:min(mold,m)) - deallocate(p) - END FUNCTION reallocate_rm -!BL - FUNCTION reallocate_im(p,n,m) - INTEGER(I4B), DIMENSION(:,:), POINTER :: p, reallocate_im - INTEGER(I4B), INTENT(IN) :: n,m - INTEGER(I4B) :: nold,mold,ierr - allocate(reallocate_im(n,m),stat=ierr) - if (ierr /= 0) call & - nrerror('reallocate_im: problem in attempt to allocate memory') - if (.not. associated(p)) RETURN - nold=size(p,1) - mold=size(p,2) - reallocate_im(1:min(nold,n),1:min(mold,m))=& - p(1:min(nold,n),1:min(mold,m)) - deallocate(p) - END FUNCTION reallocate_im -!BL - FUNCTION ifirstloc(mask) - LOGICAL(LGT), DIMENSION(:), INTENT(IN) :: mask - INTEGER(I4B) :: ifirstloc - INTEGER(I4B), DIMENSION(1) :: loc - loc=maxloc(merge(1,0,mask)) - ifirstloc=loc(1) - if (.not. mask(ifirstloc)) ifirstloc=size(mask)+1 - END FUNCTION ifirstloc -!BL - FUNCTION imaxloc_r(arr) - REAL(SP), DIMENSION(:), INTENT(IN) :: arr - INTEGER(I4B) :: imaxloc_r - INTEGER(I4B), DIMENSION(1) :: imax - imax=maxloc(arr(:)) - imaxloc_r=imax(1) - END FUNCTION imaxloc_r -!BL - FUNCTION imaxloc_i(iarr) - INTEGER(I4B), DIMENSION(:), INTENT(IN) :: iarr - INTEGER(I4B), DIMENSION(1) :: imax - INTEGER(I4B) :: imaxloc_i - imax=maxloc(iarr(:)) - imaxloc_i=imax(1) - END FUNCTION imaxloc_i -!BL - FUNCTION iminloc(arr) - REAL(SP), DIMENSION(:), INTENT(IN) :: arr - INTEGER(I4B), DIMENSION(1) :: imin - INTEGER(I4B) :: iminloc - imin=minloc(arr(:)) - iminloc=imin(1) - END FUNCTION iminloc -!BL - SUBROUTINE assert1(n1,string) - CHARACTER(LEN=*), INTENT(IN) :: string - LOGICAL, INTENT(IN) :: n1 - if (.not. n1) then - write (*,*) 'nrerror: an assertion failed with this tag:', & - string - STOP 'program terminated by assert1' - end if - END SUBROUTINE assert1 -!BL - SUBROUTINE assert2(n1,n2,string) - CHARACTER(LEN=*), INTENT(IN) :: string - LOGICAL, INTENT(IN) :: n1,n2 - if (.not. (n1 .and. n2)) then - write (*,*) 'nrerror: an assertion failed with this tag:', & - string - STOP 'program terminated by assert2' - end if - END SUBROUTINE assert2 -!BL - SUBROUTINE assert3(n1,n2,n3,string) - CHARACTER(LEN=*), INTENT(IN) :: string - LOGICAL, INTENT(IN) :: n1,n2,n3 - if (.not. (n1 .and. n2 .and. n3)) then - write (*,*) 'nrerror: an assertion failed with this tag:', & - string - STOP 'program terminated by assert3' - end if - END SUBROUTINE assert3 -!BL - SUBROUTINE assert4(n1,n2,n3,n4,string) - CHARACTER(LEN=*), INTENT(IN) :: string - LOGICAL, INTENT(IN) :: n1,n2,n3,n4 - if (.not. (n1 .and. n2 .and. n3 .and. n4)) then - write (*,*) 'nrerror: an assertion failed with this tag:', & - string - STOP 'program terminated by assert4' - end if - END SUBROUTINE assert4 -!BL - SUBROUTINE assert_v(n,string) - CHARACTER(LEN=*), INTENT(IN) :: string - LOGICAL, DIMENSION(:), INTENT(IN) :: n - if (.not. all(n)) then - write (*,*) 'nrerror: an assertion failed with this tag:', & - string - STOP 'program terminated by assert_v' - end if - END SUBROUTINE assert_v -!BL - FUNCTION assert_eq2(n1,n2,string) - CHARACTER(LEN=*), INTENT(IN) :: string - INTEGER, INTENT(IN) :: n1,n2 - INTEGER :: assert_eq2 - if (n1 == n2) then - assert_eq2=n1 - else - write (*,*) 'nrerror: an assert_eq failed with this tag:', & - string - STOP 'program terminated by assert_eq2' - end if - END FUNCTION assert_eq2 -!BL - FUNCTION assert_eq3(n1,n2,n3,string) - CHARACTER(LEN=*), INTENT(IN) :: string - INTEGER, INTENT(IN) :: n1,n2,n3 - INTEGER :: assert_eq3 - if (n1 == n2 .and. n2 == n3) then - assert_eq3=n1 - else - write (*,*) 'nrerror: an assert_eq failed with this tag:', & - string - STOP 'program terminated by assert_eq3' - end if - END FUNCTION assert_eq3 -!BL - FUNCTION assert_eq4(n1,n2,n3,n4,string) - CHARACTER(LEN=*), INTENT(IN) :: string - INTEGER, INTENT(IN) :: n1,n2,n3,n4 - INTEGER :: assert_eq4 - if (n1 == n2 .and. n2 == n3 .and. n3 == n4) then - assert_eq4=n1 - else - write (*,*) 'nrerror: an assert_eq failed with this tag:', & - string - STOP 'program terminated by assert_eq4' - end if - END FUNCTION assert_eq4 -!BL - FUNCTION assert_eqn(nn,string) - CHARACTER(LEN=*), INTENT(IN) :: string - INTEGER, DIMENSION(:), INTENT(IN) :: nn - INTEGER :: assert_eqn - if (all(nn(2:) == nn(1))) then - assert_eqn=nn(1) - else - write (*,*) 'nrerror: an assert_eq failed with this tag:', & - string - STOP 'program terminated by assert_eqn' - end if - END FUNCTION assert_eqn -!BL - SUBROUTINE nrerror(string) - CHARACTER(LEN=*), INTENT(IN) :: string - write (*,*) 'nrerror: ',string - STOP 'program terminated by nrerror' - END SUBROUTINE nrerror -!BL - FUNCTION arth_r(first,increment,n) - REAL(SP), INTENT(IN) :: first,increment - INTEGER(I4B), INTENT(IN) :: n - REAL(SP), DIMENSION(n) :: arth_r - INTEGER(I4B) :: k,k2 - REAL(SP) :: temp - if (n > 0) arth_r(1)=first - if (n <= NPAR_ARTH) then - do k=2,n - arth_r(k)=arth_r(k-1)+increment - end do - else - do k=2,NPAR2_ARTH - arth_r(k)=arth_r(k-1)+increment - end do - temp=increment*NPAR2_ARTH - k=NPAR2_ARTH - do - if (k >= n) exit - k2=k+k - arth_r(k+1:min(k2,n))=temp+arth_r(1:min(k,n-k)) - temp=temp+temp - k=k2 - end do - end if - END FUNCTION arth_r -!BL - FUNCTION arth_d(first,increment,n) - REAL(DP), INTENT(IN) :: first,increment - INTEGER(I4B), INTENT(IN) :: n - REAL(DP), DIMENSION(n) :: arth_d - INTEGER(I4B) :: k,k2 - REAL(DP) :: temp - if (n > 0) arth_d(1)=first - if (n <= NPAR_ARTH) then - do k=2,n - arth_d(k)=arth_d(k-1)+increment - end do - else - do k=2,NPAR2_ARTH - arth_d(k)=arth_d(k-1)+increment - end do - temp=increment*NPAR2_ARTH - k=NPAR2_ARTH - do - if (k >= n) exit - k2=k+k - arth_d(k+1:min(k2,n))=temp+arth_d(1:min(k,n-k)) - temp=temp+temp - k=k2 - end do - end if - END FUNCTION arth_d -!BL - FUNCTION arth_i(first,increment,n) - INTEGER(I4B), INTENT(IN) :: first,increment,n - INTEGER(I4B), DIMENSION(n) :: arth_i - INTEGER(I4B) :: k,k2,temp - if (n > 0) arth_i(1)=first - if (n <= NPAR_ARTH) then - do k=2,n - arth_i(k)=arth_i(k-1)+increment - end do - else - do k=2,NPAR2_ARTH - arth_i(k)=arth_i(k-1)+increment - end do - temp=increment*NPAR2_ARTH - k=NPAR2_ARTH - do - if (k >= n) exit - k2=k+k - arth_i(k+1:min(k2,n))=temp+arth_i(1:min(k,n-k)) - temp=temp+temp - k=k2 - end do - end if - END FUNCTION arth_i -!BL -!BL - FUNCTION geop_r(first,factor,n) - REAL(SP), INTENT(IN) :: first,factor - INTEGER(I4B), INTENT(IN) :: n - REAL(SP), DIMENSION(n) :: geop_r - INTEGER(I4B) :: k,k2 - REAL(SP) :: temp - if (n > 0) geop_r(1)=first - if (n <= NPAR_GEOP) then - do k=2,n - geop_r(k)=geop_r(k-1)*factor - end do - else - do k=2,NPAR2_GEOP - geop_r(k)=geop_r(k-1)*factor - end do - temp=factor**NPAR2_GEOP - k=NPAR2_GEOP - do - if (k >= n) exit - k2=k+k - geop_r(k+1:min(k2,n))=temp*geop_r(1:min(k,n-k)) - temp=temp*temp - k=k2 - end do - end if - END FUNCTION geop_r -!BL - FUNCTION geop_d(first,factor,n) - REAL(DP), INTENT(IN) :: first,factor - INTEGER(I4B), INTENT(IN) :: n - REAL(DP), DIMENSION(n) :: geop_d - INTEGER(I4B) :: k,k2 - REAL(DP) :: temp - if (n > 0) geop_d(1)=first - if (n <= NPAR_GEOP) then - do k=2,n - geop_d(k)=geop_d(k-1)*factor - end do - else - do k=2,NPAR2_GEOP - geop_d(k)=geop_d(k-1)*factor - end do - temp=factor**NPAR2_GEOP - k=NPAR2_GEOP - do - if (k >= n) exit - k2=k+k - geop_d(k+1:min(k2,n))=temp*geop_d(1:min(k,n-k)) - temp=temp*temp - k=k2 - end do - end if - END FUNCTION geop_d -!BL - FUNCTION geop_i(first,factor,n) - INTEGER(I4B), INTENT(IN) :: first,factor,n - INTEGER(I4B), DIMENSION(n) :: geop_i - INTEGER(I4B) :: k,k2,temp - if (n > 0) geop_i(1)=first - if (n <= NPAR_GEOP) then - do k=2,n - geop_i(k)=geop_i(k-1)*factor - end do - else - do k=2,NPAR2_GEOP - geop_i(k)=geop_i(k-1)*factor - end do - temp=factor**NPAR2_GEOP - k=NPAR2_GEOP - do - if (k >= n) exit - k2=k+k - geop_i(k+1:min(k2,n))=temp*geop_i(1:min(k,n-k)) - temp=temp*temp - k=k2 - end do - end if - END FUNCTION geop_i -!BL - FUNCTION geop_c(first,factor,n) - COMPLEX(SP), INTENT(IN) :: first,factor - INTEGER(I4B), INTENT(IN) :: n - COMPLEX(SP), DIMENSION(n) :: geop_c - INTEGER(I4B) :: k,k2 - COMPLEX(SP) :: temp - if (n > 0) geop_c(1)=first - if (n <= NPAR_GEOP) then - do k=2,n - geop_c(k)=geop_c(k-1)*factor - end do - else - do k=2,NPAR2_GEOP - geop_c(k)=geop_c(k-1)*factor - end do - temp=factor**NPAR2_GEOP - k=NPAR2_GEOP - do - if (k >= n) exit - k2=k+k - geop_c(k+1:min(k2,n))=temp*geop_c(1:min(k,n-k)) - temp=temp*temp - k=k2 - end do - end if - END FUNCTION geop_c -!BL - FUNCTION geop_dv(first,factor,n) - REAL(DP), DIMENSION(:), INTENT(IN) :: first,factor - INTEGER(I4B), INTENT(IN) :: n - REAL(DP), DIMENSION(size(first),n) :: geop_dv - INTEGER(I4B) :: k,k2 - REAL(DP), DIMENSION(size(first)) :: temp - if (n > 0) geop_dv(:,1)=first(:) - if (n <= NPAR_GEOP) then - do k=2,n - geop_dv(:,k)=geop_dv(:,k-1)*factor(:) - end do - else - do k=2,NPAR2_GEOP - geop_dv(:,k)=geop_dv(:,k-1)*factor(:) - end do - temp=factor**NPAR2_GEOP - k=NPAR2_GEOP - do - if (k >= n) exit - k2=k+k - geop_dv(:,k+1:min(k2,n))=geop_dv(:,1:min(k,n-k))*& - spread(temp,2,size(geop_dv(:,1:min(k,n-k)),2)) - temp=temp*temp - k=k2 - end do - end if - END FUNCTION geop_dv -!BL -!BL - RECURSIVE FUNCTION cumsum_r(arr,seed) RESULT(ans) - REAL(SP), DIMENSION(:), INTENT(IN) :: arr - REAL(SP), OPTIONAL, INTENT(IN) :: seed - REAL(SP), DIMENSION(size(arr)) :: ans - INTEGER(I4B) :: n,j - REAL(SP) :: sd - n=size(arr) - if (n == 0_i4b) RETURN - sd=0.0_sp - if (present(seed)) sd=seed - ans(1)=arr(1)+sd - if (n < NPAR_CUMSUM) then - do j=2,n - ans(j)=ans(j-1)+arr(j) - end do - else - ans(2:n:2)=cumsum_r(arr(2:n:2)+arr(1:n-1:2),sd) - ans(3:n:2)=ans(2:n-1:2)+arr(3:n:2) - end if - END FUNCTION cumsum_r -!BL - RECURSIVE FUNCTION cumsum_i(arr,seed) RESULT(ans) - INTEGER(I4B), DIMENSION(:), INTENT(IN) :: arr - INTEGER(I4B), OPTIONAL, INTENT(IN) :: seed - INTEGER(I4B), DIMENSION(size(arr)) :: ans - INTEGER(I4B) :: n,j,sd - n=size(arr) - if (n == 0_i4b) RETURN - sd=0_i4b - if (present(seed)) sd=seed - ans(1)=arr(1)+sd - if (n < NPAR_CUMSUM) then - do j=2,n - ans(j)=ans(j-1)+arr(j) - end do - else - ans(2:n:2)=cumsum_i(arr(2:n:2)+arr(1:n-1:2),sd) - ans(3:n:2)=ans(2:n-1:2)+arr(3:n:2) - end if - END FUNCTION cumsum_i -!BL -!BL - RECURSIVE FUNCTION cumprod(arr,seed) RESULT(ans) - REAL(SP), DIMENSION(:), INTENT(IN) :: arr - REAL(SP), OPTIONAL, INTENT(IN) :: seed - REAL(SP), DIMENSION(size(arr)) :: ans - INTEGER(I4B) :: n,j - REAL(SP) :: sd - n=size(arr) - if (n == 0_i4b) RETURN - sd=1.0_sp - if (present(seed)) sd=seed - ans(1)=arr(1)*sd - if (n < NPAR_CUMPROD) then - do j=2,n - ans(j)=ans(j-1)*arr(j) - end do - else - ans(2:n:2)=cumprod(arr(2:n:2)*arr(1:n-1:2),sd) - ans(3:n:2)=ans(2:n-1:2)*arr(3:n:2) - end if - END FUNCTION cumprod -!BL -!BL - FUNCTION poly_rr(x,coeffs) - REAL(SP), INTENT(IN) :: x - REAL(SP), DIMENSION(:), INTENT(IN) :: coeffs - REAL(SP) :: poly_rr - REAL(SP) :: pow - REAL(SP), DIMENSION(:), ALLOCATABLE :: vec - INTEGER(I4B) :: i,n,nn - n=size(coeffs) - if (n <= 0) then - poly_rr=0.0_sp - else if (n < NPAR_POLY) then - poly_rr=coeffs(n) - do i=n-1,1,-1 - poly_rr=x*poly_rr+coeffs(i) - end do - else - allocate(vec(n+1)) - pow=x - vec(1:n)=coeffs - do - vec(n+1)=0.0_sp - nn=ishft(n+1,-1) - vec(1:nn)=vec(1:n:2)+pow*vec(2:n+1:2) - if (nn == 1) exit - pow=pow*pow - n=nn - end do - poly_rr=vec(1) - deallocate(vec) - end if - END FUNCTION poly_rr -!BL - FUNCTION poly_dd(x,coeffs) - REAL(DP), INTENT(IN) :: x - REAL(DP), DIMENSION(:), INTENT(IN) :: coeffs - REAL(DP) :: poly_dd - REAL(DP) :: pow - REAL(DP), DIMENSION(:), ALLOCATABLE :: vec - INTEGER(I4B) :: i,n,nn - n=size(coeffs) - if (n <= 0) then - poly_dd=0.0_dp - else if (n < NPAR_POLY) then - poly_dd=coeffs(n) - do i=n-1,1,-1 - poly_dd=x*poly_dd+coeffs(i) - end do - else - allocate(vec(n+1)) - pow=x - vec(1:n)=coeffs - do - vec(n+1)=0.0_dp - nn=ishft(n+1,-1) - vec(1:nn)=vec(1:n:2)+pow*vec(2:n+1:2) - if (nn == 1) exit - pow=pow*pow - n=nn - end do - poly_dd=vec(1) - deallocate(vec) - end if - END FUNCTION poly_dd -!BL - FUNCTION poly_rc(x,coeffs) - COMPLEX(SPC), INTENT(IN) :: x - REAL(SP), DIMENSION(:), INTENT(IN) :: coeffs - COMPLEX(SPC) :: poly_rc - COMPLEX(SPC) :: pow - COMPLEX(SPC), DIMENSION(:), ALLOCATABLE :: vec - INTEGER(I4B) :: i,n,nn - n=size(coeffs) - if (n <= 0) then - poly_rc=0.0_sp - else if (n < NPAR_POLY) then - poly_rc=coeffs(n) - do i=n-1,1,-1 - poly_rc=x*poly_rc+coeffs(i) - end do - else - allocate(vec(n+1)) - pow=x - vec(1:n)=coeffs - do - vec(n+1)=0.0_sp - nn=ishft(n+1,-1) - vec(1:nn)=vec(1:n:2)+pow*vec(2:n+1:2) - if (nn == 1) exit - pow=pow*pow - n=nn - end do - poly_rc=vec(1) - deallocate(vec) - end if - END FUNCTION poly_rc -!BL - FUNCTION poly_cc(x,coeffs) - COMPLEX(SPC), INTENT(IN) :: x - COMPLEX(SPC), DIMENSION(:), INTENT(IN) :: coeffs - COMPLEX(SPC) :: poly_cc - COMPLEX(SPC) :: pow - COMPLEX(SPC), DIMENSION(:), ALLOCATABLE :: vec - INTEGER(I4B) :: i,n,nn - n=size(coeffs) - if (n <= 0) then - poly_cc=0.0_sp - else if (n < NPAR_POLY) then - poly_cc=coeffs(n) - do i=n-1,1,-1 - poly_cc=x*poly_cc+coeffs(i) - end do - else - allocate(vec(n+1)) - pow=x - vec(1:n)=coeffs - do - vec(n+1)=0.0_sp - nn=ishft(n+1,-1) - vec(1:nn)=vec(1:n:2)+pow*vec(2:n+1:2) - if (nn == 1) exit - pow=pow*pow - n=nn - end do - poly_cc=vec(1) - deallocate(vec) - end if - END FUNCTION poly_cc -!BL - FUNCTION poly_rrv(x,coeffs) - REAL(SP), DIMENSION(:), INTENT(IN) :: coeffs,x - REAL(SP), DIMENSION(size(x)) :: poly_rrv - INTEGER(I4B) :: i,n,m - m=size(coeffs) - n=size(x) - if (m <= 0) then - poly_rrv=0.0_sp - else if (m < n .or. m < NPAR_POLY) then - poly_rrv=coeffs(m) - do i=m-1,1,-1 - poly_rrv=x*poly_rrv+coeffs(i) - end do - else - do i=1,n - poly_rrv(i)=poly_rr(x(i),coeffs) - end do - end if - END FUNCTION poly_rrv -!BL - FUNCTION poly_ddv(x,coeffs) - REAL(DP), DIMENSION(:), INTENT(IN) :: coeffs,x - REAL(DP), DIMENSION(size(x)) :: poly_ddv - INTEGER(I4B) :: i,n,m - m=size(coeffs) - n=size(x) - if (m <= 0) then - poly_ddv=0.0_dp - else if (m < n .or. m < NPAR_POLY) then - poly_ddv=coeffs(m) - do i=m-1,1,-1 - poly_ddv=x*poly_ddv+coeffs(i) - end do - else - do i=1,n - poly_ddv(i)=poly_dd(x(i),coeffs) - end do - end if - END FUNCTION poly_ddv -!BL - FUNCTION poly_msk_rrv(x,coeffs,mask) - REAL(SP), DIMENSION(:), INTENT(IN) :: coeffs,x - LOGICAL(LGT), DIMENSION(:), INTENT(IN) :: mask - REAL(SP), DIMENSION(size(x)) :: poly_msk_rrv - poly_msk_rrv=unpack(poly_rrv(pack(x,mask),coeffs),mask,0.0_sp) - END FUNCTION poly_msk_rrv -!BL - FUNCTION poly_msk_ddv(x,coeffs,mask) - REAL(DP), DIMENSION(:), INTENT(IN) :: coeffs,x - LOGICAL(LGT), DIMENSION(:), INTENT(IN) :: mask - REAL(DP), DIMENSION(size(x)) :: poly_msk_ddv - poly_msk_ddv=unpack(poly_ddv(pack(x,mask),coeffs),mask,0.0_dp) - END FUNCTION poly_msk_ddv -!BL -!BL - RECURSIVE FUNCTION poly_term_rr(a,b) RESULT(u) - REAL(SP), DIMENSION(:), INTENT(IN) :: a - REAL(SP), INTENT(IN) :: b - REAL(SP), DIMENSION(size(a)) :: u - INTEGER(I4B) :: n,j - n=size(a) - if (n <= 0) RETURN - u(1)=a(1) - if (n < NPAR_POLYTERM) then - do j=2,n - u(j)=a(j)+b*u(j-1) - end do - else - u(2:n:2)=poly_term_rr(a(2:n:2)+a(1:n-1:2)*b,b*b) - u(3:n:2)=a(3:n:2)+b*u(2:n-1:2) - end if - END FUNCTION poly_term_rr -!BL - RECURSIVE FUNCTION poly_term_cc(a,b) RESULT(u) - COMPLEX(SPC), DIMENSION(:), INTENT(IN) :: a - COMPLEX(SPC), INTENT(IN) :: b - COMPLEX(SPC), DIMENSION(size(a)) :: u - INTEGER(I4B) :: n,j - n=size(a) - if (n <= 0) RETURN - u(1)=a(1) - if (n < NPAR_POLYTERM) then - do j=2,n - u(j)=a(j)+b*u(j-1) - end do - else - u(2:n:2)=poly_term_cc(a(2:n:2)+a(1:n-1:2)*b,b*b) - u(3:n:2)=a(3:n:2)+b*u(2:n-1:2) - end if - END FUNCTION poly_term_cc -!BL -!BL - FUNCTION zroots_unity(n,nn) - INTEGER(I4B), INTENT(IN) :: n,nn - COMPLEX(SPC), DIMENSION(nn) :: zroots_unity - INTEGER(I4B) :: k - REAL(SP) :: theta - zroots_unity(1)=1.0 - theta=TWOPI/n - k=1 - do - if (k >= nn) exit - zroots_unity(k+1)=cmplx(cos(k*theta),sin(k*theta),SPC) - zroots_unity(k+2:min(2*k,nn))=zroots_unity(k+1)*& - zroots_unity(2:min(k,nn-k)) - k=2*k - end do - END FUNCTION zroots_unity -!BL - FUNCTION outerprod_r(a,b) - REAL(SP), DIMENSION(:), INTENT(IN) :: a,b - REAL(SP), DIMENSION(size(a),size(b)) :: outerprod_r - outerprod_r = spread(a,dim=2,ncopies=size(b)) * & - spread(b,dim=1,ncopies=size(a)) - END FUNCTION outerprod_r -!BL - FUNCTION outerprod_d(a,b) - REAL(DP), DIMENSION(:), INTENT(IN) :: a,b - REAL(DP), DIMENSION(size(a),size(b)) :: outerprod_d - outerprod_d = spread(a,dim=2,ncopies=size(b)) * & - spread(b,dim=1,ncopies=size(a)) - END FUNCTION outerprod_d -!BL - FUNCTION outerdiv(a,b) - REAL(SP), DIMENSION(:), INTENT(IN) :: a,b - REAL(SP), DIMENSION(size(a),size(b)) :: outerdiv - outerdiv = spread(a,dim=2,ncopies=size(b)) / & - spread(b,dim=1,ncopies=size(a)) - END FUNCTION outerdiv -!BL - FUNCTION outersum(a,b) - REAL(SP), DIMENSION(:), INTENT(IN) :: a,b - REAL(SP), DIMENSION(size(a),size(b)) :: outersum - outersum = spread(a,dim=2,ncopies=size(b)) + & - spread(b,dim=1,ncopies=size(a)) - END FUNCTION outersum -!BL - FUNCTION outerdiff_r(a,b) - REAL(SP), DIMENSION(:), INTENT(IN) :: a,b - REAL(SP), DIMENSION(size(a),size(b)) :: outerdiff_r - outerdiff_r = spread(a,dim=2,ncopies=size(b)) - & - spread(b,dim=1,ncopies=size(a)) - END FUNCTION outerdiff_r -!BL - FUNCTION outerdiff_d(a,b) - REAL(DP), DIMENSION(:), INTENT(IN) :: a,b - REAL(DP), DIMENSION(size(a),size(b)) :: outerdiff_d - outerdiff_d = spread(a,dim=2,ncopies=size(b)) - & - spread(b,dim=1,ncopies=size(a)) - END FUNCTION outerdiff_d -!BL - FUNCTION outerdiff_i(a,b) - INTEGER(I4B), DIMENSION(:), INTENT(IN) :: a,b - INTEGER(I4B), DIMENSION(size(a),size(b)) :: outerdiff_i - outerdiff_i = spread(a,dim=2,ncopies=size(b)) - & - spread(b,dim=1,ncopies=size(a)) - END FUNCTION outerdiff_i -!BL - FUNCTION outerand(a,b) - LOGICAL(LGT), DIMENSION(:), INTENT(IN) :: a,b - LOGICAL(LGT), DIMENSION(size(a),size(b)) :: outerand - outerand = spread(a,dim=2,ncopies=size(b)) .and. & - spread(b,dim=1,ncopies=size(a)) - END FUNCTION outerand -!BL - SUBROUTINE scatter_add_r(dest,source,dest_index) - REAL(SP), DIMENSION(:), INTENT(OUT) :: dest - REAL(SP), DIMENSION(:), INTENT(IN) :: source - INTEGER(I4B), DIMENSION(:), INTENT(IN) :: dest_index - INTEGER(I4B) :: m,n,j,i - n=assert_eq2(size(source),size(dest_index),'scatter_add_r') - m=size(dest) - do j=1,n - i=dest_index(j) - if (i > 0 .and. i <= m) dest(i)=dest(i)+source(j) - end do - END SUBROUTINE scatter_add_r - SUBROUTINE scatter_add_d(dest,source,dest_index) - REAL(DP), DIMENSION(:), INTENT(OUT) :: dest - REAL(DP), DIMENSION(:), INTENT(IN) :: source - INTEGER(I4B), DIMENSION(:), INTENT(IN) :: dest_index - INTEGER(I4B) :: m,n,j,i - n=assert_eq2(size(source),size(dest_index),'scatter_add_d') - m=size(dest) - do j=1,n - i=dest_index(j) - if (i > 0 .and. i <= m) dest(i)=dest(i)+source(j) - end do - END SUBROUTINE scatter_add_d - SUBROUTINE scatter_max_r(dest,source,dest_index) - REAL(SP), DIMENSION(:), INTENT(OUT) :: dest - REAL(SP), DIMENSION(:), INTENT(IN) :: source - INTEGER(I4B), DIMENSION(:), INTENT(IN) :: dest_index - INTEGER(I4B) :: m,n,j,i - n=assert_eq2(size(source),size(dest_index),'scatter_max_r') - m=size(dest) - do j=1,n - i=dest_index(j) - if (i > 0 .and. i <= m) dest(i)=max(dest(i),source(j)) - end do - END SUBROUTINE scatter_max_r - SUBROUTINE scatter_max_d(dest,source,dest_index) - REAL(DP), DIMENSION(:), INTENT(OUT) :: dest - REAL(DP), DIMENSION(:), INTENT(IN) :: source - INTEGER(I4B), DIMENSION(:), INTENT(IN) :: dest_index - INTEGER(I4B) :: m,n,j,i - n=assert_eq2(size(source),size(dest_index),'scatter_max_d') - m=size(dest) - do j=1,n - i=dest_index(j) - if (i > 0 .and. i <= m) dest(i)=max(dest(i),source(j)) - end do - END SUBROUTINE scatter_max_d -!BL - SUBROUTINE diagadd_rv(mat,diag) - REAL(SP), DIMENSION(:,:), INTENT(INOUT) :: mat - REAL(SP), DIMENSION(:), INTENT(IN) :: diag - INTEGER(I4B) :: j,n - n = assert_eq2(size(diag),min(size(mat,1),size(mat,2)),'diagadd_rv') - do j=1,n - mat(j,j)=mat(j,j)+diag(j) - end do - END SUBROUTINE diagadd_rv -!BL - SUBROUTINE diagadd_r(mat,diag) - REAL(SP), DIMENSION(:,:), INTENT(INOUT) :: mat - REAL(SP), INTENT(IN) :: diag - INTEGER(I4B) :: j,n - n = min(size(mat,1),size(mat,2)) - do j=1,n - mat(j,j)=mat(j,j)+diag - end do - END SUBROUTINE diagadd_r -!BL - SUBROUTINE diagmult_rv(mat,diag) - REAL(SP), DIMENSION(:,:), INTENT(INOUT) :: mat - REAL(SP), DIMENSION(:), INTENT(IN) :: diag - INTEGER(I4B) :: j,n - n = assert_eq2(size(diag),min(size(mat,1),size(mat,2)),'diagmult_rv') - do j=1,n - mat(j,j)=mat(j,j)*diag(j) - end do - END SUBROUTINE diagmult_rv -!BL - SUBROUTINE diagmult_r(mat,diag) - REAL(SP), DIMENSION(:,:), INTENT(INOUT) :: mat - REAL(SP), INTENT(IN) :: diag - INTEGER(I4B) :: j,n - n = min(size(mat,1),size(mat,2)) - do j=1,n - mat(j,j)=mat(j,j)*diag - end do - END SUBROUTINE diagmult_r -!BL - FUNCTION get_diag_rv(mat) - REAL(SP), DIMENSION(:,:), INTENT(IN) :: mat - REAL(SP), DIMENSION(size(mat,1)) :: get_diag_rv - INTEGER(I4B) :: j - j=assert_eq2(size(mat,1),size(mat,2),'get_diag_rv') - do j=1,size(mat,1) - get_diag_rv(j)=mat(j,j) - end do - END FUNCTION get_diag_rv -!BL - FUNCTION get_diag_dv(mat) - REAL(DP), DIMENSION(:,:), INTENT(IN) :: mat - REAL(DP), DIMENSION(size(mat,1)) :: get_diag_dv - INTEGER(I4B) :: j - j=assert_eq2(size(mat,1),size(mat,2),'get_diag_dv') - do j=1,size(mat,1) - get_diag_dv(j)=mat(j,j) - end do - END FUNCTION get_diag_dv -!BL - SUBROUTINE put_diag_rv(diagv,mat) - REAL(SP), DIMENSION(:), INTENT(IN) :: diagv - REAL(SP), DIMENSION(:,:), INTENT(INOUT) :: mat - INTEGER(I4B) :: j,n - n=assert_eq2(size(diagv),min(size(mat,1),size(mat,2)),'put_diag_rv') - do j=1,n - mat(j,j)=diagv(j) - end do - END SUBROUTINE put_diag_rv -!BL - SUBROUTINE put_diag_r(scal,mat) - REAL(SP), INTENT(IN) :: scal - REAL(SP), DIMENSION(:,:), INTENT(INOUT) :: mat - INTEGER(I4B) :: j,n - n = min(size(mat,1),size(mat,2)) - do j=1,n - mat(j,j)=scal - end do - END SUBROUTINE put_diag_r -!BL - SUBROUTINE unit_matrix(mat) - REAL(SP), DIMENSION(:,:), INTENT(OUT) :: mat - INTEGER(I4B) :: i,n - n=min(size(mat,1),size(mat,2)) - mat(:,:)=0.0_sp - do i=1,n - mat(i,i)=1.0_sp - end do - END SUBROUTINE unit_matrix -!BL - FUNCTION upper_triangle(j,k,extra) - INTEGER(I4B), INTENT(IN) :: j,k - INTEGER(I4B), OPTIONAL, INTENT(IN) :: extra - LOGICAL(LGT), DIMENSION(j,k) :: upper_triangle - INTEGER(I4B) :: n - n=0 - if (present(extra)) n=extra - upper_triangle=(outerdiff(arth_i(1,1,j),arth_i(1,1,k)) < n) - END FUNCTION upper_triangle -!BL - FUNCTION lower_triangle(j,k,extra) - INTEGER(I4B), INTENT(IN) :: j,k - INTEGER(I4B), OPTIONAL, INTENT(IN) :: extra - LOGICAL(LGT), DIMENSION(j,k) :: lower_triangle - INTEGER(I4B) :: n - n=0 - if (present(extra)) n=extra - lower_triangle=(outerdiff(arth_i(1,1,j),arth_i(1,1,k)) > -n) - END FUNCTION lower_triangle -!BL - FUNCTION vabs(v) - REAL(SP), DIMENSION(:), INTENT(IN) :: v - REAL(SP) :: vabs - vabs=sqrt(dot_product(v,v)) - END FUNCTION vabs -!BL -END MODULE nrutil diff --git a/odeint.f90 b/odeint.f90 deleted file mode 100755 index e07d854..0000000 --- a/odeint.f90 +++ /dev/null @@ -1,95 +0,0 @@ -MODULE ode_path - USE nrtype - INTEGER(I4B) :: nok,nbad,kount - LOGICAL(LGT), SAVE :: save_steps=.false. - REAL(SP) :: dxsav - REAL(SP), DIMENSION(:), POINTER :: xp - REAL(SP), DIMENSION(:,:), POINTER :: yp -END MODULE ode_path - - SUBROUTINE odeint(ystart,x1,x2,eps,h1,hmin,derivs,rkqs) - USE nrtype; USE nrutil, ONLY : nrerror,reallocate - USE ode_path - IMPLICIT NONE - REAL(SP), DIMENSION(:), INTENT(INOUT) :: ystart - REAL(SP), INTENT(IN) :: x1,x2,eps,h1,hmin - INTERFACE - SUBROUTINE derivs(x,y,dydx) - USE nrtype - IMPLICIT NONE - REAL(SP), INTENT(IN) :: x - REAL(SP), DIMENSION(:), INTENT(IN) :: y - REAL(SP), DIMENSION(:), INTENT(OUT) :: dydx - END SUBROUTINE derivs -!BL - SUBROUTINE rkqs(y,dydx,x,htry,eps,yscal,hdid,hnext,derivs) - USE nrtype - IMPLICIT NONE - REAL(SP), DIMENSION(:), INTENT(INOUT) :: y - REAL(SP), DIMENSION(:), INTENT(IN) :: dydx,yscal - REAL(SP), INTENT(INOUT) :: x - REAL(SP), INTENT(IN) :: htry,eps - REAL(SP), INTENT(OUT) :: hdid,hnext - INTERFACE - SUBROUTINE derivs(x,y,dydx) - USE nrtype - IMPLICIT NONE - REAL(SP), INTENT(IN) :: x - REAL(SP), DIMENSION(:), INTENT(IN) :: y - REAL(SP), DIMENSION(:), INTENT(OUT) :: dydx - END SUBROUTINE derivs - END INTERFACE - END SUBROUTINE rkqs - END INTERFACE - REAL(SP), PARAMETER :: TINY=1.0e-30_sp - INTEGER(I4B), PARAMETER :: MAXSTP=10000 - INTEGER(I4B) :: nstp - REAL(SP) :: h,hdid,hnext,x,xsav - REAL(SP), DIMENSION(size(ystart)) :: dydx,y,yscal - x=x1 - h=sign(h1,x2-x1) - nok=0 - nbad=0 - kount=0 - y(:)=ystart(:) - nullify(xp,yp) - if (save_steps) then - xsav=x-2.0_sp*dxsav - allocate(xp(256)) - allocate(yp(size(ystart),size(xp))) - end if - do nstp=1,MAXSTP - call derivs(x,y,dydx) - yscal(:)=abs(y(:))+abs(h*dydx(:))+TINY - if (save_steps .and. (abs(x-xsav) > abs(dxsav))) & - call save_a_step - if ((x+h-x2)*(x+h-x1) > 0.0) h=x2-x - call rkqs(y,dydx,x,h,eps,yscal,hdid,hnext,derivs) - if (hdid == h) then - nok=nok+1 - else - nbad=nbad+1 - end if - if ((x-x2)*(x2-x1) >= 0.0) then - ystart(:)=y(:) - if (save_steps) call save_a_step - RETURN - end if - if (abs(hnext) < hmin)& - call nrerror('stepsize smaller than minimum in odeint') - h=hnext - end do - call nrerror('too many steps in odeint') - CONTAINS -!BL - SUBROUTINE save_a_step - kount=kount+1 - if (kount > size(xp)) then - xp=>reallocate(xp,2*size(xp)) - yp=>reallocate(yp,size(yp,1),size(xp)) - end if - xp(kount)=x - yp(:,kount)=y(:) - xsav=x - END SUBROUTINE save_a_step - END SUBROUTINE odeint diff --git a/polint.f90 b/polint.f90 deleted file mode 100755 index cbb17fb..0000000 --- a/polint.f90 +++ /dev/null @@ -1,31 +0,0 @@ - SUBROUTINE polint(xa,ya,x,y,dy) - USE nrtype; USE nrutil, ONLY : assert_eq,iminloc,nrerror - IMPLICIT NONE - REAL(SP), DIMENSION(:), INTENT(IN) :: xa,ya - REAL(SP), INTENT(IN) :: x - REAL(SP), INTENT(OUT) :: y,dy - INTEGER(I4B) :: m,n,ns - REAL(SP), DIMENSION(size(xa)) :: c,d,den,ho - n=assert_eq(size(xa),size(ya),'polint') - c=ya - d=ya - ho=xa-x - ns=iminloc(abs(x-xa)) - y=ya(ns) - ns=ns-1 - do m=1,n-1 - den(1:n-m)=ho(1:n-m)-ho(1+m:n) - if (any(den(1:n-m) == 0.0)) & - call nrerror('polint: calculation failure') - den(1:n-m)=(c(2:n-m+1)-d(1:n-m))/den(1:n-m) - d(1:n-m)=ho(1+m:n)*den(1:n-m) - c(1:n-m)=ho(1:n-m)*den(1:n-m) - if (2*ns < n-m) then - dy=c(ns+1) - else - dy=d(ns) - ns=ns-1 - end if - y=y+dy - end do - END SUBROUTINE polint diff --git a/random.f90 b/random.f90 deleted file mode 100644 index 0b5f278..0000000 --- a/random.f90 +++ /dev/null @@ -1,1595 +0,0 @@ -MODULE random -! A module for random number generation from the following distributions: -! -! Distribution Function/subroutine name -! -! Normal (Gaussian) random_normal -! Gamma random_gamma -! Chi-squared random_chisq -! Exponential random_exponential -! Weibull random_Weibull -! Beta random_beta -! t random_t -! Multivariate normal random_mvnorm -! Generalized inverse Gaussian random_inv_gauss -! Poisson random_Poisson -! Binomial random_binomial1 * -! random_binomial2 * -! Negative binomial random_neg_binomial -! von Mises random_von_Mises -! Cauchy random_Cauchy -! -! Generate a random ordering of the integers 1 .. N -! random_order -! Initialize (seed) the uniform random number generator for ANY compiler -! seed_random_number - -! Lognormal - see note below. - -! ** Two functions are provided for the binomial distribution. -! If the parameter values remain constant, it is recommended that the -! first function is used (random_binomial1). If one or both of the -! parameters change, use the second function (random_binomial2). - -! The compilers own random number generator, SUBROUTINE RANDOM_NUMBER(r), -! is used to provide a source of uniformly distributed random numbers. - -! N.B. At this stage, only one random number is generated at each call to -! one of the functions above. - -! The module uses the following functions which are included here: -! bin_prob to calculate a single binomial probability -! lngamma to calculate the logarithm to base e of the gamma function - -! Some of the code is adapted from Dagpunar's book: -! Dagpunar, J. 'Principles of random variate generation' -! Clarendon Press, Oxford, 1988. ISBN 0-19-852202-9 -! -! In most of Dagpunar's routines, there is a test to see whether the value -! of one or two floating-point parameters has changed since the last call. -! These tests have been replaced by using a logical variable FIRST. -! This should be set to .TRUE. on the first call using new values of the -! parameters, and .FALSE. if the parameter values are the same as for the -! previous call. - -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -! Lognormal distribution -! If X has a lognormal distribution, then log(X) is normally distributed. -! Here the logarithm is the natural logarithm, that is to base e, sometimes -! denoted as ln. To generate random variates from this distribution, generate -! a random deviate from the normal distribution with mean and variance equal -! to the mean and variance of the logarithms of X, then take its exponential. - -! Relationship between the mean & variance of log(X) and the mean & variance -! of X, when X has a lognormal distribution. -! Let m = mean of log(X), and s^2 = variance of log(X) -! Then -! mean of X = exp(m + 0.5s^2) -! variance of X = (mean(X))^2.[exp(s^2) - 1] - -! In the reverse direction (rarely used) -! variance of log(X) = log[1 + var(X)/(mean(X))^2] -! mean of log(X) = log(mean(X) - 0.5var(log(X)) - -! N.B. The above formulae relate to population parameters; they will only be -! approximate if applied to sample values. -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - -! Version 1.13, 2 October 2000 -! Changes from version 1.01 -! 1. The random_order, random_Poisson & random_binomial routines have been -! replaced with more efficient routines. -! 2. A routine, seed_random_number, has been added to seed the uniform random -! number generator. This requires input of the required number of seeds -! for the particular compiler from a specified I/O unit such as a keyboard. -! 3. Made compatible with Lahey's ELF90. -! 4. Marsaglia & Tsang algorithm used for random_gamma when shape parameter > 1. -! 5. INTENT for array f corrected in random_mvnorm. - -! Author: Alan Miller -! CSIRO Division of Mathematical & Information Sciences -! Private Bag 10, Clayton South MDC -! Clayton 3169, Victoria, Australia -! Phone: (+61) 3 9545-8016 Fax: (+61) 3 9545-8080 -! e-mail: amiller @ bigpond.net.au - -IMPLICIT NONE -REAL, PRIVATE :: zero = 0.0, half = 0.5, one = 1.0, two = 2.0, & - vsmall = TINY(1.0), vlarge = HUGE(1.0) -PRIVATE :: integral -INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(12, 60) - - -CONTAINS - - -FUNCTION random_normal() RESULT(fn_val) - -! Adapted from the following Fortran 77 code -! ALGORITHM 712, COLLECTED ALGORITHMS FROM ACM. -! THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE, -! VOL. 18, NO. 4, DECEMBER, 1992, PP. 434-435. - -! The function random_normal() returns a normally distributed pseudo-random -! number with zero mean and unit variance. - -! The algorithm uses the ratio of uniforms method of A.J. Kinderman -! and J.F. Monahan augmented with quadratic bounding curves. - -REAL :: fn_val - -! Local variables -REAL :: s = 0.449871, t = -0.386595, a = 0.19600, b = 0.25472, & - r1 = 0.27597, r2 = 0.27846, u, v, x, y, q - -! Generate P = (u,v) uniform in rectangle enclosing acceptance region - -DO - CALL RANDOM_NUMBER(u) - CALL RANDOM_NUMBER(v) - v = 1.7156 * (v - half) - -! Evaluate the quadratic form - x = u - s - y = ABS(v) - t - q = x**2 + y*(a*y - b*x) - -! Accept P if inside inner ellipse - IF (q < r1) EXIT -! Reject P if outside outer ellipse - IF (q > r2) CYCLE -! Reject P if outside acceptance region - IF (v**2 < -4.0*LOG(u)*u**2) EXIT -END DO - -! Return ratio of P's coordinates as the normal deviate -fn_val = v/u -RETURN - -END FUNCTION random_normal - - - -FUNCTION random_gamma(s, first) RESULT(fn_val) - -! Adapted from Fortran 77 code from the book: -! Dagpunar, J. 'Principles of random variate generation' -! Clarendon Press, Oxford, 1988. ISBN 0-19-852202-9 - -! FUNCTION GENERATES A RANDOM GAMMA VARIATE. -! CALLS EITHER random_gamma1 (S > 1.0) -! OR random_exponential (S = 1.0) -! OR random_gamma2 (S < 1.0). - -! S = SHAPE PARAMETER OF DISTRIBUTION (0 < REAL). - -REAL, INTENT(IN) :: s -LOGICAL, INTENT(IN) :: first -REAL :: fn_val - -IF (s <= zero) THEN - WRITE(*, *) 'SHAPE PARAMETER VALUE MUST BE POSITIVE' - STOP -END IF - -IF (s > one) THEN - fn_val = random_gamma1(s, first) -ELSE IF (s < one) THEN - fn_val = random_gamma2(s, first) -ELSE - fn_val = random_exponential() -END IF - -RETURN -END FUNCTION random_gamma - - - -FUNCTION random_gamma1(s, first) RESULT(fn_val) - -! Uses the algorithm in -! Marsaglia, G. and Tsang, W.W. (2000) `A simple method for generating -! gamma variables', Trans. om Math. Software (TOMS), vol.26(3), pp.363-372. - -! Generates a random gamma deviate for shape parameter s >= 1. - -REAL, INTENT(IN) :: s -LOGICAL, INTENT(IN) :: first -REAL :: fn_val - -! Local variables -REAL, SAVE :: c, d -REAL :: u, v, x - -IF (first) THEN - d = s - one/3. - c = one/SQRT(9.0*d) -END IF - -! Start of main loop -DO - -! Generate v = (1+cx)^3 where x is random normal; repeat if v <= 0. - - DO - x = random_normal() - v = (one + c*x)**3 - IF (v > zero) EXIT - END DO - -! Generate uniform variable U - - CALL RANDOM_NUMBER(u) - IF (u < one - 0.0331*x**4) THEN - fn_val = d*v - EXIT - ELSE IF (LOG(u) < half*x**2 + d*(one - v + LOG(v))) THEN - fn_val = d*v - EXIT - END IF -END DO - -RETURN -END FUNCTION random_gamma1 - - - -FUNCTION random_gamma2(s, first) RESULT(fn_val) - -! Adapted from Fortran 77 code from the book: -! Dagpunar, J. 'Principles of random variate generation' -! Clarendon Press, Oxford, 1988. ISBN 0-19-852202-9 - -! FUNCTION GENERATES A RANDOM VARIATE IN [0,INFINITY) FROM -! A GAMMA DISTRIBUTION WITH DENSITY PROPORTIONAL TO -! GAMMA2**(S-1) * EXP(-GAMMA2), -! USING A SWITCHING METHOD. - -! S = SHAPE PARAMETER OF DISTRIBUTION -! (REAL < 1.0) - -REAL, INTENT(IN) :: s -LOGICAL, INTENT(IN) :: first -REAL :: fn_val - -! Local variables -REAL :: r, x, w -REAL, SAVE :: a, p, c, uf, vr, d - -IF (s <= zero .OR. s >= one) THEN - WRITE(*, *) 'SHAPE PARAMETER VALUE OUTSIDE PERMITTED RANGE' - STOP -END IF - -IF (first) THEN ! Initialization, if necessary - a = one - s - p = a/(a + s*EXP(-a)) - IF (s < vsmall) THEN - WRITE(*, *) 'SHAPE PARAMETER VALUE TOO SMALL' - STOP - END IF - c = one/s - uf = p*(vsmall/a)**s - vr = one - vsmall - d = a*LOG(a) -END IF - -DO - CALL RANDOM_NUMBER(r) - IF (r >= vr) THEN - CYCLE - ELSE IF (r > p) THEN - x = a - LOG((one - r)/(one - p)) - w = a*LOG(x)-d - ELSE IF (r > uf) THEN - x = a*(r/p)**c - w = x - ELSE - fn_val = zero - RETURN - END IF - - CALL RANDOM_NUMBER(r) - IF (one-r <= w .AND. r > zero) THEN - IF (r*(w + one) >= one) CYCLE - IF (-LOG(r) <= w) CYCLE - END IF - EXIT -END DO - -fn_val = x -RETURN - -END FUNCTION random_gamma2 - - - -FUNCTION random_chisq(ndf, first) RESULT(fn_val) - -! Generates a random variate from the chi-squared distribution with -! ndf degrees of freedom - -INTEGER, INTENT(IN) :: ndf -LOGICAL, INTENT(IN) :: first -REAL :: fn_val - -fn_val = two * random_gamma(half*ndf, first) -RETURN - -END FUNCTION random_chisq - - - -FUNCTION random_exponential() RESULT(fn_val) - -! Adapted from Fortran 77 code from the book: -! Dagpunar, J. 'Principles of random variate generation' -! Clarendon Press, Oxford, 1988. ISBN 0-19-852202-9 - -! FUNCTION GENERATES A RANDOM VARIATE IN [0,INFINITY) FROM -! A NEGATIVE EXPONENTIAL DlSTRIBUTION WlTH DENSITY PROPORTIONAL -! TO EXP(-random_exponential), USING INVERSION. - -REAL :: fn_val - -! Local variable -REAL :: r - -DO - CALL RANDOM_NUMBER(r) - IF (r > zero) EXIT -END DO - -fn_val = -LOG(r) -RETURN - -END FUNCTION random_exponential - - - -FUNCTION random_Weibull(a) RESULT(fn_val) - -! Generates a random variate from the Weibull distribution with -! probability density: -! a -! a-1 -x -! f(x) = a.x e - -REAL, INTENT(IN) :: a -REAL :: fn_val - -! For speed, there is no checking that a is not zero or very small. - -fn_val = random_exponential() ** (one/a) -RETURN - -END FUNCTION random_Weibull - - - -FUNCTION random_beta(aa, bb, first) RESULT(fn_val) - -! Adapted from Fortran 77 code from the book: -! Dagpunar, J. 'Principles of random variate generation' -! Clarendon Press, Oxford, 1988. ISBN 0-19-852202-9 - -! FUNCTION GENERATES A RANDOM VARIATE IN [0,1] -! FROM A BETA DISTRIBUTION WITH DENSITY -! PROPORTIONAL TO BETA**(AA-1) * (1-BETA)**(BB-1). -! USING CHENG'S LOG LOGISTIC METHOD. - -! AA = SHAPE PARAMETER FROM DISTRIBUTION (0 < REAL) -! BB = SHAPE PARAMETER FROM DISTRIBUTION (0 < REAL) - -REAL, INTENT(IN) :: aa, bb -LOGICAL, INTENT(IN) :: first -REAL :: fn_val - -! Local variables -REAL, PARAMETER :: aln4 = 1.3862944 -REAL :: a, b, g, r, s, x, y, z -REAL, SAVE :: d, f, h, t, c -LOGICAL, SAVE :: swap - -IF (aa <= zero .OR. bb <= zero) THEN - WRITE(*, *) 'IMPERMISSIBLE SHAPE PARAMETER VALUE(S)' - STOP -END IF - -IF (first) THEN ! Initialization, if necessary - a = aa - b = bb - swap = b > a - IF (swap) THEN - g = b - b = a - a = g - END IF - d = a/b - f = a+b - IF (b > one) THEN - h = SQRT((two*a*b - f)/(f - two)) - t = one - ELSE - h = b - t = one/(one + (a/(vlarge*b))**b) - END IF - c = a+h -END IF - -DO - CALL RANDOM_NUMBER(r) - CALL RANDOM_NUMBER(x) - s = r*r*x - IF (r < vsmall .OR. s <= zero) CYCLE - IF (r < t) THEN - x = LOG(r/(one - r))/h - y = d*EXP(x) - z = c*x + f*LOG((one + d)/(one + y)) - aln4 - IF (s - one > z) THEN - IF (s - s*z > one) CYCLE - IF (LOG(s) > z) CYCLE - END IF - fn_val = y/(one + y) - ELSE - IF (4.0*s > (one + one/d)**f) CYCLE - fn_val = one - END IF - EXIT -END DO - -IF (swap) fn_val = one - fn_val -RETURN -END FUNCTION random_beta - - - -FUNCTION random_t(m) RESULT(fn_val) - -! Adapted from Fortran 77 code from the book: -! Dagpunar, J. 'Principles of random variate generation' -! Clarendon Press, Oxford, 1988. ISBN 0-19-852202-9 - -! FUNCTION GENERATES A RANDOM VARIATE FROM A -! T DISTRIBUTION USING KINDERMAN AND MONAHAN'S RATIO METHOD. - -! M = DEGREES OF FREEDOM OF DISTRIBUTION -! (1 <= 1NTEGER) - -INTEGER, INTENT(IN) :: m -REAL :: fn_val - -! Local variables -REAL, SAVE :: s, c, a, f, g -REAL :: r, x, v - -REAL, PARAMETER :: three = 3.0, four = 4.0, quart = 0.25, & - five = 5.0, sixteen = 16.0 -INTEGER :: mm = 0 - -IF (m < 1) THEN - WRITE(*, *) 'IMPERMISSIBLE DEGREES OF FREEDOM' - STOP -END IF - -IF (m /= mm) THEN ! Initialization, if necessary - s = m - c = -quart*(s + one) - a = four/(one + one/s)**c - f = sixteen/a - IF (m > 1) THEN - g = s - one - g = ((s + one)/g)**c*SQRT((s+s)/g) - ELSE - g = one - END IF - mm = m -END IF - -DO - CALL RANDOM_NUMBER(r) - IF (r <= zero) CYCLE - CALL RANDOM_NUMBER(v) - x = (two*v - one)*g/r - v = x*x - IF (v > five - a*r) THEN - IF (m >= 1 .AND. r*(v + three) > f) CYCLE - IF (r > (one + v/s)**c) CYCLE - END IF - EXIT -END DO - -fn_val = x -RETURN -END FUNCTION random_t - - - -SUBROUTINE random_mvnorm(n, h, d, f, first, x, ier) - -! Adapted from Fortran 77 code from the book: -! Dagpunar, J. 'Principles of random variate generation' -! Clarendon Press, Oxford, 1988. ISBN 0-19-852202-9 - -! N.B. An extra argument, ier, has been added to Dagpunar's routine - -! SUBROUTINE GENERATES AN N VARIATE RANDOM NORMAL -! VECTOR USING A CHOLESKY DECOMPOSITION. - -! ARGUMENTS: -! N = NUMBER OF VARIATES IN VECTOR -! (INPUT,INTEGER >= 1) -! H(J) = J'TH ELEMENT OF VECTOR OF MEANS -! (INPUT,REAL) -! X(J) = J'TH ELEMENT OF DELIVERED VECTOR -! (OUTPUT,REAL) -! -! D(J*(J-1)/2+I) = (I,J)'TH ELEMENT OF VARIANCE MATRIX (J> = I) -! (INPUT,REAL) -! F((J-1)*(2*N-J)/2+I) = (I,J)'TH ELEMENT OF LOWER TRIANGULAR -! DECOMPOSITION OF VARIANCE MATRIX (J <= I) -! (OUTPUT,REAL) - -! FIRST = .TRUE. IF THIS IS THE FIRST CALL OF THE ROUTINE -! OR IF THE DISTRIBUTION HAS CHANGED SINCE THE LAST CALL OF THE ROUTINE. -! OTHERWISE SET TO .FALSE. -! (INPUT,LOGICAL) - -! ier = 1 if the input covariance matrix is not +ve definite -! = 0 otherwise - -INTEGER, INTENT(IN) :: n -REAL, INTENT(IN) :: h(:), d(:) ! d(n*(n+1)/2) -REAL, INTENT(IN OUT) :: f(:) ! f(n*(n+1)/2) -REAL, INTENT(OUT) :: x(:) -LOGICAL, INTENT(IN) :: first -INTEGER, INTENT(OUT) :: ier - -! Local variables -INTEGER :: j, i, m -REAL :: y, v -INTEGER, SAVE :: n2 - -IF (n < 1) THEN - WRITE(*, *) 'SIZE OF VECTOR IS NON POSITIVE' - STOP -END IF - -ier = 0 -IF (first) THEN ! Initialization, if necessary - n2 = 2*n - IF (d(1) < zero) THEN - ier = 1 - RETURN - END IF - - f(1) = SQRT(d(1)) - y = one/f(1) - DO j = 2,n - f(j) = d(1+j*(j-1)/2) * y - END DO - - DO i = 2,n - v = d(i*(i-1)/2+i) - DO m = 1,i-1 - v = v - f((m-1)*(n2-m)/2+i)**2 - END DO - - IF (v < zero) THEN - ier = 1 - RETURN - END IF - - v = SQRT(v) - y = one/v - f((i-1)*(n2-i)/2+i) = v - DO j = i+1,n - v = d(j*(j-1)/2+i) - DO m = 1,i-1 - v = v - f((m-1)*(n2-m)/2+i)*f((m-1)*(n2-m)/2 + j) - END DO ! m = 1,i-1 - f((i-1)*(n2-i)/2 + j) = v*y - END DO ! j = i+1,n - END DO ! i = 2,n -END IF - -x(1:n) = h(1:n) -DO j = 1,n - y = random_normal() - DO i = j,n - x(i) = x(i) + f((j-1)*(n2-j)/2 + i) * y - END DO ! i = j,n -END DO ! j = 1,n - -RETURN -END SUBROUTINE random_mvnorm - - - -FUNCTION random_inv_gauss(h, b, first) RESULT(fn_val) - -! Adapted from Fortran 77 code from the book: -! Dagpunar, J. 'Principles of random variate generation' -! Clarendon Press, Oxford, 1988. ISBN 0-19-852202-9 - -! FUNCTION GENERATES A RANDOM VARIATE IN [0,INFINITY] FROM -! A REPARAMETERISED GENERALISED INVERSE GAUSSIAN (GIG) DISTRIBUTION -! WITH DENSITY PROPORTIONAL TO GIG**(H-1) * EXP(-0.5*B*(GIG+1/GIG)) -! USING A RATIO METHOD. - -! H = PARAMETER OF DISTRIBUTION (0 <= REAL) -! B = PARAMETER OF DISTRIBUTION (0 < REAL) - -REAL, INTENT(IN) :: h, b -LOGICAL, INTENT(IN) :: first -REAL :: fn_val - -! Local variables -REAL :: ym, xm, r, w, r1, r2, x -REAL, SAVE :: a, c, d, e -REAL, PARAMETER :: quart = 0.25 - -IF (h < zero .OR. b <= zero) THEN - WRITE(*, *) 'IMPERMISSIBLE DISTRIBUTION PARAMETER VALUES' - STOP -END IF - -IF (first) THEN ! Initialization, if necessary - IF (h > quart*b*SQRT(vlarge)) THEN - WRITE(*, *) 'THE RATIO H:B IS TOO SMALL' - STOP - END IF - e = b*b - d = h + one - ym = (-d + SQRT(d*d + e))/b - IF (ym < vsmall) THEN - WRITE(*, *) 'THE VALUE OF B IS TOO SMALL' - STOP - END IF - - d = h - one - xm = (d + SQRT(d*d + e))/b - d = half*d - e = -quart*b - r = xm + one/xm - w = xm*ym - a = w**(-half*h) * SQRT(xm/ym) * EXP(-e*(r - ym - one/ym)) - IF (a < vsmall) THEN - WRITE(*, *) 'THE VALUE OF H IS TOO LARGE' - STOP - END IF - c = -d*LOG(xm) - e*r -END IF - -DO - CALL RANDOM_NUMBER(r1) - IF (r1 <= zero) CYCLE - CALL RANDOM_NUMBER(r2) - x = a*r2/r1 - IF (x <= zero) CYCLE - IF (LOG(r1) < d*LOG(x) + e*(x + one/x) + c) EXIT -END DO - -fn_val = x - -RETURN -END FUNCTION random_inv_gauss - - - -FUNCTION random_Poisson(mu, first) RESULT(ival) -!********************************************************************** -! Translated to Fortran 90 by Alan Miller from: -! RANLIB -! -! Library of Fortran Routines for Random Number Generation -! -! Compiled and Written by: -! -! Barry W. Brown -! James Lovato -! -! Department of Biomathematics, Box 237 -! The University of Texas, M.D. Anderson Cancer Center -! 1515 Holcombe Boulevard -! Houston, TX 77030 -! -! This work was supported by grant CA-16672 from the National Cancer Institute. - -! GENerate POIsson random deviate - -! Function - -! Generates a single random deviate from a Poisson distribution with mean mu. - -! Arguments - -! mu --> The mean of the Poisson distribution from which -! a random deviate is to be generated. -! REAL mu - -! Method - -! For details see: - -! Ahrens, J.H. and Dieter, U. -! Computer Generation of Poisson Deviates -! From Modified Normal Distributions. -! ACM Trans. Math. Software, 8, 2 -! (June 1982),163-179 - -! TABLES: COEFFICIENTS A0-A7 FOR STEP F. FACTORIALS FACT -! COEFFICIENTS A(K) - FOR PX = FK*V*V*SUM(A(K)*V**K)-DEL - -! SEPARATION OF CASES A AND B - -! .. Scalar Arguments .. -REAL, INTENT(IN) :: mu -LOGICAL, INTENT(IN) :: first -INTEGER :: ival -! .. -! .. Local Scalars .. -REAL :: b1, b2, c, c0, c1, c2, c3, del, difmuk, e, fk, fx, fy, g, & - omega, px, py, t, u, v, x, xx -REAL, SAVE :: s, d, p, q, p0 -INTEGER :: j, k, kflag -LOGICAL, SAVE :: full_init -INTEGER, SAVE :: l, m -! .. -! .. Local Arrays .. -REAL, SAVE :: pp(35) -! .. -! .. Data statements .. -REAL, PARAMETER :: a0 = -.5, a1 = .3333333, a2 = -.2500068, a3 = .2000118, & - a4 = -.1661269, a5 = .1421878, a6 = -.1384794, & - a7 = .1250060 - -REAL, PARAMETER :: fact(10) = (/ 1., 1., 2., 6., 24., 120., 720., 5040., & - 40320., 362880. /) - -! .. -! .. Executable Statements .. -IF (mu > 10.0) THEN -! C A S E A. (RECALCULATION OF S, D, L IF MU HAS CHANGED) - - IF (first) THEN - s = SQRT(mu) - d = 6.0*mu*mu - -! THE POISSON PROBABILITIES PK EXCEED THE DISCRETE NORMAL -! PROBABILITIES FK WHENEVER K >= M(MU). L=IFIX(MU-1.1484) -! IS AN UPPER BOUND TO M(MU) FOR ALL MU >= 10 . - - l = mu - 1.1484 - full_init = .false. - END IF - - -! STEP N. NORMAL SAMPLE - random_normal() FOR STANDARD NORMAL DEVIATE - - g = mu + s*random_normal() - IF (g > 0.0) THEN - ival = g - -! STEP I. IMMEDIATE ACCEPTANCE IF ival IS LARGE ENOUGH - - IF (ival>=l) RETURN - -! STEP S. SQUEEZE ACCEPTANCE - SAMPLE U - - fk = ival - difmuk = mu - fk - CALL RANDOM_NUMBER(u) - IF (d*u >= difmuk*difmuk*difmuk) RETURN - END IF - -! STEP P. PREPARATIONS FOR STEPS Q AND H. -! (RECALCULATIONS OF PARAMETERS IF NECESSARY) -! .3989423=(2*PI)**(-.5) .416667E-1=1./24. .1428571=1./7. -! THE QUANTITIES B1, B2, C3, C2, C1, C0 ARE FOR THE HERMITE -! APPROXIMATIONS TO THE DISCRETE NORMAL PROBABILITIES FK. -! C=.1069/MU GUARANTEES MAJORIZATION BY THE 'HAT'-FUNCTION. - - IF (.NOT. full_init) THEN - omega = .3989423/s - b1 = .4166667E-1/mu - b2 = .3*b1*b1 - c3 = .1428571*b1*b2 - c2 = b2 - 15.*c3 - c1 = b1 - 6.*b2 + 45.*c3 - c0 = 1. - b1 + 3.*b2 - 15.*c3 - c = .1069/mu - full_init = .true. - END IF - - IF (g < 0.0) GO TO 50 - -! 'SUBROUTINE' F IS CALLED (KFLAG=0 FOR CORRECT RETURN) - - kflag = 0 - GO TO 70 - -! STEP Q. QUOTIENT ACCEPTANCE (RARE CASE) - - 40 IF (fy-u*fy <= py*EXP(px-fx)) RETURN - -! STEP E. EXPONENTIAL SAMPLE - random_exponential() FOR STANDARD EXPONENTIAL -! DEVIATE E AND SAMPLE T FROM THE LAPLACE 'HAT' -! (IF T <= -.6744 THEN PK < FK FOR ALL MU >= 10.) - - 50 e = random_exponential() - CALL RANDOM_NUMBER(u) - u = u + u - one - t = 1.8 + SIGN(e, u) - IF (t <= (-.6744)) GO TO 50 - ival = mu + s*t - fk = ival - difmuk = mu - fk - -! 'SUBROUTINE' F IS CALLED (KFLAG=1 FOR CORRECT RETURN) - - kflag = 1 - GO TO 70 - -! STEP H. HAT ACCEPTANCE (E IS REPEATED ON REJECTION) - - 60 IF (c*ABS(u) > py*EXP(px+e) - fy*EXP(fx+e)) GO TO 50 - RETURN - -! STEP F. 'SUBROUTINE' F. CALCULATION OF PX, PY, FX, FY. -! CASE ival < 10 USES FACTORIALS FROM TABLE FACT - - 70 IF (ival>=10) GO TO 80 - px = -mu - py = mu**ival/fact(ival+1) - GO TO 110 - -! CASE ival >= 10 USES POLYNOMIAL APPROXIMATION -! A0-A7 FOR ACCURACY WHEN ADVISABLE -! .8333333E-1=1./12. .3989423=(2*PI)**(-.5) - - 80 del = .8333333E-1/fk - del = del - 4.8*del*del*del - v = difmuk/fk - IF (ABS(v)>0.25) THEN - px = fk*LOG(one + v) - difmuk - del - ELSE - px = fk*v*v* (((((((a7*v+a6)*v+a5)*v+a4)*v+a3)*v+a2)*v+a1)*v+a0) - del - END IF - py = .3989423/SQRT(fk) - 110 x = (half - difmuk)/s - xx = x*x - fx = -half*xx - fy = omega* (((c3*xx + c2)*xx + c1)*xx + c0) - IF (kflag <= 0) GO TO 40 - GO TO 60 - -!--------------------------------------------------------------------------- -! C A S E B. mu < 10 -! START NEW TABLE AND CALCULATE P0 IF NECESSARY - -ELSE - IF (first) THEN - m = MAX(1, INT(mu)) - l = 0 - p = EXP(-mu) - q = p - p0 = p - END IF - -! STEP U. UNIFORM SAMPLE FOR INVERSION METHOD - - DO - CALL RANDOM_NUMBER(u) - ival = 0 - IF (u <= p0) RETURN - -! STEP T. TABLE COMPARISON UNTIL THE END PP(L) OF THE -! PP-TABLE OF CUMULATIVE POISSON PROBABILITIES -! (0.458=PP(9) FOR MU=10) - - IF (l == 0) GO TO 150 - j = 1 - IF (u > 0.458) j = MIN(l, m) - DO k = j, l - IF (u <= pp(k)) GO TO 180 - END DO - IF (l == 35) CYCLE - -! STEP C. CREATION OF NEW POISSON PROBABILITIES P -! AND THEIR CUMULATIVES Q=PP(K) - - 150 l = l + 1 - DO k = l, 35 - p = p*mu / k - q = q + p - pp(k) = q - IF (u <= q) GO TO 170 - END DO - l = 35 - END DO - - 170 l = k - 180 ival = k - RETURN -END IF - -RETURN -END FUNCTION random_Poisson - - - -FUNCTION random_binomial1(n, p, first) RESULT(ival) - -! FUNCTION GENERATES A RANDOM BINOMIAL VARIATE USING C.D.Kemp's method. -! This algorithm is suitable when many random variates are required -! with the SAME parameter values for n & p. - -! P = BERNOULLI SUCCESS PROBABILITY -! (0 <= REAL <= 1) -! N = NUMBER OF BERNOULLI TRIALS -! (1 <= INTEGER) -! FIRST = .TRUE. for the first call using the current parameter values -! = .FALSE. if the values of (n,p) are unchanged from last call - -! Reference: Kemp, C.D. (1986). `A modal method for generating binomial -! variables', Commun. Statist. - Theor. Meth. 15(3), 805-813. - -INTEGER, INTENT(IN) :: n -REAL, INTENT(IN) :: p -LOGICAL, INTENT(IN) :: first -INTEGER :: ival - -! Local variables - -INTEGER :: ru, rd -INTEGER, SAVE :: r0 -REAL :: u, pd, pu -REAL, SAVE :: odds_ratio, p_r -REAL, PARAMETER :: zero = 0.0, one = 1.0 - -IF (first) THEN - r0 = (n+1)*p - p_r = bin_prob(n, p, r0) - odds_ratio = p / (one - p) -END IF - -CALL RANDOM_NUMBER(u) -u = u - p_r -IF (u < zero) THEN - ival = r0 - RETURN -END IF - -pu = p_r -ru = r0 -pd = p_r -rd = r0 -DO - rd = rd - 1 - IF (rd >= 0) THEN - pd = pd * (rd+1) / (odds_ratio * (n-rd)) - u = u - pd - IF (u < zero) THEN - ival = rd - RETURN - END IF - END IF - - ru = ru + 1 - IF (ru <= n) THEN - pu = pu * (n-ru+1) * odds_ratio / ru - u = u - pu - IF (u < zero) THEN - ival = ru - RETURN - END IF - END IF -END DO - -! This point should not be reached, but just in case: - -ival = r0 -RETURN - -END FUNCTION random_binomial1 - - - -FUNCTION bin_prob(n, p, r) RESULT(fn_val) -! Calculate a binomial probability - -INTEGER, INTENT(IN) :: n, r -REAL, INTENT(IN) :: p -REAL :: fn_val - -! Local variable -REAL :: one = 1.0 - -fn_val = EXP( lngamma(DBLE(n+1)) - lngamma(DBLE(r+1)) - lngamma(DBLE(n-r+1)) & - + r*LOG(p) + (n-r)*LOG(one - p) ) -RETURN - -END FUNCTION bin_prob - - - -FUNCTION lngamma(x) RESULT(fn_val) - -! Logarithm to base e of the gamma function. -! -! Accurate to about 1.e-14. -! Programmer: Alan Miller - -! Latest revision of Fortran 77 version - 28 February 1988 - -REAL (dp), INTENT(IN) :: x -REAL (dp) :: fn_val - -! Local variables - -REAL (dp) :: a1 = -4.166666666554424D-02, a2 = 2.430554511376954D-03, & - a3 = -7.685928044064347D-04, a4 = 5.660478426014386D-04, & - temp, arg, product, lnrt2pi = 9.189385332046727D-1, & - pi = 3.141592653589793D0 -LOGICAL :: reflect - -! lngamma is not defined if x = 0 or a negative integer. - -IF (x > 0.d0) GO TO 10 -IF (x /= INT(x)) GO TO 10 -fn_val = 0.d0 -RETURN - -! If x < 0, use the reflection formula: -! gamma(x) * gamma(1-x) = pi * cosec(pi.x) - -10 reflect = (x < 0.d0) -IF (reflect) THEN - arg = 1.d0 - x -ELSE - arg = x -END IF - -! Increase the argument, if necessary, to make it > 10. - -product = 1.d0 -20 IF (arg <= 10.d0) THEN - product = product * arg - arg = arg + 1.d0 - GO TO 20 -END IF - -! Use a polynomial approximation to Stirling's formula. -! N.B. The real Stirling's formula is used here, not the simpler, but less -! accurate formula given by De Moivre in a letter to Stirling, which -! is the one usually quoted. - -arg = arg - 0.5D0 -temp = 1.d0/arg**2 -fn_val = lnrt2pi + arg * (LOG(arg) - 1.d0 + & - (((a4*temp + a3)*temp + a2)*temp + a1)*temp) - LOG(product) -IF (reflect) THEN - temp = SIN(pi * x) - fn_val = LOG(pi/temp) - fn_val -END IF -RETURN -END FUNCTION lngamma - - - -FUNCTION random_binomial2(n, pp, first) RESULT(ival) -!********************************************************************** -! Translated to Fortran 90 by Alan Miller from: -! RANLIB -! -! Library of Fortran Routines for Random Number Generation -! -! Compiled and Written by: -! -! Barry W. Brown -! James Lovato -! -! Department of Biomathematics, Box 237 -! The University of Texas, M.D. Anderson Cancer Center -! 1515 Holcombe Boulevard -! Houston, TX 77030 -! -! This work was supported by grant CA-16672 from the National Cancer Institute. - -! GENerate BINomial random deviate - -! Function - -! Generates a single random deviate from a binomial -! distribution whose number of trials is N and whose -! probability of an event in each trial is P. - -! Arguments - -! N --> The number of trials in the binomial distribution -! from which a random deviate is to be generated. -! INTEGER N - -! P --> The probability of an event in each trial of the -! binomial distribution from which a random deviate -! is to be generated. -! REAL P - -! FIRST --> Set FIRST = .TRUE. for the first call to perform initialization -! the set FIRST = .FALSE. for further calls using the same pair -! of parameter values (N, P). -! LOGICAL FIRST - -! random_binomial2 <-- A random deviate yielding the number of events -! from N independent trials, each of which has -! a probability of event P. -! INTEGER random_binomial - -! Method - -! This is algorithm BTPE from: - -! Kachitvichyanukul, V. and Schmeiser, B. W. -! Binomial Random Variate Generation. -! Communications of the ACM, 31, 2 (February, 1988) 216. - -!********************************************************************** - -!*****DETERMINE APPROPRIATE ALGORITHM AND WHETHER SETUP IS NECESSARY - -! .. -! .. Scalar Arguments .. -REAL, INTENT(IN) :: pp -INTEGER, INTENT(IN) :: n -LOGICAL, INTENT(IN) :: first -INTEGER :: ival -! .. -! .. Local Scalars .. -REAL :: alv, amaxp, f, f1, f2, u, v, w, w2, x, x1, x2, ynorm, z, z2 -REAL, PARAMETER :: zero = 0.0, half = 0.5, one = 1.0 -INTEGER :: i, ix, ix1, k, mp -INTEGER, SAVE :: m -REAL, SAVE :: p, q, xnp, ffm, fm, xnpq, p1, xm, xl, xr, c, al, xll, & - xlr, p2, p3, p4, qn, r, g - -! .. -! .. Executable Statements .. - -!*****SETUP, PERFORM ONLY WHEN PARAMETERS CHANGE - -IF (first) THEN - p = MIN(pp, one-pp) - q = one - p - xnp = n * p -END IF - -IF (xnp > 30.) THEN - IF (first) THEN - ffm = xnp + p - m = ffm - fm = m - xnpq = xnp * q - p1 = INT(2.195*SQRT(xnpq) - 4.6*q) + half - xm = fm + half - xl = xm - p1 - xr = xm + p1 - c = 0.134 + 20.5 / (15.3 + fm) - al = (ffm-xl) / (ffm - xl*p) - xll = al * (one + half*al) - al = (xr - ffm) / (xr*q) - xlr = al * (one + half*al) - p2 = p1 * (one + c + c) - p3 = p2 + c / xll - p4 = p3 + c / xlr - END IF - -!*****GENERATE VARIATE, Binomial mean at least 30. - - 20 CALL RANDOM_NUMBER(u) - u = u * p4 - CALL RANDOM_NUMBER(v) - -! TRIANGULAR REGION - - IF (u <= p1) THEN - ix = xm - p1 * v + u - GO TO 110 - END IF - -! PARALLELOGRAM REGION - - IF (u <= p2) THEN - x = xl + (u-p1) / c - v = v * c + one - ABS(xm-x) / p1 - IF (v > one .OR. v <= zero) GO TO 20 - ix = x - ELSE - -! LEFT TAIL - - IF (u <= p3) THEN - ix = xl + LOG(v) / xll - IF (ix < 0) GO TO 20 - v = v * (u-p2) * xll - ELSE - -! RIGHT TAIL - - ix = xr - LOG(v) / xlr - IF (ix > n) GO TO 20 - v = v * (u-p3) * xlr - END IF - END IF - -!*****DETERMINE APPROPRIATE WAY TO PERFORM ACCEPT/REJECT TEST - - k = ABS(ix-m) - IF (k <= 20 .OR. k >= xnpq/2-1) THEN - -! EXPLICIT EVALUATION - - f = one - r = p / q - g = (n+1) * r - IF (m < ix) THEN - mp = m + 1 - DO i = mp, ix - f = f * (g/i-r) - END DO - - ELSE IF (m > ix) THEN - ix1 = ix + 1 - DO i = ix1, m - f = f / (g/i-r) - END DO - END IF - - IF (v > f) THEN - GO TO 20 - ELSE - GO TO 110 - END IF - END IF - -! SQUEEZING USING UPPER AND LOWER BOUNDS ON LOG(F(X)) - - amaxp = (k/xnpq) * ((k*(k/3. + .625) + .1666666666666)/xnpq + half) - ynorm = -k * k / (2.*xnpq) - alv = LOG(v) - IF (alvynorm + amaxp) GO TO 20 - -! STIRLING'S (actually de Moivre's) FORMULA TO MACHINE ACCURACY FOR -! THE FINAL ACCEPTANCE/REJECTION TEST - - x1 = ix + 1 - f1 = fm + one - z = n + 1 - fm - w = n - ix + one - z2 = z * z - x2 = x1 * x1 - f2 = f1 * f1 - w2 = w * w - IF (alv - (xm*LOG(f1/x1) + (n-m+half)*LOG(z/w) + (ix-m)*LOG(w*p/(x1*q)) + & - (13860.-(462.-(132.-(99.-140./f2)/f2)/f2)/f2)/f1/166320. + & - (13860.-(462.-(132.-(99.-140./z2)/z2)/z2)/z2)/z/166320. + & - (13860.-(462.-(132.-(99.-140./x2)/x2)/x2)/x2)/x1/166320. + & - (13860.-(462.-(132.-(99.-140./w2)/w2)/w2)/w2)/w/166320.) > zero) THEN - GO TO 20 - ELSE - GO TO 110 - END IF - -ELSE -! INVERSE CDF LOGIC FOR MEAN LESS THAN 30 - IF (first) THEN - qn = q ** n - r = p / q - g = r * (n+1) - END IF - - 90 ix = 0 - f = qn - CALL RANDOM_NUMBER(u) - 100 IF (u >= f) THEN - IF (ix > 110) GO TO 90 - u = u - f - ix = ix + 1 - f = f * (g/ix - r) - GO TO 100 - END IF -END IF - -110 IF (pp > half) ix = n - ix -ival = ix -RETURN - -END FUNCTION random_binomial2 - - - - -FUNCTION random_neg_binomial(sk, p) RESULT(ival) - -! Adapted from Fortran 77 code from the book: -! Dagpunar, J. 'Principles of random variate generation' -! Clarendon Press, Oxford, 1988. ISBN 0-19-852202-9 - -! FUNCTION GENERATES A RANDOM NEGATIVE BINOMIAL VARIATE USING UNSTORED -! INVERSION AND/OR THE REPRODUCTIVE PROPERTY. - -! SK = NUMBER OF FAILURES REQUIRED (Dagpunar's words!) -! = the `power' parameter of the negative binomial -! (0 < REAL) -! P = BERNOULLI SUCCESS PROBABILITY -! (0 < REAL < 1) - -! THE PARAMETER H IS SET SO THAT UNSTORED INVERSION ONLY IS USED WHEN P <= H, -! OTHERWISE A COMBINATION OF UNSTORED INVERSION AND -! THE REPRODUCTIVE PROPERTY IS USED. - -REAL, INTENT(IN) :: sk, p -INTEGER :: ival - -! Local variables -! THE PARAMETER ULN = -LOG(MACHINE'S SMALLEST REAL NUMBER). - -REAL, PARAMETER :: h = 0.7 -REAL :: q, x, st, uln, v, r, s, y, g -INTEGER :: k, i, n - -IF (sk <= zero .OR. p <= zero .OR. p >= one) THEN - WRITE(*, *) 'IMPERMISSIBLE DISTRIBUTION PARAMETER VALUES' - STOP -END IF - -q = one - p -x = zero -st = sk -IF (p > h) THEN - v = one/LOG(p) - k = st - DO i = 1,k - DO - CALL RANDOM_NUMBER(r) - IF (r > zero) EXIT - END DO - n = v*LOG(r) - x = x + n - END DO - st = st - k -END IF - -s = zero -uln = -LOG(vsmall) -IF (st > -uln/LOG(q)) THEN - WRITE(*, *) ' P IS TOO LARGE FOR THIS VALUE OF SK' - STOP -END IF - -y = q**st -g = st -CALL RANDOM_NUMBER(r) -DO - IF (y > r) EXIT - r = r - y - s = s + one - y = y*p*g/s - g = g + one -END DO - -ival = x + s + half -RETURN -END FUNCTION random_neg_binomial - - - -FUNCTION random_von_Mises(k, first) RESULT(fn_val) - -! Algorithm VMD from: -! Dagpunar, J.S. (1990) `Sampling from the von Mises distribution via a -! comparison of random numbers', J. of Appl. Statist., 17, 165-168. - -! Fortran 90 code by Alan Miller -! CSIRO Division of Mathematical & Information Sciences - -! Arguments: -! k (real) parameter of the von Mises distribution. -! first (logical) set to .TRUE. the first time that the function -! is called, or the first time with a new value -! for k. When first = .TRUE., the function sets -! up starting values and may be very much slower. - -REAL, INTENT(IN) :: k -LOGICAL, INTENT(IN) :: first -REAL :: fn_val - -! Local variables - -INTEGER :: j, n -INTEGER, SAVE :: nk -REAL, PARAMETER :: pi = 3.14159265 -REAL, SAVE :: p(20), theta(0:20) -REAL :: sump, r, th, lambda, rlast -REAL (dp) :: dk - -IF (first) THEN ! Initialization, if necessary - IF (k < zero) THEN - WRITE(*, *) '** Error: argument k for random_von_Mises = ', k - RETURN - END IF - - nk = k + k + one - IF (nk > 20) THEN - WRITE(*, *) '** Error: argument k for random_von_Mises = ', k - RETURN - END IF - - dk = k - theta(0) = zero - IF (k > half) THEN - -! Set up array p of probabilities. - - sump = zero - DO j = 1, nk - IF (j < nk) THEN - theta(j) = ACOS(one - j/k) - ELSE - theta(nk) = pi - END IF - -! Numerical integration of e^[k.cos(x)] from theta(j-1) to theta(j) - - CALL integral(theta(j-1), theta(j), p(j), dk) - sump = sump + p(j) - END DO - p(1:nk) = p(1:nk) / sump - ELSE - p(1) = one - theta(1) = pi - END IF ! if k > 0.5 -END IF ! if first - -CALL RANDOM_NUMBER(r) -DO j = 1, nk - r = r - p(j) - IF (r < zero) EXIT -END DO -r = -r/p(j) - -DO - th = theta(j-1) + r*(theta(j) - theta(j-1)) - lambda = k - j + one - k*COS(th) - n = 1 - rlast = lambda - - DO - CALL RANDOM_NUMBER(r) - IF (r > rlast) EXIT - n = n + 1 - rlast = r - END DO - - IF (n .NE. 2*(n/2)) EXIT ! is n even? - CALL RANDOM_NUMBER(r) -END DO - -fn_val = SIGN(th, (r - rlast)/(one - rlast) - half) -RETURN -END FUNCTION random_von_Mises - - - -SUBROUTINE integral(a, b, result, dk) - -! Gaussian integration of exp(k.cosx) from a to b. - -REAL (dp), INTENT(IN) :: dk -REAL, INTENT(IN) :: a, b -REAL, INTENT(OUT) :: result - -! Local variables - -REAL (dp) :: xmid, range, x1, x2, & - x(3) = (/0.238619186083197_dp, 0.661209386466265_dp, 0.932469514203152_dp/), & - w(3) = (/0.467913934572691_dp, 0.360761573048139_dp, 0.171324492379170_dp/) -INTEGER :: i - -xmid = (a + b)/2._dp -range = (b - a)/2._dp - -result = 0._dp -DO i = 1, 3 - x1 = xmid + x(i)*range - x2 = xmid - x(i)*range - result = result + w(i)*(EXP(dk*COS(x1)) + EXP(dk*COS(x2))) -END DO - -result = result * range -RETURN -END SUBROUTINE integral - - - -FUNCTION random_Cauchy() RESULT(fn_val) - -! Generate a random deviate from the standard Cauchy distribution - -REAL :: fn_val - -! Local variables -REAL :: v(2) - -DO - CALL RANDOM_NUMBER(v) - v = two*(v - half) - IF (ABS(v(2)) < vsmall) CYCLE ! Test for zero - IF (v(1)**2 + v(2)**2 < one) EXIT -END DO -fn_val = v(1) / v(2) - -RETURN -END FUNCTION random_Cauchy - - - -SUBROUTINE random_order(order, n) - -! Generate a random ordering of the integers 1 ... n. - -INTEGER, INTENT(IN) :: n -INTEGER, INTENT(OUT) :: order(n) - -! Local variables - -INTEGER :: i, j, k -REAL :: wk - -DO i = 1, n - order(i) = i -END DO - -! Starting at the end, swap the current last indicator with one -! randomly chosen from those preceeding it. - -DO i = n, 2, -1 - CALL RANDOM_NUMBER(wk) - j = 1 + i * wk - IF (j < i) THEN - k = order(i) - order(i) = order(j) - order(j) = k - END IF -END DO - -RETURN -END SUBROUTINE random_order - - - -SUBROUTINE seed_random_number(iounit) - -INTEGER, INTENT(IN) :: iounit - -! Local variables - -INTEGER :: k -INTEGER, ALLOCATABLE :: seed(:) - -CALL RANDOM_SEED(SIZE=k) -ALLOCATE( seed(k) ) - -WRITE(*, '(a, i2, a)')' Enter ', k, ' integers for random no. seeds: ' -READ(*, *) seed -WRITE(iounit, '(a, (7i10))') ' Random no. seeds: ', seed -CALL RANDOM_SEED(PUT=seed) - -DEALLOCATE( seed ) - -RETURN -END SUBROUTINE seed_random_number - - -END MODULE random diff --git a/rkck.f90 b/rkck.f90 deleted file mode 100755 index f83ff85..0000000 --- a/rkck.f90 +++ /dev/null @@ -1,42 +0,0 @@ - SUBROUTINE rkck(y,dydx,x,h,yout,yerr,derivs) - USE nrtype; USE nrutil, ONLY : assert_eq - IMPLICIT NONE - REAL(SP), DIMENSION(:), INTENT(IN) :: y,dydx - REAL(SP), INTENT(IN) :: x,h - REAL(SP), DIMENSION(:), INTENT(OUT) :: yout,yerr - INTERFACE - SUBROUTINE derivs(x,y,dydx) - USE nrtype - IMPLICIT NONE - REAL(SP), INTENT(IN) :: x - REAL(SP), DIMENSION(:), INTENT(IN) :: y - REAL(SP), DIMENSION(:), INTENT(OUT) :: dydx - END SUBROUTINE derivs - END INTERFACE - INTEGER(I4B) :: ndum - REAL(SP), DIMENSION(size(y)) :: ak2,ak3,ak4,ak5,ak6,ytemp - REAL(SP), PARAMETER :: A2=0.2_sp,A3=0.3_sp,A4=0.6_sp,A5=1.0_sp,& - A6=0.875_sp,B21=0.2_sp,B31=3.0_sp/40.0_sp,B32=9.0_sp/40.0_sp,& - B41=0.3_sp,B42=-0.9_sp,B43=1.2_sp,B51=-11.0_sp/54.0_sp,& - B52=2.5_sp,B53=-70.0_sp/27.0_sp,B54=35.0_sp/27.0_sp,& - B61=1631.0_sp/55296.0_sp,B62=175.0_sp/512.0_sp,& - B63=575.0_sp/13824.0_sp,B64=44275.0_sp/110592.0_sp,& - B65=253.0_sp/4096.0_sp,C1=37.0_sp/378.0_sp,& - C3=250.0_sp/621.0_sp,C4=125.0_sp/594.0_sp,& - C6=512.0_sp/1771.0_sp,DC1=C1-2825.0_sp/27648.0_sp,& - DC3=C3-18575.0_sp/48384.0_sp,DC4=C4-13525.0_sp/55296.0_sp,& - DC5=-277.0_sp/14336.0_sp,DC6=C6-0.25_sp - ndum=assert_eq(size(y),size(dydx),size(yout),size(yerr),'rkck') - ytemp=y+B21*h*dydx - call derivs(x+A2*h,ytemp,ak2) - ytemp=y+h*(B31*dydx+B32*ak2) - call derivs(x+A3*h,ytemp,ak3) - ytemp=y+h*(B41*dydx+B42*ak2+B43*ak3) - call derivs(x+A4*h,ytemp,ak4) - ytemp=y+h*(B51*dydx+B52*ak2+B53*ak3+B54*ak4) - call derivs(x+A5*h,ytemp,ak5) - ytemp=y+h*(B61*dydx+B62*ak2+B63*ak3+B64*ak4+B65*ak5) - call derivs(x+A6*h,ytemp,ak6) - yout=y+h*(C1*dydx+C3*ak3+C4*ak4+C6*ak6) - yerr=h*(DC1*dydx+DC3*ak3+DC4*ak4+DC5*ak5+DC6*ak6) - END SUBROUTINE rkck diff --git a/rkqs.f90 b/rkqs.f90 deleted file mode 100755 index 6202256..0000000 --- a/rkqs.f90 +++ /dev/null @@ -1,43 +0,0 @@ - SUBROUTINE rkqs(y,dydx,x,htry,eps,yscal,hdid,hnext,derivs) - USE nrtype; USE nrutil, ONLY : assert_eq,nrerror - USE nr, ONLY : rkck - IMPLICIT NONE - REAL(SP), DIMENSION(:), INTENT(INOUT) :: y - REAL(SP), DIMENSION(:), INTENT(IN) :: dydx,yscal - REAL(SP), INTENT(INOUT) :: x - REAL(SP), INTENT(IN) :: htry,eps - REAL(SP), INTENT(OUT) :: hdid,hnext - INTERFACE - SUBROUTINE derivs(x,y,dydx) - USE nrtype - IMPLICIT NONE - REAL(SP), INTENT(IN) :: x - REAL(SP), DIMENSION(:), INTENT(IN) :: y - REAL(SP), DIMENSION(:), INTENT(OUT) :: dydx - END SUBROUTINE derivs - END INTERFACE - INTEGER(I4B) :: ndum - REAL(SP) :: errmax,h,htemp,xnew - REAL(SP), DIMENSION(size(y)) :: yerr,ytemp - REAL(SP), PARAMETER :: SAFETY=0.9_sp,PGROW=-0.2_sp,PSHRNK=-0.25_sp,& - ERRCON=1.89e-4 - ndum=assert_eq(size(y),size(dydx),size(yscal),'rkqs') - h=htry - do - call rkck(y,dydx,x,h,ytemp,yerr,derivs) - errmax=maxval(abs(yerr(:)/yscal(:)))/eps - if (errmax <= 1.0) exit - htemp=SAFETY*h*(errmax**PSHRNK) - h=sign(max(abs(htemp),0.1_sp*abs(h)),h) - xnew=x+h - if (xnew == x) call nrerror('stepsize underflow in rkqs') - end do - if (errmax > ERRCON) then - hnext=SAFETY*h*(errmax**PGROW) - else - hnext=5.0_sp*h - end if - hdid=h - x=x+h - y(:)=ytemp(:) - END SUBROUTINE rkqs diff --git a/variables.f90 b/variables.f90 index 99b4615..5c833a4 100644 --- a/variables.f90 +++ b/variables.f90 @@ -3,7 +3,7 @@ !>@brief !>variables for the shallow water model on sphere module variables - use nrtype + use numerics_type !>@author !>Paul J. Connolly, The University of Manchester !>@brief @@ -17,13 +17,13 @@ module variables ! variables for grid integer(i4b) :: ip, jp, ntim, o_halo, ipstart, jpstart integer(i4b), dimension(2) :: coords - real(sp) :: f, re, g, rho, dt - real(sp), dimension(:,:), allocatable :: f_cor, h, hs, u,v, height, & + real(wp) :: f, re, g, rho, dt + real(wp), dimension(:,:), allocatable :: f_cor, h, hs, u,v, height, & dx, dy, x, y, & recqdp, recqdp_s, recqdq_s, redq_s, redq, & recq, cq_s, cq, dp1, dq, recqdq - real(sp), dimension(:), allocatable :: phi, theta, phin, thetan, u_nudge, & + real(wp), dimension(:), allocatable :: phi, theta, phin, thetan, u_nudge, & dphi, dtheta, dphin, dthetan end type grid @@ -43,7 +43,7 @@ module variables viscous_dissipation, & dissipate_h, nudge, restart integer(i4b) :: initial_winds, ip, jp, subgrid_model - real(sp) :: wind_factor, wind_shift, wind_reduce, vis, & + real(wp) :: wind_factor, wind_shift, wind_reduce, vis, & runtime, dt, output_interval, & grav, rho, Re, rotation_period_hours, scale_height, & slat, nlat, slat_thresh, nlat_thresh, nudge_timescale, & @@ -56,7 +56,7 @@ module variables !>variables for mpi type mpi_vars integer(i4b) :: rank, id, error, dx, dy - real(sp) :: wtime + real(wp) :: wtime logical, dimension(2) :: periods logical :: reorder=.true. integer(i4b) :: ndim=2 diff --git a/zbrent.f90 b/zbrent.f90 deleted file mode 100755 index 1a74d9a..0000000 --- a/zbrent.f90 +++ /dev/null @@ -1,78 +0,0 @@ - FUNCTION zbrent(func,x1,x2,tol) - USE nrtype; USE nrutil, ONLY : nrerror - IMPLICIT NONE - REAL(SP), INTENT(IN) :: x1,x2,tol - REAL(SP) :: zbrent - INTERFACE - FUNCTION func(x) - USE nrtype - IMPLICIT NONE - REAL(SP), INTENT(IN) :: x - REAL(SP) :: func - END FUNCTION func - END INTERFACE - INTEGER(I4B), PARAMETER :: ITMAX=100 - REAL(SP), PARAMETER :: EPS=epsilon(x1) - INTEGER(I4B) :: iter - REAL(SP) :: a,b,c,d,e,fa,fb,fc,p,q,r,s,tol1,xm - a=x1 - b=x2 - fa=func(a) - fb=func(b) - if ((fa > 0.0 .and. fb > 0.0) .or. (fa < 0.0 .and. fb < 0.0)) & - call nrerror('root must be bracketed for zbrent') - c=b - fc=fb - do iter=1,ITMAX - if ((fb > 0.0 .and. fc > 0.0) .or. (fb < 0.0 .and. fc < 0.0)) then - c=a - fc=fa - d=b-a - e=d - end if - if (abs(fc) < abs(fb)) then - a=b - b=c - c=a - fa=fb - fb=fc - fc=fa - end if - tol1=2.0_sp*EPS*abs(b)+0.5_sp*tol - xm=0.5_sp*(c-b) - if (abs(xm) <= tol1 .or. fb == 0.0) then - zbrent=b - RETURN - end if - if (abs(e) >= tol1 .and. abs(fa) > abs(fb)) then - s=fb/fa - if (a == c) then - p=2.0_sp*xm*s - q=1.0_sp-s - else - q=fa/fc - r=fb/fc - p=s*(2.0_sp*xm*q*(q-r)-(b-a)*(r-1.0_sp)) - q=(q-1.0_sp)*(r-1.0_sp)*(s-1.0_sp) - end if - if (p > 0.0) q=-q - p=abs(p) - if (2.0_sp*p < min(3.0_sp*xm*q-abs(tol1*q),abs(e*q))) then - e=d - d=p/q - else - d=xm - e=d - end if - else - d=xm - e=d - end if - a=b - fa=fb - b=b+merge(d,sign(tol1,xm), abs(d) > tol1 ) - fb=func(b) - end do - call nrerror('zbrent: exceeded maximum iterations') - zbrent=b - END FUNCTION zbrent