fix logic in subduction flux calculation

src/define_bc_segmented_s_line_parabola.f90

-!------------------------------------------------------------------------------|
-!------------------------------------------------------------------------------|
-!                                                                              |
-!              ||===\\                                                         |
-!              ||    \\                                                        |
-!              ||     ||   //==\\   ||  ||   //==||  ||/==\\                   |
-!              ||     ||  ||    ||  ||  ||  ||   ||  ||    ||                  |
-!              ||    //   ||    ||  ||  ||  ||   ||  ||                        |
-!              ||===//     \\==//    \\==\\  \\==\\  ||                        |
-!                                                                              |
-!------------------------------------------------------------------------------|
-!------------------------------------------------------------------------------|
-!                                                                              |
-!              DEFINE_BC_SEGMENTED_S_LINE_PARABOLA   August 2015               |
-!                                                                              |
-!------------------------------------------------------------------------------|
-!------------------------------------------------------------------------------|
-
-subroutine define_bc_segmented_s_line_parabola(params,osolve,vo,bcdef,nest)
-
-!------------------------------------------------------------------------------|
-!((((((((((((((((  Purpose of the routine  )))))))))))))))))))))))))))))))))))))
-!------------------------------------------------------------------------------|
-! This routine assigns the velocity boundary conditions for the segmented s-line
-! geometry
-
-!------------------------------------------------------------------------------|
-!((((((((((((((((  declaration of the subroutine arguments  ))))))))))))))))))))
-!------------------------------------------------------------------------------|
-
-use definitions
-!use mpi
-
-implicit none
-
-include 'mpif.h'
-
-type (parameters) params
-type (octreesolve) osolve
-type (void) vo
-type (bc_definition) bcdef
-type (nest_info) nest
-
-!------------------------------------------------------------------------------|
-!(((((((((((((((( declaration of the subroutine internal variables )))))))))))))
-!------------------------------------------------------------------------------|
-
-integer i,iproc,nproc,ierr
-double precision :: eps,lsf0,pi,lorig,h,x1,x2,phi,yend,cper,cscl,xstart,ystart
-double precision :: theta,l,vin,vzfluxscl,cntvel,dxy,xend,xsym,ymax,xwidth
-double precision :: ywidth,xdisp,ydisp,nb
-double precision,dimension(:),allocatable :: x0,ldisp
-integer ie,ij,j,jp,nelemx,nelemz
-double precision :: base,startp,endp,alpha,kink1,kink2,flytt,a,d
-double precision :: wmax,uend,vmag,dipangle
-double precision e(1),f(1),y(1)
-
-!------------------------------------------------------------------------------|
-!------------------------------------------------------------------------------|
-
-call mpi_comm_size (mpi_comm_world,nproc,ierr)
-call mpi_comm_rank (mpi_comm_world,iproc,ierr)
-
-eps=1.d-10
-osolve%kfix=0
-osolve%kfixt=0
-pi=atan(1.d0)*4.d0
-
-l=bcdef%bc_parameters(1)
-base=bcdef%bc_parameters(2)
-startp=bcdef%bc_parameters(3)
-endp=bcdef%bc_parameters(4)
-vin=bcdef%bc_parameters(5)
-alpha=bcdef%bc_parameters(6)
-kink1=bcdef%bc_parameters(7)
-kink2=bcdef%bc_parameters(8)
-
-nelemx=idint(bcdef%bc_parameters(9))
-nelemz=idint(bcdef%bc_parameters(10))
-
-nb=2**params%levelmax_oct
-dxy=1.d0/2**(params%levelmax_oct+1.d0)
-
-alpha=alpha*pi/180
-
-flytt=(kink2-kink1) * tan(alpha)
-
-a = (base-l) /((endp-startp)*(endp-startp))
-d = startp
-
-wmax=(2*vin*l)/(endp-startp)
-uend=wmax/(2*a*(endp-startp))
-vmag=sqrt(uend**2.d0+wmax**2.d0)
-
-do i=1,osolve%nnode
-
-   if (osolve%x(i).lt.eps) then
-      osolve%kfix((i-1)*3+1)=1 ; osolve%u(i)=1.d0
-      osolve%kfix((i-1)*3+2)=1 ; osolve%v(i)=0.d0
-      osolve%kfix((i-1)*3+3)=1 ; osolve%w(i)=0.d0
-   endif
-
-   if (osolve%x(i).gt.1.d0-eps) then
-      osolve%kfix((i-1)*3+1)=1 ; osolve%u(i)=0.d0
-   endif
-
-   if (osolve%y(i).lt.eps) then
-      osolve%kfix((i-1)*3+2)=1 ; osolve%v(i)=0.d0
-   endif
-
-   if (osolve%y(i).gt.1.d0-eps) then
-      osolve%kfix((i-1)*3+2)=1 ; osolve%v(i)=0.d0
-   endif
-
-   if (osolve%z(i).lt.eps) then
-      osolve%kfix((i-1)*3+1)=1 ; osolve%u(i)=0.d0
-      osolve%kfix((i-1)*3+2)=1 ; osolve%v(i)=0.d0
-      osolve%kfix((i-1)*3+3)=1 ; osolve%w(i)=0.d0
-      osolve%kfixt(i)=1 ;  osolve%temp(i)=1.d0
-
-      ! Velocity in first margin-normal convergence segment
-      if (osolve%y(i) .lt. kink1) then
-         ! Slab dip angle
-         dipangle = atan(2*a*(osolve%x(i)-startp))
-         ! Velocity before transition to subduction
-         if (osolve%x(i) .le. startp-real(nelemz)*dxy) then
-            osolve%u(i) = vin
-         ! Velocity in transition to subduction
-         elseif (osolve%x(i) .le. startp+real(nelemz)*dxy) then
-            ! Velocity contribution from incoming horizontal velocity
-            ! (diminishes to zero across transition)
-            osolve%u(i) = vin * (1.d0 - (osolve%x(i)-(startp-real(nelemz)*dxy))/(2.0*real(nelemz)*dxy))
-            ! Velocity contribution from subduction region
-            ! (increases from zero to full magnitude across transition)
-            osolve%u(i) = osolve%u(i) + vmag*cos(dipangle) * (osolve%x(i)-(startp-real(nelemz)*dxy))/(2.0*real(nelemz)*dxy)
-            osolve%w(i) = vmag*sin(dipangle) * (osolve%x(i)-(startp-real(nelemz)*dxy))/(2.0*real(nelemz)*dxy)
-         ! Velocity in subduction region
-         elseif (osolve%x(i) .le. endp-real(nelemz)*dxy) then
-            osolve%u(i) = vmag*cos(dipangle)
-            osolve%w(i) = vmag*sin(dipangle)
-         ! Velocity in transition from subduction region
-         elseif (osolve%x(i) .le. endp+real(nelemz)*dxy) then
-            ! Velocity contribution from subduction region
-            ! (decreases to zero across transition)
-            osolve%u(i) = vmag*cos(dipangle) * (1.d0 - (osolve%x(i)-(endp-real(nelemz)*dxy))/(2.0*real(nelemz)*dxy))
-            osolve%w(i) = vmag*sin(dipangle) * (1.d0 - (osolve%x(i)-(endp-real(nelemz)*dxy))/(2.0*real(nelemz)*dxy))
-         endif
-      ! Velocity in oblique convergence segment
-      elseif (osolve%y(i) .lt. kink2) then
-         ! Slab dip angle
-         dipangle = atan(2*a*(osolve%x(i)-((osolve%y(i)-kink1)*tan(alpha)+startp)))
-         ! Velocity before transition to subduction
-         if (osolve%x(i) .le. (osolve%y(i)-kink1)*tan(alpha)+startp-real(nelemz)*dxy) then
-           osolve%u(i) = vin
-         ! Velocity in transition to subduction
-         elseif (osolve%x(i) .le. (osolve%y(i)-kink1)*tan(alpha)+startp+real(nelemz)*dxy) then
-           ! Velocity contribution from incoming horizontal velocity
-           ! (diminishes to zero across transition)
-           osolve%u(i)= vin * (1.d0 - (osolve%x(i)-((osolve%y(i)-kink1)*tan(alpha)+startp-real(nelemz)*dxy))/(2.0*real(nelemz)*dxy))
-           ! Velocity contribution from subduction region
-           ! (increases from zero to full magnitude across transition)
-           osolve%u(i)= osolve%u(i) + vmag*cos(dipangle) * (osolve%x(i)-((osolve%y(i)-kink1)*tan(alpha)+startp-real(nelemz)*dxy))/(2.0*real(nelemz)*dxy)
-           osolve%w(i) = vmag*sin(dipangle) * (osolve%x(i)-((osolve%y(i)-kink1)*tan(alpha)+startp-real(nelemz)*dxy))/(2.0*real(nelemz)*dxy)
-         ! Velocity in subduction region
-         elseif (osolve%x(i) .le. (osolve%y(i)-kink1)*tan(alpha)+endp-real(nelemz)*dxy) then
-           osolve%u(i) = vmag*cos(dipangle)
-           osolve%w(i) = vmag*sin(dipangle)
-         ! Velocity in transition from subduction region
-         elseif (osolve%x(i) .le. (osolve%y(i)-kink1)*tan(alpha)+endp+real(nelemz)*dxy) then
-           ! Velocity contribution from subduction region
-           ! (decreases to zero across transition)
-           osolve%u(i) = vmag*cos(dipangle) * (1.d0 - (osolve%x(i)-((osolve%y(i)-kink1)*tan(alpha)+endp-real(nelemz)*dxy))/(2.0*real(nelemz)*dxy))
-           osolve%w(i) = vmag*sin(dipangle) * (1.d0 - (osolve%x(i)-((osolve%y(i)-kink1)*tan(alpha)+endp-real(nelemz)*dxy))/(2.0*real(nelemz)*dxy))
-         endif
-      ! Velocity in second margin-normal convergence segment
-      else
-        ! Slab dip angle
-        dipangle = atan(2*a*(osolve%x(i)-(flytt+startp)))
-        ! Velocity before transition to subduction
-        if (osolve%x(i) .le. flytt+startp-real(nelemz)*dxy) then
-          osolve%u(i) = vin
-        ! Velocity in transition to subduction
-        elseif (osolve%x(i) .le. flytt+startp+real(nelemz)*dxy) then
-          ! Velocity contribution from incoming horizontal velocity
-          ! (diminishes to zero across transition)
-          osolve%u(i) = vin * (1.d0 - (osolve%x(i)-((flytt+startp)-real(nelemz)*dxy))/(2.0*real(nelemz)*dxy))
-          ! Velocity contribution from subduction region
-          ! (increases from zero to full magnitude across transition)
-          osolve%u(i) = osolve%u(i) + vmag * cos(dipangle) * (osolve%x(i)-((flytt+startp)-real(nelemz)*dxy))/(2.0*real(nelemz)*dxy)
-          osolve%w(i) = vmag * sin(dipangle) * (osolve%x(i)-((flytt+startp)-real(nelemz)*dxy))/(2.0*real(nelemz)*dxy)
-        ! Velocity in subduction region
-        elseif (osolve%x(i) .le. flytt+endp-real(nelemz)*dxy) then
-          osolve%u(i) = vmag * cos(dipangle)
-          osolve%w(i) = vmag * sin(dipangle)
-        ! Velocity in transition from subduction region
-        elseif (osolve%x(i) .le. flytt+endp+real(nelemz)*dxy) then
-          ! Velocity contribution from subduction region
-          ! (decreases to zero across transition)
-          osolve%u(i) = vmag * cos(dipangle) * (1.d0 - (osolve%x(i)-((flytt+endp)-real(nelemz)*dxy))/(2.0*real(nelemz)*dxy))
-          osolve%w(i) = vmag * sin(dipangle) * (1.d0 - (osolve%x(i)-((flytt+endp)-real(nelemz)*dxy))/(2.0*real(nelemz)*dxy))
-        endif
-      endif
-   endif
-
-   if (osolve%z(i).gt.1.d0-eps) then
-      osolve%kfix((i-1)*3+1)=1 ; osolve%u(i)=0.d0
-      osolve%kfix((i-1)*3+2)=1 ; osolve%v(i)=0.d0
-      osolve%kfix((i-1)*3+3)=1 ; osolve%w(i)=0.d0
-      osolve%kfixt(i)=1 ;  osolve%temp(i)=0.d0
-   endif
-
-   if (.not.vo%influid(i)) then
-      osolve%kfixt(i)=1
-      osolve%temp(i)=0.d0
-   endif
-   if (abs(bcdef%utrans).gt.eps) then
-      if (osolve%kfix((i-1)*3+1)==1) osolve%u(i)=osolve%u(i)+bcdef%utrans
-   endif
-   if (abs(bcdef%vtrans).gt.eps) then
-      if (osolve%kfix((i-1)*3+2)==1) osolve%v(i)=osolve%v(i)+bcdef%vtrans
-   endif
-   if (osolve%kfix((i-1)*3+1)==1) osolve%u(i)=osolve%u(i)*vin
-   if (osolve%kfix((i-1)*3+2)==1) osolve%v(i)=osolve%v(i)*vin
-   if (osolve%kfix((i-1)*3+3)==1) osolve%w(i)=osolve%w(i)*vin
-enddo
-
-if (params%isobc) then
-  call define_isostasy_bc(params,osolve,bcdef)
-endif
-
-end
-!------------------------------------------------------------------------------|
+!------------------------------------------------------------------------------|
+!------------------------------------------------------------------------------|
+!                                                                              |
+!              ||===\\                                                         |
+!              ||    \\                                                        |
+!              ||     ||   //==\\   ||  ||   //==||  ||/==\\                   |
+!              ||     ||  ||    ||  ||  ||  ||   ||  ||    ||                  |
+!              ||    //   ||    ||  ||  ||  ||   ||  ||                        |
+!              ||===//     \\==//    \\==\\  \\==\\  ||                        |
+!                                                                              |
+!------------------------------------------------------------------------------|
+!------------------------------------------------------------------------------|
+!                                                                              |
+!              DEFINE_BC_SEGMENTED_S_LINE_PARABOLA   August 2015               |
+!                                                                              |
+!------------------------------------------------------------------------------|
+!------------------------------------------------------------------------------|
+
+subroutine define_bc_segmented_s_line_parabola(params,osolve,vo,bcdef,nest)
+
+!------------------------------------------------------------------------------|
+!((((((((((((((((  Purpose of the routine  )))))))))))))))))))))))))))))))))))))
+!------------------------------------------------------------------------------|
+! This routine assigns the velocity boundary conditions for the segmented s-line
+! geometry
+
+!------------------------------------------------------------------------------|
+!((((((((((((((((  declaration of the subroutine arguments  ))))))))))))))))))))
+!------------------------------------------------------------------------------|
+
+use definitions
+!use mpi
+
+implicit none
+
+include 'mpif.h'
+
+type (parameters) params
+type (octreesolve) osolve
+type (void) vo
+type (bc_definition) bcdef
+type (nest_info) nest
+
+!------------------------------------------------------------------------------|
+!(((((((((((((((( declaration of the subroutine internal variables )))))))))))))
+!------------------------------------------------------------------------------|
+
+integer i,iproc,nproc,ierr
+double precision :: eps,lsf0,pi,lorig,h,x1,x2,phi,yend,cper,cscl,xstart,ystart
+double precision :: theta,l,vin,vzfluxscl,cntvel,dxy,xend,xsym,ymax,xwidth
+double precision :: ywidth,xdisp,ydisp,nb
+double precision,dimension(:),allocatable :: x0,ldisp
+integer ie,ij,j,jp,nelemx,nelemz
+double precision :: base,startp,endp,alpha,kink1,kink2,flytt,a,d
+double precision :: wmax,uend,vmag,dipangle
+double precision e(1),f(1),y(1)
+
+!------------------------------------------------------------------------------|
+!------------------------------------------------------------------------------|
+
+call mpi_comm_size (mpi_comm_world,nproc,ierr)
+call mpi_comm_rank (mpi_comm_world,iproc,ierr)
+
+eps=1.d-10
+osolve%kfix=0
+osolve%kfixt=0
+pi=atan(1.d0)*4.d0
+
+l=bcdef%bc_parameters(1)
+base=bcdef%bc_parameters(2)
+startp=bcdef%bc_parameters(3)
+endp=bcdef%bc_parameters(4)
+vin=bcdef%bc_parameters(5)
+alpha=bcdef%bc_parameters(6)
+kink1=bcdef%bc_parameters(7)
+kink2=bcdef%bc_parameters(8)
+
+nelemx=idint(bcdef%bc_parameters(9))
+nelemz=idint(bcdef%bc_parameters(10))
+
+nb=2**params%levelmax_oct
+dxy=1.d0/2**(params%levelmax_oct+1.d0)
+
+alpha=alpha*pi/180
+
+flytt=(kink2-kink1) * tan(alpha)
+
+a = (base-l) /((endp-startp)*(endp-startp))
+d = startp
+
+wmax=(2*(vin+bcdef%utrans)*l)/(endp-startp)
+uend=wmax/(2*a*(endp-startp))
+vmag=sqrt(uend**2.d0+wmax**2.d0)
+
+do i=1,osolve%nnode
+
+   if (osolve%x(i).lt.eps) then
+      osolve%kfix((i-1)*3+1)=1 ; osolve%u(i)=1.d0
+      osolve%kfix((i-1)*3+2)=1 ; osolve%v(i)=0.d0
+      osolve%kfix((i-1)*3+3)=1 ; osolve%w(i)=0.d0
+   endif
+
+   if (osolve%x(i).gt.1.d0-eps) then
+      osolve%kfix((i-1)*3+1)=1 ; osolve%u(i)=0.d0
+   endif
+
+   if (osolve%y(i).lt.eps) then
+      osolve%kfix((i-1)*3+2)=1 ; osolve%v(i)=0.d0
+   endif
+
+   if (osolve%y(i).gt.1.d0-eps) then
+      osolve%kfix((i-1)*3+2)=1 ; osolve%v(i)=0.d0
+   endif
+
+   if (osolve%z(i).lt.eps) then
+      osolve%kfix((i-1)*3+1)=1 ; osolve%u(i)=0.d0
+      osolve%kfix((i-1)*3+2)=1 ; osolve%v(i)=0.d0
+      osolve%kfix((i-1)*3+3)=1 ; osolve%w(i)=0.d0
+      osolve%kfixt(i)=1 ;  osolve%temp(i)=1.d0
+
+      ! Velocity in first margin-normal convergence segment
+      if (osolve%y(i) .lt. kink1) then
+         ! Slab dip angle
+         dipangle = atan(2*a*(osolve%x(i)-startp))
+         ! Velocity before transition to subduction
+         if (osolve%x(i) .le. startp-real(nelemz)*dxy) then
+            osolve%u(i) = vin
+         ! Velocity in transition to subduction
+         elseif (osolve%x(i) .le. startp+real(nelemz)*dxy) then
+            ! Velocity contribution from incoming horizontal velocity
+            ! (diminishes to zero across transition)
+            osolve%u(i) = vin * (1.d0 - (osolve%x(i)-(startp-real(nelemz)*dxy))/(2.0*real(nelemz)*dxy))
+            ! Velocity contribution from subduction region
+            ! (increases from zero to full magnitude across transition)
+            osolve%u(i) = osolve%u(i) + vmag*cos(dipangle) * (osolve%x(i)-(startp-real(nelemz)*dxy))/(2.0*real(nelemz)*dxy)
+            osolve%w(i) = vmag*sin(dipangle) * (osolve%x(i)-(startp-real(nelemz)*dxy))/(2.0*real(nelemz)*dxy)
+         ! Velocity in subduction region
+         elseif (osolve%x(i) .le. endp-real(nelemz)*dxy) then
+            osolve%u(i) = vmag*cos(dipangle)
+            osolve%w(i) = vmag*sin(dipangle)
+         ! Velocity in transition from subduction region
+         elseif (osolve%x(i) .le. endp+real(nelemz)*dxy) then
+            ! Velocity contribution from subduction region
+            ! (decreases to zero across transition)
+            osolve%u(i) = vmag*cos(dipangle) * (1.d0 - (osolve%x(i)-(endp-real(nelemz)*dxy))/(2.0*real(nelemz)*dxy))
+            osolve%w(i) = vmag*sin(dipangle) * (1.d0 - (osolve%x(i)-(endp-real(nelemz)*dxy))/(2.0*real(nelemz)*dxy))
+         endif
+      ! Velocity in oblique convergence segment
+      elseif (osolve%y(i) .lt. kink2) then
+         ! Slab dip angle
+         dipangle = atan(2*a*(osolve%x(i)-((osolve%y(i)-kink1)*tan(alpha)+startp)))
+         ! Velocity before transition to subduction
+         if (osolve%x(i) .le. (osolve%y(i)-kink1)*tan(alpha)+startp-real(nelemz)*dxy) then
+           osolve%u(i) = vin
+         ! Velocity in transition to subduction
+         elseif (osolve%x(i) .le. (osolve%y(i)-kink1)*tan(alpha)+startp+real(nelemz)*dxy) then
+           ! Velocity contribution from incoming horizontal velocity
+           ! (diminishes to zero across transition)
+           osolve%u(i)= vin * (1.d0 - (osolve%x(i)-((osolve%y(i)-kink1)*tan(alpha)+startp-real(nelemz)*dxy))/(2.0*real(nelemz)*dxy))
+           ! Velocity contribution from subduction region
+           ! (increases from zero to full magnitude across transition)
+           osolve%u(i)= osolve%u(i) + vmag*cos(dipangle) * (osolve%x(i)-((osolve%y(i)-kink1)*tan(alpha)+startp-real(nelemz)*dxy))/(2.0*real(nelemz)*dxy)
+           osolve%w(i) = vmag*sin(dipangle) * (osolve%x(i)-((osolve%y(i)-kink1)*tan(alpha)+startp-real(nelemz)*dxy))/(2.0*real(nelemz)*dxy)
+         ! Velocity in subduction region
+         elseif (osolve%x(i) .le. (osolve%y(i)-kink1)*tan(alpha)+endp-real(nelemz)*dxy) then
+           osolve%u(i) = vmag*cos(dipangle)
+           osolve%w(i) = vmag*sin(dipangle)
+         ! Velocity in transition from subduction region
+         elseif (osolve%x(i) .le. (osolve%y(i)-kink1)*tan(alpha)+endp+real(nelemz)*dxy) then
+           ! Velocity contribution from subduction region
+           ! (decreases to zero across transition)
+           osolve%u(i) = vmag*cos(dipangle) * (1.d0 - (osolve%x(i)-((osolve%y(i)-kink1)*tan(alpha)+endp-real(nelemz)*dxy))/(2.0*real(nelemz)*dxy))
+           osolve%w(i) = vmag*sin(dipangle) * (1.d0 - (osolve%x(i)-((osolve%y(i)-kink1)*tan(alpha)+endp-real(nelemz)*dxy))/(2.0*real(nelemz)*dxy))
+         endif
+      ! Velocity in second margin-normal convergence segment
+      else
+        ! Slab dip angle
+        dipangle = atan(2*a*(osolve%x(i)-(flytt+startp)))
+        ! Velocity before transition to subduction
+        if (osolve%x(i) .le. flytt+startp-real(nelemz)*dxy) then
+          osolve%u(i) = vin
+        ! Velocity in transition to subduction
+        elseif (osolve%x(i) .le. flytt+startp+real(nelemz)*dxy) then
+          ! Velocity contribution from incoming horizontal velocity
+          ! (diminishes to zero across transition)
+          osolve%u(i) = vin * (1.d0 - (osolve%x(i)-((flytt+startp)-real(nelemz)*dxy))/(2.0*real(nelemz)*dxy))
+          ! Velocity contribution from subduction region
+          ! (increases from zero to full magnitude across transition)
+          osolve%u(i) = osolve%u(i) + vmag * cos(dipangle) * (osolve%x(i)-((flytt+startp)-real(nelemz)*dxy))/(2.0*real(nelemz)*dxy)
+          osolve%w(i) = vmag * sin(dipangle) * (osolve%x(i)-((flytt+startp)-real(nelemz)*dxy))/(2.0*real(nelemz)*dxy)
+        ! Velocity in subduction region
+        elseif (osolve%x(i) .le. flytt+endp-real(nelemz)*dxy) then
+          osolve%u(i) = vmag * cos(dipangle)
+          osolve%w(i) = vmag * sin(dipangle)
+        ! Velocity in transition from subduction region
+        elseif (osolve%x(i) .le. flytt+endp+real(nelemz)*dxy) then
+          ! Velocity contribution from subduction region
+          ! (decreases to zero across transition)
+          osolve%u(i) = vmag * cos(dipangle) * (1.d0 - (osolve%x(i)-((flytt+endp)-real(nelemz)*dxy))/(2.0*real(nelemz)*dxy))
+          osolve%w(i) = vmag * sin(dipangle) * (1.d0 - (osolve%x(i)-((flytt+endp)-real(nelemz)*dxy))/(2.0*real(nelemz)*dxy))
+        endif
+      endif
+   endif
+
+   if (osolve%z(i).gt.1.d0-eps) then
+      osolve%kfix((i-1)*3+1)=1 ; osolve%u(i)=0.d0
+      osolve%kfix((i-1)*3+2)=1 ; osolve%v(i)=0.d0
+      osolve%kfix((i-1)*3+3)=1 ; osolve%w(i)=0.d0
+      osolve%kfixt(i)=1 ;  osolve%temp(i)=0.d0
+   endif
+
+   if (.not.vo%influid(i)) then
+      osolve%kfixt(i)=1
+      osolve%temp(i)=0.d0
+   endif
+   if (abs(bcdef%utrans).gt.eps) then
+      if (osolve%kfix((i-1)*3+1)==1) osolve%u(i)=osolve%u(i)+bcdef%utrans
+   endif
+   if (abs(bcdef%vtrans).gt.eps) then
+      if (osolve%kfix((i-1)*3+2)==1) osolve%v(i)=osolve%v(i)+bcdef%vtrans
+   endif
+   if (osolve%kfix((i-1)*3+1)==1) osolve%u(i)=osolve%u(i)*vin
+   if (osolve%kfix((i-1)*3+2)==1) osolve%v(i)=osolve%v(i)*vin
+   if (osolve%kfix((i-1)*3+3)==1) osolve%w(i)=osolve%w(i)*vin
+enddo
+
+if (params%isobc) then
+  call define_isostasy_bc(params,osolve,bcdef)
+endif
+
+end
+!------------------------------------------------------------------------------|
\ No newline at end of file