anchor.diff:       &          ( 
        &            + AiNa0 
        &            - KijNaj0
  -c     &            + pt50 * A0_0(:,iflow,jflow)*alpha_LF*shp0(:,j)
  +     &            + KijNbjCn0 * shp0(:,j)
  +CCc     &            + pt50 * A0_0(:,iflow,jflow)*alpha_LF*shp0(:,j)
        &            + pt50 * (A0_0(:,iflow,jflow)+A0_1(:,iflow,jflow))*alpha_LF*shp0(:,j)
  -c     &            +   s * KijNajC0(:,isd,iflow,jflow))*n0(:,isd)
        &            + e*this_mu(:,iflow,jflow)/length_h * ctc(:,iflow,jflow) * shp0(:,j)
        &          ) * shp0(:,i) * WdetJ0
                 enddo
  
positive.code: c
                  AiNa0 = zero
                  KijNaj0 = zero
  c
                  do isd = 1,nsd
                    AiNa0 = AiNa0 + pt50 * Ai0(:,isd,iflow,jflow) * n0(:,isd) * shp0(:,j)
                    do jsd = 1,nsd
                      KijNaj0 = KijNaj0 + pt50 * Kij0(:,isd,jsd,iflow,jflow) * n0(:,isd) * shg0(:,j,jsd)
                    enddo
                  enddo
  c
  
negative.code: c
                  AiNa1 = zero
                  KijNaj1 = zero
                  KijNbjCn0 = zero
  c
                  do isd = 1,nsd
                    AiNa1 = AiNa1 + pt50 * Ai1(:,isd,iflow,jflow) * n0(:,isd) * shp1(:,j)
                    KijNbjCn0 = KijNbjCn0 + pt50 * s * KijNajC0(:,isd,iflow,jflow) * n0(:,isd)
                    do jsd = 1,nsd
                      KijNaj1 = KijNaj1 + pt50 * Kij1(:,isd,jsd,iflow,jflow) * shg1(:,j,jsd) * n0(:,isd)
                    enddo
                  enddo
  c
  