USE mg3d.cpl
!  USE rtchecks.cpl
! stok3d
FILE ALTEZZE
ALTEZZE=CREATE('sto')
nx=2∗8
ny=nx∗6
!     nz=nx∗3 DIV 2 
nz=2∗nx

hmax=1.
! valore massimo di h compatibile con la memoria prevista

INTEGER ny0
ny0=FLOOR(2∗nx∗hmax) +5
!   ny0=nx∗2

pi=ATAN(1)∗4 
pfo=0.5
! parametro formula
pbc=2
protra = 2
! protra1 =lamina; protra2 = e5; protra3 = profilo a v;protra4=coseno 
prolo = 3
! prolo1 =2d;  prolo2 = profilo iperbolico; prolo3=BE 

inco=1
! inco0 =senza correzione; inco1 = con correzione

REAL fx0

IF protra=1 THEN
  fx0 =0
ELSE IF protra=2 THEN
  fx0= 0.05
ELSE IF protra=3 THEN
  fx0= 0.25
ELSE IF protra=4 THEN
  fx0= pi/2
END IF 

INTEGER it,js,wp,btm(-1..nx+1,-1..nz+1),btf(0..nx,0..nz)
REAL ad(0..nx,0..ny0,0..nz),al(0..nx,0..ny0,0..nz),ar(0..nx,0..ny0,0..nz), \
ab(0..nx,0..ny0,0..nz),af(0..nx,0..ny0,0..nz),rhs(0..nx,0..ny,0..nz)
REAL up(0..nx,0..ny0,0..nz),right(0..nx,0..ny0,0..nz),down(0..nx,0..ny0,0..nz),\
left(0..nx,0..ny0,0..nz),back(0..nx,0..ny0,0..nz),forw(0..nx,0..ny0,0..nz)
REAL bcu(0..nx,0..ny0,0..nz),bcr(0..nx,0..ny0,0..nz),bcd(0..nx,0..ny0,0..nz), \
bcl(0..nx,0..ny0,0..nz),bcb(0..nx,0..ny0,0..nz),bcf(0..nx,0..ny0,0..nz)
REAL coenoto(0..nx,0..ny0,0..nz)
REAL m1(-1..nx+1,0..ny,-1..nz+1),m2(-1..nx+1,0..ny,-1..nz+1), \
m3(-1..nx+1,0..ny,-1..nz+1),fi(-1..nx+1,0..ny,-1..nz+1)


REAL vzn(1..1500),vsalfa(1..1500),vcalfa(1..1500),vdel(1..1500),vk1(1..1500),vk2(1..1500)
REAL  zn,salfa,calfa,del,k1,k2,c2,c4,c6,c2l
INTEGER nid(0..nx,0..ny0,0..nz)
REAL fill,fillm,fils,deno,rapp,hlon,htra,dh,prec,ka 
REAL denx,deny,denz,ris(1..50,1..50)
REAL ax0,axl,axr,ay0,ayd,ayu,az0,azf,azb
INTEGER datof    ! datofalso per interrompere con READ
jsta=9
numh=1
numk=1
numf=16
k0=1
delth=0.200
deltf=0.050
fil0=0.10-deltf
rap0=0.300- delth
lx=1.

h=lx/nx/2
ly=ny∗h
lz=nz∗h∗4.
eps=h/100
dcor=2.05∗h


REAL FUNCTION tanh( REAL xx)
  RESULT=(EXP(xx)-EXP(-xx))/(EXP(xx)+EXP(-xx))
END tanh 

REAL FUNCTION fx( REAL xx)
  REAL xs
  xs=0.
  IF protra = 1 THEN
    !   Lamina
    RESULT=1
    IF xx > xs  THEN RESULT=0
  ELSE IF protra = 2 THEN 
    !   e5
    RESULT=(1-xx∗4)**5
    IF xx>0.25 THEN RESULT= 0
  ELSE IF protra = 3 THEN 
    !   v
    RESULT=1-xx∗4
    IF xx>0.25 THEN RESULT= 0
  ELSE IF protra = 4 THEN 
    !   coseno
    RESULT=COS(2∗pi∗xx)
    IF xx>0.25 THEN RESULT= 0
  END IF 
  IF RESULT <0 THEN RESULT=0
END fx

REAL FUNCTION gz(REAL zz)
  REAL l1,l2,l3,l4,zr,a1,a2,a3,b1,b2,y1
  nr=4
  IF prolo =1 THEN
    !  2d
    RESULT=rapp
  ELSE IF prolo=2 THEN
    !  Profilo iperbolico
    RESULT=rapp∗(1.+tanh(ka∗2./rapp∗(fill∗lz/2.-zz)))/2
  ELSE IF prolo=3 THEN
    
    !  BE 
    
    tal=1
    cal=1/(1+tal**2)**0.5
    sal=tal∗cal
    r0=0.347826
    !   r0=0.7
    !     l0=0.5108695 
    l0=fill
    l2=(rapp-r0+tal∗l0+r0∗(1+tal**2)**0.5)/tal
    l1=l2-(rapp-r0+r0/(1+tal**2)**0.5)/tal
    y1=rapp-r0+r0∗cal
    l1=r0∗sal+l0
    l2=l1+y1/tal
    !     WRITE 'z,y1,l0,l1,l2',zz,y1,l0,l1,l2 
    l3=l1+(rapp-2∗r0∗(1-cal))/tal
    l4=l3+r0∗sal
    b1=(l2-fill+0.125-(nr+1)∗rapp)/(l2-fill+0.125)**nr
    b2=(nr∗rapp-l2+fill-0.125)/(l2-fill+0.125)**(nr+1)
    
    IF zz <l0 THEN
      RESULT =rapp
    ELSE IF zz <l1 THEN 
      RESULT=rapp-r0+(r0**2-(zz-l0)**2)**0.5
    ELSE IF zz<l2  THEN 
      RESULT = y1+tal∗(l1-zz)
    ELSE IF zz<lz-l2  THEN 
      RESULT=0
    ELSE IF zz<lz-l1  THEN 
      RESULT = y1+tal∗(zz-lz+l1)
    ELSE IF zz<lz-l0  THEN 
      RESULT=rapp-r0+(r0**2-(zz-lz+l0)**2)**0.5
    ELSE IF zz <lz  THEN
      RESULT =rapp
    END IF 
    
    IF RESULT<0  THEN RESULT=0
  END IF 
  
  !(
    ! Profilo arrotondato superiormente e con retta
    IF zz< l0 THEN
      RESULT=rapp
    ELSE IF zz<l1 THEN
      zr=(zz-l0)/(l1-l0)
      a1=(rapp-y1)∗nr/(l1-l0) -tal
      a2=nr∗(nr-1)∗(rapp-y1)/2/(l1-l0)∗∗2 -a1∗nr/(l1-l0)
      RESULT= rapp + zr∗∗nr∗(y1-rapp+a1∗(zz-l1)+ a2∗(zz-l1)∗∗2)
    ELSE IF zz < l2 THENJ    RESULT=tal∗(l2-zz)
  ELSE IF zz> l2 THEN
    RESULT=0
  END IF   
  !) 

!(    profilo continuo
  
  IF zz<= fill-.125  THEN RESULT=rapp
  IF zz>fill-.125  THEN 
    zr=(zz-fill+0.125)
    RESULT=rapp+b1∗zr∗∗nr+b2∗zr∗∗(nr+1)
  END IF
  !)
!(   Raccordo inferiore o
  IF zz>l1 THEN
    a3=(-y1+l2-l1-a2∗(l2-l1)∗∗2)/(l2-l1)∗∗3
    RESULT=y1-(zz-l1)+a2∗(zz-l1)∗∗2+a3∗(zz-l1)∗∗3
  END IF
  !) 
!     IF zz>l3 THEN RESULT=r0-(r0∗∗2-(zz-l4)∗∗2)∗∗.5
!     IF zz>l4 THEN RESULT =0

END    gz


REAL FUNCTION fondo(REAL xx,zz)
  RESULT=fx(xx)∗gz(zz)+fx(0.5-xx)∗gz(lz/2.-zz)
END fondo



REAL FUNCTION vert(REAL yy,x1,x2,zz)
  REAL f1,f2,xm,fm,x1c,x2c
  f1=fondo(x1,zz)-yy
  f2=fondo(x2,zz)-yy
  IF f1∗f2 > 0 THEN
    IF ABS(f1)  <  ABS(f2)  THEN
      RESULT=x1
    ELSE
      RESULT=x2
    END IF
    EXIT vert
  END IF
  x1c=x1
  x2c=x2
  DO
    xm=(f2∗x1c-f1∗x2c)/(f2-f1)   
    fm=fondo(xm,zz)-yy
    IF fm∗f1 > 0 THEN
      x1c=xm
      f1=fm
    ELSE
      x2c=xm
      f2=fm
    END IF
  UNTIL  (ABS(fm) < 0.0001) OR (ABS(x2c-x1c)) < 0.00005 
  RESULT=xm
END vert

REAL FUNCTION vertz(REAL  yy,z1,z2,xx)
  REAL f1,f2,zm,fm,z1c,z2c
  f1=fondo(xx,z1)-yy
  f2=fondo(xx,z2)-yy
  IF(f1∗f2 > 0) THEN
    IF(ABS(f1) < ABS(f2)) THEN
      RESULT =z1
    ELSE
      RESULT =z2
    END IF
    EXIT vertz
  END IF
  z1c=z1
  z2c=z2
  DO
    zm=(f2∗z1c-f1∗z2c)/(f2-f1)   
    fm=fondo(xx,zm)-yy
    IF (fm∗f1 > 0) THEN
      z1c=zm
      f1=fm
    ELSE
      z2c=zm
      f2=fm
    END IF
  UNTIL (ABS(fm) < 0.0001) OR (ABS(z2c-z1c)) <0.00005 
  RESULT=zm
END vertz



REAL FUNCTION delta(REAL xx,zz,calf,salf)
  REAL dz,dzp,xv,talfa,fp,gp,pv
  IF xx <0.25 THEN
    xv=0
    pv=1
  ELSE 
    xv=0.5
    pv=-1
  END IF
  IF calf = 0   THEN
    talfa=100000000000000.
  ELSE
    talfa=salf/calf
  END IF
  fp=(fondo(xv,zz)-fondo(xv+pv∗eps,zz))/eps
  IF xx=0 THEN
    gp=0
  ELSE
    gp=(fondo(xv,zz+eps)-fondo(xv,zz))/eps
  END IF
  IF rapp=0 THEN
    RESULT=pi/2
  ELSE
    RESULT=ATAN(fx0/rapp∗(1-gp∗talfa)**0.5)
  END IF
  !(   caso lamina     
    RESULT=0.
    !) 
END delta


REAL FUNCTION kappa(REAL xx, zz,calf,salf) 
  REAL cot,s,e,kv,kn
  IF ABS(delta(xx,zz,calf,salf))  = 0 THEN
    RESULT=0.5
  ELSE
    kn=0.5
    s=SIN(2.∗delta(xx,zz,calf,salf))
    DO
      kv=kn
      IF wp = 1 THEN
        cot=((1.-s**2∗kv**2)**0.5-1.)/s/kv 
        kn=(ATAN(1./cot)+pi)/(pi-delta(xx,zz,calf,salf))
      ELSE
        cot=(1-(1.-s**2∗kv**2)**0.5)/s/kv 
        kn=ATAN(1./cot)/(pi-delta(xx,zz,calf,salf))
      END IF
    WHILE  ABS(kv-kn) > 0.0001 
    RESULT=kn
  END IF
END kappa


SUBROUTINE catalogo
  
  REAL xx,yy,zz,dis,zve,znu,zno,z1,z2,d1,d2,dm,xv,yv,calf,salf,dpm,z3,norma
  INTEGER kv,kvm,kid
  kid=0
  LOOP FOR  ica=0 TO nid.HI
    LOOP FOR  jca=0 TO nid(0).HI
      LOOP FOR  kca=0 TO nid(0,0).HI
        nid(ica,jca,kca)=0
      REPEAT LOOP
    REPEAT LOOP
  REPEAT LOOP
  LOOP FOR  ica=0 TO nx
    LOOP FOR  kca=0 TO  nz
      
      IF ica < 2   OR  ica > nx-2   THEN
        ! ∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗ btm(0,0) max btm ∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗
        
        LOOP FOR  jca=btm(ica,kca) TO  btm(0,0)+2
          
          xx=ica∗h
          yy=jca∗h
          zz=kca∗h
          IF ica∗h < 0.25 THEN
            xv=0.
          ELSE
            xv=0.5
          END IF
          IF kca > 0   THEN
            z1=0.
            z2=2.∗nz∗h
            d1=(xx-xv)**2+(yy-fondo(xv,z1))**2+(zz-z1)**2
            d2=(xx-xv)**2+(yy-fondo(xv,z2))**2+(zz-z2)**2
            DO
              !   WRITE 'z1,z2',z1,z2
              zno=0.5∗z1+0.5∗z2
              dm=(xx-xv)**2+(yy-fondo(xv,zno))**2+(zz-zno)**2
              dpm=(xx-xv)**2+(yy-fondo(xv,zno+eps))**2+(zz-zno-eps)**2-
      ((xx-xv)**2+(yy-fondo(xv,zno-eps))**2+(zz-zno+eps)**2)
              !   WRITE 'xx,yy,zz,d1,d2,dm,zno,dpm',xx,yy,zz,d1,d2,dm,zno,dpm
              IF dpm > 0 THEN
                z2=zno 
                d2=dm
              ELSE
                z1=zno 
                d1=dm
              END IF
              !  WRITE 'z1,z2,d1,d2',z1,z2,d1,d2
            WHILE  ABS(z1-z2) > h/1000000
          ELSE
            
            zno=0.
            calf=1.0
            salf=0.
            dis=((xx-xv)**2+(yy-fondo(xv,0.))**2)**0.5
          END IF
          yv=fondo(xv,zno)
          dis=((xx-xv)**2+(yy-yv)**2+(zz-zno)**2)**0.5
          
          IF yv > 1.05∗h AND jca>1 THEN 
            !    IF  dis < dcor    THEN
            IF  dis < dcor    THEN
              norma=((yy-yv)**2+(zz-zno)**2)**0.5
              IF yy < yv THEN norma=-norma
              IF ABS(norma) < 0.001 THEN
                calf=1.
                salf=0.
              ELSE
                calf=ABS((yy-yv)/norma)
                salf=(zz-zno)/norma
              END IF
              z3=(zz-zno)∗calf-(yy-yv)∗salf 
              !   IF ABS(z3) <  0.9∗h  THEN  
              kid=kid+1
              vzn(kid)=zno
              vsalfa(kid)=salf
              vcalfa(kid)=calf
              nid(ica,jca,kca)=kid
              vdel(kid)=delta(xx,zno,calf,salf)
              vk1(kid)=kappa(xx,zno,calf,salf)
              vk2(kid)=pi/(pi-vdel(kid))∗(5.+wp)/4.
              
              !   END IF
            END IF
          END IF
        REPEAT LOOP
      END IF
    REPEAT LOOP
  REPEAT LOOP
  WRITE 'punti da correggere',kid                     
END catalogo


REAL FUNCTION derc(REAL fc,fs,fd,ds,dd,den)
  RESULT=(fd∗ds**2-fs∗dd**2+fc∗(dd**2-ds**2))/den/h
END derc

REAL FUNCTION ders(REAL f0,f1,f2,d1,d2,den)
  RESULT=(f1∗d2**2-f2∗d1**2-f0∗(d2**2-d1**2))/den/h
END ders

REAL FUNCTION des(REAL fl,fc,fr,frr,le,ri,rr)
  REAL den,alr,alrr
  alr=le+ri
  alrr=alr+rr
  IF le< pfo AND rr>0  THEN
    den =alr∗alrr∗rr
    RESULT=(fr∗alrr**2-frr∗alr**2-fl∗(alrr**2-alr**2))/h/den
  ELSE
    den=le∗ri∗alr
    RESULT=(fc∗alr**2-fr∗le**2-fl∗(alr**2-le**2))/h/den
  END IF
END des



! ∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗



!     funzioni locali

REAL FUNCTION teta(REAL xx,yy)
  REAL tet,rq
  INTEGER si
  IF  xx = 0 THEN si=0
  IF xx < 0  THEN si=-1
  IF xx > 0  THEN si=1
  rq=xx∗xx+yy∗yy     
  IF  yy = 0   THEN
    tet=pi/2∗si
  ELSE
    tet=ATAN(xx/yy)
    IF yy < 0 THEN tet=tet+pi∗si
  END IF
  RESULT=tet
END teta


SUBROUTINE fun1(REAL fi1^,m11^,m21^,xx,yy)
  REAL tet,e,a,d,b,fil,fis,k3,es,cot,s,c,r
  r=(xx∗xx+yy∗yy)**0.5
  tet=teta(xx,yy)
  IF del > 0 THEN
    s=SIN(2∗del)
    c=COS(2∗del)
    k3=(pi-del)∗k1
    IF ABS(pi/2-k3)  < 0.005  THEN
      cot=0
    ELSE
      cot=1/TAN(k3)
    END IF
    IF wp = 1  THEN
      e=1-2∗cot/(cot+TAN(del))
      !    d=1./(1.-e∗(1.+k1)-k1∗(s∗cot+c)) 
      d=cot/(cot∗(1+e∗(1+k1)+k1∗c)-s∗k1)
      b=1
      a=4∗d-1
    ELSE
      e=2.∗TAN(del)∗cot/(1.-cot∗TAN(del))+1.
      d=cot/(cot∗(1-k1∗c+e∗(1+k1))+k1∗s)
      a=1
      b=1-4∗d
    END IF
    es=1.+k1
    IF (wp = 1) THEN
      fis=r**es∗COS((1.-k1)∗tet)
      fil=r**es∗COS(es∗tet)
      fi1=d∗(fis+e∗fil)
      m11=a∗r**k1∗SIN(k1∗tet)
      m21=b∗r**k1∗COS(k1∗tet)
    ELSE
      fis=r**es∗SIN((1.-k1)∗tet)
      fil=r**es∗SIN(es∗tet)
      fi1=d∗(fis+e∗fil)
      m11=a∗r**k1∗COS(k1∗tet)
      m21=b∗r**k1∗SIN(k1∗tet)
    END IF
  ELSE
    ! Lamina
    IF wp = -1 THEN
      m11=r**0.5∗COS(tet/2.)
      m21=0.
      fis=r**1.5∗SIN(tet∗0.5)
      fil=r**1.5∗SIN(tet∗1.5)
      fi1=0.25∗(fis+fil)
    ELSE
      !    m11=0.
      m11=-r**0.5∗SIN(0.5∗tet)/2
      m21=r**0.5∗COS(0.5∗tet)
      fis=r**1.5∗COS(tet∗0.5)
      fil=r**1.5∗COS(tet∗1.5)
      !    fi1=.25∗(fis+fil/3)
      fi1=(fis+fil∗3)/8
    END IF
  END IF
END  fun1 



SUBROUTINE fun2(REAL fi2^,m12^,m22^,m32^,xx,yy)
  REAL tet,es,k23,r
  k23=pi/2./(pi-del)
  r=(xx∗xx+yy∗yy)**0.5
  tet=teta(xx,yy)
  IF del > 0 THEN
    es=k2-1.
    IF (wp = 1) THEN
      fi2=-r**k2∗COS(k2∗tet)/k2 
      m12=r**es∗SIN(es∗tet)
      m22=-r**es∗COS(es∗tet)
      m32=r**k23∗COS(k23∗tet)
    ELSE
      fi2=r**k2∗SIN(k2∗tet)/k2 
      m12=r**es∗COS(es∗tet)
      m22=r**es∗SIN(es∗tet)
      m32=r**k2∗SIN(k2∗tet)
    END IF
    ! Caso lamina
  ELSE
    IF  wp = -1  THEN
      fi2=xx
      m12=1.
      m22=xx
      m32=xx
    ELSE
      fi2=-2./3.∗r**1.5∗COS(1.5∗tet)
      m12=r**0.5∗SIN(0.5∗tet)
      m22=-r**0.5∗COS(0.5∗tet)
      m32=r**0.5∗COS(0.5∗tet)
    END IF
  END IF
END fun2




SUBROUTINE costanti( INTEGER icor,jcor,kcor)
  INTEGER ne,ive,kve,ke,n6,k6
  neq=120
  REAL xv,yv,den1,den2,denc6,dc6
  REAL m1a,m2a,m3a,den,xc,yc,zc,lun,er
  REAL fia,fib,fic,m11,m12,m21,m22,m32,fi1,fi2
  REAL b1,b2,b3,ai(1..neq,1..2),a6(1..neq),bi(1..neq)
  REAL b6(1..neq),bl(1..neq),ail(1..neq),a11,a12,a21,a22
  REAL f1l,f2l,fl,fip,fip1,fip2,xp,yp,erv
  IF icor∗h < 0.25  THEN
    xv=0
    ive=0
  ELSE
    xv=0.5
    ive=nx
  END IF
  yv=fondo(xv,zn)
  jve=btm(ive,kcor)
  LOOP FOR  iv=1 TO  neq
    bi(iv)=0.
    LOOP FOR  jv=1 TO 2
      ai(iv,jv)=0.
    REPEAT LOOP
  REPEAT LOOP
  LOOP FOR  iv=1 TO neq
    a6(iv)=0.
    ail(iv)=0.
    b6(iv)=0.
    bl(iv)=0.
  REPEAT LOOP
  ke=0
  k6=0
  
  LOOP FOR  iv=-1 TO 1
    LOOP FOR  kv=-1 TO  1
      LOOP FOR  jv=-1 TO 1
        
        IF  icor+iv > -1 AND icor+iv < nx+1     THEN
          IF   kcor+kv < nz+1 AND  kcor+kv > -1    THEN
            IF   jcor+jv > btm(icor+iv,kcor+kv) -1 THEN
              xc=(icor+iv)∗h-xv
              yc=((jcor+jv)∗h-yv)∗calfa+(h∗kcor+kv∗h-zn)∗salfa
              zc=-((jcor+jv)∗h-yv)∗salfa+(h∗kcor+kv∗h-zn)∗calfa
              lun=1
              m1a=m1(icor+iv,jcor+jv,kcor+kv)
              m2a=m2(icor+iv,jcor+jv,kcor+kv)
              m3a=m3(icor+iv,jcor+jv,kcor+kv)
              fic=fi(icor+iv,jcor+jv,kcor+kv)
              fun1(fi1,m11,m21,xc,yc)
              fun2(fi2,m12,m22,m32,xc,yc)
              ke=ke+1
              k6=k6+1
              ai(ke,1)=m11
              ai(ke,2)=m12
              bi(ke)=m1a
              ail(k6)=m22
              a6(k6)=m32
              b6(k6)=m3a∗calfa-m2a∗salfa
              bl(k6)=m2a∗calfa+m3a∗salfa
              ke=ke+1
              ai(ke,1)=fi1/lun
              ai(ke,2)=fi2/lun
              bi(ke)=fic/lun
              ke=ke+1
              ai(ke,1)=m21
              ai(ke,2)=m22
              bi(ke)=m2a∗calfa+m3a∗salfa
            END IF
          END IF
        END IF
      REPEAT LOOP
    REPEAT LOOP
  REPEAT LOOP
  ne=ke
  n6=k6
  a11=0.
  a12=0.
  a21=0.
  a22=0.
  b1=0.
  b2=0.
  c6=0.
  denc6=0.
  dc6=0.
  LOOP FOR  ke=1 TO n6
    denc6=denc6+a6(ke)**2
    dc6=dc6+a6(ke)∗b6(ke)
  REPEAT LOOP
  c6=dc6/denc6
  
  denc6=0.
  dc6=0.
  LOOP FOR  ke=1 TO n6
    denc6=denc6+ail(ke)**2
    dc6=dc6+ail(ke)∗bl(ke)
  REPEAT LOOP
  c2l=dc6/denc6
  
  LOOP FOR  iv=1 TO  ne
    a11=ai(iv,1)∗ai(iv,1)+a11
    a12=ai(iv,2)∗ai(iv,1)+a12
    a21=ai(iv,1)∗ai(iv,2)+a21
    a22=ai(iv,2)∗ai(iv,2)+a22
    b1=bi(iv)∗ai(iv,1)+b1
    b2=bi(iv)∗ai(iv,2)+b2
  REPEAT LOOP
  den=a11∗a22-a12∗a21
  den1=b1∗a22-b2∗a12
  den2=a11∗b2-b1∗a21
  c2=den1/den
  c4=den2/den
END costanti


REAL FUNCTION m1loc(REAL xx,yy)
  REAL fi1,fi2,m11,m12,m21,m22,m32
  fun1(fi1,m11,m21,xx,yy)
  fun2(fi2,m12,m22,m32,xx,yy)
  RESULT=c2∗m11+c4∗m12
END m1loc



REAL FUNCTION m2loc(REAL xx,yy)
  REAL fi1,fi2,m11,m12,m21,m22,m32
  fun1(fi1,m11,m21,xx,yy)
  fun2(fi2,m12,m22,m32,xx,yy)
  RESULT=c2∗m21+c4∗m22
END m2loc



REAL FUNCTION m3loc(REAL xx,yy)
  REAL fi2,m12,m22,m32
  fun2(fi2,m12,m22,m32,xx,yy)
  RESULT=c6∗m32
END m3loc



REAL FUNCTION filoc(REAL xx,yy)
  REAL fi1,fi2,m11,m12,m21,m22,m32
  fun1(fi1,m11,m21,xx,yy)
  fun2(fi2,m12,m22,m32,xx,yy)
  RESULT=c2∗fi1+c4∗fi2
END filoc


REAL FUNCTION fico(INTEGER ic,jc,kc)
  REAL xx,yy,xv
  IF ic∗h < 0.25 THEN
    xv=0
  ELSE
    xv=0.5
  END IF
  xx=ic∗h-xv
  yy=(jc∗h-fondo(xv,zn))∗calfa+(kc∗h-zn)∗salfa
  RESULT=filoc(xx,yy)
END fico




REAL FUNCTION m1co(INTEGER ic,jc,kc)
  REAL xx,yy,xv
  IF ic∗h < 0.25 THEN
    xv=0
  ELSE
    xv=0.5
  END IF
  xx=ic∗h-xv
  yy=(jc∗h-fondo(xv,zn))∗calfa+(kc∗h-zn)∗salfa
  RESULT=m1loc(xx,yy)
END m1co




REAL FUNCTION m2co(INTEGER ic,jc,kc)
  REAL xx,yy,xv
  IF ic∗h < 0.25 THEN
    xv=0
  ELSE
    xv=0.5
  END IF
  xx=ic∗h-xv
  yy=(jc∗h-fondo(xv,zn))∗calfa+(kc∗h-zn)∗salfa
  RESULT=m2loc(xx,yy)∗calfa-m3loc(xx,yy)∗salfa
END m2co




REAL FUNCTION m3co( INTEGER ic,jc,kc)
  REAL xx,yy,xv
  IF ic∗h < 0.25 THEN
    xv=0
  ELSE
    xv=0.5
  END IF
  xx=ic∗h-xv
  yy=(jc∗h-fondo(xv,zn))∗calfa+(kc∗h-zn)∗salfa
  RESULT=m2loc(xx,yy)∗salfa+m3loc(xx,yy)∗calfa
  
END m3co




REAL FUNCTION fxd(REAL tang; INTEGER i,j,k,tra,co)
  INTEGER ile,ifu,kfu
  REAL fr,fd,fu,fuu,ld
  IF tang>0 THEN
    ile=1
  ELSE
    ile=-1
  END IF
  ifu=i+tra∗ile
  kfu=k+(1-tra)∗ile
  ld=ad(ifu,j,kfu)-ad(i,j,k)
  IF co=0 THEN
    fr=fi(ifu,j,kfu)
    IF btm(ifu,kfu)> j THEN fr=0
    fd=fi(ifu,j-1,kfu)
    IF btm(ifu,kfu)> j-1  THEN fd=0
    fu=fi(ifu,j+1,kfu)
    fuu=fi(ifu,j+2,kfu)
  ELSE
    fr=fico(ifu,j,kfu)
    IF btm(ifu,kfu)> j THEN fr=0
    fd=fico(ifu,j-1,kfu)
    IF btm(ifu,kfu)> j-1  THEN fd=0
    fu=fico(ifu,j+1,kfu)
    fuu=fico(ifu,j+2,kfu)
  END IF
  RESULT =(fd/h+des(fd,fr,fu,fuu,ad(ifu,j,kfu),1,1)∗ld)∗ile
END fxd



REAL FUNCTION fyl( INTEGER i,j,k,tra,co,leftside)
  INTEGER ile,ifu,kfu
  REAL fl,fc,fr,frr,ld,lrd
  ile=leftside
  IF ile=1 THEN
    ld=al(i,j+1,k)
    lrd=ld-al(i,j,k)
  ELSE
    ld=ar(i,j+1,k)
    lrd=ld-ar(i,j,k)
  END IF
  ifu=i-tra∗ile
  kfu=k-(1-tra)∗ile
  IF co=0 THEN
    fr=fi(i+tra∗ile,j+1,k+(1-tra)∗ile)
    IF btm(i+tra∗ile,k+(1-tra∗ile))> j+1 THEN fr=0
    frr=fi(i+2∗tra∗ile,j+1,k+2∗(1-tra)∗ile)
    IF btm(i+2∗tra∗ile,k+2∗(1-tra∗ile))> j+1 THEN frr=0
    fl=fi(ifu,j+1,kfu)
    IF btm(ifu,kfu)> j+1  THEN fl=0
    fc=fi(i,j+1,k)
  ELSE
    fr=fico(i+tra∗ile,j+1,k+(1-tra)∗ile)
    IF btm(i+tra∗ile,k+(1-tra∗ile))> j+1 THEN fr=0
    frr=fico(i+2∗tra∗ile,j+1,k+2∗(1-tra)∗ile)
    IF btm(i+2∗tra∗ile,k+2∗(1-tra∗ile))> j+1 THEN frr=0
    fl=fico(ifu,j+1,kfu)
    IF btm(ifu,kfu)> j+1  THEN fl=0
    fc=fico(i,j+1,k)
  END IF
  RESULT =fl/h+des(fl,fc,fr,frr,ld,1,1)∗lrd
END fyl

SUBROUTINE notom1
  REAL den,cl,cr,cd,cf,cb,ml,mr,md,mf,mb,fix,alr,rr,ll,
  tang,ficen,fl,fll,fr,frr,fb,fff,fd,fu,flad
  LOOP FOR  ino=0 TO nx
    LOOP FOR  kno=0 TO nz
      LOOP FOR  jno=btm(ino,kno) TO  btf(ino,kno)-1
        alr=al(ino,jno,kno)+ar(ino,jno,kno)
        den=al(ino,jno,kno)∗ar(ino,jno,kno)∗alr
        cl=0.
        cd=0.
        cr=0.
        cb=0.
        cf=0.
        mr=0.
        md=0.
        ml=0.
        mf=0.
        mb=0.
        fu=fi(ino,jno+1,kno)
        ficen=fi(ino,jno,kno)
        
        IF btm(ino-1,kno) > jno  THEN 
          fl=0
          rr=0 
        ELSE
          fl=fi(ino-1,jno,kno)
          rr=al(ino-1,jno,kno)
        END IF
        
        
        IF ino=0 THEN
          fll=wp∗fi(ino+2,jno,kno)
          IF btm(ino+2,kno)> jno THEN fll=0
        ELSE
          fll=fi(ino-2,jno,kno)
          IF btm(ino-2,kno)> jno THEN fll=0
        END IF
        
        IF btm(ino+1,kno) > jno  THEN 
          fr=0
          ll=0
        ELSE
          fr=fi(ino+1,jno,kno)
          ll=ar(ino+1,jno,kno)
        END IF
        
        IF ino=nx THEN
          frr=wp∗fi(ino-2,jno,kno)
          IF btm(ino-2,kno)> jno THEN frr=0
        ELSE
          frr=fi(ino+2,jno,kno)
          IF btm(ino+2,kno)> jno THEN frr=0
        END IF
        
        IF btm(ino,kno-1) > jno  THEN 
          fb=0
        ELSE
          fb=fi(ino,jno,kno-1)
        END IF
        
        IF btm(ino,kno+1) > jno  THEN 
          fff=0
        ELSE
          fff=fi(ino,jno,kno+1)
        END IF
        
        IF  al(ino,jno,kno) < pfo   THEN
          fix=ders(ficen,fr,frr,1,2,2)
        ELSE IF  ar(ino,jno,kno) < pfo   THEN
          fix=-ders(ficen,fl,fll,1,2,2)
        ELSE
          fix=derc(ficen,fl,fr,al(ino,jno,kno),ar(ino,jno,kno),den) 
        END IF
        
        IF btm(ino-1,kno) > jno THEN
          cl=bcl(ino,jno,kno)
          IF  al(ino,jno,kno) < pfo   THEN
            ml=ders(0,fr,frr,alr,alr+1,alr∗(alr+1))
          ELSE
            ml=ders(0,ficen,fr,al(ino,jno,kno),alr,den)
          END IF
          !    ml=des(fl,ficen,fr,frr,al(ino,jno,kno),ar(ino,jno,kno),rr)
        END IF
        IF  btm(ino+1,kno) > jno  THEN
          cr=bcr(ino,jno,kno)
          IF  ar(ino,jno,kno) < pfo   THEN
            mr=-ders(0,fl,fll,alr,alr+1,alr∗(alr+1))
          ELSE
            mr=-ders(0,ficen,fl,ar(ino,jno,kno),alr,den)
          END IF
          
        END IF
        
        ! calcolo md
        IF jno = btm(ino,kno)  THEN
          cd=bcd(ino,jno,kno)
          tang=(fondo(ino∗h-eps,kno∗h)-fondo(ino∗h+eps,h∗kno))/2/eps
          IF pbc=1 THEN
            md=fix
          ELSE IF pbc=2 THEN
            IF ino=0 OR ino= nx THEN
              md=fix
            ELSE
              IF tang=0 THEN 
                md=0
              ELSE
                md=fxd(tang,ino,jno,kno,1,0)    
              END IF
            END IF
          END IF
        END IF
        ! fine calcolo md
        
        IF btm(ino,kno-1) > jno  THEN
          cb=bcb(ino,jno,kno)
          
          mb=fix        
          
        END IF
        
        IF btm(ino,kno+1) > jno  THEN
          cf=bcf(ino,jno,kno)
          mf=fix        
        END IF
        
        rhs(ino,jno,kno)=cl∗ml+cr∗mr+cd∗md+cf∗mf+cb∗mb
      REPEAT LOOP
      LOOP FOR  jno=btf(ino,kno) TO  ny-1
        rhs(ino,jno,kno)=0.
      REPEAT LOOP
      rhs(ino,ny,kno)=h/3∗(1-wp)/2
    REPEAT LOOP
  REPEAT LOOP
END notom1




SUBROUTINE notom2
  REAL cl,cr,cf,cd,cb,ml,mr,md,mb,mf,den,fiy,fd,ad1,tang
  REAL fir,ficen,fu,fl,fr,fb,fff
  LOOP FOR  ino=0 TO nx
    LOOP FOR  kno=0 TO nz
      LOOP FOR  jno=btm(ino,kno) TO  btf(ino,kno)-1
        cl=0.
        cd=0.
        cr=0.
        cf=0.
        cb=0.
        ml=0.
        md=0.
        mr=0.
        mf=0.
        mb=0.
        fu=fi(ino,jno+1,kno)
        ficen=fi(ino,jno,kno)
        
        IF btm(ino-1,kno) > jno  THEN 
          fl=0
        ELSE
          fl=fi(ino-1,jno,kno)
        END IF
        
        IF btm(ino+1,kno) > jno  THEN 
          fr=0
        ELSE
          fr=fi(ino+1,jno,kno)
        END IF
        
        IF btm(ino,kno-1) > jno  THEN 
          fb=0
        ELSE
          fb=fi(ino,jno,kno-1)
        END IF
        
        IF btm(ino,kno+1) > jno  THEN 
          fff=0
        ELSE
          fff=fi(ino,jno,kno+1)
        END IF
        
        ad1=ad(ino,jno,kno)+1.
        den=ad(ino,jno,kno)∗ad1
        IF ad(ino,jno,kno) < pfo THEN
          fiy=ders(ficen,fu,fi(ino,jno+2,kno),1,2,2)
        ELSE
          fd=fi(ino,jno-1,kno)
          IF jno = btm(ino,kno) THEN fd=0.
          fiy=derc(ficen,fd,fu,ad(ino,jno,kno),1,den)
        END IF
        tang=(fondo(ino∗h-eps,kno∗h)-fondo(ino∗h+eps,h∗kno))/2/eps
        IF  btm(ino-1,kno) > jno  THEN
          cl=bcl(ino,jno,kno)
          IF pbc=1 THEN
            ml =fiy
          ELSE IF pbc=2 THEN
            IF ino =1 AND jno = btm(ino-1,kno) THEN
              !   ml =fiy
              ml=0
              
            ELSE
              ml =fyl(ino,jno,kno,1,0,1)
            END IF
          END IF
        END IF
        IF  btm(ino+1,kno) > jno  THEN
          cr=bcr(ino,jno,kno)
          IF pbc=1 THEN
            mr =fiy
          ELSE IF pbc=2 THEN
            IF ino =nx-1 AND jno = btm(ino+1,kno) THEN
              !   mr =fiy
              mr=0
            ELSE
              mr =fyl(ino,jno,kno,1,0,-1)
            END IF
          END IF
        END IF
        
        IF  jno = btm(ino,kno)  THEN
          cd=bcd(ino,jno,kno)
          IF( ino=0 OR ino=nx) AND jno >1 THEN
            md=0
          ELSE 
            IF ad(ino,jno,kno) < pfo THEN
              md=ders(0,fu,fi(ino,jno+2,kno),ad1,ad1+1,ad1∗(ad1+1))
            ELSE
              md=ders(0,ficen,fu,ad(ino,jno,kno),ad1,den)
            END IF
          END IF
        END IF
        
        IF btm(ino,kno-1) > jno  THEN
          cb=bcb(ino,jno,kno)
          IF pbc=1 THEN
            mb=fiy
          ELSE IF pbc=2 THEN
            !  mb=fiy
            mb =fyl(ino,jno,kno,0,0,1)
          END IF
        END IF
        
        IF btm(ino,kno+1) > jno  THEN
          cf=bcf(ino,jno,kno)
          IF pbc=1 THEN
            mf=fiy
          ELSE  IF pbc=2 THEN
            !  mf=fiy
            mf =fyl(ino,jno,kno,0,0,-1)
          END IF
        END IF
        
        rhs(ino,jno,kno)=cl∗ml+cr∗mr+cd∗md+mb∗cb+mf∗cf
      REPEAT LOOP
      LOOP FOR  jno=btf(ino,kno) TO  ny-1
        rhs(ino,jno,kno)=0.
      REPEAT LOOP
      rhs(ino,ny,kno)=0
    REPEAT LOOP
  REPEAT LOOP
END notom2




SUBROUTINE notom3
  REAL den,fiz,cf,cb,cl,cr,cd,ml,mr,md,mf,mb,abf,fiq,fit,
  tangz,fu,ficen,fl,fr,fb,fff,ffff,fbb
  LOOP FOR  ino=0 TO nx
    LOOP FOR  kno=0 TO nz
      LOOP FOR  jno=btm(ino,kno) TO btf(ino,kno)-1
        abf=ab(ino,jno,kno)+af(ino,jno,kno)
        den=ab(ino,jno,kno)∗af(ino,jno,kno)∗abf
        cb=0.
        cf=0.
        cl=0.
        cr=0.
        mf=0.
        mb=0.
        cd=0.
        md=0.
        ml=0.
        mr=0.
        fu=fi(ino,jno+1,kno)
        ficen=fi(ino,jno,kno)
        
        IF btm(ino-1,kno) > jno  THEN 
          fl=0
        ELSE
          fl=fi(ino-1,jno,kno)
        END IF
        
        IF btm(ino+1,kno) > jno  THEN 
          fr=0
        ELSE
          fr=fi(ino+1,jno,kno)
        END IF
        
        IF btm(ino,kno-1) > jno  THEN 
          fb=0
        ELSE
          fb=fi(ino,jno,kno-1)
        END IF
        
        IF btm(ino,kno+1) > jno  THEN 
          fff=0
        ELSE
          fff=fi(ino,jno,kno+1)
        END IF
        
        ad1=ad(ino,jno,kno)+1.
        IF  ab(ino,jno,kno) < pfo  AND  kno <  nz  THEN
          fiz=ders(ficen,fff,fi(ino,jno,kno+2),1,2,2)
        ELSE IF  af(ino,jno,kno) < pfo AND  kno > 0  THEN
          fiz=-ders(ficen,fb,fi(ino,jno,kno-2),1,2,2)
        ELSE
          fiz=derc(ficen,fb,fff,ab(ino,jno,kno),af(ino,jno,kno),den)
        END IF
        IF (ab(ino,jno,kno) < pfo AND kno=nz) OR (af(ino,jno,kno) < pfo AND kno =0) THEN \
        fiz=(fff- fb)/abf/h
        
        IF btm(ino,kno-1)   > jno THEN
          cb=bcb(ino,jno,kno)
          IF  ab(ino,jno,kno) < pfo AND  kno < nz THEN
            mb=ders(0,fff,fi(ino,jno,kno+2),abf,abf+1,abf∗(abf+1))
          ELSE
            mb=ders(0,ficen,fff,ab(ino,jno,kno),abf,den)
          END IF
          IF    ab(ino,jno,kno) < pfo AND   kno=nz THEN mb=fiz 
          IF ino=0 OR ino=nx THEN mb=0
        END IF
        
        IF btm(ino,kno+1)   > jno THEN
          cf=bcf(ino,jno,kno)
          IF  af(ino,jno,kno) < pfo AND  kno > 0  THEN
            mf=-ders(0,fb,fi(ino,jno,kno-2),abf,abf+1,abf∗(abf+1))
          ELSE
            mf=-ders(0,ficen,fb,af(ino,jno,kno),abf,den)
          END IF
          IF ino=0 OR ino=nx THEN mf=0
        END IF
        IF    af(ino,jno,kno) < pfo AND   kno=0  THEN mf=fiz 
        
        
        IF  jno = btm(ino,kno)  THEN
          cd=bcd(ino,jno,kno)
          tangz=(fondo(ino∗h,kno∗h-eps)-fondo(ino∗h,h∗kno+eps))/2/eps
          IF pbc =1 THEN
            md=fiz
          ELSE IF pbc =2 THEN
            !  md=fiz
            md=fxd(tangz,ino,jno,kno,0,0)
            IF ino=0 OR ino=nx THEN md=0
            !(   IF ino=0 OR ino=nx THEN
                IF tangz>0 THEN
                  md=(tangz/ad(ino,jno,kno+1))∗∗1.5/h∗fff
                ELSE IF tangz<0 THEN
                  md=-(-tangz/ad(ino,jno,kno-1))∗∗1.5/h∗fb
                ELSE
                  md=0
                END IF
              END IF
              !)
          END IF
          IF tangz=0 THEN md=0
          IF kno=0 AND wp=-1 THEN md=0
        END IF
        ! fine calcolo md
        
        IF btm(ino-1,kno) > jno  THEN
          cl=bcl(ino,jno,kno)
          IF pbc =1 THEN
            ml=fiz
            !   ml=0
          ELSE IF pbc =2 THEN
            !  ml=fiz
            ml=0
          END IF
        END IF
        
        IF btm(ino+1,kno) > jno  THEN
          cr=bcr(ino,jno,kno)
          IF pbc =1 THEN
            mr=fiz
            !    mr=0
          ELSE IF pbc =2 THEN
            !  mr=fiz
            mr=0
          END IF
        END IF
        
        rhs(ino,jno,kno)=cd∗md+cb∗mb+cf∗mf+cl∗ml+cr∗mr
      REPEAT LOOP
      LOOP FOR  jno=btf(ino,kno) TO ny-1
        rhs(ino,jno,kno)=0.
      REPEAT LOOP
      rhs(ino,ny,kno)=h/3∗(wp+1)/2
    REPEAT LOOP
  REPEAT LOOP
END notom3





SUBROUTINE notofia
  REAL den1,den2,den3,m1x,m2y,m3z,alr,abf,fff,fb,md,ml,mr,mb,mf,ad1
  REAL ficen,fu,fl,fr,fll,frr
  LOOP FOR  ino=0 TO nx
    LOOP FOR  kno=0 TO nz
      LOOP FOR  jno=btm(ino,kno) TO btf(ino,kno)-1
        alr=al(ino,jno,kno)+ar(ino,jno,kno)
        abf=ab(ino,jno,kno)+af(ino,jno,kno)
        ad1=ad(ino,jno,kno)+1.
        den1=al(ino,jno,kno)∗ar(ino,jno,kno)∗alr
        den2=ad(ino,jno,kno)∗ad1
        den3=af(ino,jno,kno)∗ab(ino,jno,kno)∗abf
        ficen=fi(ino,jno,kno)
        fu=fi(ino,jno+1,kno)
        ml=m1(ino-1,jno,kno)
        mr=m1(ino+1,jno,kno)
        md=m2(ino,jno-1,kno)
        mb=m3(ino,jno,kno-1)
        mf=m3(ino,jno,kno+1)
        fb=fi(ino,jno,kno-1)
        fff=fi(ino,jno,kno+1)
        IF btm(ino,kno-1) > jno  THEN fb=0. 
        IF btm(ino,kno+1) > jno  THEN fff=0. 
        
        IF btm(ino-1,kno) > jno  THEN 
          fl=0
        ELSE
          fl=fi(ino-1,jno,kno)
        END IF
        
        
        IF ino=0 THEN
          fll=wp∗fi(ino+2,jno,kno)
          IF btm(ino+2,kno)> jno THEN fll=0
        ELSE
          fll=fi(ino-2,jno,kno)
          IF btm(ino-2,kno)> jno THEN fll=0
        END IF
        
        IF btm(ino+1,kno) > jno  THEN 
          fr=0
        ELSE
          fr=fi(ino+1,jno,kno)
        END IF
        
        IF ino=nx THEN
          frr=wp∗fi(ino-2,jno,kno)
          IF btm(ino-2,kno)> jno THEN frr=0
        ELSE
          frr=fi(ino+2,jno,kno)
          IF btm(ino+2,kno)> jno THEN frr=0
        END IF
        
        IF   jno = btm(ino,kno) THEN    
          IF ad(ino,jno,kno) < pfo THEN
            md=ders(0,fu,fi(ino,jno+2,kno),ad1,ad1+1,ad1∗(ad1+1))
          ELSE
            md=ders(0,ficen,fu,ad(ino,jno,kno),ad1,den2)
          END IF
          IF (ino=0 OR ino= nx) AND jno >1 THEN md=0
        END IF
        IF  btm(ino+1,kno) > jno  THEN
          IF  ar(ino,jno,kno) < pfo    THEN
            mr=-ders(0,fl,fll,alr,alr+1,alr∗(alr+1))
          ELSE
            mr=-ders(0,ficen,fl,ar(ino,jno,kno),alr,den1)
          END IF
        END IF
        IF  btm(ino-1,kno) > jno  THEN
          IF al(ino,jno,kno) < pfo   THEN
            ml=ders(0,fr,frr,alr,alr+1,alr∗(alr+1))
          ELSE
            ml=ders(0,ficen,fr,al(ino,jno,kno),alr,den1)
          END IF
        END IF
        IF  btm(ino,kno+1) > jno  THEN
          IF  af(ino,jno,kno) < pfo AND  kno > 0  THEN
            mf=-ders(0,fi(ino,jno,kno-1),fi(ino,jno,kno-2),abf,abf+1,abf∗(abf+1))
          ELSE
            mf=-ders(0,fi(ino,jno,kno),fi(ino,jno,kno-1),af(ino,jno,kno),abf,den3)
            IF af(ino,jno,kno) < pfo AND  kno=0  THEN mf=(fff-fb)/abf/h
          END IF
          IF ino=0 OR ino= nx  THEN mf=0
        END IF
        IF  btm(ino,kno-1) > jno  THEN
          IF  ab(ino,jno,kno) < pfo AND  kno <  nz  THEN
            mb=ders(0,fi(ino,jno,kno+1),fi(ino,jno,kno+2),abf,abf+1,abf∗(abf+1))
          ELSE
            mb=ders(0,fi(ino,jno,kno),fi(ino,jno,kno+1),ab(ino,jno,kno),abf,den3)
            IF  ab(ino,jno,kno) < pfo AND  kno =nz  THEN mb=(fff-fb)/abf/h
          END IF
          IF ino=0 OR ino= nx  THEN mb=0
        END IF
        IF  al(ino,jno,kno) < pfo AND  ino <  nx  THEN
          m1x=ders(m1(ino,jno,kno),mr,m1(ino+2,jno,kno),1,2,2)
        ELSE IF  ar(ino,jno,kno) < pfo AND  ino > 0  THEN
          m1x=-ders(m1(ino,jno,kno),ml,m1(ino-2,jno,kno),1,2,2)
        ELSE
          m1x=derc(m1(ino,jno,kno),ml,mr,al(ino,jno,kno),ar(ino,jno,kno),den1)
        END IF
        IF ad(ino,jno,kno) < pfo THEN 
          m2y=ders(m2(ino,jno,kno),m2(ino,jno+1,kno),m2(ino,jno+2,kno),1,2,2)
        ELSE
          m2y=derc(m2(ino,jno,kno),md,m2(ino,jno+1,kno),ad(ino,jno,kno),1,den2)
        END IF
        IF  ab(ino,jno,kno) < pfo AND  kno <  nz  THEN
          m3z=ders(m3(ino,jno,kno),mf,m3(ino,jno,kno+2),1,2,2)
        ELSE IF  af(ino,jno,kno) < pfo AND  kno > 0  THEN
          m3z=-ders(m3(ino,jno,kno),mb,m3(ino,jno,kno-2),1,2,2)
        ELSE
          m3z=derc(m3(ino,jno,kno),mb,mf,ab(ino,jno,kno),af(ino,jno,kno),den3)
          IF(af(ino,jno,kno) < pfo AND kno=0) OR (ab(ino,jno,kno) < pfo  AND kno=nz ) THEN  m3z=(mf-mb)/abf/h  
          
        END IF
        rhs(ino,jno,kno)=-(m1x+m2y+m3z)∗coenoto(ino,jno,kno)∗h∗h
      REPEAT LOOP
      LOOP FOR  jno=btf(ino,kno) TO  ny-1
        rhs(ino,jno,kno)=-(m1(ino+1,jno,kno)-m1(ino-1,jno,kno)+m2(ino,jno+1,kno)-m2(ino,jno-1,kno)+m3(ino,jno,kno+1)-m3(ino,jno,kno-1))∗h/12
      REPEAT LOOP
      rhs(ino,ny,kno)=-(m1(ino+1,ny,kno)-m1(ino-1,ny,kno)-2.∗m2(ino,ny-1,kno)+m3(ino,ny,kno+1)-m3(ino,ny,kno-1))∗h/12
    REPEAT LOOP
  REPEAT LOOP
END  notofia

! ∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗





SUBROUTINE corm1
  REAL ficen,ml,mr,md,mf,mb,fix,alr,rhs1,cen,lef,rig,u,d,ba,fo,xv
  REAL den,cl,cr,cd,cf,cb,fl,fll,fr,fb,fff,fu,fd,frr,tang,ll,rr
  INTEGER ndis
  LOOP FOR  k=0 TO nz
    LOOP FOR  i=0 TO nx
      IF  i < 2  OR  i > nx-2 THEN
        IF i∗h < 0.25 THEN
          xv=0.
        ELSE
          xv=0.5
        END IF
        
        IF  fondo(xv,k∗h)>2∗h THEN  
          LOOP FOR  j=btm(i,k) TO btm(0,0)+2
            ndis=nid(i,j,k)
            IF  ndis > 0  THEN 
              del=vdel(ndis)
              k1=vk1(ndis)
              k2=vk2(ndis)
              zn=vzn(ndis)
              calfa=vcalfa(ndis)
              salfa=vsalfa(ndis)
              costanti(i,j,k)
              alr=al(i,j,k)+ar(i,j,k)
              den=al(i,j,k)∗ar(i,j,k)∗alr
              u=m1co(i,j+1,k)
              cen=m1co(i,j,k)
              ba=0.
              fo=0.
              lef=0.
              rig=0.
              mr=0.
              md=0.
              ml=0.
              mf=0.
              mb=0.
              cr=0.
              cl=0.
              cd=0.
              cb=0.
              cf=0.
              d=0
              ficen=fico(i,j,k)
              fu=fico(i,j+1,k)
              IF btm(i-1,k) > j THEN
                fl=0
                rr=0
              ELSE
                fl=fico(i-1,j,k)
                rr=al(i-1,j,k)
              END IF
              
              IF i=0 THEN
                fll=wp∗fico(i+2,j,k)
                IF btm(i+2,k)> j THEN fll=0
              ELSE
                fll=fico(i-2,j,k)
                IF btm(i-2,k)> j THEN fll=0
              END IF
              
              IF btm(i+1,k) > j THEN
                fr=0
                ll=0
              ELSE
                fr=fico(i+1,j,k)
                ll=al(i+1,j,k)
              END IF
              
              IF i=nx THEN
                frr=wp∗fico(i-2,j,k)
                IF btm(i-2,k)> j THEN frr=0
              ELSE
                frr=fico(i+2,j,k)
                IF btm(i+2,k)> j THEN frr=0
              END IF
              
              IF btm(i,k-1) > j THEN
                fb=0
              ELSE
                fb=fico(i,j,k-1)
              END IF
              IF btm(i,k+1) > j THEN
                fff=0
              ELSE
                fff=fico(i,j,k+1)
              END IF
              IF btm(i,k)=j THEN
                fd=0
              ELSE
                fd=fico(i,j-1,k)
              END IF
              
              IF  al(i,j,k) < pfo   THEN
                fix=ders(ficen,fr,frr,1,2,2)
              ELSE IF  ar(i,j,k) < pfo   THEN
                fix=-ders(ficen,fl,fll,1,2,2)
              ELSE
                fix=derc(ficen,fl,fr,al(i,j,k),ar(i,j,k),den) 
              END IF
              IF  btm(i-1,k) >  j  THEN
                cl=bcl(i,j,k)
                lef=0.
                IF  al(i,j,k)  < pfo   THEN
                  ml=ders(0,fr,frr,alr,alr+1,alr∗(alr+1))
                ELSE
                  ml=ders(0,ficen,fr,al(i,j,k),alr,den)
                END IF
              ELSE
                lef=m1co(i-1,j,k)
              END IF
              
              IF  btm(i+1,k) >  j  THEN
                cr=bcr(i,j,k)
                rig=0.
                IF  ar(i,j,k) < pfo   THEN
                  mr=-ders(0,fl,fll,alr,alr+1.,alr∗(alr+1.))
                ELSE
                  mr=-ders(0,ficen,fl,ar(i,j,k),alr,den)
                END IF
              ELSE
                rig=m1co(i+1,j,k)
              END IF
              
              IF j > btm(i,k)  THEN
                d=m1co(i,j-1,k)
              ELSE
                ! calcolo md  ∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗
                cd=bcd(i,j,k)
                d=0.
                tang=(fondo(i∗h-eps,k∗h)-fondo(i∗h+eps,h∗k))/2/eps
                IF pbc=1 THEN
                  md=fix
                ELSE IF pbc=2 THEN
                  IF i=0 OR i= nx THEN
                    md=fix
                  ELSE
                    IF tang=0 THEN 
                      md=0
                    ELSE
                      md=fxd(tang,i,j,k,1,1)    
                    END IF
                  END IF
                END IF
              END IF
              
              IF btm(i,k-1) > j  THEN
                cb=bcb(i,j,k)
                ba=0.
                IF pbc =1 THEN
                  mb=fix
                ELSE IF pbc =2 THEN
                  mb=fix
                END IF
              ELSE
                ba=m1co(i,j,k-1)
              END IF
              
              IF btm(i,k+1) > j  THEN
                cf=bcf(i,j,k)
                fo=0
                IF pbc =1 THEN
                  mf=fix
                ELSE IF pbc =2 THEN
                  mf=fix
                END IF
                
              ELSE
                fo=m1co(i,j,k+1)
              END IF
              
              rhs1=cl∗ml+cd∗md+cb∗mb+cf∗mf+cr∗mr
              rhs(i,j,k)=rhs(i,j,k)-left(i,j,k)∗lef-right(i,j,k)∗rig-
      down(i,j,k)∗d-up(i,j,k)∗u-back(i,j,k)∗ba-forw(i,j,k)∗fo+cen-rhs1
            END IF
          REPEAT LOOP
        END IF
      END IF
    REPEAT LOOP
  REPEAT LOOP
END corm1



SUBROUTINE corm2
  REAL ficen,ml,mr,md,mf,mb,fiy,fd,rhs1,cen,lef,rig,u,d,ba,fo,xv
  REAL den,cl,cr,cd,cf,cb,tang,fiy0,fir,fu,fl,fr,fb,fff
  INTEGER ndis
  LOOP FOR  k=0 TO nz
    LOOP FOR  i=0 TO nx
      IF  i < 2  OR  i > nx-2 THEN
        IF i∗h < 0.25 THEN
          xv=0.
        ELSE
          xv=0.5
        END IF
        IF fondo(xv,k∗h)>2∗h THEN 
          LOOP FOR  j=btm(i,k) TO btm(0,0)+2
            ndis=nid(i,j,k)
            IF  ndis > 0  THEN 
              del=vdel(ndis)
              k1=vk1(ndis)
              k2=vk2(ndis)
              zn=vzn(ndis)
              calfa=vcalfa(ndis)
              salfa=vsalfa(ndis)
              costanti(i,j,k)
              ad1=ad(i,j,k)+1.
              den=ad(i,j,k)∗ad1
              cl=0.
              cd=0.
              cr=0.
              cf=0.
              cb=0.
              ml=0.
              md=0.
              mr=0.
              mf=0.
              mb=0.
              ba=0
              fo=0
              lef=0
              rig=0
              d=0
              u=m2co(i,j+1,k)
              cen=m2co(i,j,k)
              ficen=fico(i,j,k)
              fu=fico(i,j+1,k)
              IF btm(i-1,k) > j THEN
                fl=0
              ELSE
                fl=fico(i-1,j,k)
              END IF
              IF btm(i+1,k) > j THEN
                fr=0
              ELSE
                fr=fico(i+1,j,k)
              END IF
              
              IF btm(i,k-1) > j THEN
                fb=0
              ELSE
                fb=fico(i,j,k-1)
              END IF
              IF btm(i,k+1) > j THEN
                fff=0
              ELSE
                fff=fico(i,j,k+1)
              END IF
              IF btm(i,k)=j THEN
                fd=0
              ELSE
                fd=fico(i,j-1,k)
              END IF
              
              IF ad(i,j,k) < pfo THEN
                fiy=ders(ficen,fu,fico(i,j+2,k),1,2,2)
              ELSE
                fiy=derc(ficen,fd,fu,ad(i,j,k),1,den)
              END IF
              
              IF j=btm(i,k) THEN
                IF ad(i,j,k) < pfo THEN
                  fiy0=ders(0,fu,fico(i,j+2,k),ad1,ad1+1,ad1∗(ad1+1))
                ELSE
                  fiy0=ders(0,ficen,fu,ad(i,j,k),ad1,den)
                END IF
              END IF
              
              tang=(fondo(i∗h-eps,k∗h)-fondo(i∗h+eps,h∗k))/2/eps
              
              IF  btm(i-1,k) > j  THEN
                cl=bcl(i,j,k)
                lef=0.
                IF pbc=1 THEN
                  ml=fiy
                ELSE IF pbc=2 THEN
                  IF i=1 AND j= btm(i-1,k) THEN
                    !   ml=fiy
                    ml=0
                  ELSE
                    ml=fyl(i,j,k,1,1,1)
                  END IF
                END IF
              ELSE
                lef=m2co(i-1,j,k)
              END IF
              
              IF  btm(i+1,k) >  j  THEN
                cr=bcr(i,j,k)
                rig=0.
                IF pbc=1 THEN
                  mr=fiy
                ELSE  IF pbc=2 THEN
                  IF i=nx-1 AND j= btm(i+1,k) THEN
                    !   mr=fiy
                    mr=0
                  ELSE
                    mr=fyl(i,j,k,1,1,-1)
                  END IF
                END IF
              ELSE
                rig=m2co(i+1,j,k)
              END IF
              
              IF  j > btm(i,k)  THEN
                d=m2co(i,j-1,k)
              ELSE
                cd=bcd(i,j,k)
                d=0.
                IF (i=0 OR i=nx) AND j >1 THEN
                  md=0
                ELSE
                  md=fiy0
                END IF
              END IF
              
              IF btm(i,k-1) >  j  THEN
                cb=bcb(i,j,k)
                ba=0
                IF pbc=1 THEN
                  mb=fiy
                ELSE IF pbc=2 THEN
                  !  mb=fiy
                  mb=fyl(i,j,k,0,1,1)
                END IF
              ELSE
                ba=m2co(i,j,k-1)
              END IF
              IF btm(i,k+1) >  j  THEN
                cf=bcf(i,j,k)
                fo=0
                IF pbc=1 THEN
                  mf=fiy
                ELSE IF pbc=2 THEN
                  !  mf=fiy
                  mf=fyl(i,j,k,0,1,-1)
                END IF
              ELSE
                fo=m2co(i,j,k+1)
              END IF
              
              rhs1=cl∗ml+cd∗md+cb∗mb+cf∗mf+cr∗mr
              rhs(i,j,k)=rhs(i,j,k)-left(i,j,k)∗lef-right(i,j,k)∗rig-down(i,j,k) \
              ∗d-up(i,j,k)∗u-back(i,j,k)∗ba-forw(i,j,k)∗fo+cen-rhs1
            END IF
          REPEAT LOOP
        END IF
      END IF
    REPEAT LOOP
  REPEAT LOOP
END corm2

SUBROUTINE corm3
  REAL ficen,ml,mr,md,mf,mb,fiz,rhs1,cen,lef,rig,u,d,ba,fo,abf,xv
  REAL den,cl,cr,cd,cf,cb,fiq,fit,tangz,fl,fr,fb,fff,fu,ffff,fbb
  INTEGER ndis
  LOOP FOR  k=0 TO nz
    LOOP FOR  i=0 TO nx
      IF  i < 2  OR  i > nx-2 THEN
        IF i∗h < 0.25 THEN
          xv=0.
        ELSE
          xv=0.5
        END IF
        IF fondo(xv,h∗k)>2∗h THEN  
          LOOP FOR  j=btm(i,k) TO btm(0,0)+2
            ndis=nid(i,j,k)
            
            IF  ndis > 0  THEN 
              del=vdel(ndis)
              k1=vk1(ndis)
              k2=vk2(ndis)
              zn=vzn(ndis)
              calfa=vcalfa(ndis)
              salfa=vsalfa(ndis)
              costanti(i,j,k)
              !   WRITE 'c6,i,j,k',c6,i,j,k
              !    IF (i=0 AND j=btm(0,0)) OR (i=nx AND j=btm(nx,0))THEN WRITE i,j,c6
              abf=ab(i,j,k)+af(i,j,k)
              den=ab(i,j,k)∗af(i,j,k)∗abf
              ad1=ad(i,j,k)+1
              cb=0.
              cf=0.
              cl=0.
              cr=0.
              cd=0.
              mf=0.
              mb=0.
              md=0.
              ml=0.
              mr=0.
              u=m3co(i,j+1,k)
              cen=m3co(i,j,k)
              ficen=fico(i,j,k)
              fu=fico(i,j+1,k)
              IF btm(i-1,k) > j THEN
                fl=0
              ELSE
                fl=fico(i-1,j,k)
              END IF
              IF btm(i+1,k) > j THEN
                fr=0
              ELSE
                fr=fico(i+1,j,k)
              END IF
              
              IF btm(i,k-1) > j THEN
                fb=0
              ELSE
                fb=fico(i,j,k-1)
              END IF
              IF btm(i,k+1) > j THEN
                fff=0
              ELSE
                fff=fico(i,j,k+1)
              END IF
              ! calcolo fiz
              IF  ab(i,j,k) < pfo  AND  k <  nz  THEN
                fiz=ders(ficen,fff,fico(i,j,k+2),1,2,2)
              ELSE IF  af(i,j,k) < pfo AND  k > 0  THEN
                fiz=-ders(ficen,fb,fico(i,j,k-2),1,2,2)
              ELSE
                fiz=derc(ficen,fb,fff,ab(i,j,k),af(i,j,k),den)
              END IF
              ! fine calcolo fiz
              IF(ab(i,j,k) < pfo AND k=nz) OR (af(i,j,k)<pfo AND k=0) THEN fiz=(fff-fb)/abf/h
              
              IF  btm(i,k-1) >  j  THEN
                ba=0.
                cb=bcb(i,j,k)
                IF  ab(i,j,k) < pfo AND  k < nz  THEN
                  mb=ders(0,fff,fico(i,j,k+2),abf,abf+1,abf∗(abf+1))
                ELSE
                  mb=ders(0,ficen,fff,ab(i,j,k),abf,den)
                END IF
                IF i=0 OR i=nx THEN mb=0
              ELSE
                ba=m3co(i,j,k-1)
              END IF
              
              !   IF  ab(i,j,k) < pfo  AND  k = nz  THEN  mb=fiz
              
              IF  btm(i,k+1) >  j  THEN
                fo=0.
                cf=bcf(i,j,k)
                IF  af(i,j,k) < pfo AND  k > 0  THEN
                  mf=-ders(0,fb,fico(i,j,k-2),abf,abf+1,abf∗(abf+1))
                ELSE
                  mf=-ders(0,ficen,fb,af(i,j,k),abf,den)
                END IF
                IF i=0 OR i=nx THEN mf=0
              ELSE
                fo=m3co(i,j,k+1)
              END IF
              
              !  IF      af(i,j,k) < pfo AND   k=0  THEN mf=fiz
              
              ! calcolo md  ∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗
              IF j > btm(i,k)  THEN
                d=m3co(i,j-1,k)
              ELSE
                d=0
                cd=bcd(i,j,k)
                tangz=(fondo(i∗h,k∗h-eps)-fondo(i∗h,h∗k+eps))/2/eps
                IF pbc=1 THEN
                  md=fiz
                ELSE IF pbc=2 THEN
                  !   md=fiz
                  md=fxd(tangz,i,j,k,0,1)
                  IF i=0 OR i=nx THEN md=0
                  !(    IF i=0 OR i=nx THEN
                      IF tangz>0 THEN
                        md=(tangz/ad(i,j,k+1))∗∗1.5/h∗fff
                      ELSE IF tangz<0 THEN
                        md=-(-tangz/ad(i,j,k-1))∗∗1.5/h∗fb
                      ELSE
                        md=0
                      END IF
                    END IF
                    !)
                END IF
                IF tangz=0 THEN md=0
                IF k=0 AND wp=-1 THEN md=0
              END IF
              !   fine calcolo md ∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗
              
              IF btm(i-1,k) > j  THEN
                cl=bcl(i,j,k)
                lef=0
                IF pbc=1 THEN
                  ml=fiz
                  !   ml=0
                ELSE IF pbc=2 THEN
                  !   ml=fiz
                  ml=0
                END IF
              ELSE
                lef=m3co(i-1,j,k)
              END IF
              
              IF btm(i+1,k) > j  THEN
                cr=bcr(i,j,k)
                rig=0
                IF pbc=1 THEN
                  mr=fiz
                  !   mr=0
                ELSE  IF pbc=2 THEN
                  !  mr=fiz
                  mr=0
                END IF
              ELSE
                rig=m3co(i+1,j,k)
              END IF
              rhs1=cl∗ml+cd∗md+cb∗mb+cf∗mf+mr∗cr
              rhs(i,j,k)=rhs(i,j,k)-left(i,j,k)∗lef-right(i,j,k)∗rig-down(i,j,k) \
              ∗d-up(i,j,k)∗u-back(i,j,k)∗ba-forw(i,j,k)∗fo+cen-rhs1
            END IF
          REPEAT LOOP
        END IF
      END IF
    REPEAT LOOP
  REPEAT LOOP 
END corm3




! ∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗


SUBROUTINE corfi
  
  
  REAL den1,den2,den3
  REAL m1x,m2y,m3z,alr,abf,ml,mr,mb,mf,ad1,md,ficen,lef,rig,d,fu,ba,xv,fo
  REAL fl,fr,fb,fff,fll,frr
  INTEGER ndis
  LOOP FOR  k=0 TO nz
    LOOP FOR  i=0 TO nx
      IF  i < 2  OR  i > nx-2 THEN
        IF i∗h < 0.25 THEN
          xv=0.
        ELSE
          xv=0.5
        END IF
        IF fondo(xv,h∗k)>2∗h THEN
          LOOP FOR  j=btm(i,k) TO btm(0,0)+2
            ndis=nid(i,j,k)
            IF  ndis > 0  THEN 
              del=vdel(ndis)
              k1=vk1(ndis)
              k2=vk2(ndis)
              zn=vzn(ndis)
              calfa=vcalfa(ndis)
              salfa=vsalfa(ndis)
              costanti(i,j,k)
              alr=al(i,j,k)+ar(i,j,k)
              abf=ab(i,j,k)+af(i,j,k)
              ad1=ad(i,j,k)+1.
              den1=al(i,j,k)∗ar(i,j,k)∗alr
              den2=ad(i,j,k)∗ad1
              den3=af(i,j,k)∗ab(i,j,k)∗abf
              mr=m1co(i+1,j,k)
              rig=fico(i+1,j,k)
              fu=fico(i,j+1,k)
              fff=fico(i,j,k+1)
              IF btm(i,k+1) > j  THEN  fff=0.
              fb=fico(i,j,k-1)
              IF btm(i,k-1) > j  THEN  fb=0.
              ficen=fico(i,j,k)
              IF btm(i-1,k) > j THEN
                fl=0
              ELSE
                fl=fico(i-1,j,k)
              END IF
              
              IF i=0 THEN
                fll=wp∗fico(i+2,j,k)
                IF btm(i+2,k)> j THEN fll=0
              ELSE
                fll=fico(i-2,j,k)
                IF btm(i-2,k)> j THEN fll=0
              END IF
              
              IF btm(i+1,k) > j THEN
                fr=0
              ELSE
                fr=fico(i+1,j,k)
              END IF
              
              IF i=nx THEN
                frr=wp∗fico(i-2,j,k)
                IF btm(i-2,k)> j THEN frr=0
              ELSE
                frr=fico(i+2,j,k)
                IF btm(i+2,k)> j THEN frr=0
              END IF
              
              IF   j > btm(i,k)  THEN    
                d=fico(i,j-1,k)
                md=m2co(i,j-1,k)
              ELSE
                d=0.
                IF ad(i,j,k) < pfo THEN
                  md=ders(0,fu,fico(i,j+2,k),ad1,ad1+1,ad1∗(ad1+1))
                ELSE
                  md=ders(0,ficen,fu,ad(i,j,k),ad1,den2)
                END IF
                IF (i=0 OR i=nx) AND j>1 THEN md=0
              END IF
              IF  btm(i+1,k) > j  THEN
                rig=0.
                IF ar(i,j,k) < pfo THEN
                  mr=-ders(0,fl,fll,alr,alr+1,alr∗(alr+1))
                ELSE
                  mr=-ders(0,ficen,fl,ar(i,j,k),alr,den1)
                END IF
              ELSE
                rig=fico(i+1,j,k)
                mr=m1co(i+1,j,k)
              END IF
              IF  btm(i-1,k) > j THEN
                lef=0.
                IF  al(i,j,k) < pfo   THEN
                  ml=ders(0,fr,frr,alr,alr+1,alr∗(alr+1))
                ELSE
                  ml=ders(0,ficen,fr,al(i,j,k),alr,den1)
                END IF
              ELSE
                lef=fico(i-1,j,k)
                ml=m1co(i-1,j,k)
              END IF
              IF  btm(i,k+1) > j THEN
                fo=0.
                IF  af(i,j,k) < pfo AND  k > 0  THEN
                  mf=-ders(0.,fico(i,j,k-1),fico(i,j,k-2),abf,abf+1,abf∗(abf+1))
                ELSE
                  mf=-ders(0.,fico(i,j,k),fico(i,j,k-1),af(i,j,k),abf,den3)
                END IF
                IF  (af(i,j,k) < pfo OR  ab(i,j,k) < pfo)  AND (k=nz OR k=0)  THEN   mf=(fff-fb)/abf/h
                IF i=0 OR i=nx  THEN mf=0
              ELSE
                fo=fico(i,j,k+1)
                mf=m3co(i,j,k+1)
              END IF
              IF btm(i,k-1) >  j  THEN
                ba=0.
                IF  ab(i,j,k) < pfo AND  k < nz  THEN
                  mb=ders(0.,fico(i,j,k+1),fico(i,j,k+2),abf,abf+1,abf∗(abf+1))
                ELSE
                  mb=ders(0.,fico(i,j,k),fico(i,j,k+1),ab(i,j,k),abf,den3)
                END IF
                IF(af(i,j,k) < pfo OR  ab(i,j,k) < pfo)  AND (k=nz OR k=0 ) THEN mb=(fff-fb)/abf/h
                IF i=0 OR i=nx  THEN mb=0
              ELSE
                ba=fico(i,j,k-1)
                mb=m3co(i,j,k-1)
              END IF
              IF  al(i,j,k) < pfo AND  i <  nx  THEN
                m1x=ders(m1co(i,j,k),mr,m1co(i+2,j,k),1.,2.,2.)
              ELSE IF  ar(i,j,k) < pfo AND  i > 0  THEN
                m1x=-ders(m1co(i,j,k),ml,m1co(i-2,j,k),1.,2.,2.)
              ELSE
                m1x=derc(m1co(i,j,k),ml,mr,al(i,j,k),ar(i,j,k),den1)
              END IF
              IF ad(i,j,k) < pfo THEN 
                m2y=ders(m2co(i,j,k),m2co(i,j+1,k),m2co(i,j+2,k),1,2,2)
              ELSE
                m2y=derc(m2co(i,j,k),md,m2co(i,j+1,k),ad(i,j,k),1,den2)
              END IF
              IF  ab(i,j,k) < pfo  AND  k  < nz  THEN
                m3z=ders(m3co(i,j,k),mf,m3co(i,j,k+2),1.,2.,2.)
              ELSE IF  af(i,j,k) < pfo AND  k > 0  THEN
                m3z=-ders(m3co(i,j,k),mb,m3co(i,j,k-2),1.,2.,2.)
              ELSE
                m3z=derc(m3co(i,j,k),mb,mf,ab(i,j,k),af(i,j,k),den3)
                IF  ( af(i,j,k) < pfo OR ab(i,j,k) < pfo)  AND (k = nz OR k=0) THEN m3z=(mf-mb)/abf/h
              END IF
              rhs(i,j,k)=rhs(i,j,k)+(m1x+m2y+m3z)∗coenoto(i,j,k)∗h∗h+ficen-lef∗ \
              left(i,j,k)-rig∗right(i,j,k)-d∗down(i,j,k)-fu∗up(i,j,k)-fo∗ \
              forw(i,j,k)-ba∗back(i,j,k)
            END IF
            
          REPEAT LOOP
        END IF
      END IF
    REPEAT LOOP
  REPEAT LOOP
END corfi

! ∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗  FINE  SUBROUTINES  ∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗

LOOP kapp FOR  fk=k0 TO numk 
  ka=fk
  rapp=rap0
  LOOP FOR  fr=1 TO numh
    rapp=rapp+delth
    LOOP  FOR   ff=1 TO numf
      fill=fil0+deltf∗ff
      
      WRITE 'ly,lz.fill,rapp,rapp/h',ly,lz,fill,rapp,rapp/h
      WRITE'kappa',ka
      
      LOOP FOR   i=-1 TO nx+1
        LOOP FOR  k=-1 TO nz+1
          btm(i,k)=0
        REPEAT LOOP
      REPEAT LOOP
      
      LOOP FOR   i=0 TO nx
        LOOP FOR  k=0 TO nz
          btf(i,k)=0
        REPEAT LOOP
      REPEAT LOOP
      LOOP FOR   i=-1 TO nx+1
        LOOP FOR   j=0 TO ny
          LOOP FOR  k=-1 TO nz+1
            fi(i,j,k)=0.
            m1(i,j,k)=0.
            m2(i,j,k)=0.
            m3(i,j,k)=0.
          REPEAT LOOP
        REPEAT LOOP
      REPEAT LOOP
      LOOP FOR   i=0 TO nx
        LOOP FOR   j=0 TO ny0  
          LOOP FOR  k=0 TO nz
            ad(i,j,k)=1.
            al(i,j,k)=1.
            ar(i,j,k)=1.
            af(i,j,k)=1.
            ab(i,j,k)=1.
            left(i,j,k)=0.
            right(i,j,k)=0.
            up(i,j,k)=0.
            down(i,j,k)=0.
            back(i,j,k)=0.
            forw(i,j,k)=0.
            bcl(i,j,k)=0.
            bcr(i,j,k)=0.
            bcd(i,j,k)=0. 
            bcu(i,j,k)=0.
            bcb(i,j,k)=0.
            bcf(i,j,k)=0.
            coenoto(i,j,k)=0.
          REPEAT LOOP
        REPEAT LOOP
      REPEAT LOOP
      LOOP FOR   i=0 TO nx
        LOOP FOR   j=0 TO ny  
          LOOP FOR  k=0 TO nz
            rhs(i,j,k)=0
          REPEAT LOOP
        REPEAT LOOP
      REPEAT LOOP
      LOOP FOR   i=0 TO nx
        x=i∗h
        LOOP FOR  k=0 TO nz
          z=k∗h
          btm(i,k)=1+FLOOR(fondo(x,z)/h)
          
          IF ABS(1+fondo(x,z)/h- btm(i,k))> 1-h/100 THEN btm(i,k)=btm(i,k)+1
          
        REPEAT LOOP
      REPEAT LOOP
      LOOP FOR  k=-1 TO nz+1
        btm(nx+1,k)=btm(nx-1,k) 
        btm(-1,k)=btm(1,k)
      REPEAT LOOP
      
      LOOP FOR  i=-1 TO nx+1
        btm(i,-1)=btm(i,1)
        btm(i,nz+1)=btm(nx-i,nz-1)
      REPEAT LOOP
      !      calcolo di btf                       
      ! ∗∗∗∗∗o∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗     
      LOOP FOR  k=0 TO nz
        LOOP FOR  i=0 TO nx
          btf(i,k)=MAX(btm(i,k)+1,btm(i-1,k),btm(i+1,k))
          btf(i,k)=MAX(btf(i,k),btm(i,k-1),btm(i,k+1))
        REPEAT LOOP
        
      REPEAT LOOP 
      WRITE 'fine calcolo btf '  
      LOOP FOR  k=0 TO nz
        z=k∗h
        LOOP FOR   i=0 TO nx
          x=h∗i
          LOOP FOR   j=btm(i,k) TO btf(i,k)-1
            y=h∗j
            IF j =btm(i,k) THEN  ad(i,j,k)=j-fondo(x,z)/h
            IF btm(i-1,k) > j THEN al(i,j,k)=(x-vert(y,x,x-h,z))/h
            IF btm(i+1,k)  > j THEN ar(i,j,k)=(vert(y,x,x+h,z)-x)/h
            IF  i = nx  THEN ar(i,j,k)=al(i,j,k)
            !       WRITE'condizione di simmmetria rispetto ad x',i,j,k
            IF  btm(i,k-1) > j  THEN ab(i,j,k)=(z-vertz(y,z,z-h,x))/h
            IF btm(i,k+1) > j THEN af(i,j,k)=(vertz(y,z,z+h,x)-z)/h
            ! condizione di simmetria rispetto all'asse x
            IF z=0 THEN ab(i,j,k)=af(i,j,k)
            
            deny=ad(i,j,k)∗(1.+ad(i,j,k))
            denx=al(i,j,k)∗ar(i,j,k)∗(al(i,j,k)+ar(i,j,k))
            denz=af(i,j,k)∗ab(i,j,k)∗(af(i,j,k)+ab(i,j,k))
            dud=deny/2.
            dlr=denx/2.
            dbf=denz/2.
            axr=al(i,j,k)
            axl=ar(i,j,k)
            ax0=axl+axr
            ayu=ad(i,j,k)
            ayd=1.
            ay0=ayu+ayd 
            azf=ab(i,j,k)
            azb=af(i,j,k)
            az0=azf+azb
            deno=ax0∗dbf∗dud+ay0∗dbf∗dlr+az0∗dlr∗dud
            bcl(i,j,k)=axl∗dud∗dbf/deno
            bcd(i,j,k)=ayd∗dlr∗dbf/deno
            bcr(i,j,k)=axr∗dud∗dbf/deno
            bcu(i,j,k)=ayu∗dlr∗dbf/deno
            bcb(i,j,k)=azb∗dud∗dlr/deno
            bcf(i,j,k)=azf∗dlr∗dud/deno
            coenoto(i,j,k)=dud∗dlr∗dbf/deno
          REPEAT LOOP
          LOOP FOR  j= btf(i,k) TO ny0
            bcl(i,j,k)=1./6.
            bcd(i,j,k)=1./6.
            bcr(i,j,k)=1./6.
            bcu(i,j,k)=1./6.
            bcf(i,j,k)=1./6.
            bcb(i,j,k)=1./6.
            coenoto(i,j,k)=1./6.
          REPEAT LOOP
        REPEAT LOOP
      REPEAT LOOP
      LOOP FOR  i=0 TO nx
        LOOP FOR  k=0 TO nz
          LOOP FOR  j=btm(i,k) TO ny0
            left(i,j,k)=bcl(i,j,k)
            IF btm(i-1,k) > j THEN left(i,j,k)=0.
            right(i,j,k)=bcr(i,j,k)
            IF btm(i+1,k) > j THEN right(i,j,k)=0.
            back(i,j,k)=bcb(i,j,k)
            IF btm(i,k-1) > j THEN back(i,j,k)=0.
            forw(i,j,k)=bcf(i,j,k)
            IF btm(i,k+1) > j THEN forw(i,j,k)=0.
            down(i,j,k)=bcd(i,j,k)
            IF btm(i,k)= j THEN down(i,j,k)=0.
            up(i,j,k)=bcu(i,j,k)
            
          REPEAT LOOP
        REPEAT LOOP
      REPEAT LOOP
      !   WRITE 'al1,ar,ad1,ad,,bcl,bcr',al(1,btm(0,0)-1,0),ar(nx-1,btm(nx,0)-1,0),ad(0,
      !    btm(0,0),0),ad(nx,btm(nx,0),0),bcl(1,btm(0,0)-1,0),bcr(nx-1,btm(nx,0)-1,0)
      !    READ datof
      
      !(        calcolo diagramma
        LOOP FOR k=0 TO 32
          !    LOOP FOR k=0 TO 100
          !   ALTEZZE k/100∗lz/2∗5/8, fondo(0,k/100∗lz/2∗5/8),btm(0,k).
          ALTEZZE k, fondo(0,k∗h)/h,btm(0,k).
        REPEAT LOOP
        STOP
        !)
      !∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗∗
      LOOP FOR   fw=0 TO 1
        wp=-1+fw∗2 
        LOOP FOR   i=-1 TO nx+1
          LOOP FOR   j=0 TO ny
            LOOP FOR  k=-1 TO nz+1
              fi(i,j,k)=0.
              m1(i,j,k)=0.
              m2(i,j,k)=0.
              m3(i,j,k)=0.
            REPEAT LOOP
          REPEAT LOOP
        REPEAT LOOP
        !     calcolo di dati per correzioni
        IF inco=1 THEN  catalogo
        !( 
          LOOP FOR k=0 TO nz-3
            it=nid(0,btm(0,k),k)
            WRITE 'it,z,zn',it,k∗h,vzn(it),
            'sa,ca,de,k1,k2',vsalfa(it),vcalfa(it),vdel(it),vk1(it),vk2(it)
            ALTEZZE k,vzn(it),vsalfa(it)
          REPEAT LOOP
          STOP
          !)
        WRITE'inizio ciclo,btm(0,0),btm(0,1),btm(0,-1)',btm(0,0),btm(0,1),btm(0,-1)
        js=0
        it=1
        DO
          
          IF js=0 THEN prec=((wp+1)∗m3(0,ny,0)+(1-wp)∗m1(0,ny,0))/2
          IF it<jsta THEN prec= 1000
          
          notom1
          IF it >  10 AND inco=1  THEN    corm1
          vbmg(m1,btm,btf,left,down,right,up,back,forw,rhs,nx,ny,nz,-wp,-wp)
          
          notofia
          IF it >  10 AND inco=1 THEN   corfi
          vbmg(fi,btm,btf,left,down,right,up,back,forw,rhs,nx,ny,nz,wp,-wp)
          notom2
          IF it > 10 AND inco=1  THEN  corm2
          vbmg(m2,btm,btf,left,down,right,up,back,forw,rhs,nx,ny,nz,wp,-wp)
          
          notom3
          IF it >  10 AND inco=1   THEN    corm3
          vbmg(m3,btm,btf,left,down,right,up,back,forw,rhs,nx,ny,nz,wp,wp)
          !    vbrb(m3,btm,btf,left,down,right,up,back,forw,rhs,nx,ny,nz,wp,wp)
          
          IF  js =jsta   THEN
            WRITE'm1c,m2c,m3c,fi',m1(0,ny,0):9.9,m2(0,ny,0),m3(0,ny,0):9.11,fi(0,ny,0),it
            js=0
          ELSE
            js=js+1
          END IF
          it=it+1
          IF wp  = 1 THEN hlon=m3(0,ny,0)-ny∗h+fondo(0.,0.)
          IF wp =  -1 THEN htra=m1(0,ny,0)-ny∗h+fondo(0.,0.)
        UNTIL ABS(prec-((wp+1)∗m3(0,ny,0)+(1-wp)∗m1(0,ny,0))/2) < 0.000001
      REPEAT LOOP
      ! ∗∗∗∗∗∗∗∗∗∗∗ fine ciclo parita' ∗∗∗∗∗∗∗∗∗∗
      
      dh=hlon-htra
      WRITE 'hlon,htra,dh',hlon,htra,dh
      fillm=fill
      fils=fill
      IF prolo =3 THEN 
        fillm=(fill +rapp/2 +0.144074)∗2/lz
        fils=(fill+0.144074)∗2/lz 
      END IF
      
      ! ALTEZZE rapp:15.10,fillm:15.10,dh:15.10,hlon:15.10,htra:15.10
      ALTEZZE rapp:10.6,fils:10.6,fillm:10.6,dh:15.10,hlon:15.10,htra:15.10
      <∗ fflush(ALTEZZE_); ∗>
      
    REPEAT LOOP
    ! ∗∗∗∗∗∗∗∗∗∗   fine ciclo fill  ∗∗∗∗∗∗∗∗∗∗∗ 
  REPEAT LOOP
  ! ∗∗∗∗∗∗∗∗∗∗∗∗  fine ciclo rapporto ∗∗∗∗∗∗∗      
REPEAT kapp 
! ∗∗∗∗∗∗∗∗∗∗∗∗  fine ciclo kappa ∗∗∗∗∗∗∗