text:
!
!     CalculiX - A 3-dimensional finite element program
!              Copyright (C) 1998-2019 Guido Dhondt
!
!     This program is free software; you can redistribute it and/or
!     modify it under the terms of the GNU General Public License as
!     published by the Free Software Foundation(version 2);
!
!
!     This program is distributed in the hope that it will be useful,
!     but WITHOUT ANY WARRANTY; without even the implied warranty of
!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!     GNU General Public License for more details.
!
!     You should have received a copy of the GNU General Public License
!     along with this program; if not, write to the Free Software
!     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
!
!
c>   \brief Calculating the slave-master triangulation
c>    and generating the integration points needed for the calculation of the coupling matrices
c>    for details see phd-thesis Sitzmann Appendix A
c>
c> @param   [in]     tieset      name and dependent surface of tie set
c> @param   [in]     ntie        number of contraints
c> @param   [in]     itietri     (1,i)pointer to node where trangulation starts for i (2,i) pointer to end
c> @param   [in]     ipkon       pointer into field kon
c> @param   [in]     kon         Field containing the connectivity of the elements in succesive order
c> @param   [in]     lakon	 (i) label for element i
c> @param   [in]     set         (i)name of set_i
c> @param   [in]     cg          field containing centers of gravity
c> @param   [in]     straight    (1:4 5:8 9:13,i)coeffs of plane equation for edges of triagle_i (13:16,i) coeffs of plane containing triagle
c> @param   [out]    nintpoint   number of generated integration points
c> @param   [in]     koncont     (1:3,i) nodes of triagle_i (4,i) element face
c> @param   [in]     co          field containing the coordinates of all nodes
c> @param   [in]     vold        field containing the displacements
c> @param   [in,out] x,y,z       ONLY HELP FIELD
c> @param   [in,out] xo,yo,zo    ONLY HELP FIELD
c> @param   [in,out] nx,ny,nz    ONLY HELP FIELD
c> @param   [in]     nset	 number of sets
c> @param   [in]     iinc        current increment
c> @param   [in]     iit         current iteration
c> @param   [in,out] islavsurf   islavsurf(1,i) slaveface i islavsurf(2,i) # integration points generated before looking at face i
c> @param   [out]    imastsurf   pointer into pmastsurf
c> @param   [out]    pmastsurf   field storing position and weights for integration points on master side
c> @param   [in]     itiefac     pointer into field islavsurf: (1,i) beginning slave_i (2,i) end of slave_i
C> @param   [in]     islavnode   fields containing nodes of slace surfaces
c> @param   [in]     nslavnode   (i) for contraint i pointer into field islavnode
C> @param   [in]     slavnor     normal vektors in the nods of slave surface
C> @param   [in]     slavtan      tangetial vektors in the nodes of slave surface NOT USED!
c> @param   [in]     imastop     (l,i) for edge l in triagle i neightbouring triangle
c> @param   [out]    gapmints    stores gaps between master and slave side
c> @param   [in,out] islavact	 (i) indicates, if slave node i is active (=-3 no-slave-node, =-2 no-LM-node, =-1 no-gap-node, =0 inactive node, =1 sticky node, =2 slipping/active node)
c> @param   [in]     mi		 (1) max # of integration points per element (2) max degree of freedom per element
c> @param   [in]     ncont	 number of triangles in triagulation of master side
c> @param   [in]     ipe	 (i) pointer to ime for node i
c> @param   [in]     ime         ... cataloging the edges with node i
c> @param   [out]    pslavsurf   field storing position xil, etal and weight for integration point on slave side
c> @param   [in]     i           current tie
c> @param   [in]     l           current face
c> @param   [in]     ntri        # triangles
c> @param   [in]     tietol      field containing the clearance definitions
c> @param   [in]     reltime     relative time within step, needed for shrink
c> @param   [in]     nmethod     analysis method, needed for shrink
c>
      subroutine slavintmortar(tieset,ntie,itietri,ipkon,kon,
     &     lakon,set,cg,straight,nintpoint,
     &     koncont,co,vold,xo,yo,zo,x,y,z,nx,ny,nz,nset,
     &     iinc,iit,
     &     islavsurf,imastsurf,pmastsurf,itiefac,
     &     islavnode,nslavnode,slavnor,slavtan,imastop,gapmints,
     &     islavact,mi,ncont,ipe,ime,pslavsurf,i,l,ntri,tietol,
     &     reltime,nmethod)

!     Author: Sitzmann, Saskia
!
      implicit none
!
      logical debug,nogap,shrink
!
      character*8 lakon(*)
      character*81 tieset(3,*),set(*)
!
      integer ntie,nset,nintpoint,imastop(3,*),ncont,
     &     itietri(2,ntie),ipkon(*),kon(*),koncont(4,*),node,
     &     neigh(10),iflag,kneigh,i,j,k,l,ii,
     &     itri,ipos,nx(*),ny(*),
     &     nz(*),index1,nmethod,
     &     nelemm,jfacem,indexe,iit,iinc,
     &     nnodelem,nope,m1,
     &     islavsurf(2,*),islavnode(*),nslavnode(ntie+1),
     &     imastsurf(*),itiefac(2,*),ifaces,nelems,jfaces,mi(*),
     &     m,nopes,konl(26),id,islavact(*),indexnode(8),
     &     imface(8),ntria,imfacecorner(8,8),line,
     &     iactiveline(3,3*ncont),icoveredmelem(3*ncont),
     &     nactiveline,ipe(*),ime(4,*),k1,j1,ncoveredmelem,
     &     info,ntri,nintpfirst,nodem(8),nodem2(8),ithree,
     &     il,ifac,getiface,ifacem,idummy,nopemm,nmp,k2,j2
!
      real*8 cg(3,*),straight(16,*),co(3,*),vold(0:mi(2),*),
     &     xo(*),yo(*),zo(*),x(*),y(*),z(*),
     &     pmastsurf(2,*),xl2m(3,8),xl2s(3,8),xl2m2(3,8),
     &     pmiddle(3),xl2sr(3,8),xl2sp(3,8),xl2s2(3,8),
     &     slavnor(3,*),slavtan(6,*),dd,xns(3,8),areaslav,
     &     al,xn(3),gapmints(*),slavstraight(36),
     &     err2,dist,distmin,tietol(3,*),clearance,reltime,
     &     pslavsurf(3,*),err,xquad(2,8),cl2s(3,8),cl2m(3,8),
     &     xtri(2,6),xi,et,xsj2(3),xs2(3,2),shp2(7,8),anglesm
!
      include "gauss.f"
!
c      debug=.true.
      debug=.false.
      data iflag /2/
      data ithree /3/
!
      data xquad /-1, -1,
     &     1, -1,
     &     1, 1,
     &     -1, 1,
     &     0, -1,
     &     1, 0,
     &     0, 1,
     &     -1, 0/
!
      data xtri /0, 0,
     &     1, 0,
     &     0, 1,
     &     0.5, 0,
     &     0.5, 0.5,
     &     0, 0.5/
!
      kneigh=1
      err=0.1
      areaslav=0.0
      nintpfirst=nintpoint
      islavsurf(2,l)=nintpoint
      if(debug)WRITE(30,*) '#SLAVINTMORTAR iinc',iinc, 'face',l
      if(debug)WRITE(20,*) '#SLAVINTMORTAR iinc',iinc, 'face',l
c     WRITE(*,*) '#SLAVINTMORTAR iit',iit, 'face',l
!
!     get clearance and shrink
!
      clearance=tietol(3,i)
      if(nmethod.eq.4)then
         shrink=.false.
      else
         shrink=.true.
      endif
c      shrink=.false.

c      write(*,*)'tietol',tietol(1,i),tietol(2,i),tietol(3,i)
!
!     Research of the contact integration points
!
      ifaces = islavsurf(1,l)
      nelems = int(ifaces/10)
      jfaces = ifaces - nelems*10
c      if(nelems.eq.39816 .or.nelems.eq.40364)debug=.true.
!
!     get nope,nopes
!
      call getnumberofnodes(nelems,jfaces,lakon,nope,nopes,idummy)
!
!     actual position of the nodes belonging to the
!     slave surface
!
      do j=1,nope
         konl(j)=kon(ipkon(nelems)+j)
      enddo
!
      do m=1,nopes
         do j=1,3
            ifac=getiface(m,jfaces,nope)
            xl2s(j,m)=co(j,konl(ifac))+
     &           vold(j,konl(ifac))
            cl2s(j,m)=co(j,konl(ifac))
         enddo
c         write(*,*)'vold',(vold(k,konl(ifac)),k=1,3)
      enddo
!
!     slightly reducing the size of the slave surface in
!     an aleatoric way needed for contact search
!
      do j=1,3
         pmiddle(j)=0.d0
         do m=1,nopes
            pmiddle(j)=pmiddle(j)+xl2s(j,m)
         enddo
         pmiddle(j)=pmiddle(j)/nopes
      enddo
      do j=1,3
         do m=1,nopes
            xl2sr(j,m)=xl2s(j,m)-0.05*err*(xl2s(j,m)-pmiddle(j))
c     xl2sr(j,m)=xl2s(j,m)
         enddo
      enddo

c     write(*,*)'slavm xn',(xn(k),k=1,3)
!
!
!     Resort vertices for quadratic elements
!
      if(nopes.eq.3 .or. nopes.eq.4)then
         do j=1, nopes
            do k=1,3
               xl2s2(k,j)= xl2s(k,j)
            enddo
         enddo
      else
         do j=1,int(nopes/2.0)
            do k=1,3
               xl2s2(k,2*j-1)= xl2s(k,j)
               xl2s2(k,2*j)= xl2s(k,(int(nopes/2.0))+j)
            enddo
         enddo
      endif
c      do j=1,nopes
c       write(*,*)'xs',(xl2s(k,j),k=1,3)
c       write(*,*)'cs',(cl2s(k,j),k=1,3)
c      enddo
c     write(20,*)'*******************************'
c     do j=1,nopes
c     write(20,*)'xs2',(xl2s2(k,j),k=1,3)
c     enddo
!
!     calculate the mean normal vector on the Slave Surface
!     +
!     determine the equations of the triangle/quadrilateral
!     (mean)plane and of the planes boardering the
!     triangle/quadrilateral
!
      if(nopes.eq.3) then
         call straighteq3d(xl2s2,slavstraight)
         do k=1,3
            xn(k)= slavstraight(4*nopes+k)
         enddo
      else
         do k=1,3
            xn(k)=0.d0
         enddo
!
         do m = 1, nopes
            if(nopes.eq.4 .or. nopes.eq.8)then
               xi = xquad(1,m)
               et = xquad(2,m)
            else
               xi = xtri(1,m)
               et = xtri(2,m)
            endif
            if(nopes.eq.8)then
               call shape8q(xi,et,xl2s,xsj2,xs2,shp2,iflag)
            elseif(nopes.eq.4)then
               call shape4q(xi,et,xl2s,xsj2,xs2,shp2,iflag)
            elseif(nopes.eq.6)then
               call shape6tri(xi,et,xl2s,xsj2,xs2,shp2,iflag)
            else
               call shape3tri(xi,et,xl2s,xsj2,xs2,shp2,iflag)
            endif
            dd = dsqrt(xsj2(1)*xsj2(1) + xsj2(2)*xsj2(2)
     &           + xsj2(3)*xsj2(3))
            xsj2(1) = xsj2(1)/dd
            xsj2(2) = xsj2(2)/dd
            xsj2(3) = xsj2(3)/dd
!
            do k=1,3
               xn(k) = xn(k)
     &              +xsj2(k)
            enddo
         enddo
!
!     normalizing the mean normal on the Slave surface
!
         dd=dsqrt(xn(1)**2+xn(2)**2+xn(3)**2)
         do k=1,3
            xn(k)=xn(k)/dd
         enddo
         call approxplane(xl2s2,slavstraight,xn,nopes)
      endif
c      call exit(201)
!
!     Project slave nodes to meanplane, needed for Sutherland-Hodgman
!
      do j=1, nopes
         al=-xn(1)*xl2s2(1,j)-xn(2)*
     &        xl2s2(2,j)-xn(3)*xl2s2(3,j)-
     &        slavstraight(nopes*4+4)
         if(nopes.eq.3)then
            do k=1,3
               xl2sp(k,j)= xl2s2(k,j)
            enddo
         else
            do k=1,3
               xl2sp(k,j)= xl2s2(k,j)+al*xn(k)
            enddo
         endif
      enddo
c      do j=1,nopes
c      write(20,*)'SIM:xsp',(xl2sp(k,j),k=1,3)
c      enddo
!
!     determine the triangles corresponding to the corner
!     nodes
!
      ntria=0
      do j=1,8
         imface(j)=0
         do k=1,8
            imfacecorner(j,k)=0
         enddo
      enddo
!     check for nogap-nodes
      do j=1,nopes
!        here unreduced coordinates are used
         call neartriangle(xl2s(1,j),xn,xo,yo,zo,x,y,z,nx,ny,nz,
     &        ntri,neigh,kneigh,itietri,ntie,straight,imastop,itri,i)
         ifac= getiface(j,jfaces,nope)
         node= konl(ifac)
         if(debug) then
            write(20,*) neigh(1),neigh(1)+itietri(1,i)-1
            write(20,*) 'itri',itri,'node',node
         endif
         call nident(islavnode(nslavnode(i)+1), node,
     &        nslavnode(i+1)-nslavnode(i), id)
         if(itri.eq.0) then
!     no triangle found
c            if(islavact(nslavnode(i)+id).gt.-1)then
!     node was possibly active in last increment
               islavact(nslavnode(i)+id)=-1
c            endif
         else
!     triangle found
      if(debug)write(20,*) ' node',node,islavnode(nslavnode(i)+id),
     &            id,islavact(nslavnode(i)+id)
               if(islavact(nslavnode(i)+id).lt.0) then
                  islavact(nslavnode(i)+id)=0
               endif
      if(debug)write(20,*) ' node',node,islavnode(nslavnode(i)+id),
     &            id,islavact(nslavnode(i)+id)
         endif
      enddo
c
      do j=1,3
         do m=1,nopes
            xl2sr(j,m)=xl2s(j,m)-2*err*(xl2s(j,m)-pmiddle(j))
         enddo
      enddo
      distmin=1.1
      do j=1,nopes
!        for initial master elements reduced coordinates are used
         call neartriangle(xl2sr(1,j),xn,xo,yo,zo,x,y,z,nx,ny,nz,
     &        ntri,neigh,kneigh,itietri,ntie,straight,imastop,itri,i)
         ifac= getiface(j,jfaces,nope)
         node= konl(ifac)
         if(itri.eq.0) then
            cycle
         endif
         dist= -(straight(13,itri)*xl2sr(1,j)+
     &        straight(14,itri)*xl2sr(2,j)+
     &        straight(15,itri)*xl2sr(3,j)+
     &        straight(16,itri))/
     &        (straight(13,itri)*xn(1)+
     &        straight(14,itri)*xn(2)+
     &        straight(15,itri)*xn(3))
         if(dist.lt.distmin)distmin=dist
         if(debug)write(20,*) 'j',j,'dist',dist,distmin
         ifacem=koncont(4,itri)
         if(debug)write(20,*)'noder ',node, 'itri',itri,
     &        'ifacem',ifacem
!
         call nident(imface,ifacem,ntria,id)
         if(id.gt.0) then
            if(imface(id).eq.ifacem) then
               imfacecorner(j,id)=1
               cycle
            endif
         endif
!
!     triangle was not covered yet: add to stack
!
!     angle criteria
!
         anglesm=xn(1)*straight(13,itri)
     &        +xn(2)*straight(14,itri)
     &        +xn(3)*straight(15,itri)
         if(debug)write(20,*)'cos alpa',anglesm
!
!     check angle to avoid problem with the dual gap
         if(anglesm.lt.-0.7)then
!     angle between surface normals between 135 and 225 degree
!     angle between surfaces between 0 and 45 degree
            ntria=ntria+1
            do k=ntria,id+2,-1
               imface(k)=imface(k-1)
               do m=1,j-1
                  imfacecorner(m,k)=imfacecorner(m,k-1)
               enddo
            enddo
            imface(id+1)=ifacem
            imfacecorner(j,id+1)=1
            do m=1,j-1
               imfacecorner(m,id+1)=0
            enddo
         endif
      enddo
      if(debug)then
         write(20,*)'ifacem n1:n8 face', l,ntria
         do k=1,ntria
            write(20,*) imface (k),imfacecorner(1:nopes,k)
         enddo
      endif
!
      nactiveline=0
!
!     treating the corner elements first
!
      ncoveredmelem=0
      do j=1,ntria
         if(debug)write(20,*) 'corner triangle j',j
         ifacem=imface(j)
         nelemm=int(ifacem/10.d0)
         jfacem=ifacem-10*nelemm
!
!        add master element to covered stack
!
         call nident(icoveredmelem,nelemm, ncoveredmelem,id)
         if(id.ne.0 .and. icoveredmelem(id).eq.nelemm)then
!     master elment was already threated
            cycle
         else
!     add master element to covered elements
            ncoveredmelem=ncoveredmelem+1
            do ii=ncoveredmelem,id+2,-1
               icoveredmelem(ii)=icoveredmelem(ii-1)
            enddo
            icoveredmelem(id+1)=nelemm
         endif
         if(debug)write(20,*)itri,imface(j), ifacem,nelemm,jfacem
         call getnumberofnodes(nelemm,jfacem,lakon,nopemm,
     &        nnodelem,idummy)
!
!     determining the nodes of the face
!
         do j1=1,nopemm
            konl(j1)=kon(ipkon(nelemm)+j1)
         enddo
         do k1=1,nnodelem
            ifac=getiface(k1,jfacem,nopemm)
            nodem(k1)=konl(ifac)
            do j1=1,3
               xl2m(j1,k1)=co(j1,konl(ifac))+
     &              vold(j1,konl(ifac))
               cl2m(j1,k1)=co(j1,konl(ifac))
            enddo
         enddo
         dd=dsqrt(xn(1)**2+xn(2)**2+xn(3)**2)
         if(debug)then
            write(20,*) 'dd',dd
            write(20,100)(xn(k),k=1,3)
            write(20,100)(slavstraight(nopes*4+k),k=1,3)
            write(20,*) 'SIM: xl2s'
            do j1=1,nopes
               write(20,*)(xl2s(k,j1),k=1,3)
            enddo
            write(20,*) 'SIM: xl2m'
            do j1=1,nnodelem
               write(20,*)(xl2m(k,j1),k=1,3)
            enddo
         endif
 100     format('SIM: xns',3(3x,e15.8))
 101     format(3(e15.8))
!
         if(debug) write(20,*) 'TT: itri',nelemm
!
!     divide master element into konvex subelements
!
         if(nnodelem.eq.3 .or.nnodelem.eq.4)then
!
!     no loop needed
!
            nmp=nnodelem
            do k2=1,nnodelem
               nodem2(k2)=nodem(k2)
               do j2=1,3
                  xl2m2(j2,k2)=xl2m(j2,k2)
               enddo
            enddo
            call treatmasterface_mortar(
     &           nopes,slavstraight,xn,xns,co,xl2s,xl2sp,
     &           ipe,ime,iactiveline,nactiveline,
     &           ifacem,
     &           nintpoint,pslavsurf,imastsurf,pmastsurf,
     &           xl2m,nnodelem,xl2m2,nmp,nodem2,mi,
     &           gapmints,l,areaslav,debug,clearance,
     &           cl2s,cl2m,shrink,reltime)
!
         elseif(nnodelem.eq.6)then
!
!     tri6 surface is divided into 4 tri3 surfaces
!
!     1. triangle
c            write(20,*)'*********tria 1**************'
            nmp=3
            nodem2(1)=nodem(1)
            nodem2(2)=nodem(4)
            nodem2(3)=nodem(6)
            do j2=1,3
               xl2m2(j2,1)=xl2m(j2,1)
            enddo
            do j2=1,3
               xl2m2(j2,2)=xl2m(j2,4)
            enddo
            do j2=1,3
               xl2m2(j2,3)=xl2m(j2,6)
            enddo
            call treatmasterface_mortar(
     &           nopes,slavstraight,xn,xns,co,xl2s,xl2sp,
     &           ipe,ime,iactiveline,nactiveline,
     &           ifacem,
     &           nintpoint,pslavsurf,imastsurf,pmastsurf,
     &           xl2m,nnodelem,xl2m2,nmp,nodem2,mi,
     &           gapmints,l,areaslav,debug,clearance,
     &           cl2s,cl2m,shrink,reltime)
!
!     2. triangle
!
c            write(20,*)'*********tria 2**************'
            nmp=3
            nodem2(1)=nodem(4)
            nodem2(2)=nodem(2)
            nodem2(3)=nodem(5)
            do j2=1,3
               xl2m2(j2,1)=xl2m(j2,4)
            enddo
            do j2=1,3
               xl2m2(j2,2)=xl2m(j2,2)
            enddo
            do j2=1,3
               xl2m2(j2,3)=xl2m(j2,5)
            enddo
            call treatmasterface_mortar(
     &           nopes,slavstraight,xn,xns,co,xl2s,xl2sp,
     &           ipe,ime,iactiveline,nactiveline,
     &           ifacem,
     &           nintpoint,pslavsurf,imastsurf,pmastsurf,
     &           xl2m,nnodelem,xl2m2,nmp,nodem2,mi,
     &           gapmints,l,areaslav,debug,clearance,
     &           cl2s,cl2m,shrink,reltime)
!
!     3. triangle
!
c            write(20,*)'*********tria 3**************'
            nmp=3
            nodem2(1)=nodem(5)
            nodem2(2)=nodem(3)
            nodem2(3)=nodem(6)
            do j2=1,3
               xl2m2(j2,1)=xl2m(j2,5)
            enddo
            do j2=1,3
               xl2m2(j2,2)=xl2m(j2,3)
            enddo
            do j2=1,3
               xl2m2(j2,3)=xl2m(j2,6)
            enddo
            call treatmasterface_mortar(
     &           nopes,slavstraight,xn,xns,co,xl2s,xl2sp,
     &           ipe,ime,iactiveline,nactiveline,
     &           ifacem,
     &           nintpoint,pslavsurf,imastsurf,pmastsurf,
     &           xl2m,nnodelem,xl2m2,nmp,nodem2,mi,
     &           gapmints,l,areaslav,debug,clearance,
     &           cl2s,cl2m,shrink,reltime)
!
!     4. triangle
!
c            write(20,*)'*********tria 4**************'
            nmp=3
            nodem2(1)=nodem(4)
            nodem2(2)=nodem(5)
            nodem2(3)=nodem(6)
            do j2=1,3
               xl2m2(j2,1)=xl2m(j2,4)
            enddo
            do j2=1,3
               xl2m2(j2,2)=xl2m(j2,5)
            enddo
            do j2=1,3
               xl2m2(j2,3)=xl2m(j2,6)
            enddo
            call treatmasterface_mortar(
     &           nopes,slavstraight,xn,xns,co,xl2s,xl2sp,
     &           ipe,ime,iactiveline,nactiveline,
     &           ifacem,
     &           nintpoint,pslavsurf,imastsurf,pmastsurf,
     &           xl2m,nnodelem,xl2m2,nmp,nodem2,mi,
     &           gapmints,l,areaslav,debug,clearance,
     &           cl2s,cl2m,shrink,reltime)
!
         elseif(nnodelem.eq.8)then
!
!     quad8 surface is divided into 4 tri3 + 1 quad4 surfaces
!
!     1. triangle
            nmp=3
            nodem2(1)=nodem(1)
            nodem2(2)=nodem(5)
            nodem2(3)=nodem(8)
            do j2=1,3
               xl2m2(j2,1)=xl2m(j2,1)
            enddo
            do j2=1,3
               xl2m2(j2,2)=xl2m(j2,5)
            enddo
            do j2=1,3
               xl2m2(j2,3)=xl2m(j2,8)
            enddo
            call treatmasterface_mortar(
     &           nopes,slavstraight,xn,xns,co,xl2s,xl2sp,
     &           ipe,ime,iactiveline,nactiveline,
     &           ifacem,
     &           nintpoint,pslavsurf,imastsurf,pmastsurf,
     &           xl2m,nnodelem,xl2m2,nmp,nodem2,mi,
     &           gapmints,l,areaslav,debug,clearance,
     &           cl2s,cl2m,shrink,reltime)
!
!     2. triangle
!
            nmp=3
            nodem2(1)=nodem(5)
            nodem2(2)=nodem(2)
            nodem2(3)=nodem(6)
            do j2=1,3
               xl2m2(j2,1)=xl2m(j2,5)
            enddo
            do j2=1,3
               xl2m2(j2,2)=xl2m(j2,2)
            enddo
            do j2=1,3
               xl2m2(j2,3)=xl2m(j2,6)
            enddo
            call treatmasterface_mortar(
     &           nopes,slavstraight,xn,xns,co,xl2s,xl2sp,
     &           ipe,ime,iactiveline,nactiveline,
     &           ifacem,
     &           nintpoint,pslavsurf,imastsurf,pmastsurf,
     &           xl2m,nnodelem,xl2m2,nmp,nodem2,mi,
     &           gapmints,l,areaslav,debug,clearance,
     &           cl2s,cl2m,shrink,reltime)
!
!     3. triangle
!
            nmp=3
            nodem2(1)=nodem(6)
            nodem2(2)=nodem(3)
            nodem2(3)=nodem(7)
            do j2=1,3
               xl2m2(j2,1)=xl2m(j2,6)
            enddo
            do j2=1,3
               xl2m2(j2,2)=xl2m(j2,3)
            enddo
            do j2=1,3
               xl2m2(j2,3)=xl2m(j2,7)
            enddo
            call treatmasterface_mortar(
     &           nopes,slavstraight,xn,xns,co,xl2s,xl2sp,
     &           ipe,ime,iactiveline,nactiveline,
     &           ifacem,
     &           nintpoint,pslavsurf,imastsurf,pmastsurf,
     &           xl2m,nnodelem,xl2m2,nmp,nodem2,mi,
     &           gapmints,l,areaslav,debug,clearance,
     &           cl2s,cl2m,shrink,reltime)
!
!     4. triangle
!
            nmp=3
            nodem2(1)=nodem(7)
            nodem2(2)=nodem(4)
            nodem2(3)=nodem(8)
            do j2=1,3
               xl2m2(j2,1)=xl2m(j2,7)
            enddo
            do j2=1,3
               xl2m2(j2,2)=xl2m(j2,4)
            enddo
            do j2=1,3
               xl2m2(j2,3)=xl2m(j2,8)
            enddo
            call treatmasterface_mortar(
     &           nopes,slavstraight,xn,xns,co,xl2s,xl2sp,
     &           ipe,ime,iactiveline,nactiveline,
     &           ifacem,
     &           nintpoint,pslavsurf,imastsurf,pmastsurf,
     &           xl2m,nnodelem,xl2m2,nmp,nodem2,mi,
     &           gapmints,l,areaslav,debug,clearance,
     &           cl2s,cl2m,shrink,reltime)
!
!     quad
!
            nmp=4
            nodem2(1)=nodem(5)
            nodem2(2)=nodem(6)
            nodem2(3)=nodem(7)
            nodem2(4)=nodem(8)
            do j2=1,3
               xl2m2(j2,1)=xl2m(j2,5)
            enddo
            do j2=1,3
               xl2m2(j2,2)=xl2m(j2,6)
            enddo
            do j2=1,3
              xl2m2(j2,3)=xl2m(j2,7)
           enddo
           do j2=1,3
              xl2m2(j2,4)=xl2m(j2,8)
           enddo
           call treatmasterface_mortar(
     &          nopes,slavstraight,xn,xns,co,xl2s,xl2sp,
     &          ipe,ime,iactiveline,nactiveline,
     &          ifacem,
     &          nintpoint,pslavsurf,imastsurf,pmastsurf,
     &          xl2m,nnodelem,xl2m2,nmp,nodem2,mi,
     &          gapmints,l,areaslav,debug,clearance,
     &           cl2s,cl2m,shrink,reltime)
        endif
      enddo
!
!     corners of the Slave surface have already been treated
!
      do j=1,nopes
         imfacecorner(j,1)=0
      enddo
!
!     retrieving all triangles by neighborhood search
!
      do
         line=iactiveline(1,1)
         if(nactiveline.eq.0) exit
         if(koncont(4,ime(2,line)).eq.iactiveline(2,1)) then
            itri=imastop(ime(3,line),ime(2,line))
         else
            itri=ime(2,line)
         endif
!     check whether still in contact tie
         if(itri.gt.itietri(2,i) .or. itri.lt.itietri(1,i))then
c            if(itri.ne.0)then
c     write(*,*)'sim:tie',i,'face',l,'itiri',itri
c     write(*,*)' mintri',itietri(1,i),'maxtri',itietri(2,i)
c            endif
            itri=0
         endif
!
         if(itri.eq.0) then
            nactiveline=nactiveline-1
            do il=1,nactiveline
               do k=1,3
                  iactiveline(k,il)=iactiveline(k,il+1)
               enddo
            enddo
            cycle
         endif
!
         ifacem=koncont(4,itri)
         nelemm=int(koncont(4,itri)/10.d0)
         jfacem=koncont(4,itri)-10*nelemm
!
!     add master element to covered stack
!
         call nident(icoveredmelem,nelemm, ncoveredmelem,id)
         if(id .gt. 0 .and. icoveredmelem(id).eq.nelemm)then
!     master elment was already threated
            nactiveline=nactiveline-1
            do il=1,nactiveline
               do k=1,3
                  iactiveline(k,il)=iactiveline(k,il+1)
               enddo
            enddo
            cycle
         else
!     add master element to covered elements
            ncoveredmelem=ncoveredmelem+1
            do ii=ncoveredmelem,id+2,-1
               icoveredmelem(ii)=icoveredmelem(ii-1)
            enddo
            icoveredmelem(id+1)=nelemm
         endif
         indexe=ipkon(nelemm)
         call getnumberofnodes(nelemm,jfacem,lakon,nopemm,
     &        nnodelem,idummy)

!
!     determining the nodes of the face
!
         do j1=1,nopemm
            konl(j1)=kon(ipkon(nelemm)+j1)
         enddo
         do k1=1,nnodelem
            ifac=getiface(k1,jfacem,nopemm)
            nodem(k1)=konl(ifac)
            do j1=1,3
               xl2m(j1,k1)=co(j1,konl(ifac))+
     &              vold(j1,konl(ifac))
               cl2m(j1,k1)=co(j1,konl(ifac))
            enddo
         enddo
         if(debug)then
            write(20,*) 'dd',dd ,'ifacem', ifacem
            write(20,100)(xn(k),k=1,3)
            write(20,*) 'SIM: xl2s'
            do j1=1,nopes
               write(20,*)(xl2s(k,j1),k=1,3)
            enddo
            write(20,*) 'SIM: xl2m'
            do j1=1,nnodelem
               write(20,*)(xl2m(k,j1),k=1,3)
            enddo
         endif
!
!     divide master element into konvex subelements
!
         if(nnodelem.eq.3 .or.nnodelem.eq.4)then
!
!     no loop needed
!
            nmp=nnodelem
            do k2=1,nnodelem
               nodem2(k2)=nodem(k2)
               do j2=1,3
                  xl2m2(j2,k2)=xl2m(j2,k2)
               enddo
            enddo
            call treatmasterface_mortar(
     &           nopes,slavstraight,xn,xns,co,xl2s,xl2sp,
     &           ipe,ime,iactiveline,nactiveline,
     &           ifacem,
     &           nintpoint,pslavsurf,imastsurf,pmastsurf,
     &           xl2m,nnodelem,xl2m2,nmp,nodem2,mi,
     &           gapmints,l,areaslav,debug,clearance,
     &           cl2s,cl2m,shrink,reltime)
!
         elseif(nnodelem.eq.6)then
!
!     tri6 surface is divided into 4 tri3 surfaces
!
!     1. triangle
c            write(20,*)'*********tria 1**************'
            nmp=3
            nodem2(1)=nodem(1)
            nodem2(2)=nodem(4)
            nodem2(3)=nodem(6)
            do j2=1,3
               xl2m2(j2,1)=xl2m(j2,1)
            enddo
            do j2=1,3
               xl2m2(j2,2)=xl2m(j2,4)
            enddo
            do j2=1,3
               xl2m2(j2,3)=xl2m(j2,6)
            enddo
            call treatmasterface_mortar(
     &           nopes,slavstraight,xn,xns,co,xl2s,xl2sp,
     &           ipe,ime,iactiveline,nactiveline,
     &           ifacem,
     &           nintpoint,pslavsurf,imastsurf,pmastsurf,
     &           xl2m,nnodelem,xl2m2,nmp,nodem2,mi,
     &           gapmints,l,areaslav,debug,clearance,
     &           cl2s,cl2m,shrink,reltime)
!
!     2. triangle
!
c            write(20,*)'*********tria 2**************'
            nmp=3
            nodem2(1)=nodem(4)
            nodem2(2)=nodem(2)
            nodem2(3)=nodem(5)
            do j2=1,3
               xl2m2(j2,1)=xl2m(j2,4)
            enddo
            do j2=1,3
               xl2m2(j2,2)=xl2m(j2,2)
            enddo
            do j2=1,3
               xl2m2(j2,3)=xl2m(j2,5)
            enddo
            call treatmasterface_mortar(
     &           nopes,slavstraight,xn,xns,co,xl2s,xl2sp,
     &           ipe,ime,iactiveline,nactiveline,
     &           ifacem,
     &           nintpoint,pslavsurf,imastsurf,pmastsurf,
     &           xl2m,nnodelem,xl2m2,nmp,nodem2,mi,
     &           gapmints,l,areaslav,debug,clearance,
     &           cl2s,cl2m,shrink,reltime)
!
!     3. triangle
!
c            write(20,*)'*********tria 3**************'
            nmp=3
            nodem2(1)=nodem(5)
            nodem2(2)=nodem(3)
            nodem2(3)=nodem(6)
            do j2=1,3
               xl2m2(j2,1)=xl2m(j2,5)
            enddo
            do j2=1,3
               xl2m2(j2,2)=xl2m(j2,3)
            enddo
            do j2=1,3
               xl2m2(j2,3)=xl2m(j2,6)
            enddo
            call treatmasterface_mortar(
     &           nopes,slavstraight,xn,xns,co,xl2s,xl2sp,
     &           ipe,ime,iactiveline,nactiveline,
     &           ifacem,
     &           nintpoint,pslavsurf,imastsurf,pmastsurf,
     &           xl2m,nnodelem,xl2m2,nmp,nodem2,mi,
     &           gapmints,l,areaslav,debug,clearance,
     &           cl2s,cl2m,shrink,reltime)
!
!     4. triangle
!
c            write(20,*)'*********tria 4**************'
            nmp=3
            nodem2(1)=nodem(4)
            nodem2(2)=nodem(5)
            nodem2(3)=nodem(6)
            do j2=1,3
               xl2m2(j2,1)=xl2m(j2,4)
            enddo
            do j2=1,3
               xl2m2(j2,2)=xl2m(j2,5)
            enddo
            do j2=1,3
               xl2m2(j2,3)=xl2m(j2,6)
            enddo
            call treatmasterface_mortar(
     &           nopes,slavstraight,xn,xns,co,xl2s,xl2sp,
     &           ipe,ime,iactiveline,nactiveline,
     &           ifacem,
     &           nintpoint,pslavsurf,imastsurf,pmastsurf,
     &           xl2m,nnodelem,xl2m2,nmp,nodem2,mi,
     &           gapmints,l,areaslav,debug,clearance,
     &           cl2s,cl2m,shrink,reltime)
!
         elseif(nnodelem.eq.8)then
!
!     quad8 surface is divided into 4 tri3 + 1 quad4 surfaces
!
!     1. triangle
            nmp=3
            nodem2(1)=nodem(1)
            nodem2(2)=nodem(5)
            nodem2(3)=nodem(8)
            do j2=1,3
               xl2m2(j2,1)=xl2m(j2,1)
            enddo
            do j2=1,3
               xl2m2(j2,2)=xl2m(j2,5)
            enddo
            do j2=1,3
               xl2m2(j2,3)=xl2m(j2,8)
            enddo
            call treatmasterface_mortar(
     &           nopes,slavstraight,xn,xns,co,xl2s,xl2sp,
     &           ipe,ime,iactiveline,nactiveline,
     &           ifacem,
     &           nintpoint,pslavsurf,imastsurf,pmastsurf,
     &           xl2m,nnodelem,xl2m2,nmp,nodem2,mi,
     &           gapmints,l,areaslav,debug,clearance,
     &           cl2s,cl2m,shrink,reltime)
!
!     2. triangle
!
            nmp=3
            nodem2(1)=nodem(5)
            nodem2(2)=nodem(2)
            nodem2(3)=nodem(6)
            do j2=1,3
               xl2m2(j2,1)=xl2m(j2,5)
            enddo
            do j2=1,3
               xl2m2(j2,2)=xl2m(j2,2)
            enddo
            do j2=1,3
               xl2m2(j2,3)=xl2m(j2,6)
            enddo
            call treatmasterface_mortar(
     &           nopes,slavstraight,xn,xns,co,xl2s,xl2sp,
     &           ipe,ime,iactiveline,nactiveline,
     &           ifacem,
     &           nintpoint,pslavsurf,imastsurf,pmastsurf,
     &           xl2m,nnodelem,xl2m2,nmp,nodem2,mi,
     &           gapmints,l,areaslav,debug,clearance,
     &           cl2s,cl2m,shrink,reltime)
!
!     3. triangle
!
            nmp=3
            nodem2(1)=nodem(6)
            nodem2(2)=nodem(3)
            nodem2(3)=nodem(7)
            do j2=1,3
               xl2m2(j2,1)=xl2m(j2,6)
            enddo
            do j2=1,3
               xl2m2(j2,2)=xl2m(j2,3)
            enddo
            do j2=1,3
               xl2m2(j2,3)=xl2m(j2,7)
            enddo
            call treatmasterface_mortar(
     &           nopes,slavstraight,xn,xns,co,xl2s,xl2sp,
     &           ipe,ime,iactiveline,nactiveline,
     &           ifacem,
     &           nintpoint,pslavsurf,imastsurf,pmastsurf,
     &           xl2m,nnodelem,xl2m2,nmp,nodem2,mi,
     &           gapmints,l,areaslav,debug,clearance,
     &           cl2s,cl2m,shrink,reltime)
!
!     4. triangle
!
            nmp=3
            nodem2(1)=nodem(7)
            nodem2(2)=nodem(4)
            nodem2(3)=nodem(8)
            do j2=1,3
               xl2m2(j2,1)=xl2m(j2,7)
            enddo
            do j2=1,3
               xl2m2(j2,2)=xl2m(j2,4)
            enddo
            do j2=1,3
               xl2m2(j2,3)=xl2m(j2,8)
            enddo
            call treatmasterface_mortar(
     &           nopes,slavstraight,xn,xns,co,xl2s,xl2sp,
     &           ipe,ime,iactiveline,nactiveline,
     &           ifacem,
     &           nintpoint,pslavsurf,imastsurf,pmastsurf,
     &           xl2m,nnodelem,xl2m2,nmp,nodem2,mi,
     &           gapmints,l,areaslav,debug,clearance,
     &           cl2s,cl2m,shrink,reltime)
!
!     quad
!
            nmp=4
            nodem2(1)=nodem(5)
            nodem2(2)=nodem(6)
            nodem2(3)=nodem(7)
            nodem2(4)=nodem(8)
            do j2=1,3
               xl2m2(j2,1)=xl2m(j2,5)
            enddo
            do j2=1,3
               xl2m2(j2,2)=xl2m(j2,6)
            enddo
            do j2=1,3
              xl2m2(j2,3)=xl2m(j2,7)
           enddo
           do j2=1,3
              xl2m2(j2,4)=xl2m(j2,8)
           enddo
           call treatmasterface_mortar(
     &          nopes,slavstraight,xn,xns,co,xl2s,xl2sp,
     &          ipe,ime,iactiveline,nactiveline,
     &          ifacem,
     &          nintpoint,pslavsurf,imastsurf,pmastsurf,
     &          xl2m,nnodelem,xl2m2,nmp,nodem2,mi,
     &          gapmints,l,areaslav,debug,clearance,
     &           cl2s,cl2m,shrink,reltime)
        endif
      enddo
      islavsurf(2,l+1)=nintpoint
      if(debug) then
         if(areaslav.lt.1.e-12)write(*,*)'areaslav(',l,')=',
     &        areaslav
      endif
!
      if(debug)then
         write(20,*) 'mint2d', (nintpoint-islavsurf(2,l))
         write(20,*)'areaslav(',l,')=',areaslav
         write(20,*)'*********************************'
      endif
!
      return
      end

max_stars_repo_name: alleindrach/calculix-desktop

max_issues_repo_name: alleindrach/calculix-desktop

max_forks_repo_name: alleindrach/calculix-desktop


MODULE module_sf_noahmpdrv

  USE module_sf_noahmplsm
  USE module_sf_urban
  USE module_sf_noahdrv, ONLY : SOIL_VEG_GEN_PARM
  USE module_sf_noah_seaice
  USE module_sf_noahmp_glacier
  USE MODULE_RA_GFDLETA, ONLY: CAL_MON_DAY

  USE module_data_gocart_dust

CONTAINS

  SUBROUTINE noahmplsm(ITIMESTEP, YR, JULIAN, COSZIN, XLATIN, &
                  DZ8W, DT, DZS, NSOIL, DX, &
         IVGTYP, ISLTYP, VEGFRA, VEGMAX, TMN, &
   XLAND, XICE,XICE_THRES, ISICE, ISURBAN, &
                 IDVEG, IOPT_CRS, IOPT_BTR, IOPT_RUN, IOPT_SFC, IOPT_FRZ, &
              IOPT_INF, IOPT_RAD, IOPT_ALB, IOPT_SNF,IOPT_TBOT, IOPT_STC, &
               IZ0TLND, &
                   T3D, QV3D, U_PHY, V_PHY, SWDOWN, GLW, &
                 P8W3D, RAINBL, &
                   TSK, HFX, QFX, LH, GRDFLX, SMSTAV, &
                SMSTOT,SFCRUNOFF, UDRUNOFF, ALBEDO, SNOWC, SMOIS, &
    SH2O, TSLB, SNOW, SNOWH, CANWAT, ACSNOM, &
  ACSNOW, EMISS, QSFC, &
               ISNOWXY, TVXY, TGXY, CANICEXY, CANLIQXY, EAHXY, &
          TAHXY, CMXY, CHXY, FWETXY, SNEQVOXY, ALBOLDXY, &
               QSNOWXY, WSLAKEXY, ZWTXY, WAXY, WTXY, TSNOXY, &
        ZSNSOXY, SNICEXY, SNLIQXY, LFMASSXY, RTMASSXY, STMASSXY, &
         WOODXY, STBLCPXY, FASTCPXY, XLAIXY, XSAIXY, TAUSSXY, &
         T2MVXY, T2MBXY, Q2MVXY, Q2MBXY, &
         TRADXY, NEEXY, GPPXY, NPPXY, FVEGXY, RUNSFXY, &
        RUNSBXY, ECANXY, EDIRXY, ETRANXY, FSAXY, FIRAXY, &
         APARXY, PSNXY, SAVXY, SAGXY, RSSUNXY, RSSHAXY, &
  BGAPXY, WGAPXY, TGVXY, TGBXY, CHVXY, CHBXY, &
   SHGXY, SHCXY, SHBXY, EVGXY, EVBXY, GHVXY, &
   GHBXY, IRGXY, IRCXY, IRBXY, TRXY, EVCXY, &
              CHLEAFXY, CHUCXY, CHV2XY, CHB2XY, &
               ids,ide, jds,jde, kds,kde, &
               ims,ime, jms,jme, kms,kme, &
               its,ite, jts,jte, kts,kte )

    IMPLICIT NONE

    INTEGER, INTENT(IN ) :: ITIMESTEP
    INTEGER, INTENT(IN ) :: YR
    REAL, INTENT(IN ) :: JULIAN
    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: COSZIN
    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: XLATIN
    REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) :: DZ8W
    REAL, INTENT(IN ) :: DT
    REAL, DIMENSION(1:nsoil), INTENT(IN ) :: DZS
    INTEGER, INTENT(IN ) :: NSOIL
    REAL, INTENT(IN ) :: DX
    INTEGER, DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: IVGTYP
    INTEGER, DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: ISLTYP
    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: VEGFRA
    REAL, DIMENSION( ims:ime , jms:jme ), INTENT(IN ) :: VEGMAX
    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: TMN
    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: XLAND
    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: XICE
    REAL, INTENT(IN ) :: XICE_THRES
    INTEGER, INTENT(IN ) :: ISICE
    INTEGER, INTENT(IN ) :: ISURBAN
    INTEGER, INTENT(IN ) :: IDVEG
    INTEGER, INTENT(IN ) :: IOPT_CRS
    INTEGER, INTENT(IN ) :: IOPT_BTR
    INTEGER, INTENT(IN ) :: IOPT_RUN
    INTEGER, INTENT(IN ) :: IOPT_SFC
    INTEGER, INTENT(IN ) :: IOPT_FRZ
    INTEGER, INTENT(IN ) :: IOPT_INF
    INTEGER, INTENT(IN ) :: IOPT_RAD
    INTEGER, INTENT(IN ) :: IOPT_ALB
    INTEGER, INTENT(IN ) :: IOPT_SNF
    INTEGER, INTENT(IN ) :: IOPT_TBOT
    INTEGER, INTENT(IN ) :: IOPT_STC
    INTEGER, INTENT(IN ) :: IZ0TLND
    REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) :: T3D
    REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) :: QV3D
    REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) :: U_PHY
    REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) :: V_PHY
    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: SWDOWN
    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: GLW
    REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) :: P8W3D
    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: RAINBL

    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: TSK
    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: HFX
    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: QFX
    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: LH
    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: GRDFLX
    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: SMSTAV
    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: SMSTOT
    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: SFCRUNOFF
    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: UDRUNOFF
    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: ALBEDO
    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: SNOWC
    REAL, DIMENSION( ims:ime, 1:nsoil, jms:jme ), INTENT(INOUT) :: SMOIS
    REAL, DIMENSION( ims:ime, 1:nsoil, jms:jme ), INTENT(INOUT) :: SH2O
    REAL, DIMENSION( ims:ime, 1:nsoil, jms:jme ), INTENT(INOUT) :: TSLB
    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: SNOW
    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: SNOWH
    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: CANWAT
    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: ACSNOM
    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: ACSNOW
    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: EMISS
    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: QSFC

    INTEGER, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: ISNOWXY
    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: TVXY
    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: TGXY
    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: CANICEXY
    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: CANLIQXY
    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: EAHXY
    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: TAHXY
    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: CMXY
    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: CHXY
    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: FWETXY
    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: SNEQVOXY
    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: ALBOLDXY
    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: QSNOWXY
    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: WSLAKEXY
    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: ZWTXY
    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: WAXY
    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: WTXY
    REAL, DIMENSION( ims:ime,-2:0, jms:jme ), INTENT(INOUT) :: TSNOXY
    REAL, DIMENSION( ims:ime,-2:NSOIL, jms:jme ), INTENT(INOUT) :: ZSNSOXY
    REAL, DIMENSION( ims:ime,-2:0, jms:jme ), INTENT(INOUT) :: SNICEXY
    REAL, DIMENSION( ims:ime,-2:0, jms:jme ), INTENT(INOUT) :: SNLIQXY
    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: LFMASSXY
    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: RTMASSXY
    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: STMASSXY
    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: WOODXY
    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: STBLCPXY
    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: FASTCPXY
    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: XLAIXY
    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: XSAIXY
    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: TAUSSXY

    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT ) :: T2MVXY
    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT ) :: T2MBXY
    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT ) :: Q2MVXY
    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT ) :: Q2MBXY
    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT ) :: TRADXY
    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT ) :: NEEXY
    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT ) :: GPPXY
    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT ) :: NPPXY
    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT ) :: FVEGXY
    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT ) :: RUNSFXY
    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT ) :: RUNSBXY
    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT ) :: ECANXY
    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT ) :: EDIRXY
    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT ) :: ETRANXY
    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT ) :: FSAXY
    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT ) :: FIRAXY
    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT ) :: APARXY
    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT ) :: PSNXY
    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT ) :: SAVXY
    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT ) :: SAGXY
    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT ) :: RSSUNXY
    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT ) :: RSSHAXY
    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT ) :: BGAPXY
    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT ) :: WGAPXY
    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT ) :: TGVXY
    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT ) :: TGBXY
    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT ) :: CHVXY
    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT ) :: CHBXY
    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT ) :: SHGXY
    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT ) :: SHCXY
    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT ) :: SHBXY
    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT ) :: EVGXY
    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT ) :: EVBXY
    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT ) :: GHVXY
    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT ) :: GHBXY
    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT ) :: IRGXY
    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT ) :: IRCXY
    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT ) :: IRBXY
    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT ) :: TRXY
    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT ) :: EVCXY
    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT ) :: CHLEAFXY
    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT ) :: CHUCXY
    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT ) :: CHV2XY
    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT ) :: CHB2XY
    INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
         & ims,ime, jms,jme, kms,kme, &
         & its,ite, jts,jte, kts,kte

    REAL :: COSZ
    REAL :: LAT
    REAL :: Z_ML
    INTEGER :: VEGTYP
    INTEGER :: SOILTYP
    REAL :: FVEG
    REAL :: FVGMAX
    REAL :: TBOT
    REAL :: T_ML
    REAL :: Q_ML
    REAL :: U_ML
    REAL :: V_ML
    REAL :: SWDN
    REAL :: LWDN
    REAL :: P_ML
    REAL :: PSFC
    REAL :: PRCP

    REAL :: FSH
    REAL :: SSOIL
    REAL :: SALB
    REAL :: FSNO
    REAL, DIMENSION( 1:NSOIL) :: SMC
    REAL, DIMENSION( 1:NSOIL) :: SMH2O
    REAL, DIMENSION(-2:NSOIL) :: STC
    REAL :: SWE
    REAL :: SNDPTH
    REAL :: EMISSI
    REAL :: QSFC1D

    INTEGER :: ISNOW
    REAL :: TV
    REAL :: TG
    REAL :: CANICE
    REAL :: CANLIQ
    REAL :: EAH
    REAL :: TAH
    REAL :: CM
    REAL :: CH
    REAL :: FWET
    REAL :: SNEQVO
    REAL :: ALBOLD
    REAL :: QSNOW
    REAL :: WSLAKE
    REAL :: ZWT
    REAL :: WA
    REAL :: WT
    REAL, DIMENSION(-2:NSOIL) :: ZSNSO
    REAL, DIMENSION(-2: 0) :: SNICE
    REAL, DIMENSION(-2: 0) :: SNLIQ
    REAL :: LFMASS
    REAL :: RTMASS
    REAL :: STMASS
    REAL :: WOOD
    REAL :: STBLCP
    REAL :: FASTCP
    REAL :: PLAI
    REAL :: PSAI
    REAL :: TAUSS

    REAL :: T2MV
    REAL :: T2MB
    REAL :: Q2MV
    REAL :: Q2MB
    REAL :: TRAD
    REAL :: NEE
    REAL :: GPP
    REAL :: NPP
    REAL :: FVEGMP
    REAL :: RUNSF
    REAL :: RUNSB
    REAL :: ECAN
    REAL :: ETRAN
    REAL :: ESOIL
    REAL :: FSA
    REAL :: FIRA
    REAL :: APAR
    REAL :: PSN
    REAL :: SAV
    REAL :: SAG
    REAL :: RSSUN
    REAL :: RSSHA
    REAL :: BGAP
    REAL :: WGAP
    REAL :: TGV
    REAL :: TGB
    REAL :: CHV
    REAL :: CHB
    REAL :: IRC
    REAL :: IRG
    REAL :: SHC
    REAL :: SHG
    REAL :: EVG
    REAL :: GHV
    REAL :: IRB
    REAL :: SHB
    REAL :: EVB
    REAL :: GHB
    REAL :: TR
    REAL :: EVC
    REAL :: CHLEAF
    REAL :: CHUC
    REAL :: CHV2
    REAL :: CHB2

    REAL :: FPICE
    REAL :: FCEV
    REAL :: FGEV
    REAL :: FCTR
    REAL :: QSNBOT
    REAL :: PONDING
    REAL :: PONDING1
    REAL :: PONDING2

    REAL :: FSR
    REAL, DIMENSION(-2:0) :: FICEOLD
    REAL :: CO2PP
    REAL :: O2PP
    REAL, DIMENSION(1:NSOIL) :: ZSOIL
    REAL :: FOLN

    REAL :: QC
    REAL :: PBLH
    REAL :: DZ8W1D

    INTEGER :: I
    INTEGER :: J
    INTEGER :: K
    INTEGER :: ICE
    INTEGER :: SLOPETYP
    LOGICAL :: IPRINT

    INTEGER :: ISC
    INTEGER :: IST
    INTEGER :: YEARLEN

    INTEGER, PARAMETER :: NSNOW = 3
    REAL, PARAMETER :: CO2 = 395.e-06
    REAL, PARAMETER :: O2 = 0.209
    REAL, PARAMETER :: undefined_value = -1.E36

    CALL NOAHMP_OPTIONS(IDVEG ,IOPT_CRS ,IOPT_BTR ,IOPT_RUN ,IOPT_SFC ,IOPT_FRZ , &
                     IOPT_INF ,IOPT_RAD ,IOPT_ALB ,IOPT_SNF ,IOPT_TBOT, IOPT_STC )

    IPRINT = .false.

    YEARLEN = 365
    if (mod(YR,4) == 0) then
       YEARLEN = 366
       if (mod(YR,100) == 0) then
          YEARLEN = 365
          if (mod(YR,400) == 0) then
             YEARLEN = 366
          endif
       endif
    endif

    ZSOIL(1) = -DZS(1)
    DO K = 2, NSOIL
       ZSOIL(K) = -DZS(K) + ZSOIL(K-1)
    END DO

    JLOOP : DO J=jts,jte

       IF(ITIMESTEP == 1)THEN
          DO I=its,ite
             IF((XLAND(I,J)-1.5) >= 0.) THEN
                IF(XICE(I,J) == 1. .AND. IPRINT) PRINT *,' sea-ice at water point, I=',I,'J=',J
                SMSTAV(I,J) = 1.0
                SMSTOT(I,J) = 1.0
                DO K = 1, NSOIL
                   SMOIS(I,K,J) = 1.0
                    TSLB(I,K,J) = 273.16
                ENDDO
             ELSE
                IF(XICE(I,J) == 1.) THEN
                   SMSTAV(I,J) = 1.0
                   SMSTOT(I,J) = 1.0
                   DO K = 1, NSOIL
                      SMOIS(I,K,J) = 1.0
                   ENDDO
                ENDIF
             ENDIF
          ENDDO
       ENDIF

   ILOOP : DO I = its, ite

    IF (XICE(I,J) >= XICE_THRES) THEN
       ICE = 1
    ELSE IF ( IVGTYP(I,J) == ISICE ) THEN
       ICE = -1
    ELSE
       ICE=0
    ENDIF

    IF((XLAND(I,J)-1.5) >= 0.) CYCLE ILOOP

    IF ( ICE == 1) THEN

       SH2O (i,1:NSOIL,j) = 1.0
       XLAIXY(i,j) = 0.01

       CYCLE ILOOP

    ELSE

       COSZ = COSZIN (I,J)
       LAT = XLATIN (I,J)
       Z_ML = 0.5*DZ8W(I,1,J)
       VEGTYP = IVGTYP(I,J)
       SOILTYP= ISLTYP(I,J)
       FVEG = VEGFRA(I,J)/100.
       FVGMAX = VEGMAX (I,J)/100.
       TBOT = TMN(I,J)
       T_ML = T3D(I,1,J)
       Q_ML = QV3D(I,1,J)/(1.0+QV3D(I,1,J))
       U_ML = U_PHY(I,1,J)
       V_ML = V_PHY(I,1,J)
       SWDN = SWDOWN(I,J)
       LWDN = GLW(I,J)
       P_ML =(P8W3D(I,KTS+1,J)+P8W3D(I,KTS,J))*0.5

       PSFC = P8W3D(I,1,J)
       PRCP = RAINBL(I,J)/DT

       ISNOW = ISNOWXY (I,J)
       SMC ( 1:NSOIL) = SMOIS (I, 1:NSOIL,J)
       SMH2O( 1:NSOIL) = SH2O (I, 1:NSOIL,J)
       STC (-NSNOW+1: 0) = TSNOXY (I,-NSNOW+1: 0,J)
       STC ( 1:NSOIL) = TSLB (I, 1:NSOIL,J)
       SWE = SNOW (I,J)
       SNDPTH = SNOWH (I,J)
       QSFC1D = QSFC (I,J)

       TV = TVXY (I,J)
       TG = TGXY (I,J)
       CANLIQ = CANLIQXY(I,J)
       CANICE = CANICEXY(I,J)
       EAH = EAHXY (I,J)
       TAH = TAHXY (I,J)
       CM = CMXY (I,J)
       CH = CHXY (I,J)
       FWET = FWETXY (I,J)
       SNEQVO = SNEQVOXY(I,J)
       ALBOLD = ALBOLDXY(I,J)
       QSNOW = QSNOWXY (I,J)
       WSLAKE = WSLAKEXY(I,J)
       ZWT = ZWTXY (I,J)
       WA = WAXY (I,J)
       WT = WTXY (I,J)
       ZSNSO(-NSNOW+1:NSOIL) = ZSNSOXY (I,-NSNOW+1:NSOIL,J)
       SNICE(-NSNOW+1: 0) = SNICEXY (I,-NSNOW+1: 0,J)
       SNLIQ(-NSNOW+1: 0) = SNLIQXY (I,-NSNOW+1: 0,J)
       LFMASS = LFMASSXY(I,J)
       RTMASS = RTMASSXY(I,J)
       STMASS = STMASSXY(I,J)
       WOOD = WOODXY (I,J)
       STBLCP = STBLCPXY(I,J)
       FASTCP = FASTCPXY(I,J)
       PLAI = XLAIXY (I,J)
       PSAI = XSAIXY (I,J)
       TAUSS = TAUSSXY (I,J)

       FICEOLD = 0.0
       FICEOLD(ISNOW+1:0) = SNICEXY(I,ISNOW+1:0,J) &
           /(SNICEXY(I,ISNOW+1:0,J)+SNLIQXY(I,ISNOW+1:0,J))
       CO2PP = CO2 * P_ML
       O2PP = O2 * P_ML
       FOLN = 1.0
       QC = undefined_value
       PBLH = undefined_value
       DZ8W1D = DZ8W (I,1,J)
       SLOPETYP = 1
       IST = 1
       ISC = 4

       IF(SOILTYP == 14 .AND. XICE(I,J) == 0.) THEN
          IF(IPRINT) PRINT *, ' SOIL TYPE FOUND TO BE WATER AT A LAND-POINT'
          IF(IPRINT) PRINT *, i,j,'RESET SOIL in surfce.F'
          SOILTYP = 7
       ENDIF

       IF( IVGTYP(I,J) == ISURBAN .or. IVGTYP(I,J) == 31 .or. &
            IVGTYP(I,J) == 32 .or. IVGTYP(I,J) == 33) THEN
          VEGTYP = ISURBAN
       ENDIF
       IF(VEGTYP == 25) FVEG = 0.0
       IF(VEGTYP == 25) PLAI = 0.0
       IF(VEGTYP == 26) FVEG = 0.0
       IF(VEGTYP == 26) PLAI = 0.0
       IF(VEGTYP == 27) FVEG = 0.0
       IF(VEGTYP == 27) PLAI = 0.0

       CALL REDPRM (VEGTYP,SOILTYP,SLOPETYP,ZSOIL,NSOIL,ISURBAN)

    IF ( ICE == -1 ) THEN

      CALL NOAHMP_OPTIONS_GLACIER(IDVEG ,IOPT_CRS ,IOPT_BTR ,IOPT_RUN ,IOPT_SFC ,IOPT_FRZ , &
                      IOPT_INF ,IOPT_RAD ,IOPT_ALB ,IOPT_SNF ,IOPT_TBOT, IOPT_STC )

      TBOT = MIN(TBOT,263.15)
      CALL NOAHMP_GLACIER( I, J, COSZ, NSNOW, NSOIL, DT, &
                               T_ML, P_ML, U_ML, V_ML, Q_ML, SWDN, &
                               PRCP, LWDN, TBOT, Z_ML, FICEOLD, ZSOIL, &
                              QSNOW, SNEQVO, ALBOLD, CM, CH, ISNOW, &
                                SWE, SMC, ZSNSO, SNDPTH, SNICE, SNLIQ, &
                                TGB, STC, SMH2O, TAUSS, QSFC1D, &
                                FSA, FSR, FIRA, FSH, FGEV, SSOIL, &
                               TRAD, ESOIL, RUNSF, RUNSB, SAG, SALB, &
                              QSNBOT,PONDING,PONDING1,PONDING2, T2MB, Q2MB, &
         EMISSI, FPICE, CHB2 )

       FSNO = 1.0
       TV = undefined_value
       TG = TGB
       CANICE = undefined_value
       CANLIQ = undefined_value
       EAH = undefined_value
       TAH = undefined_value
       FWET = undefined_value
       WSLAKE = undefined_value
       ZWT = undefined_value
       WA = undefined_value
       WT = undefined_value
       LFMASS = undefined_value
       RTMASS = undefined_value
       STMASS = undefined_value
       WOOD = undefined_value
       STBLCP = undefined_value
       FASTCP = undefined_value
       PLAI = undefined_value
       PSAI = undefined_value
       T2MV = undefined_value
       Q2MV = undefined_value
       NEE = undefined_value
       GPP = undefined_value
       NPP = undefined_value
       FVEGMP = 0.0
       ECAN = undefined_value
       ETRAN = undefined_value
       APAR = undefined_value
       PSN = undefined_value
       SAV = undefined_value
       RSSUN = undefined_value
       RSSHA = undefined_value
       BGAP = undefined_value
       WGAP = undefined_value
       TGV = undefined_value
       CHV = undefined_value
       CHB = CH
       IRC = undefined_value
       IRG = undefined_value
       SHC = undefined_value
       SHG = undefined_value
       EVG = undefined_value
       GHV = undefined_value
       IRB = FIRA
       SHB = FSH
       EVB = FGEV
       GHB = SSOIL
       TR = undefined_value
       EVC = undefined_value
       CHLEAF = undefined_value
       CHUC = undefined_value
       CHV2 = undefined_value
       FCEV = undefined_value
       FCTR = undefined_value

       QFX(I,J) = ESOIL

    ELSE
    goto 1000
if(i==1.and.j==8) then
    print*,I , J , LAT , YEARLEN , JULIAN , COSZ
    print*,'DT'
    print*,DT , DX , DZ8W1D , NSOIL , ZSOIL , 3
    print*,'FVEG'
    print*,FVEG , FVGMAX , VEGTYP , ISURBAN , ICE , IST
    print*,ISC
    print*,IZ0TLND
    print*,'T_ML'
    print*,T_ML , P_ML , PSFC , U_ML , V_ML , Q_ML
    print*,'QC'
    print*,QC , SWDN , LWDN , PRCP , TBOT , CO2PP
    print*,'O2PP'
    print*,O2PP , FOLN , FICEOLD , PBLH , Z_ML
    print*,'ALBOLD'
    print*,ALBOLD , SNEQVO
    print*,'STC'
    print*,STC , SMH2O , SMC , TAH , EAH , FWET
    print*,'CANLIQ'
    print*,CANLIQ , CANICE , TV , TG , QSFC1D , QSNOW
    print*,'ISNOW'
    print*,ISNOW , ZSNSO , SNDPTH , SWE , SNICE , SNLIQ
    print*,'ZWT'
    print*,ZWT , WA , WT , WSLAKE , LFMASS , RTMASS
    print*,'STMASS'
    print*,STMASS , WOOD , STBLCP , FASTCP , PLAI , PSAI
    print*,'CM'
    print*,CM , CH , TAUSS
    print*,'FSA'
    print*,FSA , FSR , FIRA , FSH , SSOIL , FCEV
    print*,'FGEV'
    print*,FGEV , FCTR , ECAN , ETRAN , ESOIL , TRAD
    print*,'TGB'
    print*, TGB , TGV , T2MV , T2MB
    print*,'Q2MV'
    print*, Q2MV , Q2MB , RUNSF , RUNSB , APAR
    print*,'PSN'
    print*,PSN , SAV , SAG , FSNO , NEE , GPP
    print*,'NPP'
    print*,NPP , FVEGMP , SALB , QSNBOT , PONDING , PONDING1
    print*,'PONDING2'
    print*,PONDING2, RSSUN , RSSHA , BGAP , WGAP
    print*,'CHV'
    print*, CHV , CHB , EMISSI
end if

1000 continue

       CALL NOAHMP_SFLX (&
            I , J , LAT , YEARLEN , JULIAN , COSZ , &
            DT , DX , DZ8W1D , NSOIL , ZSOIL , NSNOW , &
            FVEG , FVGMAX , VEGTYP , ISURBAN , ICE , IST , &
            ISC , &
            IZ0TLND , &
            T_ML , P_ML , PSFC , U_ML , V_ML , Q_ML , &
            QC , SWDN , LWDN , PRCP , TBOT , CO2PP , &
            O2PP , FOLN , FICEOLD , PBLH , Z_ML , &
            ALBOLD , SNEQVO , &
            STC , SMH2O , SMC , TAH , EAH , FWET , &
            CANLIQ , CANICE , TV , TG , QSFC1D , QSNOW , &
            ISNOW , ZSNSO , SNDPTH , SWE , SNICE , SNLIQ , &
            ZWT , WA , WT , WSLAKE , LFMASS , RTMASS , &
            STMASS , WOOD , STBLCP , FASTCP , PLAI , PSAI , &
            CM , CH , TAUSS , &
            FSA , FSR , FIRA , FSH , SSOIL , FCEV , &
            FGEV , FCTR , ECAN , ETRAN , ESOIL , TRAD , &
            TGB , TGV , T2MV , T2MB , Q2MV , Q2MB , &
            RUNSF , RUNSB , APAR , PSN , SAV , SAG , &
            FSNO , NEE , GPP , NPP , FVEGMP , SALB , &
            QSNBOT , PONDING , PONDING1, PONDING2, RSSUN , RSSHA , &
            BGAP , WGAP , CHV , CHB , EMISSI , &
            SHG , SHC , SHB , EVG , EVB , GHV , &
     GHB , IRG , IRC , IRB , TR , EVC , &
     CHLEAF , CHUC , CHV2 , CHB2 , FPICE )

            QFX(I,J) = ECAN + ESOIL + ETRAN

   ENDIF

             TSK (I,J) = TRAD
             HFX (I,J) = FSH
             LH (I,J) = FCEV + FGEV + FCTR
             GRDFLX (I,J) = SSOIL
      SMSTAV (I,J) = 0.0
             SMSTOT (I,J) = 0.0
             SFCRUNOFF(I,J) = SFCRUNOFF(I,J) + RUNSF * DT
             UDRUNOFF (I,J) = UDRUNOFF(I,J) + RUNSB * DT
             IF ( SALB > -999 ) THEN
                ALBEDO(I,J) = SALB
             ENDIF
             SNOWC (I,J) = FSNO
             SMOIS (I, 1:NSOIL,J) = SMC ( 1:NSOIL)
             SH2O (I, 1:NSOIL,J) = SMH2O ( 1:NSOIL)
             TSLB (I, 1:NSOIL,J) = STC ( 1:NSOIL)
             SNOW (I,J) = SWE
             SNOWH (I,J) = SNDPTH
             CANWAT (I,J) = CANLIQ + CANICE
             ACSNOW (I,J) = ACSNOW(I,J) + PRCP * FPICE
             ACSNOM (I,J) = ACSNOM(I,J) + QSNBOT*DT + PONDING + PONDING1 + PONDING2
             EMISS (I,J) = EMISSI
             QSFC (I,J) = QSFC1D

             ISNOWXY (I,J) = ISNOW
             TVXY (I,J) = TV
             TGXY (I,J) = TG
             CANLIQXY (I,J) = CANLIQ
             CANICEXY (I,J) = CANICE
             EAHXY (I,J) = EAH
             TAHXY (I,J) = TAH
             CMXY (I,J) = CM
             CHXY (I,J) = CH
             FWETXY (I,J) = FWET
             SNEQVOXY (I,J) = SNEQVO
             ALBOLDXY (I,J) = ALBOLD
             QSNOWXY (I,J) = QSNOW
             WSLAKEXY (I,J) = WSLAKE
             ZWTXY (I,J) = ZWT
             WAXY (I,J) = WA
             WTXY (I,J) = WT
             TSNOXY (I,-NSNOW+1: 0,J) = STC (-NSNOW+1: 0)
             ZSNSOXY (I,-NSNOW+1:NSOIL,J) = ZSNSO (-NSNOW+1:NSOIL)
             SNICEXY (I,-NSNOW+1: 0,J) = SNICE (-NSNOW+1: 0)
             SNLIQXY (I,-NSNOW+1: 0,J) = SNLIQ (-NSNOW+1: 0)
             LFMASSXY (I,J) = LFMASS
             RTMASSXY (I,J) = RTMASS
             STMASSXY (I,J) = STMASS
             WOODXY (I,J) = WOOD
             STBLCPXY (I,J) = STBLCP
             FASTCPXY (I,J) = FASTCP
             XLAIXY (I,J) = PLAI
             XSAIXY (I,J) = PSAI
             TAUSSXY (I,J) = TAUSS

             T2MVXY (I,J) = T2MV
             T2MBXY (I,J) = T2MB
             Q2MVXY (I,J) = Q2MV/(1.0 - Q2MV)
             Q2MBXY (I,J) = Q2MB/(1.0 - Q2MB)
             TRADXY (I,J) = TRAD
             NEEXY (I,J) = NEE
             GPPXY (I,J) = GPP
             NPPXY (I,J) = NPP
             FVEGXY (I,J) = FVEGMP
             RUNSFXY (I,J) = RUNSF
             RUNSBXY (I,J) = RUNSB
             ECANXY (I,J) = ECAN
             EDIRXY (I,J) = ESOIL
             ETRANXY (I,J) = ETRAN
             FSAXY (I,J) = FSA
             FIRAXY (I,J) = FIRA
             APARXY (I,J) = APAR
             PSNXY (I,J) = PSN
             SAVXY (I,J) = SAV
             SAGXY (I,J) = SAG
             RSSUNXY (I,J) = RSSUN
             RSSHAXY (I,J) = RSSHA
             BGAPXY (I,J) = BGAP
             WGAPXY (I,J) = WGAP
             TGVXY (I,J) = TGV
             TGBXY (I,J) = TGB
             CHVXY (I,J) = CHV
             CHBXY (I,J) = CHB
             IRCXY (I,J) = IRC
             IRGXY (I,J) = IRG
             SHCXY (I,J) = SHC
             SHGXY (I,J) = SHG
             EVGXY (I,J) = EVG
             GHVXY (I,J) = GHV
             IRBXY (I,J) = IRB
             SHBXY (I,J) = SHB
             EVBXY (I,J) = EVB
             GHBXY (I,J) = GHB
             TRXY (I,J) = TR
             EVCXY (I,J) = EVC
             CHLEAFXY (I,J) = CHLEAF
             CHUCXY (I,J) = CHUC
             CHV2XY (I,J) = CHV2
             CHB2XY (I,J) = CHB2

          ENDIF

      ENDDO ILOOP
   ENDDO JLOOP

  END SUBROUTINE noahmplsm

  SUBROUTINE NOAHMP_INIT ( MMINLU, SNOW , SNOWH , CANWAT , ISLTYP , IVGTYP, &
       TSLB , SMOIS , SH2O , DZS , FNDSOILW , FNDSNOWH , ISICE, &
       TSK, isnowxy , tvxy ,tgxy ,canicexy , TMN, XICE, &
       canliqxy ,eahxy ,tahxy ,cmxy ,chxy , &
       fwetxy ,sneqvoxy ,alboldxy ,qsnowxy ,wslakexy ,zwtxy ,waxy , &
       wtxy ,tsnoxy ,zsnsoxy ,snicexy ,snliqxy ,lfmassxy ,rtmassxy , &
       stmassxy ,woodxy ,stblcpxy ,fastcpxy ,xsaixy , &

       t2mvxy ,t2mbxy ,chstarxy, &

       NSOIL, restart, &
       allowed_to_read , &
       ids,ide, jds,jde, kds,kde, &
       ims,ime, jms,jme, kms,kme, &
       its,ite, jts,jte, kts,kte )

    INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
         & ims,ime, jms,jme, kms,kme, &
         & its,ite, jts,jte, kts,kte

    INTEGER, INTENT(IN) :: NSOIL, ISICE

    LOGICAL, INTENT(IN) :: restart, &
         & allowed_to_read

    REAL, DIMENSION( NSOIL), INTENT(IN) :: DZS

    REAL, DIMENSION( ims:ime, NSOIL, jms:jme ) , &
         & INTENT(INOUT) :: SMOIS, &
         & SH2O, &
         & TSLB

    REAL, DIMENSION( ims:ime, jms:jme ) , &
         & INTENT(INOUT) :: SNOW, &
         & SNOWH, &
         & CANWAT

    INTEGER, DIMENSION( ims:ime, jms:jme ), &
         & INTENT(IN) :: ISLTYP, &
                                     IVGTYP

    LOGICAL, INTENT(IN) :: FNDSOILW, &
         & FNDSNOWH

    REAL, DIMENSION(ims:ime,jms:jme), INTENT(IN) :: TSK
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: TMN
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(IN) :: XICE
    INTEGER, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: isnowxy
    REAL, DIMENSION(ims:ime,-2:NSOIL,jms:jme), INTENT(INOUT) :: zsnsoxy
    REAL, DIMENSION(ims:ime,-2: 0,jms:jme), INTENT(INOUT) :: tsnoxy
    REAL, DIMENSION(ims:ime,-2: 0,jms:jme), INTENT(INOUT) :: snicexy
    REAL, DIMENSION(ims:ime,-2: 0,jms:jme), INTENT(INOUT) :: snliqxy
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: tvxy
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: tgxy
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: canicexy
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: canliqxy
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: eahxy
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: tahxy
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: cmxy
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: chxy
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: fwetxy
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: sneqvoxy
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: alboldxy
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: qsnowxy
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: wslakexy
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: zwtxy
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: waxy
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: wtxy
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: lfmassxy
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: rtmassxy
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: stmassxy
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: woodxy
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: stblcpxy
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: fastcpxy
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: xsaixy

    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: t2mvxy
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: t2mbxy
    REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: chstarxy

    REAL, DIMENSION(1:NSOIL) :: ZSOIL

    REAL :: BX, SMCMAX, PSISAT

    REAL, PARAMETER :: BLIM = 5.5
    REAL, PARAMETER :: HLICE = 3.335E5
    REAL, PARAMETER :: GRAV = 9.81
    REAL, PARAMETER :: T0 = 273.15

    INTEGER :: errflag

    character(len=80) :: err_message
    character(len=4) :: MMINSL
    character(len=*), intent(in) :: MMINLU
    MMINSL='STAS'

    call read_mp_veg_parameters(trim(MMINLU))

    IF ( allowed_to_read ) THEN
       CALL wrf_message( 'INITIALIZE THREE Noah LSM RELATED TABLES' )
       CALL SOIL_VEG_GEN_PARM( MMINLU, MMINSL )
    ENDIF

    IF( .NOT. restart ) THEN

       itf=min0(ite,ide-1)
       jtf=min0(jte,jde-1)

       errflag = 0
       DO j = jts,jtf
          DO i = its,itf
             IF ( ISLTYP( i,j ) .LT. 1 ) THEN
                errflag = 1
                WRITE(err_message,*)"module_sf_noahlsm.F: lsminit: out of range ISLTYP ",i,j,ISLTYP( i,j )
                CALL wrf_message(err_message)
             ENDIF
          ENDDO
       ENDDO
       IF ( errflag .EQ. 1 ) THEN
          CALL wrf_error_fatal3("<stdin>",928,&
"module_sf_noahlsm.F: lsminit: out of range value "// &
               "of ISLTYP. Is this field in the input?" )
       ENDIF

       do I=1,NSLTYPE
          porosity(i)=maxsmc(i)
       enddo

       DO J = jts , jtf
          DO I = its , itf
     IF(IVGTYP(I,J)==ISICE .AND. XICE(I,J) <= 0.0) THEN
              DO NS=1, NSOIL
         SMOIS(I,NS,J) = 1.0
         SH2O(I,NS,J) = 0.0
         TSLB(I,NS,J) = MIN(TSLB(I,NS,J),263.15)
              END DO

  SNOW(I,J) = MAX(SNOW(I,J), 10.0)
     ELSE

             BX = BB(ISLTYP(I,J))
             SMCMAX = MAXSMC(ISLTYP(I,J))
             PSISAT = SATPSI(ISLTYP(I,J))
             IF ( ( BX > 0.0 ) .AND. ( SMCMAX > 0.0 ) .AND. ( PSISAT > 0.0 ) ) THEN
                DO NS=1, NSOIL
                   IF ( TSLB(I,NS,J) < 273.149 ) THEN
                      FK=(( (HLICE/(GRAV*(-PSISAT))) * &
                           ((TSLB(I,NS,J)-T0)/TSLB(I,NS,J)) )**(-1/BX) )*SMCMAX
                      FK = MAX(FK, 0.02)
                      SH2O(I,NS,J) = MIN( FK, SMOIS(I,NS,J) )
                   ELSE
                      SH2O(I,NS,J)=SMOIS(I,NS,J)
                   ENDIF
                END DO
             ELSE
                DO NS=1, NSOIL
                   SH2O(I,NS,J)=SMOIS(I,NS,J)
                END DO
             ENDIF
            ENDIF
          ENDDO
       ENDDO

       IF(.NOT.FNDSNOWH)THEN

          CALL wrf_message( 'SNOW HEIGHT NOT FOUND - VALUE DEFINED IN LSMINIT' )
          DO J = jts,jtf
             DO I = its,itf
                SNOWH(I,J)=SNOW(I,J)*0.005
             ENDDO
          ENDDO
       ENDIF

       DO J = jts,jtf
          DO I = its,itf
             tvxy (I,J) = TSK(I,J)
             tgxy (I,J) = TSK(I,J)
             CANWAT (I,J) = 0.0
             canliqxy (I,J) = CANWAT(I,J)
             canicexy (I,J) = 0.
             eahxy (I,J) = 2000.
             tahxy (I,J) = TSK(I,J)

             t2mvxy (I,J) = TSK(I,J)
             t2mbxy (I,J) = TSK(I,J)
             chstarxy (I,J) = 0.1

             cmxy (I,J) = 0.0
             chxy (I,J) = 0.0
             fwetxy (I,J) = 0.0
             sneqvoxy (I,J) = 0.0
             alboldxy (I,J) = 0.65
             qsnowxy (I,J) = 0.0
             wslakexy (I,J) = 0.0

             waxy (I,J) = 4900.
             wtxy (I,J) = waxy(i,j)
             zwtxy (I,J) = (25. + 2.0) - waxy(i,j)/1000/0.2

             lfmassxy (I,J) = 50.
             stmassxy (I,J) = 50.0
             rtmassxy (I,J) = 500.0
             woodxy (I,J) = 500.0
             stblcpxy (I,J) = 1000.0
             fastcpxy (I,J) = 1000.0
             xsaixy (I,J) = 0.1

          enddo
       enddo

       ZSOIL(1) = -DZS(1)
       DO NS=2, NSOIL
          ZSOIL(NS) = ZSOIL(NS-1) - DZS(NS)
       END DO

       CALL snow_init ( ims , ime , jms , jme , its , itf , jts , jtf , 3 , &
            & NSOIL , zsoil , snow , tgxy , snowh , &
            & zsnsoxy , tsnoxy , snicexy , snliqxy , isnowxy )

    ENDIF
  END SUBROUTINE NOAHMP_INIT

  SUBROUTINE SNOW_INIT ( ims , ime , jms , jme , its , itf , jts , jtf , &
       & NSNOW , NSOIL , ZSOIL , SWE , TGXY , SNODEP , &
       & ZSNSOXY , TSNOXY , SNICEXY ,SNLIQXY , ISNOWXY )
    IMPLICIT NONE
    INTEGER, INTENT(IN) :: ims, ime, jms, jme
    INTEGER, INTENT(IN) :: its, itf, jts, jtf
    INTEGER, INTENT(IN) :: NSNOW
    INTEGER, INTENT(IN) :: NSOIL
    REAL, INTENT(IN), DIMENSION(ims:ime, jms:jme) :: SWE
    REAL, INTENT(IN), DIMENSION(ims:ime, jms:jme) :: SNODEP
    REAL, INTENT(IN), DIMENSION(ims:ime, jms:jme) :: TGXY
    REAL, INTENT(IN), DIMENSION(1:NSOIL) :: ZSOIL
    INTEGER, INTENT(OUT), DIMENSION(ims:ime, jms:jme) :: ISNOWXY
    REAL, INTENT(OUT), DIMENSION(ims:ime, -NSNOW+1:NSOIL,jms:jme) :: ZSNSOXY
    REAL, INTENT(OUT), DIMENSION(ims:ime, -NSNOW+1: 0,jms:jme) :: TSNOXY
    REAL, INTENT(OUT), DIMENSION(ims:ime, -NSNOW+1: 0,jms:jme) :: SNICEXY
    REAL, INTENT(OUT), DIMENSION(ims:ime, -NSNOW+1: 0,jms:jme) :: SNLIQXY
    INTEGER :: I,J,IZ
    REAL, DIMENSION(-NSNOW+1: 0) :: DZSNO
    REAL, DIMENSION(-NSNOW+1:NSOIL) :: DZSNSO
    DO J = jts , jtf
       DO I = its , itf
          IF ( SNODEP(I,J) < 0.025 ) THEN
             ISNOWXY(I,J) = 0
             DZSNO(-NSNOW+1:0) = 0.
          ELSE
             IF ( ( SNODEP(I,J) >= 0.025 ) .AND. ( SNODEP(I,J) <= 0.05 ) ) THEN
                ISNOWXY(I,J) = -1
                DZSNO(0) = SNODEP(I,J)
             ELSE IF ( ( SNODEP(I,J) > 0.05 ) .AND. ( SNODEP(I,J) <= 0.10 ) ) THEN
                ISNOWXY(I,J) = -2
                DZSNO(-1) = SNODEP(I,J)/2.
                DZSNO( 0) = SNODEP(I,J)/2.
             ELSE IF ( (SNODEP(I,J) > 0.10 ) .AND. ( SNODEP(I,J) <= 0.25 ) ) THEN
                ISNOWXY(I,J) = -2
                DZSNO(-1) = 0.05
                DZSNO( 0) = SNODEP(I,J) - DZSNO(-1)
             ELSE IF ( ( SNODEP(I,J) > 0.25 ) .AND. ( SNODEP(I,J) <= 0.45 ) ) THEN
                ISNOWXY(I,J) = -3
                DZSNO(-2) = 0.05
                DZSNO(-1) = 0.5*(SNODEP(I,J)-DZSNO(-2))
                DZSNO( 0) = 0.5*(SNODEP(I,J)-DZSNO(-2))
             ELSE IF ( SNODEP(I,J) > 0.45 ) THEN
                ISNOWXY(I,J) = -3
                DZSNO(-2) = 0.05
                DZSNO(-1) = 0.20
                DZSNO( 0) = SNODEP(I,J) - DZSNO(-1) - DZSNO(-2)
             ELSE
                CALL wrf_error_fatal3("<stdin>",1118,&
"Problem with the logic assigning snow layers.")
             END IF
          END IF
          TSNOXY (I,-NSNOW+1:0,J) = 0.
          SNICEXY(I,-NSNOW+1:0,J) = 0.
          SNLIQXY(I,-NSNOW+1:0,J) = 0.
          DO IZ = ISNOWXY(I,J)+1 , 0
             TSNOXY(I,IZ,J) = TGXY(I,J)
             SNLIQXY(I,IZ,J) = 0.00
             SNICEXY(I,IZ,J) = 1.00 * DZSNO(IZ) * (SWE(I,J)/SNODEP(I,J))
          END DO
          DO IZ = ISNOWXY(I,J)+1 , 0
             DZSNSO(IZ) = -DZSNO(IZ)
          END DO
          DZSNSO(1) = ZSOIL(1)
          DO IZ = 2 , NSOIL
             DZSNSO(IZ) = (ZSOIL(IZ) - ZSOIL(IZ-1))
          END DO
          ZSNSOXY(I,ISNOWXY(I,J)+1,J) = DZSNSO(ISNOWXY(I,J)+1)
          DO IZ = ISNOWXY(I,J)+2 , NSOIL
             ZSNSOXY(I,IZ,J) = ZSNSOXY(I,IZ-1,J) + DZSNSO(IZ)
          ENDDO
       END DO
    END DO
  END SUBROUTINE SNOW_INIT
END MODULE module_sf_noahmpdrv


text:
subroutine inizializza

!--------------------------------------------------- !
!  F. Califano, 2006                                 !
!--------------------------------------------------- !

!**************************************************
!  MPI PARALLEL VERSION: VALENTINI-FAGANELLO 2009
!   3D PARALLEL VERSION: FAGANELLO 2010
!**************************************************

use box_mod
use deriv_mod
use parameter_mod
use poisson_open
use dom_distr_mod

IMPLICIT NONE

integer  :: iy, iz, ily, ilz

real(dp) :: zz, c1, p_a0, p_a1, p_b0, p_b1, p_a12
real(dp) :: alpha_1_3x3, alpha_1_3x5, alpha_1_bc

! --------------------------------- !
! *** Set up derivate Fourier       !
! --------------------------------- !

     call drffti( ny, savey )
     call drffti( nz, savez )

!---------------------------------------------- !
!   Fourier wave vector for y and z-derivative  !
!---------------------------------------------- !

     do iy = 2, ny1, 2
        ky_1(iy/2) = iy * hyl / (2.0 * ny)
        ky_4(iy/2) = iy * hyl / (2.0)
     enddo

      ky_4 = ky_4**4 / ny

     do iz = 2, nz1, 2
        kz_1(iz/2) = iz * hzl / (2.0 * nz)
        kz_4(iz/2) = iz * hzl / (2.0)
     enddo

      kz_4 = kz_4**4 / nz

! --------------------------------- !
! *** Set up derivate Compatte      !
! *** Parametri solver tridiagonale !
! --------------------------------- !

    TRANS = 'N'
    info  = 0
    nrhs  = 1
    ndb   = nx

!    L.H.S. coefficients of implicit I

        alpha_1_3x3 = 1.0d0 / 4.0d0
        alpha_1_3x5 = 1.0d0 / 3.0d0

!    R.H.S. coefficients of implicit I derivative.

        a_1_3x3 = 3.0d0 / (4.0  * dx)
        a_1_3x5 = 7.0d0 / (9.0  * dx)
        b_1_3x5 = 1.0d0 / (36.0 * dx)

!    L.H.S coefficients for FREE boundary conditions, I derivative.

        alpha_1_bc = 3.0d0

!    R.H.S coefficients boundary conditions of I derivative.

        a_1_bc = - (11.0d0 + 2.0d0 * alpha_1_bc) / (6.0d0 * dx)
        b_1_bc = (6.0d0 - alpha_1_bc) / (2.0d0 * dx)
        c_1_bc = (2.0d0 * alpha_1_bc - 3.0d0) / (2.0d0 * dx)
        d_1_bc = (2.0d0 - alpha_1_bc) / (6.0d0 * dx)

!        Differenze finite compatte 3x5
!        R.H.S arrays for LU decomposition, I derivative
!        TRE CASI, d.., w.., q.. a seconda delle cond. contorno

!             Diagonale

        dd_1 = 1.0
        wd_1 = 1.0
        qd_1 = 1.0

!        dm_1 e dp_1, wm_1 e wp_1, qm_1 e qp_1 sono
!        sotto e sopra diagonale, dim. nx - 1

        dm_1 = alpha_1_3x5
        dp_1 = alpha_1_3x5

        wm_1 = alpha_1_3x5
        wp_1 = alpha_1_3x5

        qm_1 = alpha_1_3x5
        qp_1 = alpha_1_3x5

!    Differenze finite compatte 3x3 in i = 2, nx-1

        dm_1(1)   = alpha_1_3x3
        dp_1(2)   = alpha_1_3x3
        dm_1(nx2) = alpha_1_3x3
        dp_1(nx1) = alpha_1_3x3

        wm_1(1)   = alpha_1_3x3
        wp_1(2)   = alpha_1_3x3
        wm_1(nx2) = alpha_1_3x3
        wp_1(nx1) = alpha_1_3x3

        qm_1(1)   = alpha_1_3x3
        qp_1(2)   = alpha_1_3x3
        qm_1(nx2) = alpha_1_3x3
        qp_1(nx1) = alpha_1_3x3

!     Working arrays

        dw_1 = 0.0
        ww_1 = 0.0
        qq_1 = 0.0

!!! **** CONDIZIONI AL CONTORNO *** !!!

!     Boundary condition calculated with internal points
!        (left and right), I derivative, i=1, i=nx

        dp_1(1)   = alpha_1_bc
        dm_1(nx1) = alpha_1_bc

!     FREE sleep boundary condition (left and right)
!       d/dx = 0 ai bordi

        wp_1(1)   = 0.0
        wm_1(nx1) = 0.0

!     FREE^2 sleep boundary condition (left and right)
!       d/dx = cste ai bordi (f'_1 - f'_2 = 0 e idem a destra)

        qp_1(1)   = - 1.0
        qm_1(nx1) = - 1.0

! DERIVATA SENZA COND. AL CONTORNO (si usano i punti interni al bordo)
!     Fattorizzazione LU derivata I

    CALL DGTTRF(nx, dm_1, dd_1, dp_1, dw_1, ipv_d, info)

      if (info > 0 .or. info < 0) then
      write(*,*) 'Problemi fattorizzazione LU der I, info:', info
      stop
      endif

! DERIVATA CON COND. AL CONTORNO di derivata nulla
!     Fattorizzazione LU derivata I

    CALL DGTTRF(nx, wm_1, wd_1, wp_1, ww_1, ipv_w, info)

      if (info > 0 .or. info < 0) then
      write(*,*) 'Problemi fattorizzazione LU der I, info:', info
      stop
      endif

! DERIVATA CON COND. AL CONTORNO di derivata costante
!     Fattorizzazione LU derivata I

    CALL DGTTRF(nx, qm_1, qd_1, qp_1, qq_1, ipv_q, info)

      if (info > 0 .or. info < 0) then
      write(*,*) 'Problemi fattorizzazione LU der I, info:', info
      stop
      endif

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!  SOLVER DI POISSON direzione con cond. cont. open

! ------------------------------------------------------------ !
! *** Set up coefficienti derivate Compatte per Poisson solver !
! *** Parametri solver tridiagonale                            !
! ------------------------------------------------------------ !

    TRNAS = 'N'
    ifno  = 0
    irhs  = 1
    idb   = nx

do ilz = 0, n2z-1
 do ily = 0, ny/2-1

! ===================================================================== !
! coefficients for the EQUATION:                                        !
!                                                                       !
! E_i + (a1/a0) (E_i-1 + E_i+1) = (b1/a0) (f_i-1 + f_i+1) + (b0/a0) f_i !
!                                                                       !
!         viene da (1 + de^2 nabla^2) E = RHS                           !
!     vedi Twofluids_equation                                           !
! ===================================================================== !

 c1 = rapm + ily * ily * hyl * hyl + ilz * ilz * hzl * hzl     ! rapm = 1 / de^2

p_a0 = 10.0 * c1 * dx * dx + 24.0
p_a1 =        c1 * dx * dx - 12.0
p_b0 = 10.0 * dx * dx * rapm
p_b1 =        dx * dx * rapm

! ======================================= !
! R.H.S. coefficients                     !
! ======================================= !

pss_rhs_0(ily,ilz) = p_b0 / p_a0
pss_rhs_1(ily,ilz) = p_b1 / p_a0

! ======================================= !
! L.H.S. coefficients                     !
! ======================================= !

p_a12 = p_a1 / p_a0

  poiss_d(:,ily,ilz) = 1.0
  poiss_m(:,ily,ilz) = p_a12
  poiss_p(:,ily,ilz) = p_a12

  poiss_w = 0.d0

! condizioni al contorno: E = cost a sinistra e a destra per tutti i modi

  poiss_p(1,ily,ilz)   = 0.0d0
  poiss_m(nx1,ily,ilz) = 0.0d0

! condizioni al contorno: dE/dx = 0.0 a sinistra e a destra per il modo ky=0

!  poiss_p(1,il)   = - 1.0d0
!  poiss_m(nx1,il) = - 1.0d0

 CALL DGTTRF(nx, poiss_m(:,ily,ilz), poiss_d(:,ily,ilz), poiss_p(:,ily,ilz), poiss_w, ipv_p(:,ily,ilz), ifno)

      if (ifno > 0 .or. ifno < 0) then
      write(*,*) 'Problemi fattorizzazione LU, info:', info
      stop
      endif

 enddo
enddo

end subroutine

max_stars_repo_name: sscerr/eTF

max_issues_repo_name: sscerr/eTF

max_forks_repo_name: sscerr/eTF


! { dg-do run }
! PR56786 Error on embedded spaces
integer :: i(3)
namelist /nml/ i

i = -42
open(99,status='scratch')
write(99,'(a)') '&nml i(3 ) = 5 /'
rewind(99)
read(99,nml=nml)
close(99)
if (i(1)/=-42 .or. i(2)/=-42 .or. i(3)/=5) STOP 1

! Shorten the file so the read hits EOF

open(99,status='scratch')
write(99,'(a)') '&nml i(3 ) = 5 '
rewind(99)
read(99,nml=nml, end=30)
STOP 2
! Shorten some more
 30 close(99)
open(99,status='scratch')
write(99,'(a)') '&nml i(3 ) ='
rewind(99)
read(99,nml=nml, end=40)
STOP 3
! Shorten some more
 40 close(99)
open(99,status='scratch')
write(99,'(a)') '&nml i(3 )'
rewind(99)
read(99,nml=nml, end=50)
STOP 4
! Shorten some more
 50 close(99)
open(99,status='scratch')
write(99,'(a)') '&nml i(3 '
rewind(99)
read(99,nml=nml, end=60)
STOP 5
 60 close(99)
end


!        Generated by TAPENADE     (INRIA, Ecuador team)
!  tapenade 3.x
!
!  Differentiation of rhow in forward (tangent) mode:
!   variations   of useful results: rhow
!   with respect to varying inputs: p t
! MIT License
!
! Copyright (c) 2020 SHEMAT-Suite
!
! Permission is hereby granted, free of charge, to any person obtaining a copy
! of this software and associated documentation files (the "Software"), to deal
! in the Software without restriction, including without limitation the rights
! to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
! copies of the Software, and to permit persons to whom the Software is
! furnished to do so, subject to the following conditions:
!
! The above copyright notice and this permission notice shall be included in all
! copies or substantial portions of the Software.
!
! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
! IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
! AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
! OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
! SOFTWARE.
!> @brief rhow calculates the density of pure water [kg/m3]
!> @param[in] p pressure [MPa]
!> @param[in] t temperature [degC]
!> @return rhow density of pure water [kg/m3]
!> @details
!> rhow(i,j,k,ismpl) calculates the density in (in kg/m^3) of pure
!> water, given temperature (t, in degC), and pressure (p,in Pa)\n \n
!>
!> Source: \n
!>
!> Magri, F. (2005), Mechanismus und Fluiddynamik der
!> Salzwasserzirkulation im Norddeutschen Becken: Ergebnisse
!> thermohaliner numerischer Simulationen, Doktorarbeit,
!> Geoforschungszentrum Potsdam \n \n
!>
!> Specifically: Equation (2.12) and Table 2-1
DOUBLE PRECISION FUNCTION g_RHOW(p, g_p, t, g_t, rhow)
  IMPLICIT NONE
! Local temperature (degC)
  DOUBLE PRECISION, INTENT(IN) :: t
  DOUBLE PRECISION, INTENT(IN) :: g_t
! Local pressure (MPa)
  DOUBLE PRECISION, INTENT(IN) :: p
  DOUBLE PRECISION, INTENT(IN) :: g_p
! Pressure in kPa
  DOUBLE PRECISION :: pk
  DOUBLE PRECISION :: g_pk
! Monomials in pressure and temperature
  DOUBLE PRECISION :: pk2, t2, t4, t6
  DOUBLE PRECISION :: g_pk2, g_t2, g_t4, g_t6
! Coefficients from Magri2005 Equation (2.12)
  DOUBLE PRECISION :: a, b, c, d, e, f, g
  DOUBLE PRECISION :: g_a, g_b, g_c, g_d, g_e, g_f, g_g
! Coefficients from Magri2005 Table 2-1
  DOUBLE PRECISION, PARAMETER :: a0=9.99792877961606d02
  DOUBLE PRECISION, PARAMETER :: a1=5.07605113140940d-04
  DOUBLE PRECISION, PARAMETER :: a2=-5.28425478164183d-10
  DOUBLE PRECISION, PARAMETER :: b0=5.13864847162196d-02
  DOUBLE PRECISION, PARAMETER :: b1=-3.61991396354483d-06
  DOUBLE PRECISION, PARAMETER :: b2=7.97204102509724d-12
  DOUBLE PRECISION, PARAMETER :: c0=-7.53557031774437d-03
  DOUBLE PRECISION, PARAMETER :: c1=6.32712093275576d-08
  DOUBLE PRECISION, PARAMETER :: c2=-1.66203631393248d-13
  DOUBLE PRECISION, PARAMETER :: d0=4.60380647957350d-05
  DOUBLE PRECISION, PARAMETER :: d1=-5.61299059722121d-10
  DOUBLE PRECISION, PARAMETER :: d2=1.80924436489400d-15
  DOUBLE PRECISION, PARAMETER :: e0=-2.26651454175013d-07
  DOUBLE PRECISION, PARAMETER :: e1=3.36874416675978d-12
  DOUBLE PRECISION, PARAMETER :: e2=-1.30352149261326d-17
  DOUBLE PRECISION, PARAMETER :: f0=6.14889851856743d-10
  DOUBLE PRECISION, PARAMETER :: f1=-1.06165223196756d-14
  DOUBLE PRECISION, PARAMETER :: f2=4.75014903737416d-20
  DOUBLE PRECISION, PARAMETER :: g0=-7.39221950969522d-13
  DOUBLE PRECISION, PARAMETER :: g1=1.42790422913922d-17
  DOUBLE PRECISION, PARAMETER :: g2=-7.13130230531541d-23
  DOUBLE PRECISION :: rhow
! Unit [kPa] is needed after Magri(2005); the forwarded p is in [MPa]
  g_pk = 1000*g_p
  pk = p*1000
! Monomials in pressure and temperature
  g_pk2 = 2*pk*g_pk
  pk2 = pk*pk
  g_t2 = 2*t*g_t
  t2 = t*t
  g_t4 = 2*t2*g_t2
  t4 = t2*t2
  g_t6 = t2*g_t4 + t4*g_t2
  t6 = t4*t2
! Coefficients from Magri2005, Equation (2.12)
  g_a = a1*g_pk + a2*g_pk2
  a = a0 + a1*pk + a2*pk2
  g_b = b1*g_pk + b2*g_pk2
  b = b0 + b1*pk + b2*pk2
  g_c = c1*g_pk + c2*g_pk2
  c = c0 + c1*pk + c2*pk2
  g_d = d1*g_pk + d2*g_pk2
  d = d0 + d1*pk + d2*pk2
  g_e = e1*g_pk + e2*g_pk2
  e = e0 + e1*pk + e2*pk2
  g_f = f1*g_pk + f2*g_pk2
  f = f0 + f1*pk + f2*pk2
  g_g = g1*g_pk + g2*g_pk2
  g = g0 + g1*pk + g2*pk2
! Water density [kg/m3]
  g_rhow = g_a + t*g_b + (b+d*t2+f*t4)*g_t + t2*g_c + c*g_t2&
&   + t*(t2*g_d+d*g_t2) + t4*g_e + e*g_t4 + t*(t4*g_f+f*g_t4&
&   ) + t6*g_g + g*g_t6
  rhow = a + b*t + c*t2 + d*t2*t + e*t4 + f*t4*t + g*t6
  RETURN
END FUNCTION g_RHOW


text:
C Copyright(C) 1999-2020 National Technology & Engineering Solutions
C of Sandia, LLC (NTESS).  Under the terms of Contract DE-NA0003525 with
C NTESS, the U.S. Government retains certain rights in this software.
C
C See packages/seacas/LICENSE for details

C=======================================================================
      SUBROUTINE SHODEV (SHOTYP)
C=======================================================================

C   --*** SHODEV *** (BLOT) Display device dependent options
C   --   Written by Amy Gilkey - revised 04/26/88
C   --
C   --SHODEV displays the device dependent options requested by the
C   --option type:
C   --   SOFTCHAR - software versus hardware characters
C   --   FONT - font to use
C   --   COLOR - number of standard colors to use and maximum on device
C   --   SPECTRUM - number of spectrum colors and maximum on device
C   --   SNAP - number of frames to snap
C   --   AUTO - automatic versus user-directed plotting
C   --
C   --Parameters:
C   --   SHOTYP - IN - the show option (see above)

C   --Routines Called:
C   --   GRGPAR - (GRPLIB) Get graphics device parameters
C   --   INTSTR - (STRLIB) Convert integers to strings
C   --   LENSTR - (STRLIB) Find string length
C   --   SQZSTR - (STRLIB) Delete extra blanks from string

      CHARACTER*(*) SHOTYP

      LOGICAL ONEDEV, TWODEV
      CHARACTER*3 DEVNAM(2)
      INTEGER IPARMS(3)
      CHARACTER*80 STRING
      CHARACTER*20 STR20

      CALL GRGPARD ('DEVICE', 1, ONEDEV, DEVNAM(1))
      CALL GRGPARD ('DEVICE', 2, TWODEV, DEVNAM(2))
      ISDEV = 1
      IF (.NOT. ONEDEV) ISDEV = 2
      IEDEV = 2
      IF (.NOT. TWODEV) IEDEV = 1
      TWODEV = (ONEDEV .AND. TWODEV)

      IF (SHOTYP .EQ. 'SOFTCHAR') THEN
         DO 100 I = ISDEV, IEDEV
            CALL GRGPAR (SHOTYP, I, IPARMS, STR20)
            L = LENSTR (STR20)
            STRING = 'Plot in ' // STR20(:L) // ' characters$'
            CALL PRTDEV (STRING, DEVNAM(I), TWODEV)
  100    CONTINUE

      ELSE IF (SHOTYP .EQ. 'FONT') THEN
         DO 110 I = ISDEV, IEDEV
            CALL GRGPAR (SHOTYP, I, IPARMS, STR20)
            L = LENSTR(STR20)
            STRING = 'Plot in ' // STR20(:L) // ' font$'
            CALL PRTDEV (STRING, DEVNAM(I), TWODEV)
  110    CONTINUE

      ELSE IF ((SHOTYP .EQ. 'COLOR') .OR. (SHOTYP .EQ. 'SPECTRUM')) THEN
         DO 120 I = ISDEV, IEDEV
            CALL GRGPAR ('COLOR', I, IPARMS, STR20)
            WRITE (STR20, 10000, IOSTAT=IDUM) IPARMS(1), IPARMS(2)
10000        FORMAT (I6, ' of ', I6)
            CALL SQZSTR (STR20, L)
            IF (IPARMS(2) .EQ. 0) THEN
               STRING = 'Number of colors to use$: single color device'
               CALL PRTDEV (STRING, DEVNAM(I), TWODEV)
            ELSE
               STRING = 'Number of standard colors to use$: '
     &            // STR20(:L)
               CALL PRTDEV (STRING, DEVNAM(I), TWODEV)
            END IF

            CALL GRGPAR ('SPECTRUM', I, IPARMS, STR20)
            IF (IPARMS(3) .GT. 0) THEN
               WRITE (STR20, 10000, IOSTAT=IDUM) IPARMS(1), IPARMS(2)
               CALL SQZSTR (STR20, L)
               IF (IPARMS(3) .EQ. 1) THEN
                  STRING = 'Number of spectrum colors to use$: '
     &               // STR20(:L)
                  CALL PRTDEV (STRING, DEVNAM(I), TWODEV)
               ELSE
                  STRING = '***$: ' // STR20(:L)
                  CALL PRTDEV (STRING, DEVNAM(I), TWODEV)
               END IF
            END IF
  120    CONTINUE

      ELSE IF (SHOTYP .EQ. 'SNAP') THEN
         I = ISDEV
         CALL GRGPAR (SHOTYP, I, IPARMS, STR20)
         IF (IPARMS(1) .LT. 0) THEN
            STRING =  'Cannot snap frames$'
            CALL PRTDEV (STRING, DEVNAM(I), .TRUE.)
         ELSE IF (IPARMS(1) .LE. 0) THEN
            STRING =  'Do not snap frames$'
            CALL PRTDEV (STRING, DEVNAM(I), TWODEV)
         ELSE
            CALL INTSTR (1, 0, IPARMS(1), STR20, L)
            STRING = 'Snap ' // STR20(:L) // ' frames for each plot$'
            CALL PRTDEV (STRING, DEVNAM(I), TWODEV)
         END IF

      ELSE IF (SHOTYP .EQ. 'AUTO') THEN
         I = ISDEV
         CALL GRGPAR (SHOTYP, I, IPARMS, STR20)
         L = LENSTR (STR20)
         STRING = 'Plot in ' // STR20(:L) // ' mode$'
         CALL PRTDEV (STRING, DEVNAM(I), TWODEV)
      END IF

      RETURN
      END

ext: f

max_stars_repo_name: jschueller/seacas

max_issues_repo_name: jschueller/seacas

max_forks_repo_name: jschueller/seacas


text:
subroutine findnode(gridfile,slat,slon,sdep,coords,ns,box,izone)
  implicit none
  integer ILONGLAT2UTM,IUTM2LONGLAT,izone
  parameter(ILONGLAT2UTM=0,IUTM2LONGLAT=1)
  integer, dimension(3) :: dims,coords,ns
  double precision :: slat,slon,sdep, dzmin
  double precision :: dlat,dlon,xx,yy,difx, e, n
  double precision, dimension(3000) :: xgrid,ygrid,zgrid
  double precision, dimension(4,2) :: box
  character (len=200) :: gridfile
  integer ix, iy, iz, nxs, nys, nzs, mxs, mys, mzs, nxst, nyst, nzst,nx,ny,nz
  integer ii, i, nxt, nyt, nzt

  dlat=dble(slat)
  dlon=dble(slon)
  call utm_geo(dlon,dlat,xx,yy,izone,ILONGLAT2UTM)

  !   Read info defining the FD grid.

  open(10,file=gridfile,status='old')
  read(10,*) nx,ny,nz
  read(10,*) (dims(i),i=1,3)
  read(10,*) nxt,nyt,nzt
  read(10,*) box(1,1),box(1,2)
  read(10,*) box(2,1),box(2,2)
  read(10,*) box(3,1),box(3,2)
  read(10,*) box(4,1),box(4,2)
  do iz=1,nz
     read(10,*) ii,zgrid(ii)
  end do
  close(10)

! Find out the source node numbers in the global grid.
  call utm2grd(nxs,nys,xx,yy,box,nx,ny)
  if (nxs.lt.1.or.nxs.gt.nx.or.nys.lt.1.or.nys.gt.ny.or.sdep.gt.-zgrid(1)/1e3) then
     print*, 'receiver located outside model box!'
     print*, 'receiver gird: ', nxs, nys, sdep
     print*, 'modeling grid: ', nx, ny, -zgrid(1)/1e3
  endif
  dzmin=10000
  if (sdep.le.0.0) then
     nzs=nz
  else
     do iz=1,nz
        if (abs(abs(zgrid(iz))-abs(sdep)*1e3).le.dzmin) then
           dzmin=abs(abs(zgrid(iz))-abs(sdep)*1e3)
           nzs=iz
        endif
     enddo
  endif

  call grd2utm(e,n,nxs,nys,box,nx,ny)
  call utm_geo(dlon,dlat,e,n,izone,IUTM2LONGLAT)

! calling utm_geo back and forth will cause a difference in northing of about 12 m !!!!!!!
! strange !!!!

!  print*, dlon, dlat
!  call utm_geo(dlon,dlat,e,n,izone,ILONGLAT2UTM)
!  print*, e, n
!  call utm_geo(dlon,dlat,e,n,izone,IUTM2LONGLAT)
!  print*, dlon, dlat
!  call utm_geo(dlon,dlat,e,n,izone,ILONGLAT2UTM)
!  print*, e, n
!  stop
!   Find out the source chunk numbers.

  mxs=(nxs-1)/nxt
  mys=(nys-1)/nyt
  mzs=(nzs-1)/nzt

!   Find out the source node numbers in the local grid.

  nxst=nxs-mxs*nxt
  nyst=nys-mys*nyt
  nzst=nzs-mzs*nzt

  write(*,*) ''
  write(*,*) '   Number of chunks (XYZ):    ',dims(1),dims(2),dims(3)
  write(*,*) '   Total grid numbers (XYZ):  ',nx,ny,nz
  write(*,*) '   Grids per chunk (XYZ):     ',nxt,nyt,nzt
  write(*,*) ''
  write(*,'(a,5(4x,e20.12))') '    Input receiver location:     ', xx,yy,sdep*1e3,slon,slat
  write(*,'(a,5(4x,e20.12))') '    Used receiver location:      ', e,n,-zgrid(nzs),dlon,dlat
  write(*,'(a,3(4x,e20.12))') '    Location differences (m):    ',abs(e-xx),abs(n-yy),abs(-zgrid(nzs)-sdep*1e3)
  write(*,*) '   Receiver grid (global):      ',nxs,nys,nzs
  write(*,*) '   Receiver chunk (0- ):        ',mxs,mys,mzs
  write(*,*) '   Receiver grid input (local): ',nxst,nyst,nzst

  coords(1)=mxs
  coords(2)=mys
  coords(3)=mzs

  ns(1)=nxst
  ns(2)=nyst
  ns(3)=nzst

  return
end subroutine findnode
!
!

max_stars_repo_name: SCECcode/cvmsi

max_issues_repo_name: SCECcode/cvmsi

max_forks_repo_name: SCECcode/cvmsi


text:
program lcstest
  use iso_varying_string
  implicit none

  type(varying_string) :: s1, s2

  s1 = "thisisatest"
  s2 = "testing123testing"
  print *, char(lcs(s1, s2))

  s1 = "1234"
  s2 = "1224533324"
  print *, char(lcs(s1, s2))

contains

  recursive function lcs(a, b) result(l)
    type(varying_string) :: l
    type(varying_string), intent(in) :: a, b

    type(varying_string) :: x, y

    l = ""
    if ( (len(a) == 0) .or. (len(b) == 0) ) return
    if ( extract(a, len(a), len(a)) == extract(b, len(b), len(b)) ) then
       l = lcs(extract(a, 1, len(a)-1), extract(b, 1, len(b)-1)) // extract(a, len(a), len(a))
    else
       x = lcs(a, extract(b, 1, len(b)-1))
       y = lcs(extract(a, 1, len(a)-1), b)
       if ( len(x) > len(y) ) then
          l = x
       else
          l = y
       end if
    end if
  end function lcs

end program lcstest

max_issues_repo_name: seanwallawalla-forks/RosettaCodeData

max_forks_repo_name: seanwallawalla-forks/RosettaCodeData


!----------------------------------------------------------------------
!----------------------------------------------------------------------
!   nag :    Heteroclinic orbits : A saddle-saddle connection
!----------------------------------------------------------------------
!----------------------------------------------------------------------
! Parameter assignment:
!
!           PAR(1) , PAR(2) : a      , c      (parameters)
!           PAR(3) , PAR(4) : eps-0  , eps-1  (radii)
!           PAR(5) , PAR(6) : mu-0   , mu-1   (eigenvalues)
!           PAR(7) , PAR(8) : v-0(1) , v-0(2) (eigenvector)
!           PAR(9) , PAR(10): v-1(1) , v-1(2) (eigenvector)
!           PAR(11)         : period
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------

      SUBROUTINE FUNC(NDIM,U,ICP,PAR,IJAC,F,DFDU,DFDP)
!     ---------- ----

      IMPLICIT NONE

      INTEGER, INTENT(IN) :: NDIM, ICP(*), IJAC
      DOUBLE PRECISION, INTENT(IN) :: U(NDIM), PAR(*)
      DOUBLE PRECISION, INTENT(OUT) :: F(NDIM)
      DOUBLE PRECISION, INTENT(INOUT) :: DFDU(NDIM,NDIM), DFDP(NDIM,*)

      DOUBLE PRECISION PERIOD, DUMMY(1)

       CALL FFFF(2,U,ICP,PAR,0,F,DUMMY)

       PERIOD=PAR(11)
       F(1)=PERIOD*F(1)
       F(2)=PERIOD*F(2)

      END SUBROUTINE FUNC

      SUBROUTINE FFFF(NDM,U,ICP,PAR,IJAC,F,DFDU)
!     ---------- ----

      IMPLICIT NONE
      INTEGER, INTENT(IN) :: NDM, ICP(*), IJAC
      DOUBLE PRECISION, INTENT(IN) :: U(NDM), PAR(*)
      DOUBLE PRECISION, INTENT(OUT) :: F(NDM)
      DOUBLE PRECISION, INTENT(INOUT) :: DFDU(NDM,NDM)

      DOUBLE PRECISION A,C

       A=PAR(1)
       C=PAR(2)
       F(1)= U(2)
       F(2)= C*U(2) - U(1) * (1-U(1)) * (U(1)-A)

      IF(IJAC.EQ.0)RETURN

       DFDU(1,1)= 0
       DFDU(1,2)= 1

       DFDU(2,1)= - (1-U(1))*(U(1)-A) + U(1)*(U(1)-A) - U(1)*(1-U(1))
       DFDU(2,2)= C

      END SUBROUTINE FFFF

      SUBROUTINE STPNT(NDIM,U,PAR,T)
!     ---------- -----

      IMPLICIT NONE
      INTEGER, INTENT(IN) :: NDIM
      DOUBLE PRECISION, INTENT(INOUT) :: U(NDIM),PAR(*)
      DOUBLE PRECISION, INTENT(IN) :: T

      DOUBLE PRECISION PERIOD, SIGN, R, V1, V2, X, E, U0, V0, EPS

! The following initialization is done only in the first call
       IF(T==0)THEN
!        Select period and connection (read from the constants-file)
         PERIOD=100.
         SIGN=1.0
         R = 0.5*SQRT(2.D0)
         V1= 1./SQRT(1+R**2)
         V2= R /SQRT(1+R**2)
         X  = -0.5*PERIOD
         E  = EXP(X/SQRT(2.D0))
         U0 = E/(1+E)
         V0 = ( E/(1+E)**2 ) / SQRT(2.D0)
         EPS= SQRT( U0**2 + V0**2 )
         PAR(1)= 0.5
         PAR(2)= 0
         PAR(3)= EPS
         PAR(4)= EPS
         PAR(5)= R
         PAR(6)=-R
         PAR(7)= SIGN*V1
         PAR(8)= SIGN*V2
         PAR(9)= -SIGN*V1
         PAR(10)= SIGN*V2
         PAR(11)= PERIOD
         PAR(12)= SIGN
       ENDIF

! Specify exact solution as starting point :

       PERIOD=PAR(11)
       SIGN  =PAR(12)

       X=PERIOD*SIGN*(T-0.5)
       E=DEXP(X/DSQRT(2.D0))
       U(1)=E/(1+E)
       U(2)=SIGN*E/(1+E)**2/DSQRT(2.D0)

      END SUBROUTINE STPNT

      SUBROUTINE BCND(NDIM,PAR,ICP,NBC,U0,U1,FB,IJAC,DBC)
!     ---------- ----

      IMPLICIT NONE

      INTEGER, INTENT(IN) :: NDIM, ICP(*), NBC, IJAC
      DOUBLE PRECISION, INTENT(IN) :: PAR(*), U0(NDIM), U1(NDIM)
      DOUBLE PRECISION, INTENT(OUT) :: FB(NBC)
      DOUBLE PRECISION, INTENT(INOUT) :: DBC(NBC,*)
! Local
      INTEGER, PARAMETER :: NDM=2
      DOUBLE PRECISION V0(NDM),V1(NDM),G0(NDM),G1(NDM)
      DOUBLE PRECISION DGDU0(NDM,NDM),DGDU1(NDM,NDM)

      V0(1)=U0(1) - PAR(3)*PAR(7)
      V0(2)=U0(2) - PAR(3)*PAR(8)
      V1(1)=U1(1) - PAR(4)*PAR(9)
      V1(2)=U1(2) - PAR(4)*PAR(10)

      CALL FFFF(NDM,V0,ICP,PAR,1,G0,DGDU0)
      CALL FFFF(NDM,V1,ICP,PAR,1,G1,DGDU1)

      FB(1)= DGDU0(1,1)*PAR(7) + DGDU0(1,2)*PAR(8) - PAR(5)*PAR(7)
      FB(2)= DGDU0(2,1)*PAR(7) + DGDU0(2,2)*PAR(8) - PAR(5)*PAR(8)
      FB(3)= DGDU1(1,1)*PAR(9) + DGDU1(1,2)*PAR(10)- PAR(6)*PAR(9)
      FB(4)= DGDU1(2,1)*PAR(9) + DGDU1(2,2)*PAR(10)- PAR(6)*PAR(10)
      FB(5)= PAR(7)**2 + PAR(8)**2 -1
      FB(6)= PAR(9)**2 + PAR(10)**2 -1
      FB(7)= G0(1)
      FB(8)= G0(2)
      FB(9)= G1(1)
      FB(10)=G1(2)

      END SUBROUTINE BCND

      SUBROUTINE ICND(NDIM,PAR,ICP,NINT,U,UOLD,UDOT,UPOLD,FI,IJAC,DINT)
!     ---------- ----

      IMPLICIT NONE

      INTEGER, INTENT(IN) :: NDIM, ICP(*), NINT, IJAC
      DOUBLE PRECISION, INTENT(IN) :: PAR(*)
      DOUBLE PRECISION, INTENT(IN) :: U(NDIM), UOLD(NDIM), UDOT(NDIM), UPOLD(NDIM)
      DOUBLE PRECISION, INTENT(OUT) :: FI(NINT)
      DOUBLE PRECISION, INTENT(INOUT) :: DINT(NINT,*)
! Local
      INTEGER, PARAMETER :: NDM=2
      DOUBLE PRECISION F(NDM),F0(NDM),DFDU(NDM,NDM),DUMMY(1)

      CALL FFFF(NDM,U   ,ICP,PAR,1,F ,DFDU)
      CALL FFFF(NDM,UOLD,ICP,PAR,0,F0,DUMMY)

       FI(1)= ( F(1) - F0(1) ) * ( DFDU(1,1)*F(1) + DFDU(1,2)*F(2) ) &
            + ( F(2) - F0(2) ) * ( DFDU(2,1)*F(1) + DFDU(2,2)*F(2) )

      END SUBROUTINE ICND

      SUBROUTINE FOPT
      END SUBROUTINE FOPT

      SUBROUTINE PVLS
      END SUBROUTINE PVLS