MODULE module_bl_boulac
      real ck_b,ceps_b,vk,temin
      parameter(ceps_b=1/1.4,ck_b=0.4,temin=0.0001,vk=0.4)
   CONTAINS
      subroutine boulac(frc_urb2d,idiff,flag_bep,dz8w,dt,u_phy,v_phy &
                      ,th_phy,rho,qv_curr,qc_curr,hfx &
                      ,qfx,ustar,cp,g &
                      ,rublten,rvblten,rthblten &
                      ,rqvblten,rqcblten &
                      ,tke,dlk,wu,wv,wt,wq,exch_h,exch_m,pblh &
                      ,a_u_bep,a_v_bep,a_t_bep,a_q_bep &
                      ,a_e_bep,b_u_bep,b_v_bep &
                      ,b_t_bep,b_q_bep,b_e_bep,dlg_bep &
                      ,dl_u_bep,sf_bep,vl_bep &
                      ,ids,ide, jds,jde, kds,kde &
                      ,ims,ime, jms,jme, kms,kme &
                      ,its,ite, jts,jte, kts,kte)
      implicit none
   INTEGER:: ids,ide, jds,jde, kds,kde, &
                                    ims,ime, jms,jme, kms,kme, &
                                    its,ite, jts,jte, kts,kte
   integer, INTENT(IN) :: idiff
   REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) :: DZ8W
   REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) :: qv_curr
   REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) :: qc_curr
   REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) :: RHO
   REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) :: TH_PHY
   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 ) :: ustar
   REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: hfx
   REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: qfx
   real, INTENT(IN ) :: g,cp
   REAL, INTENT(IN ):: DT
   REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: FRC_URB2D
   REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: PBLH
   REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) ::a_u_bep
   REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) ::a_v_bep
   REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) ::a_t_bep
   REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) ::a_q_bep
   REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) ::a_e_bep
   REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) ::b_u_bep
   REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) ::b_v_bep
   REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) ::b_t_bep
   REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) ::b_q_bep
   REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) ::b_e_bep
   REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::dlg_bep
   REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) ::dl_u_bep
   REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) ::sf_bep
   REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) ::vl_bep
   LOGICAL, INTENT(IN) :: flag_bep
      real, dimension (ims:ime, kms:kme, jms:jme) ::th_0
        real, dimension( ims:ime, kms:kme, jms:jme ), INTENT(OUT ) :: exch_h
        real, dimension( ims:ime, kms:kme, jms:jme ), INTENT(OUT ) :: exch_m
        real, dimension( ims:ime, kms:kme, jms:jme ), INTENT(INOUT ) :: tke
        real, dimension( ims:ime, kms:kme, jms:jme ), INTENT(OUT ) :: wu
        real, dimension( ims:ime, kms:kme, jms:jme ), INTENT(OUT ) :: wv
        real, dimension( ims:ime, kms:kme, jms:jme ), INTENT(OUT ) :: wt
        real, dimension( ims:ime, kms:kme, jms:jme ), INTENT(OUT ) :: wq
        real, dimension( ims:ime, kms:kme, jms:jme ), INTENT(OUT ) :: dlk
        REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(OUT ) :: RUBLTEN
        REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(OUT ) :: RVBLTEN
        REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(OUT ) :: RTHBLTEN
        REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(OUT ) :: RQVBLTEN
        REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(OUT ) :: RQCBLTEN
      real z1D(kms:kme)
      real dz1D(kms:kme)
      real u1D(kms:kme)
      real v1D(kms:kme)
      real th1D(kms:kme)
      real q1D(kms:kme)
      real qc1D(kms:kme)
      real rho1D(kms:kme)
      real rhoz1D(kms:kme)
      real tke1D(kms:kme)
      real th01D(kms:kme)
      real dlk1D(kms:kme)
      real dls1D(kms:kme)
      real exch1D(kms:kme)
      real sf1D(kms:kme)
      real vl1D(kms:kme)
      real a_u1D(kms:kme)
      real a_v1D(kms:kme)
      real a_t1D(kms:kme)
      real a_q1D(kms:kme)
      real a_qc1D(kms:kme)
      real a_e1D(kms:kme)
      real b_u1D(kms:kme)
      real b_v1D(kms:kme)
      real b_t1D(kms:kme)
      real b_q1D(kms:kme)
      real b_qc1D(kms:kme)
      real b_e1D(kms:kme)
      real dlg1D(kms:kme)
      real dl_u1D(kms:kme)
      real sh1D(kms:kme)
      real bu1D(kms:kme)
      real wu1D(kms:kme)
      real wv1D(kms:kme)
      real wt1D(kms:kme)
      real wq1D(kms:kme)
      real wqc1D(kms:kme)
      real gamma1D(kms:kme)
      real t2_1D(kms:kme)
      real w2_1D(kms:kme)
      real a_e(ims:ime,kms:kme,jms:jme)
      real b_e(ims:ime,kms:kme,jms:jme)
      real bu(ims:ime,kms:kme,jms:jme)
      real sh(ims:ime,kms:kme,jms:jme)
      real wrk(ims:ime)
      integer ix,iy,iz,id,iz_u,iw_u,ig,ir_u,ix1,iy1,igamma
      real ufrac_int
      real vect,time_tke,hour,zzz
      real ustarf,wstar,wts,t2,w2,tstar_w,zzi
      real summ1,summ2,summ3
      save time_tke,hour
        do ix=its,ite
        do iy=jts,jte
        do iz=kts,kte
         th_0(ix,iz,iy)=300.
        enddo
        enddo
        enddo
       z1D=0.
       dz1D=0.
       u1D =0.
       v1D =0.
       th1D=0.
       q1D=0.
       rho1D=0.
       rhoz1D=0.
       tke1D =0.
       th01D =0.
       dlk1D =0.
       dls1D =0.
       exch1D=0.
       sf1D =1.
       vl1D =1.
       a_u1D =0.
       a_v1D =0.
       a_t1D =0.
       a_q1D =0.
       a_qc1D =0.
       a_e1D =0.
       b_u1D =0.
       b_v1D =0.
       b_t1D =0.
       b_q1D =0.
       b_qc1D =0.
       b_e1D =0.
       dlg1D =0.
       dl_u1D=0.
       sh1D =0.
       bu1D =0.
       wu1D =0.
       wv1D =0.
       wt1D =0.
       wq1D =0.
      igamma=1
       do ix=its,ite
       do iy=jts,jte
         z1d(kts)=0.
         do iz= kts,kte
   u1D(iz)=u_phy(ix,iz,iy)
   v1D(iz)=v_phy(ix,iz,iy)
   th1D(iz)=th_phy(ix,iz,iy)
          q1D(iz)=qv_curr(ix,iz,iy)
          qc1D(iz)=qc_curr(ix,iz,iy)
          tke1D(iz)=tke(ix,iz,iy)
   rho1D(iz)=rho(ix,iz,iy)
   th01D(iz)=th_0(ix,iz,iy)
          dz1D(iz)=dz8w(ix,iz,iy)
          z1D(iz+1)=z1D(iz)+dz1D(iz)
         enddo
        rhoz1D(kts)=rho1D(kts)
        do iz=kts+1,kte
         rhoz1D(iz)=(rho1D(iz)*dz1D(iz-1)+rho1D(iz-1)*dz1D(iz))/(dz1D(iz-1)+dz1D(iz))
        enddo
        rhoz1D(kte+1)=rho1D(kte)
        if(flag_bep)then
         do iz=kts,kte
          a_e1D(iz)=a_e_bep(ix,iz,iy)
          b_e1D(iz)=b_e_bep(ix,iz,iy)
          dlg1D(iz)=(z1D(iz)+z1D(iz+1))/2.*(1.-frc_urb2d(ix,iy))+dlg_bep(ix,iz,iy)*frc_urb2d(ix,iy)
          dl_u1D(iz)=dl_u_bep(ix,iz,iy)
          if((1.-frc_urb2d(ix,iy)).lt.1.)dl_u1D(iz)=dl_u1D(iz)/frc_urb2d(ix,iy)
          vl1D(iz)=vl_bep(ix,iz,iy)
          sf1D(iz)=sf_bep(ix,iz,iy)
         enddo
         ufrac_int=frc_urb2d(ix,iy)
         sf1D(kte+1)=sf_bep(ix,1,iy)
        else
         do iz=kts,kte
          a_e1D(iz)=0.
          b_e1D(iz)=0.
          dlg1D(iz)=(z1D(iz)+z1D(iz+1))/2.
          dl_u1D(iz)=0.
          vl1D(iz)=1.
          sf1D(iz)=1.
         enddo
         ufrac_int=0.
         sf1D(kte+1)=1.
        endif
        call pbl_height(kms,kme,kts,kte,dz1d,z1d,th1D,q1D,pblh(ix,iy))
        wts=max(0.,hfx(ix,iy)/rho1D(1)/cp)
        wstar=(g*wts*pblh(ix,iy)/th01D(1))**(1./3.)
        if (wts .ne. 0.0) then
          tstar_w=wts/wstar
        else
          tstar_w=0.0
        endif
        t2_1D=0.
        w2_1D=0.
        gamma1D=0.
         do iz=kts+1,kte
    zzi=z1D(iz)/pblh(ix,iy)
    t2_1D(iz)=1.8*(zzi**(-2./3.))*(tstar_w**2.)
    w2_1D(iz)=1.8*(zzi**(2./3.))*((1.-0.8*zzi)**2.)*(wstar**2.)
        enddo
       if(igamma.eq.1)then
        do iz=kts+1,kte
         if(z1D(iz).le.1.0*pblh(ix,iy).and.wts.gt.0.)then
          gamma1D(iz)=10.*wts/wstar/pblh(ix,iy)
         else
          gamma1D(iz)=0.
         endif
        enddo
       elseif(igamma.eq.2)then
         do iz=kts+1,kte
   if(wts.gt.0)then
           if(z1D(iz).le.(1.0*pblh(ix,iy)).and.z1D(iz).gt.(0.1*pblh(ix,iy)))then
            gamma1D(iz)=g/th01D(iz)*t2_1D(iz)/w2_1D(iz)
           else
            gamma1D(iz)=0.
           endif
   endif
         enddo
       elseif(igamma.eq.3)then
         do iz=kts+1,kte
          if(z1D(iz).le.(1.0*pblh(ix,iy)).and.wts.gt.0)then
           gamma1D(iz)=2.*wstar*wts/w2_1D(iz)/pblh(ix,iy)
          else
           gamma1D(iz)=0.
          endif
         enddo
       endif
         call boulac1D(ix,iy,ufrac_int,kms,kme,kts,kte,dz1d,z1D,dt,u1D,v1D,th1D,rho1D,rhoz1D,q1D,th01D,&
                   tke1D,ustar(ix,iy),hfx(ix,iy),qfx(ix,iy),cp,g, &
                   a_e1D,b_e1D, &
                   dlg1D,dl_u1D,sf1D,vl1D,dlk1D,dls1D,exch1D,sh1D,bu1D,gamma1D)
         do iz= kts,kte
          a_e(ix,iz,iy)=a_e1D(iz)
          b_e(ix,iz,iy)=b_e1D(iz)
          if(flag_bep)then
          dlg_bep(ix,iz,iy)=dlg1D(iz)
          endif
          tke(ix,iz,iy)=tke1D(iz)
          dlk(ix,iz,iy)=dlk1D(iz)
          sh(ix,iz,iy)=sh1D(iz)
          bu(ix,iz,iy)=bu1D(iz)
          exch_h(ix,iz,iy)=exch1D(iz)
          exch_m(ix,iz,iy)=exch1D(iz)
         enddo
         if(idiff.ne.1)then
        if(flag_bep)then
         do iz=kts,kte
          a_t1D(iz)=a_t_bep(ix,iz,iy)
          b_t1D(iz)=b_t_bep(ix,iz,iy)
          a_u1D(iz)=a_u_bep(ix,iz,iy)
          b_u1D(iz)=b_u_bep(ix,iz,iy)
          a_v1D(iz)=a_v_bep(ix,iz,iy)
          b_v1D(iz)=b_v_bep(ix,iz,iy)
          a_q1D(iz)=a_q_bep(ix,iz,iy)
          b_q1D(iz)=b_q_bep(ix,iz,iy)
         enddo
        else
         do iz=kts,kte
          a_t1D(iz)=0.
          b_t1D(iz)=0.
          a_u1D(iz)=0.
          b_u1D(iz)=0.
          a_v1D(iz)=0.
          b_v1D(iz)=0.
          a_q1D(iz)=0.
          b_q1D(iz)=0.
         enddo
          b_t1D(1)=hfx(ix,iy)/dz1D(1)/rho1D(1)/cp
          b_q1D(1)=qfx(ix,iy)/dz1D(1)/rho1D(1)
          a_u1D(1)=(-ustar(ix,iy)*ustar(ix,iy)/dz1D(1)/((u1D(1)**2.+v1D(1)**2.)**.5))
          a_v1D(1)=(-ustar(ix,iy)*ustar(ix,iy)/dz1D(1)/((u1D(1)**2.+v1D(1)**2.)**.5))
        endif
        do iz=kts+1,kte
         if(z1D(iz).le.1.0*pblh(ix,iy).and.wts.gt.0.)then
          b_t1D(iz)=b_t1D(iz)-(exch1D(iz+1)*gamma1D(iz+1)-exch1D(iz)*gamma1D(iz))/dz1D(iz)
         endif
        enddo
       call diff(kms,kme,kts,kte,1,1,dt,u1D,rho1D,rhoz1D,exch1D,a_u1D,b_u1D,sf1D,vl1D,dz1D,wu1D)
       call diff(kms,kme,kts,kte,1,1,dt,v1D,rho1D,rhoz1D,exch1D,a_v1D,b_v1D,sf1D,vl1D,dz1D,wv1D)
       call diff(kms,kme,kts,kte,1,1,dt,th1D,rho1D,rhoz1D,exch1D,a_t1D,b_t1D,sf1D,vl1D,dz1D,wt1D)
       call diff(kms,kme,kts,kte,1,1,dt,q1D,rho1D,rhoz1D,exch1D,a_q1D,b_q1D,sf1D,vl1D,dz1D,wq1D)
       call diff(kms,kme,kts,kte,1,1,dt,qc1D,rho1D,rhoz1D,exch1D,a_qc1D,b_qc1D,sf1D,vl1D,dz1D,wqc1D)
         do iz= kts,kte
          rthblten(ix,iz,iy)=rthblten(ix,iz,iy)+(th1D(iz)-th_phy(ix,iz,iy))/dt
          rqvblten(ix,iz,iy)=rqvblten(ix,iz,iy)+(q1D(iz)-qv_curr(ix,iz,iy))/dt
          rqcblten(ix,iz,iy)=rqcblten(ix,iz,iy)+(qc1D(iz)-qc_curr(ix,iz,iy))/dt
          rublten(ix,iz,iy)=rublten(ix,iz,iy)+(u1D(iz)-u_phy(ix,iz,iy))/dt
          rvblten(ix,iz,iy)=rvblten(ix,iz,iy)+(v1D(iz)-v_phy(ix,iz,iy))/dt
          wu(ix,iz,iy)=wu1D(iz)
          wv(ix,iz,iy)=wv1D(iz)
          wt(ix,iz,iy)=wt1D(iz)
          wq(ix,iz,iy)=wq1D(iz)
        enddo
      endif
      enddo
      enddo
      return
      end subroutine boulac
      subroutine boulac1D(ix,iy,ufrac_int,kms,kme,kts,kte,dz,z,dt,u,v,th,rho,rhoz,qa,th0,te, &
                   ustar,hfx,qfx,cp,g, &
                   a_e,b_e, &
                   dlg,dl_u,sf,vl,dlk,dls,exch,sh,bu,gamma)
      implicit none
      integer iz,ix,iy
      integer kms,kme,kts,kte
      real z(kms:kme)
      real dz(kms:kme)
      real u(kms:kme)
      real v(kms:kme)
      real th(kms:kme)
      real rho(kms:kme)
      real g
      real cp
      real te(kms:kme)
      real qa(kms:kme)
      real th0(kms:kme)
      real dt
      real ustar
      real hfx
      real qfx
      real sf(kms:kme)
      real vl(kms:kme)
      real a_e(kms:kme)
      real b_e(kms:kme)
      real dlg(kms:kme)
      real dl_u(kms:kme)
      real ufrac_int
      real we(kms:kme),dwe(kms:kme)
      real sh(kms:kme)
      real bu(kms:kme)
      real gamma(kms:kme)
      real td(kms:kme)
      real exch(kms:kme)
      real dls(kms:kme)
      real dlk(kms:kme)
      real dlu(kms:kme)
      real dld(kms:kme)
      real rhoz(kms:kme)
      real tstar
      real beta
      real summ1,summ2,summ3,summ4
        tstar=-hfx/rho(1)/cp/ustar
       call dissip_bougeault(ix,iy,g,kms,kme,kts,kte,z,dz,te,dlu,dld,th,th0)
       call length_bougeault(ix,iy,kms,kme,kts,kte,dld,dlu,dlg,dl_u,dls,dlk)
       call cdtur_bougeault(ix,iy,kms,kme,kts,kte,te,z,dz,exch,dlk)
       call tke_bougeault(ix,iy,g,kms,kme,kts,kte,z,dz,vl,u,v,th,te,th0,ustar,tstar,exch,dls,td,sh,bu,gamma,b_e,a_e,sf,ufrac_int)
       call diff(kms,kme,kts,kte,1,1,dt,te,rho,rhoz,exch,a_e,b_e,sf,vl,dz,we)
       do iz=kts,kte
        if(te(iz).lt.temin) te(iz)=temin
       enddo
       return
       end subroutine boulac1d
         subroutine dissip_bougeault(ix,iy,g,kms,kme,kts,kte,z,dz,te,dlu,dld,th,th0)
         implicit none
         integer kms,kme,kts,kte,iz,izz,ix,iy
         real dzt,zup,beta,zup_inf,bbb,tl,zdo,zdo_sup,zzz,g
         real te(kms:kme),dlu(kms:kme),dld(kms:kme),dz(kms:kme)
         real z(kms:kme),th(kms:kme),th0(kms:kme)
         do iz=kts,kte
          zup=0.
          dlu(iz)=z(kte+1)-z(iz)-dz(iz)/2.
          zzz=0.
          zup_inf=0.
          beta=g/th0(iz)
          do izz=iz,kte-1
           dzt=(dz(izz+1)+dz(izz))/2.
           zup=zup-beta*th(iz)*dzt
           zup=zup+beta*(th(izz+1)+th(izz))*dzt/2.
           zzz=zzz+dzt
           if(te(iz).lt.zup.and.te(iz).ge.zup_inf)then
            bbb=(th(izz+1)-th(izz))/dzt
            if(bbb.ne.0)then
             tl=(-beta*(th(izz)-th(iz))+sqrt( max(0.,(beta*(th(izz)-th(iz)))**2.+2.*bbb*beta*(te(iz)-zup_inf))))/bbb/beta
            else
             if(th(izz).ne.th(iz))then
              tl=(te(iz)-zup_inf)/(beta*(th(izz)-th(iz)))
             else
              tl=0.
             endif
            endif
            dlu(iz)=zzz-dzt+tl
           endif
           zup_inf=zup
          enddo
          zdo=0.
          zdo_sup=0.
          dld(iz)=z(iz)+dz(iz)/2.
          zzz=0.
          do izz=iz,kts+1,-1
           dzt=(dz(izz-1)+dz(izz))/2.
           zdo=zdo+beta*th(iz)*dzt
           zdo=zdo-beta*(th(izz-1)+th(izz))*dzt/2.
           zzz=zzz+dzt
           if(te(iz).lt.zdo.and.te(iz).ge.zdo_sup)then
            bbb=(th(izz)-th(izz-1))/dzt
            if(bbb.ne.0.)then
             tl=(beta*(th(izz)-th(iz))+sqrt( max(0.,(beta*(th(izz)-th(iz)))**2.+2.*bbb*beta*(te(iz)-zdo_sup))))/bbb/beta
            else
             if(th(izz).ne.th(iz))then
              tl=(te(iz)-zdo_sup)/(beta*(th(izz)-th(iz)))
             else
              tl=0.
             endif
            endif
            dld(iz)=zzz-dzt+tl
           endif
           zdo_sup=zdo
          enddo
         enddo
         end subroutine dissip_bougeault
         subroutine length_bougeault(ix,iy,kms,kme,kts,kte,dld,dlu,dlg,dl_u,dls,dlk)
         implicit none
         integer kms,kme,kts,kte,iz,ix,iy
         real dlu(kms:kme),dld(kms:kme),dl_u(kms:kme)
         real dls(kms:kme),dlk(kms:kme),dlg(kms:kme)
         do iz=kts,kte
          dld(iz)=min(dld(iz),dlg(iz))
          dls(iz)=sqrt(dlu(iz)*dld(iz))
          dlk(iz)=min(dlu(iz),dld(iz))
         if(dl_u(iz).gt.0.)then
           dls(iz)=1./(1./dls(iz)+1./dl_u(iz))
           dlk(iz)=1./(1./dlk(iz)+1./dl_u(iz))
          endif
         enddo
         return
         end subroutine length_bougeault
       subroutine cdtur_bougeault(ix,iy,kms,kme,kts,kte,te,z,dz,exch,dlk)
       implicit none
       integer iz,kms,kme,kts,kte,ix,iy
       real te_m,dlk_m
       real te(kms:kme),exch(kms:kme)
       real dz(kms:kme),z(kms:kme)
       real dlk(kms:kme)
       real fact
       exch(kts)=0.
       do iz=kts+1,kte
        te_m=(te(iz-1)*dz(iz)+te(iz)*dz(iz-1))/(dz(iz)+dz(iz-1))
        dlk_m=(dlk(iz-1)*dz(iz)+dlk(iz)*dz(iz-1))/(dz(iz)+dz(iz-1))
        exch(iz)=ck_b*dlk_m*sqrt(te_m)
        exch(iz)=max(exch(iz),0.1)
       enddo
       exch(kte+1)=0.1
       return
       end subroutine cdtur_bougeault
       subroutine diff(kms,kme,kts,kte,iz1,izf,dt,co,rho,rhoz,cd,aa,bb,sf,vl,dz,fc)
        implicit none
        integer iz,iz1,izf
        integer kms,kme,kts,kte
        real dt,dzv
        real co(kms:kme),cd(kms:kme),dz(kms:kme)
        real rho(kms:kme),rhoz(kms:kme)
        real cddz(kms:kme+1),fc(kms:kme),df(kms:kme)
        real a(kms:kme,3),c(kms:kme)
        real sf(kms:kme),vl(kms:kme)
        real aa(kms:kme),bb(kms:kme)
        cddz(kts)=sf(kts)*rhoz(kts)*cd(kts)/dz(kts)
        do iz=kts+1,kte
         cddz(iz)=2.*sf(iz)*rhoz(iz)*cd(iz)/(dz(iz)+dz(iz-1))
        enddo
        cddz(kte+1)=sf(kte+1)*rhoz(kte+1)*cd(kte+1)/dz(kte)
         do iz=kts,iz1-1
          a(iz,1)=0.
          a(iz,2)=1.
          a(iz,3)=0.
          c(iz)=co(iz)
         enddo
          do iz=iz1,kte-izf
           dzv=vl(iz)*dz(iz)
           a(iz,1)=-cddz(iz)*dt/dzv/rho(iz)
           a(iz,2)=1+dt*(cddz(iz)+cddz(iz+1))/dzv/rho(iz)-aa(iz)*dt
           a(iz,3)=-cddz(iz+1)*dt/dzv/rho(iz)
           c(iz)=co(iz)+bb(iz)*dt
          enddo
          do iz=kte-(izf-1),kte
           a(iz,1)=0.
           a(iz,2)=1
           a(iz,3)=0.
           c(iz)=co(iz)
          enddo
          call invert (kms,kme,kts,kte,a,c,co)
          do iz=kts,iz1
           fc(iz)=0.
          enddo
          do iz=iz1+1,kte
           fc(iz)=-(cddz(iz)*(co(iz)-co(iz-1)))/rho(iz)
          enddo
       return
       end subroutine diff
       subroutine buoy(ix,iy,g,kms,kme,kts,kte,th,th0,exch,dz,bu,gamma,ustar,tstar,ufrac_int)
       implicit none
       integer kms,kme,kts,kte,iz,ix,iy
       real dtdz1,dtdz2,cdm,dtmdz,g
       real th(kms:kme),exch(kms:kme),dz(kms:kme),bu(kms:kme),gamma(kms:kme)
       real th0(kms:kme),ustar,tstar,ufrac_int,gammam
      bu(kts)=0.
       do iz=kts+1,kte-1
        dtdz1=2.*(th(iz)-th(iz-1))/(dz(iz-1)+dz(iz))
        dtdz2=2.*(th(iz+1)-th(iz))/(dz(iz+1)+dz(iz))
        dtmdz=0.5*(dtdz1+dtdz2)
        cdm=0.5*(exch(iz+1)+exch(iz))
        gammam=0.5*(gamma(iz+1)+gamma(iz))
        bu(iz)=-cdm*(dtmdz-gammam)*g/th0(iz)
       enddo
       bu(kte)=0.
       return
       end subroutine buoy
       subroutine shear(ix,iy,g,kms,kme,kts,kte,u,v,cdua,dz,sh,ustar,tstar,th,ufrac_int)
       implicit none
       integer kms,kme,kts,kte,iz,ix,iy
       real dudz1,dudz2,dvdz1,dvdz2,cdm,dumdz,ustar
       real tstar,th,al,phim,g
       real u(kms:kme),v(kms:kme),cdua(kms:kme),dz(kms:kme),sh(kms:kme)
       real u1,u2,v1,v2,ufrac_int
       sh(kts)=0.
       do iz=kts+1,kte-1
        u2=(dz(iz+1)*u(iz)+dz(iz)*u(iz+1))/(dz(iz)+dz(iz+1))
        u1=(dz(iz)*u(iz-1)+dz(iz-1)*u(iz))/(dz(iz-1)+dz(iz))
        v2=(dz(iz+1)*v(iz)+dz(iz)*v(iz+1))/(dz(iz)+dz(iz+1))
        v1=(dz(iz)*v(iz-1)+dz(iz-1)*v(iz))/(dz(iz-1)+dz(iz))
        cdm=0.5*(cdua(iz)+cdua(iz+1))
        dumdz=((u2-u1)/dz(iz))**2.+((v2-v1)/dz(iz))**2.
        sh(iz)=cdm*dumdz
       enddo
       sh(kte)=0.
       return
       end subroutine shear
       subroutine invert(kms,kme,kts,kte,a,c,x)
       implicit none
       integer in
       integer kts,kte,kms,kme
       real a(kms:kme,3),c(kms:kme),x(kms:kme)
        do in=kte-1,kts,-1
         c(in)=c(in)-a(in,3)*c(in+1)/a(in+1,2)
         a(in,2)=a(in,2)-a(in,3)*a(in+1,1)/a(in+1,2)
        enddo
        do in=kts+1,kte
         c(in)=c(in)-a(in,1)*c(in-1)/a(in-1,2)
        enddo
        do in=kts,kte
         x(in)=c(in)/a(in,2)
        enddo
        return
        end subroutine invert
       subroutine tke_bougeault(ix,iy,g,kms,kme,kts,kte,z,dz,vl,u,v,th,te,th0,ustar,tstar,exch, &
                         dls,td,sh,bu,gamma,b_e,a_e,sf,ufrac_int)
       implicit none
       integer kms,kme,kts,kte,iz,ix,iy
       real g,ustar,tstar,ufrac_int
       real z(kms:kme),dz(kms:kme),u(kms:kme),v(kms:kme),th(kms:kme),th0(kms:kme),te(kms:kme)
       real exch(kms:kme),dls(kms:kme),td(kms:kme),sh(kms:kme),bu(kms:kme),gamma(kms:kme)
       real a_e(kms:kme),b_e(kms:kme)
       real vl(kms:kme),sf(kms:kme)
       real te1,dl1
       call shear(ix,iy,g,kms,kme,kts,kte,u,v,exch,dz,sh,ustar,tstar,th(kts),ufrac_int)
       call buoy(ix,iy,g,kms,kme,kts,kte,th,th0,exch,dz,bu,gamma,ustar,tstar,ufrac_int)
       do iz=kts,kte
        te1=max(te(iz),temin)
        dl1=max(dls(iz),0.1)
        td(iz)=-ceps_b*sqrt(te1)/dl1
        sh(iz)=sh(iz)*sf(iz)
        bu(iz)=bu(iz)*sf(iz)
        a_e(iz)=a_e(iz)+td(iz)
        b_e(iz)=b_e(iz)+sh(iz)+bu(iz)
       enddo
       return
       end subroutine tke_bougeault
      SUBROUTINE BOULACINIT(RUBLTEN,RVBLTEN,RTHBLTEN,RQVBLTEN,RQCBLTEN, &
     & TKE_PBL,EXCH_H,RESTART,ALLOWED_TO_READ, &
     & IDS,IDE,JDS,JDE,KDS,KDE, &
     & IMS,IME,JMS,JME,KMS,KME, &
     & ITS,ITE,JTS,JTE,KTS,KTE )
      IMPLICIT NONE
      LOGICAL,INTENT(IN) :: ALLOWED_TO_READ,RESTART
      INTEGER,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE, &
     & IMS,IME,JMS,JME,KMS,KME, &
     & ITS,ITE,JTS,JTE,KTS,KTE
      REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(OUT) :: EXCH_H, &
     & RUBLTEN, &
     & RVBLTEN, &
     & RTHBLTEN, &
     & RQVBLTEN, &
     & RQCBLTEN, &
     & TKE_PBL
      INTEGER :: I,J,K,ITF,JTF,KTF
      JTF=MIN0(JTE,JDE-1)
      KTF=MIN0(KTE,KDE-1)
      ITF=MIN0(ITE,IDE-1)
      IF(.NOT.RESTART)THEN
        DO J=JTS,JTF
        DO K=KTS,KTF
        DO I=ITS,ITF
          TKE_PBL(I,K,J)=0.0001
          RUBLTEN(I,K,J)=0.
          RVBLTEN(I,K,J)=0.
          RTHBLTEN(I,K,J)=0.
          RQVBLTEN(I,K,J)=0.
          RQCBLTEN(I,K,J)=0.
          EXCH_H(I,K,J)=0.
        ENDDO
        ENDDO
        ENDDO
      ENDIF
      END SUBROUTINE BOULACINIT
       subroutine pbl_height(kms,kme,kts,kte,dz,z,th,q,pblh)
       implicit none
       integer kms,kme,kts,kte,iz
       real z(kms:kme),dz(kms:kme),th(kms:kme),q(kms:kme)
       real pblh
       real thv(kms:kme),zc(kms:kme)
       real thsfc
      do iz=kts,kte
       zc(iz)=z(iz)+dz(iz)/2.
      enddo
       do iz=kts,kte
        thv(iz)=th(iz)*(1.+0.61*q(iz))
       enddo
       pblh=0.
       thsfc=thv(kts)+0.5
       do iz=kts+1,kte
        if(pblh.eq.0.and.thv(iz).gt.thsfc)then
         pblh=zc(iz-1)+(thsfc-thv(iz-1))/(max(0.01,thv(iz)-thv(iz-1)))*(zc(iz)-zc(iz-1))
        endif
       enddo
       return
       end subroutine pbl_height
END MODULE module_bl_boulac


text:
program visual_tests

    use iso_fortran_env
    use sobseq

    integer, parameter :: n_dim = 15
    integer, parameter :: n_samples = 10

    type(multi_dim_sobol_state) :: md_rng
    real(real64), dimension(:,:), allocatable :: seq

    allocate(seq(n_samples, n_dim))

    md_rng = multi_dim_sobol_state(n_dim, "data/clean_new-joe-kuo-6.21201")
    call md_rng%populate(seq)

    print '(20F8.4)', transpose(seq)

end program visual_tests

max_stars_repo_name: renatomatz/sobseq

max_issues_repo_name: renatomatz/sobseq

max_forks_repo_name: renatomatz/sobseq


text:
      PROGRAM MAIN

      END

C comment
* comment

      SUBROUTINE foo( i, x, b )
      INTEGER            i
      REAL               x
      LOGICAL            b

      IF( i.NE.0 ) THEN
         CALL bar( -i )
      END IF

      RETURN
      END

      DOUBLE COMPLEX FUNCTION baz()

      baz = (0.0d0,0.0d0)

      RETURN
      END

ext: f

max_stars_repo_name: JavascriptID/sourcerer-app

max_issues_repo_name: JavascriptID/sourcerer-app

max_forks_repo_name: JavascriptID/sourcerer-app


text:
Debats du Senat (hansard)
1ere Session, 36e Legislature,
Volume 137, Numero 28
Le mardi 9 decembre 1997
L'honorable Gildas L. Molgat, President
Le programme d'echange de pages avec la Chambre des communes
Reponse a un article de journal sur les conditions de residence applicables a un senateur
L'histoire du vote au Canada
Le lancement du livre par le Gouverneur general
Loi de 1997 pour la mise en oeuvre de conventions fiscales
L'etat actuel et les perspectives d'avenir de l'agriculture
Depot du rapport du comite de l'agriculture et des forets
L'etat actuel et les perspectives d'avenir des forets
Depot du rapport du comite de l'agriculture et des forets
Regie interne, budgets et administration
Presentation du septieme rapport du comite
Assemblee internationale des parlementaires de langue francaise
Reunion tenue au Luxembourg-Depot du rapport de la section canadienne
Autorisation au comite de sieger en meme temps que le Senat
La conduite des travaux-Avis d'interpellation
Reduction du transfert des paiements a la province de Quebec-Demande de precisions
Les centres de detention des revendicateurs du statut de refugie-Demande de precisions
Reponses differees a des questions orales
La Declaration universelle des droits de l'homme
La commemoration du 50e anniversaire-Les plans du gouvernement
La demolition de laboratoires gouvernementaux-Le retablissement possible du financement-La position du gouvernement
Depot de reponses a des questions au Feuilleton
La defense-L'etat du projet Habillement du soldat
La defense-La situation du programme de remplacement des transports de troupes blindes
Les travaux du Senat
Deuxieme lecture-Suite du debat
Motion de sous-amendement
Regie interne, budgets et administration
Adoption du sixieme rapport du comite
Projet de loi modificatif - Premiere lecture
Le mardi 9 decembre 1997
La seance est ouverte a 14 heures, le President etant au fauteuil.
Le programme d'echange de pages avec la Chambre des communes
Son Honneur le President :
Il se specialise en psychologie.
Au nom de tous les senateurs et du Senat, je vous souhaite la bienvenue.
Reponse a un article de journal sur les conditions de residence applicables a un senateur
J'aimerais signaler quelques faits a ce propos:
Ce n'est pas vrai.
Je n'en ai pas parle au personnel juridique il y a cinq ans.
Je n'ai jamais eu de probleme en ce qui concerne mon lieu de residence.
Le senateur Lucier a souligne le mot dit.
Je me demande si Jack Aubry doute du fait que j'aie un cancer.
4. J'aimerais que ces points soient clarifies, mais pas dans les medias.
5. Je remercie les collegues qui m'ont exprime leur soutien a ce sujet.
L'histoire du vote au Canada
Le lancement du livre par le Gouverneur general
L'honorable Mabel M. DeWare :
Le mardi 9 decembre 1997
Son Honneur le President :
Honorables senateurs, quand lirons-nous ce projet de loi une troisieme fois?
Loi de 1997 pour la mise en oeuvre de conventions fiscales
Le mardi 9 decembre 1997
Le comite senatorial permanent des banques et du commerce a l'honneur de presenter son
Son Honneur le President :
Honorables senateurs, quand lirons-nous ce projet de loi une troisieme fois?
L'etat actuel et les perspectives d'avenir de l'agriculture
Depot du rapport du comite de l'agriculture et des forets
L'honorable Leonard J. Gustafson :
Son Honneur le President :
Est-ce d'accord, honorables senateurs?
Son Honneur le President :
Honorables senateurs, quand etudierons-nous le rapport?
L'etat actuel et les perspectives d'avenir des forets
Depot du rapport du comite de l'agriculture et des forets
L'honorable Leonard J. Gustafson :
Son Honneur le President :
Honorables senateurs, quand etudierons-nous le rapport?
Regie interne, budgets et administration
Presentation du septieme rapport du comite
Le mardi 9 decembre 1997
Son Honneur le President :
Honorables senateurs, quand etudierons-nous le rapport?
L'honorable Sharon Carstairs (leader adjoint du gouvernement) :
Son Honneur le President :
La permission est-elle accordee, honorables senateurs?
Assemblee internationale des parlementaires de langue francaise
Reunion tenue au Luxembourg-Depot du rapport de la section canadienne
L'honorable Pierre De Bane :
Autorisation au comite de sieger en meme temps que le Senat
L'honorable John B. Stewart :
Son Honneur le President :
Permission est-elle accordee, honorables senateurs?
L'honorable Noel A. Kinsella (chef adjoint suppleant de l'opposition) :
Honorables senateurs, j'ai deux observations a faire.
Attention aux jets de gaz poivre.
Pour resumer, je dirai que la question est importante pour tous les senateurs.
Le senateur Stewart nous a explique que c'etait sa principale preoccupation.
Son Honneur le President :
Permission accordee, honorables senateurs?
Cela depend de la reponse.
L'honorable Sharon Carstairs (leader adjoint du gouvernement) :
Honorables senateurs, nous devons faire un effort pour que les comites puissent sieger.
Ils invitent des temoins qui viennent parfois de l'exterieur de la ville.
Nous leur occasionnons des inconvenients si nous ne pouvons pas ecouter leur temoignage.
Ils ont essaye de faire leurs interventions le mardi et le jeudi.
Les senateurs d'en face ont tente de faire la meme chose.
L'honorable John Lynch-Staunton (chef de l'opposition) :
Cela a fonctionne pendant un certain temps.
Cela ne convient pas du tout.
Nous devons toutefois obtenir d'abord cette assurance.
Nous avons siege hier et c'etait un lundi.
Son Honneur le President :
Son Honneur le President :
Son Honneur le President :
Plait-il aux honorables senateurs d'adopter la motion?
La conduite des travaux-Avis d'interpellation
L'honorable William M. Kelly :
Honorables senateurs, je donne avis:
Reduction du transfert des paiements a la province de Quebec-Demande de precisions
Pourtant, les chiffres indiquent autre chose.
Pourtant, hier, le gouvernement a affirme qu'il appliquait un regime de transfert stable.
L'honorable B. Alasdair Graham (leader du gouvernement) :
L'honorable Mabel M. DeWare :
Quant au Nouveau-Brunswick, il perdra 11 millions de dollars.
Ce que le Canada atlantique perd, elles le gagnent.
L'honorable B. Alasdair Graham (leader du gouvernement) :
Je remercie le senateur DeWare d'avoir fait les calculs pour moi.
Je vais certainement porter ces considerations tres legitimes a l'attention du ministre responsable.
L'honorable B. Alasdair Graham (leader du gouvernement) :
Je ne sais pas quelle est la position de chaque province.
Les reunions des fonctionnaires des pays durent depuis quelques jours.
Je remercie le leader du gouvernement au Senat pour sa reponse.
L'honorable B. Alasdair Graham (leader du gouvernement) :
Je ne crois pas qu'il y ait eu d'accord final a ce sujet.
Et j'insiste sur cette derniere partie:
...compte tenu des engagements actuels.
L'honorable B. Alasdair Graham (leader du gouvernement) :
L'etude fournit des donnees tres importantes.
Je suis heureuse d'entendre cela.
Quelle est la position de la delegation canadienne au sujet de ce concept?
L'honorable B. Alasdair Graham (leader du gouvernement) :
L'honorable Pierre Claude Nolin :
Ma premiere question vise a connaitre ces conditions auxquelles le premier ministre faisait reference?
L'honorable B. Alasdair Graham (leader du gouvernement) :
Les Quebecois devront exprimer clairement leur volonte de quitter le Canada.
En outre, tout processus menant a la secession devrait evidemment respecter la loi.
Que dit la loi au sujet de la secession?
Le premier ministre nous a dit qu'il etait pret a negocier.
C'est ce que nous avons toujours dit.
Gardez bien a l'esprit que ce ministre est separatiste!
Je repete ses paroles:
C'est ce que nous avons toujours dit.
M. Chretien a admis que le gouvernement federal negocierait.
Je suis sur que le premier ministre tiendra compte de tout cela.
Pas a ma connaissance, honorables senateurs.
L'honorable Donald H. Oliver :
L'honorable B. Alasdair Graham (leader du gouvernement) :
Les centres de detention des revendicateurs du statut de refugie-Demande de precisions
L'honorable Donald H. Oliver :
L'honorable B. Alasdair Graham (leader du gouvernement) :
L'honorable Gerald J. Comeau :
L'honorable B. Alasdair Graham (leader du gouvernement) :
Honorables senateurs, je ne suis pas certain de la position actuelle.
La reponse, honorables senateurs, est oui, nous consultons l'industrie de la peche.
Je pense que sa position est exactement celle qu'a suggeree mon honorable collegue.
C'est une autre chose que je vais tirer au clair pour lui.
L'honorable J. Michael Forrestall :
L'honorable B. Alasdair Graham (leader du gouvernement) :
Le gouvernement s'occupe activement du dossier des helicopteres.
Le Cabinet s'est reuni jeudi dernier.
Le Cabinet se reunira encore, dans ce cas.
Je suis heureux de l'apprendre.
L'honorable J. Michael Forrestall :
Ou sont les six autres?
Sont-ils sur des navires en mer?
Nous en avons 30 en stock.
L'honorable B. Alasdair Graham (leader du gouvernement) :
Honorables senateurs, oui, nous avons 30 helicopteres en stock.
Comme je l'ai dit hier, 24 helicopteres ont ete autorises a voler.
Il en reste six a inspecter, et ils le seront d'ici le 18 decembre.
Je ferai part de la reponse a mon collegue.
Avons-nous ete critiques pour ne pas avoir ete en mesure de respecter nos obligations?
Pas que je sache.
Je crois savoir que ces helicopteres participent a des exercices en mer.
Je ne sais pas exactement ou, ni dans quel ocean.
Il se peut meme que ce soit un secret militaire.
L'honorable Brenda M. Robertson :
Je reviendrai sur ce point demain et apres-demain.
L'honorable B. Alasdair Graham (leader du gouvernement) :
Non, pour autant que je sache, mais la question est interessante.
Je trouverai la reponse.
Reponses differees a des questions orales
L'honorable Sharon Carstairs (leader adjoint du gouvernement) :
( Reponse a la question posee par l'honorable Donald H. Oliver le 25 novembre 1997 )
Pour une caisse de retraite, le rendement des investissements a long terme est primordial.
( Reponse a la question posee par l'honorable Donald H. Oliver le 19 novembre 1997 )
Par exemple, l'Office sera tenu de:
divulguer la remuneration versee a cinq de ces cinqcadres les mieux payes du Conseil
La Declaration universelle des droits de l'homme
La commemoration du 50e anniversaire-Les plans du gouvernement
( Reponse a la question posee par l'honorable Noel A. Kinsella le 23 octobre 1997 )
( Reponse a la question posee par l'honorable J. Michael Forrestall le 20 novembre 1997 )
La demolition de laboratoires gouvernementaux-Le retablissement possible du financement-La position du gouvernement
( Reponse a la question posee par l'honorable Mira Spivak le 20 novembre 1997 )
Depot de reponses a des questions au Feuilleton
La defense-L'etat du projet Habillement du soldat
La defense-La situation du programme de remplacement des transports de troupes blindes
Les travaux du Senat
Son Honneur le President :
Projet de loi sur l'Office d'investissement du regime de pensions du Canada
Deuxieme lecture-Suite du debat
Maintenant, il nous faut en terminer l'etude et l'adopter avant Noel.
Nous ne sommes pas des servantes.
Nous n'avons d'ordres a recevoir de personne.
Nous devons faire ce que nous croyons etre juste.
Leurs prestations seront legerement reduites et leurs cotisations legerement augmentees.
Je vais vous dire pour qui c'est important.
J'ai une fille de 24 ans et un fils de 21 ans.
Nombre d'entre vous avez aussi des enfants et meme des petits-enfants.
L'inflation aidant, les chiffres grossissent, mais pas celui-la.
Il reste le meme.
Elisez les liberaux et il sera eleve!
Si nous elisons les liberaux, l'inflation sera elevee?
Nous le savons deja!
Permettez-moi de poursuivre dans le meme ordre d'idees.
Je ne comprends pas.
Je tenterai d'expliquer cela aujourd'hui, et j'espere pouvoir vous convaincre.
Le RPC fait partie d'un programme plus vaste que Paul Martin presente.
Nous croyons fermement que cette mesure est une charge sociale qui fait perdre des emplois.
Hier, le senateur Kirby a dit que ce n'etait pas une charge sociale.
Honorables senateurs, il y a des accents de verite dans toute chose vilaine.
C'est pourquoi nous parlons de propagande.
Ses consultations avaient simplement une valeur symbolique.
Ce n'est pas une question de la fin qui justifie les moyens.
Dans ce projet de loi, les moyens et la fin sont egalement reprehensibles.
Cela devient une habitude, une habitude dont il est difficile de se defaire.
Il a utilise trois moyens pour parvenir a cette fin.
Le gouvernement a transfere la dette.
Il n'appartient pas aux administrateurs de donner ou de refuser leur accord.
Dans ce cas-ci, il echoue.
Il utilise alors l'argent a une autre fin.
Voyez le projet de loi.
Il releve simplement les cotisations et reduit les prestations.
N'importe qui peut faire ca.
Rappelons-nous que l'auteur du projet de loi C-2 est Paul Martin.
Paul Martin est un homme brillant.
Il s'est ensuite excuse aupres des Canadiens de les avoir mal renseignes.
Il leur est tous interdit de voler.
Il a ecrit le livre rouge pour pouvoir renegocier l'Accord de libre-echange.
Il a ete couvert de gloire apres les elections de 1993.
Je n'ai jamais entendu parler de cette entente.
Personne n'en a entendu parler.
Faisons un cadeau de Noel aux Canadiens en rejetant ce projet de loi.
Tuons ce projet de loi.
Nous connaissons trop bien le bilan des liberaux.
Je vous ai parle des promesses que le gouvernement a faites.
Nous connaissons tous son sens des responsabilites.
Je peux imaginer la scene.
Lorsque cette mesure legislative aura seme la pagaille, Paul Martin aura 95 ans.
Jean Chretien aura 95, 96 ou 102 ans.
Ils seront encore au pouvoir.
Il ne reconnaitra jamais sa responsabilite.
Je ne sais pas qui, mais je ne suis pas responsable.
Ils n'ont assume aucune responsabilite dans l'annulation du contrat des EH-101.
Ils en parlent toujours au Cabinet.
Faites-nous confiance, disent-ils.
Ils ont passe un contrat social avec leurs administres.
C'est eux qui l'ont fait.
Nous allons prendre votre argent et vous donner une pension, avaient-ils dit.
Nous avons l'audace d'atteindre 50 ans!
Ils ont limite leur progeniture.
Disons qu'ils batifolaient.
Ce n'etait pas leur faute, mais celle des scientifiques.
La cause, disaient-ils, c'etait un comportement social.
Voila ce qui a cause la ruine de ce regime de retraire.
Voila un autre programme social qui deraille.
Finis les beaux jours.
Nous avons un autre probleme.
Encore une fois, la science est intervenue avec un plan social.
C'etait le regime de pensions parfait.
Des le debut, le RPC etait simple, mais mal nomme.
Il etait mal nomme parce qu'il etait aussi un regime d'assurance-vie.
Cependant, c'est stupide quand cela fait partie d'un regime de retraite.
Il y a aussi un regime d'invalidite.
Je voudrais donner un apercu historique aux senateurs.
Beaucoup de modifications ont ete apportees au RPC au fil des ans.
Pourtant, le gouvernement a dit qu'il n'y toucherait pas.
Toutefois, il a modifie l'autre fonds.
Depuis 1973, il adopte de nouvelles lois et cree de nouveaux programmes.
Nous avons vu cela aujourd'hui quand nous avons etudie le budget des depenses.
Le gouvernement dit que nous l'avons fait quand nous etions au pouvoir.
Les criteres concernant la retraite et les gains des individus ont ete elimines en 1974.
Il n'y avait pas autant de personnes agees.
Il y avait enormement de gens qui travaillaient et beaucoup d'argent.
Le gouvernement s'est dit qu'il pouvait depenser allegrement ces fonds.
Maintenant, si l'on a plus de quatre enfants, c'est tres bien aussi.
Les fonds etaient la, le gouvernement a donc puise dedans.
Un projet de loi a ete adopte a propos des credits pour conjoint.
Les gens divorcent, et la loi doit evoluer.
Les gens qui divorcent peuvent partager les prestations.
Ce n'etait pas possible autrefois.
Au debut, ce n'etait pas permis.
Cela coute plus cher.
Je ne dis pas que ce n'est pas bien.
Il faut de la force et de l'intelligence pour cela.
Voulez-vous dire que nous ne sommes pas intelligents?
Pas autant que les femmes.
Les prestations ont ensuite ete elargies.
Je sais ce que cela signifie.
C'est tres bien, si le gouvernement veut faire cela.
Cela coute de l'argent.
Mais il se pourrait que deux hommes vivent ensemble un an sans etre gays.
Comment pourra-t-on trancher?
Vous voulez des cameras de surveillance?
Il y a donc une possibilite de fraude.
Ce n'est pas un regime de pensions.
Que ce soit bien clair.
Ce n'etait pas la faute des legislateurs.
Citez-vous le Parti reformiste?
Il a dit ceci:
Il parlait des ministres Martin et Pettigrew.
C'est un grand changement par rapport a la derniere campagne electorale.
[...] un programme de REER elargi; enfin, des allegements fiscaux pour les personnes agees.
Je crois que j'ai explique certaines des idees venant du gouvernement.
Ils veulent maintenant un plus petit RPC.
Toutefois, Preston Manning n'a jamais pu expliquer le premier programme.
Eh bien, on n'a jamais entendu ces explications.
Je peux comprendre pourquoi.
Les reformistes font de fausses accusations a la Chambre au sujet de diverses personnes.
C'est le programme reformiste typique auquel nous sommes habitues.
Vous vous souviendrez qu'ils ont propose un plan fiscal fonde sur un impot uniforme.
Avec un impot uniforme, on reduit les exemptions.
Les reformistes ne nous donnent pas d'idees.
Je me tourne maintenant vers les liberaux.
Le gouvernement ne veut pas que ce programme-ci soit connu.
Ce n'est pas vrai.
C'est la raison pour laquelle un accord a ete conclu.
Nous avons tente de le faire en 1986.
Vous n'avez pas reussi.
Parce que les liberaux refusaient d'accorder leur appui.
Je defendrai mon gouvernement tout comme l'honorable senateur defendra le sien.
Vous lisez le mauvais texte.
De toute evidence, l'honorable senateur n'a pas de notes.
Je crois que je commence a irriter les honorables senateurs d'en face.
Honorables senateurs, ce n'est pas nous qui paierons pour cela.
Qu'y a-t-il qui les interesse la-dedans?
Qu'y a-t-il qui interesse le ministre des Finances?
Comment savons-nous qu'ils ont raison?
Pourquoi ne pas prendre notre temps, honorables senateurs?
Maintenant que nous vous avons entendu, nous savons qu'ils ont raison.
Pourquoi ne pas prendre le temps...
Je voudrais aborder d'autres aspects du probleme.
Le gouvernement se propose, honorables senateurs, d'amender le regime de pensions du Canada.
J'aurai des choses a dire la-dessus tantot.
C'est la, me semble-t-il, une question importante a considerer.
A cette question de base, la reponse du gouvernement est la collectivisation obligatoire.
Voila un cas bien classique de mismanagement public.
Et maintenant, on n'apprend rien de cette triste histoire:
Or, comment evaluer la performance d'un monopole?
La liste est encore longue des folies des monopolistes.
Vous me direz que ces choses se produisent aussi dans le prive.
Je parle de l'Office d'investissement.
Vous rendez-vous compte des pouvoirs qu'auront ces gens?
Trois milliards de dollars, ce sont 3 000 millions de dollars.
Vous voyez cela, 600 millions de dollars theoriquement a chaque jour?
Pensez-vous que la grosse machine requise pour operer cela sera efficace?
L'entreprise privee sait cela.
Il semble qu'il n'y ait que le gouvernement qui ne le sache pas.
Je reviendrai la-dessus.
Mon quatrieme point porte sur la politique de placement .
Or, l'economie canadienne equivaut a 2,5 p. 100 de l'economie mondiale.
Nous avons connu cela a Quebec.
La question des objectifs du plan de pension est fondamentale.
A ce propos, il est important d'etre clair dans le statut.
Ce n'etait pas l'argent de M. Caouette, c'etait de vraies piastres!
A quel exercice de tordage de bras cela nous conduira-t-il?
Quand on collectivise une activite, on la politise.
Les contributions vont etre relevees de 70 p. 100 d'ici cinq ou sept ans.
Est-ce une bonne chose d'accroitre le role de l'Etat?
C'est ce que l'on fait ici.
Qui peut imaginer que cela a du sens de penser comme cela?
Dans le passe, qu'est-ce qui s'est produit?
Je constate simplement que cela se passe ainsi depuis 40 ans.
Regardez les resultats des societes d'Etat!
Autrement dit, la bonne gestion a commence quand le gouvernement n'etait plus implique.
Il faut etre prudent avant de transferer l'argent du secteur prive au secteur public.
Je me tourne maintenant vers la question de l'impact du plan sur les emplois.
Son Honneur le President :
Votre periode de 15 minutes est terminee.
La permission de continuer est-elle accordee, honorables senateurs?
Je voudrais que l'on me reponde, c'est une bonne question.
Un autre point d'interrogation concerne la concurrence interprovinciale previsible.
Elles vont emprunter du fonds, mais elles vont payer plus cher.
Or, comment choisir et suivre 800 entreprises?
Arbitraire et patronage risquent de se produire.
Au surplus, l'office sera partout et fixera inevitablement des conditions.
Cette concentration pouvoir, cette centralisation sera malsaine.
Je crains cette technocratie et les effets pervers qu'elle entrainerait.
Les reglements relatifs aux placements et a la conduite des administrateurs ne sont pas connus.
Ils ne sont pas prets.
Quelles seront les balises pour guider les administrateurs?
Une revue parlementaire s'impose donc avant leur mise en vigueur.
Voila une reponse intelligente au defi qui nous confronte.
L'honorable Philippe Deane Gigantes :
Vous parlez avec admiration de la gestion d'investissement privee.
Pour vous, ce n'est pas grave, mais pour les petits investisseurs, c'est grave.
Il y en tous les jours.
Mais c'est le secteur prive.
En deuxieme lecture, les projets de loi sont toujours beaux.
J'en ai vu suffisamment pendant 40 ans.
Les ministres nous promettent des miracles avec leur projet de loi.
Le discours du ministre n'a plus du tout le meme impact apres trois ans.
Personne n'est oblige de mettre de l'argent dans le Royal Trust.
Le public est compose d'individus.
Les individus sont des personnes privees.
C'est une question de choix.
L'honorable Gerry St. Germain :
Son experience est a la disposition de notre forum.
Je vais donc la poser au senateur Bolduc.
Il s'agit d'une question que j'ai posee hier au senateur Kirby.
C'est une question de responsabilite.
Nous en confions pourtant la responsabilite a un groupe de 12 personnes.
Honorables senateurs, je trouve que c'est une tres bonne question.
Dans le projet de loi, une chose est fort acceptable.
Le verificateur general va pouvoir etudier la gestion de l'Office d'investissement.
Par consequent, a travers la verification, il y aura un certain processus comptable.
L'honorable Jerahmiel S. Grafstein :
Honorables senateurs, je voudrais poser une breve question au senateur Bolduc.
C'est une question serieuse.
Honorables senateurs, je ne mets pas en doute la question de la reddition des comptes.
Rendre des comptes, ce n'est pas la meilleure solution.
A mon avis, c'est mieux que rien, mais ce n'est pas l'ideal.
L'honorable Michael A. Meighen :
Honorables senateurs, c'est justement la-dessus que je desirais une precision.
Vous avez parle d'une decentralisation des centres decisionnels, si j'ai bien compris.
Il est sur que pour certaines provinces, ce serait peut-etre plus difficile.
Ils ont recu de l'information qu'ils pouvaient comprendre.
Les investisseurs obtiennent de bons rendements.
Il y a deux lecons a tirer de cet exemple.
C'est une farce.
Pour moi, ce n'est pas le comportement d'un gouvernement responsable.
Par exemple, le 1er octobre, le senateur Michael Meighen posait cette question:
Ma question a trait aux cotisations a l'assurance-emploi.
Et il disait plus loin:
Le 23 octobre, j'ai pose la question suivante:
Il aura un effet multiplicateur.
Vous avez achete les personnes agees et vous faites payer les jeunes.
Ils ne l'oublieront pas.
Le 30 octobre, le senateur Tkachuk a pose la question suivante:
Autrement dit, c'est et ce sera percu comme un autre impot.
Le 5 novembre, le senateur Oliver a dit ce qui suit:
En fait, cela devrait permettre une transition ordonnee.
Comme on l'a deja dit avant moi, les gens vivent plus longtemps.
Le 6 novembre, le senateur Oliver a pose la question suivante:
Combien de Canadiens devront payer le prix de cette augmentation prevue des cotisations?
Son Honneur le President :
Je voudrais, si vous me le permettez, terminer mon discours.
Son Honneur le President :
La permission est-elle accordee a l'honorable senateur Stratton de poursuivre?
Le 18 novembre, le senateur Oliver a dit ceci:
Encore une fois, les Canadiens ne sont pas au courant de ce qui se passe.
Le 20 novembre, j'ai dit:
Je cite le senateur Meighen, qui demandait:
Cela fait 51 jours que nous l'attendons.
Le 25 novembre, j'ai dit:
Il est dit dans ce document:
Le gouvernement ecoutait-il quand j'ai pose la question?
Le senateur Graham a dit: Oui, bien sur.
Aucun amendement n'a ete adopte.
Les personnes agees, dans leur ensemble, ne sont pas au courant de ces changements.
Le senateur Graham a repondu par l'affirmative.
Il a dit que tel etait le plan de communication du gouvernement.
La question portait sur l'assurance-emploi et la reserve.
J'ai poursuivi ainsi:
Les ministres des Finances des provinces se sont reunis a Ottawa.
Je pense qu'ils ont dine hier soir avec le ministre, Paul Martin.
Dans un article du Globe and Mail d'hier, on disait:
Je predis que c'est sur le point d'arriver.
L'honorable Thelma J. Chalifoux :
Quelle est la position du senateur a ce sujet?
Meme dans le domaine des communications, il n'y a pas grand-chose.
Je n'etais pas d'accord avec le Regime de pensions du Canada.
Ils doivent faire 240 kilometres pour se rendre a la banque la plus proche.
Beaucoup ne parlent meme pas anglais.
Que feriez-vous pour les personnes qui n'ont rien?
Je remercie le senateur de sa question.
C'est a vous que je pose la question.
Il est impossible de privatiser immediatement.
Il suffit d'observer ce qui se passe au Chili.
Voila la solution que le gouvernement du Chili a choisie.
La seule facon de communiquer est le contact personnel.
Il faut sortir et aller rencontrer les gens chez eux.
Nous ne sommes pas au pouvoir maintenant.
Vous devriez insister aupres de votre parti.
(Sur la motion du senateur Gigantes, le debat est ajourne.)
(Le projet de loi est lu une premiere fois.)
Son Honneur le President :
Quand lirons-nous ce projet de loi une deuxieme fois?
L'honorable Orville H. Phillips :
Honorables senateurs, je prends la parole pour appuyer l'amendement du senateur Cools.
Il importe de s'interesser a ces preoccupations croissantes.
Je felicite madame le senateur Cools d'avoir presente cet amendement.
Les policiers ont suivi le camion jusqu'au domicile de M. Feeney.
Les policiers ont dit que la porte etait ouverte et qu'ils sont entres.
Ils ont trouve M. Feeney endormi.
Il avait encore du sang dans le visage et sur sa chemise.
Ils ont arrete M. Feeney sans mandat.
Je pense qu'elle possede toujours ce pouvoir.
Il n'est pas inscrit dans le Code criminel.
Il s'agit d'un de nos principes de droit herites du regime britannique.
Comme le senateur Cools l'a fait remarquer, aucune loi n'a ete violee.
Honorables senateurs, la Cour supreme n'a pas tenu compte des droits de M. Boyle.
D'apres eux, les mesures prises par la police ne modifiaient en rien la decision.
Mais, la encore, ce n'est que mon interpretation.
Je souligne encore une fois que M. Boyle avait le droit a sa vie.
C'est sa fonction, et le Parlement l'en a saisie.
Cependant, la promulgation des lois demeure la fonction du parlement.
Le Parlement est le plus haut tribunal du pays et doit toujours le rester.
Beaucoup d'elements dans son amendement valent la peine d'etre pris en consideration.
Elle a cependant oublie un aspect de notre projet de loi.
Voila pourquoi nous avons la sanction royale.
Motion de sous-amendement
L'honorable Orville H. Phillips :
Son Honneur le President :
Honorables senateurs, j'ai un probleme.
De prime abord, je ne sais s'il est possible de modifier un amendement motive.
Je dois consulter les precedents.
(Le debat est ajourne en attendant la decision de la presidence.)
L'honorable John Lynch-Staunton (chef de l'opposition) :
On a refuse de publiciser la tenue des audiences, contrairement a la coutume.
La liste des temoins changeait presque d'heure en heure.
Il a donc eu tendance a s'eloigner de la teneur de son mandant.
L'article 29 soustrait l'article 93 a l'application de la charte.
Il ne saurait etre plus clair:
La loi provinciale etablissant des ecoles linguistiques est deja mise en oeuvre.
Ces ecoles peuvent etre etablies sans modifier l'article 93.
Le Parti quebecois a vote contre cette disposition lorsqu'il etait dans l'opposition.
Les franco-protestants et les anglo-catholiques etaient particulierement preoccupes.
Elle est intitulee: Des ecoles confessionnelles et des services.
Je vous lis les deux paragraphes qui portent sur le sujet qui nous occupe.
Ils vont plutot parler d'ecole francaise ou anglaise.
Il y a donc au Quebec deux categories d'anglophones.
Le Quebec a souvent ete critique pour cette discrimination, et a juste titre.
C'est un sentiment que ne partagent certainement pas les membres du gouvernement du Quebec.
Il ne faut pas les traiter a la legere, surtout lorsqu'elles touchent des minorites.
Cette opinion a ete redigee et transmise apres les audiences du comite mixte.
Il importe que nous prenions tous connaissance de ce texte.
...qui est la Loi sur l'education adoptee par le Quebec en juin ...
...maintient le principe de la liberte de choix en matiere d'education.
Je pense que, compte tenu des auteurs de cet avis, nous devons y reflechir serieusement.
Il est aussi bien place que quiconque pour repondre aux preoccupations soulevees ici et ailleurs.
L'honorable Jerahmiel S. Grafstein :
Honorables senateurs, l'education, c'est l'affaire de tout le monde.
Je voudrais d'abord dire un mot au sujet des comites et de leurs travaux.
Je prefere les comites senatoriaux.
Les comites senatoriaux sont moins sectaires, plus minutieux et plus consciencieux.
Cela a rendu notre tache doublement difficile.
Nous devons examiner avec soin les temoignages presentes aux comites.
Elle differe sensiblement de celle du Quebec.
Vous les trouverez toutes deux en annexe du rapport qui vous est soumis pour etude.
Le gouvernement du Quebec le savait.
Nous, les parlementaires, le savions.
C'est desormais clair et sans equivoque.
C'est une fausse appellation.
L'article 93 est en substance une serie de privileges plutot que de reels droits.
Ce n'est manifestement pas le cas ici.
Meme sur ce point, nous avons eu des problemes de definition.
Malgre cela, nous avons quand meme traite cette majorite comme une minorite.
Je n'aborderai pas cet aspect du debat.
Nous avons entendu un vaste eventail de temoignages, concluants et clairs.
J'entends le senateur Wood.
J'aborderai ses sujets de preoccupation dans un instant.
Je crois qu'aucun point de vue n'a ete exclu de nos deliberations.
Le comite en a fait un examen approfondi.
Permettez-moi de rappeler brievement quels temoins nous avons entendus.
Il y a eu certaines opinions divergentes.
Il est indiscutable que certains membres de ces deux groupes s'estiment leses.
Son Honneur le President :
Honorables senateurs, je regrette de devoir interrompre le senateur, mais ses 15 minutes sont ecoulees.
Permission est-elle accordee au senateur Grafstein de poursuivre?
Ils ne sont pas sans recours, et j'aborderai cette question tout a l'heure.
Toutefois, c'est une question de droits linguistiques.
Nous y avons renonce en 1982.
Le senateur Pitfield le sait fort bien.
Il etait l'un des architectes de la Loi constitutionnelle de 1982.
Qui l'a mis la?
Nous l'avons dit dans le rapport.
Cela laisse deux groupes, petits, mais importants, qui affirment que la modification leur portera prejudice.
Ce sont les deux groupes dont le senateur Bolduc a parle.
Ces groupes avaient leurs porte-parole.
Les catholiques etaient representes par des pretres anglophones.
Dans son discours reflechi, le senateur Bolduc a fait allusion a ces groupes.
Son Honneur le President :
Je regrette, mais je dois interrompre l'honorable senateur a nouveau.
L'horloge indique 18 heures.
Nous nous sommes entendus pour ne pas voir l'heure, Votre Honneur.
Cela etant dit, leurs droits ne nous sont pas indifferents.
Il y a toujours du financement public pour l'enseignement confessionnel au Quebec.
Puisque nous nous preoccupons de l'avenir, voila ce que je prevois.
Nous ne leur avons pas refuse cette possibilite.
Je vais donner lecture de ces lettres, car elles sont importantes.
La situation n'est pas la meme au Quebec.
Neanmoins, notre assemblee ne s'est pas opposee au choix de modifier l'article 93.
Ils ne se sont pas opposes.
Cela ne veut pas dire que la demarche a ete parfaite.
Tout comme le senateur Pitfield et d'autres, j'ai des preoccupations a cet egard.
Mais cela fait partie du passe.
Le senateur Pitfield etait la.
Ils ont parle de pluralisme pour la toute premiere fois.
Honorables senateurs, j'appuie entierement le rapport et la resolution federale.
L'honorable Therese Lavoie-Roux :
Honorables senateurs, je renvoie ma collegue a la transcription des deliberations.
C'est tout ce qu'il y a de plus clair.
Nous differons d'opinion au sujet de l'alinea 23(1) a ) de la Charte.
Honorables senateurs, j'ai une question pour l'honorable senateur.
Je ne partage pas son point de vue au sujet du consensus.
Ses representants ont comparu devant le comite.
C'est l'un des plus imposants groupes au Quebec.
Encore une fois, honorables senateurs, c'est une question de fait.
Je renvoie tous les honorables senateurs au compte rendu.
Nous avons entendu des representants des commissions scolaires et des enseignants catholiques.
Oui, la resolution n'est pas parfaite.
Nous avons du adapter des outils imparfaits pour qu'ils correspondent a nos exigences minimales.
C'est une precision que mon honorable collegue a oublie de lire.
En fait, elle est deja bien avancee au Quebec.
Il semble que le nouveau systeme sera en place pour la prochaine annee scolaire.
De nouvelles commissions linguistiques seront creees.
Je n'en vois pas.
Honorables senateurs, il y a le droit constitutionnel.
Le Quebec veut se conformer a la primaute du droit.
Tel n'est pas le sens de ma question.
Il veut se conformer a la loi, tout comme Terre-Neuve.
Il veut se conformer a la Constitution.
Ma question n'a rien a voir avec cette reponse.
Il a ete renouvele en 1992 et il est maintenant protege par la Charte.
Les deux resolutions different sensiblement quant a leur contenu et a leur objet.
Il y aura des elections a l'echelon des commissions scolaires en juin.
Encore une fois, honorables senateurs, il y a deux solitudes ici.
Je crois que l'education est un domaine essentiellement provincial.
La n'est pas la question.
C'est inscrit dans la Constitution.
Je n'en disconviens pas.
C'est ce que dit l'article 93.
Non, l'article 93 ne dit rien de la sorte.
Je respecte la Constitution.
Honorables senateurs, il s'agit d'une requete legitime presentee de maniere legitime.
Ce fait a ete reconnu devant le comite.
Ils veulent supprimer l'enseignement religieux, qui leur est impose pour le moment.
Voila le vrai motif.
Souhaitons-nous cautionner cette demarche?
C'est pourtant ce qu'on nous demande de faire.
Je ne sais pas.
Toutefois, je n'insisterai pas la-dessus.
Honorables senateurs, j'essaie de suivre les pratiques et procedures du Senat.
Je parle du rapport du comite mixte special.
Pour moi, la question est d'elucider ce rapport pour les senateurs.
Pour moi, les temoignages et les consequences sont clairs.
C'est ce que j'ai dit dans mon discours.
Honorables senateurs, j'ai ecoute bien attentivement votre plaidoyer.
Vous avez plaide comme un bon avocat.
Il n'y a pas de doute la-dessus.
Vous dites que c'est un probleme provincial.
La question de mon collegue est en deux parties.
Il y a peu de chances que cela arrive.
C'est pourquoi j'appuie le rapport et la resolution.
Son Honneur le President :
Plait-il aux honorables senateurs d'adopter la motion?
(La motion est adoptee avec dissidence, et le rapport est adopte.)
MODIFICATION DE LA CONSTITUTION DU CANADA
Les paragraphes (1) a (4) de l'article 93 ne s'appliquent pas au Quebec.
2. Titre de la presente modification: Modification constitutionnelle de (annee de proclamation) (Quebec).
En bref, j'appuie cette resolution.
Selon le mot du temoin, ce serait la un veritable empilage des structures organisationnelles.
meme chose serait vraie de la ville de Quebec.
Manifestement, l'enjeu residait - et reside toujours - dans la legitimite de la modification.
Honorables senateurs, je vous incite fortement a appuyer cette resolution.
Quel lien y a-t-il entre les deux?
Je vais vous donner ma version des choses.
C'est pour que l'on commence la planification longtemps d'avance.
(Sur la motion du senateur Kinsella, au nom du senateur Simard, le debat est ajourne.
Regie interne, budgets et administration
Adoption du sixieme rapport du comite
L'honorable Sharon Carstairs (leader adjoint du gouvernement) :
Honorables senateurs, au nom du senateur Rompkey, je propose que le rapport soit adopte.
(La motion est adoptee et le rapport est adopte.)
Projet de loi modificatif - Premiere lecture
(Le projet de loi est lu une premiere fois.)
Son Honneur le President :
Honorables senateurs, quand lirons-nous ce projet de loi une deuxieme fois?
(Le projet de loi est lu une premiere fois.)
Son Honneur le President :
Honorables senateurs, quand lirons-nous ce projet de loi une deuxieme fois?
L'honorable Charlie Watt , conformement a l'avis donne le 8 decembre 1997, propose:
2. les mecanismes necessaires a l'implantation des nouvelles relations structurelles;
Que le comite presente son rapport au plus tard le 30 novembre 1999.
(Le Senat s'ajourne au mercredi 10 decembre 1997, a 13 h 30.)


text:
SUBROUTINE Adrianne_thread_state

   USE global_variables
   IMPLICIT NONE

   INTEGER :: i,j,istate
   INTEGER :: check_adrianne

   INTEGER :: free_cells,temp_number,icell
   INTEGER :: cell_numbers_i,cell_numbers_j
   INTEGER,DIMENSION(max_grid_size1*max_grid_size2) :: free_cells_row,free_cells_col

   nAdrianneStates = nAdrianneStates+1
   sudoku_Adrianne(nAdrianneStates,:,:) = sudoku
   sudoku_Adrianne_logic(nAdrianneStates,:,:,:) = sudoku_logic
   sudoku_Adrianne_logic_numbers(nAdrianneStates,:,:) = sudoku_logic_values

  IF (nAdrianneStates .gt. nAdrianneStates_max) THEN
    write(*,*) 'The puzzle cannot be solved by logic alone'
    write(*,*) 'It requires a lot of guessing and trial and error'
    write(*,*) 'Try increasing nAdrianneStates_max or use a brute force solver'
   STOP
 END IF

!  GO TO 110

   free_cells = 0
   DO i=1,max_grid_size1
     DO j=1,max_grid_size2
	IF (sudoku(i,j) .ne. 0) CYCLE
	  free_cells = free_cells+1
	  free_cells_row(free_cells) = i
	  free_cells_col(free_cells) = j
     END DO
   END DO

    DO i=1,free_cells
      cell_numbers_i = sudoku_logic_values(free_cells_row(i),free_cells_col(i))
      DO j=i+1,free_cells
	cell_numbers_j = sudoku_logic_values(free_cells_row(j),free_cells_col(j))
	IF (cell_numbers_j .lt. cell_numbers_i) THEN
	   temp_number = cell_numbers_j
	   cell_numbers_j = cell_numbers_i
	   cell_numbers_i = temp_number

	   temp_number = free_cells_row(j)
	   free_cells_row(j) = free_cells_row(i)
	   free_cells_row(i) = temp_number

	   temp_number = free_cells_col(j)
	   free_cells_col(j) = free_cells_col(i)
	   free_cells_col(i) = temp_number
         END IF
        END DO
      END DO

! 110 CONTINUE

   GO TO 220

   ! Find a square and a number to start Adrianne's Thread
iloop:   DO i=1,max_grid_size1
     DO j=1,max_grid_size2
        IF (sudoku(i,j) .ne. 0) CYCLE
          ! Make sure it is not one of the previous cells of the
          ! Adrianne's Thread array
          check_adrianne = 0
          DO istate=1,nAdrianneStates-1
             IF ((i .eq. Adrianne_state_row(istate)) .AND. &
                  (j .eq. Adrianne_state_col(istate))) THEN
                        check_adrianne = 1
                        EXIT
             END IF
          END DO
          IF (check_adrianne .eq. 0) THEN
                Adrianne_state_row(nAdrianneStates) = i
                Adrianne_state_col(nAdrianneStates) = j
                EXIT iloop
          END IF
     END DO
  END DO iloop

! GO TO 111

220 CONTINUE

  DO icell=1,free_cells
    i = free_cells_row(icell)
    j = free_cells_col(icell)
      check_adrianne = 0
          DO istate=1,nAdrianneStates-1
             IF ((i .eq. Adrianne_state_row(istate)) .AND. &
                  (j .eq. Adrianne_state_col(istate))) THEN
                        check_adrianne = 1
                        EXIT
             END IF
          END DO
          IF (check_adrianne .eq. 0) THEN
                Adrianne_state_row(nAdrianneStates) = i
                Adrianne_state_col(nAdrianneStates) = j
                EXIT
          END IF
  END DO

! 111 CONTINUE

END SUBROUTINE Adrianne_thread_state

!-----------------------------------------------------

SUBROUTINE check_contradiction_Adrianne(sudoku1,sudoku_logic_values1,sudoku_logic1,check_adrianne,max_grid_size11,max_grid_size21)

  USE global_variables
  IMPLICIT NONE

  INTEGER, DIMENSION(max_grid_size1,max_grid_size2) :: sudoku1
  INTEGER, DIMENSION(max_grid_size1,max_grid_size2,9) :: sudoku_logic1
  INTEGER,DIMENSION(max_grid_size1,max_grid_size2) :: sudoku_logic_values1

  INTEGER :: check_adrianne
  INTEGER :: max_grid_size11,max_grid_size21

  INTEGER,DIMENSION(9) :: row_check
  INTEGER,DIMENSION(9) :: col_check
  INTEGER :: i,j,inumber,k
  INTEGER :: small_cell_row,small_cell_column
  INTEGER :: row_k,column_k,ii_other,jj_other
  INTEGER :: ii,jj,cage_local,cage_no_of_cells_local,icell
  INTEGER :: cage_local_previous

  DO i=1,max_grid_size1
    cage_local_previous = 0
    DO j=1,max_grid_size2
      row_check = 0
      IF (sudoku1(i,j) .le. 0) CYCLE
      cage_local = cell_to_cage(i,j,1)
      IF (cage_local .eq. cage_local_previous) CYCLE
      cage_no_of_cells_local = cage_no_of_cells(cage_local)
      DO icell=1,cage_no_of_cells_local
	ii = cage_to_cells(cage_local,icell,1)
	jj = cage_to_cells(cage_local,icell,2)
	IF (sudoku1(ii,jj) .le. 0) CYCLE
        row_check(sudoku1(ii,jj)) = row_check(sudoku1(ii,jj)) + 1
      END DO
      DO inumber=1,9
       IF(row_check(inumber) .gt. 1) THEN
         check_adrianne = 1
         RETURN
       END IF
      END DO
      cage_local_previous = cage_local
    END DO
  END DO

  DO j=1,max_grid_size2
    cage_local_previous = 0
    DO i=1,max_grid_size1
      col_check = 0
      IF (sudoku1(i,j) .le. 0) CYCLE
      cage_local = cell_to_cage(i,j,2)
      IF (cage_local .eq. cage_local_previous) CYCLE
      cage_no_of_cells_local = cage_no_of_cells(cage_local)
      DO icell=1,cage_no_of_cells_local
        ii = cage_to_cells(cage_local,icell,1)
        jj = cage_to_cells(cage_local,icell,2)
        IF (sudoku1(ii,jj) .le. 0) CYCLE
        col_check(sudoku1(ii,jj)) = col_check(sudoku1(ii,jj)) + 1
      END DO
      DO inumber=1,9
       IF(col_check(inumber) .gt. 1) THEN
         check_adrianne = 1
         RETURN
       END IF
      END DO
      cage_local_previous = cage_local
    END DO
  END DO

 DO i=1,max_grid_size1
   DO j=1,max_grid_size2
     IF (sudoku_logic_values1(i,j) .gt. 0) THEN
     END IF
   END DO
 END DO

END SUBROUTINE check_contradiction_Adrianne

!-----------------------------------------------------------------------------------

max_stars_repo_name: surendra-jain/kakuro-solver

max_issues_repo_name: surendra-jain/kakuro-solver

max_forks_repo_name: surendra-jain/kakuro-solver


text:
   subroutine stat3d (title, npin, kd, nassm, power, pave, pmax, pmax2d)
   use mod_coregeom, only : icore, jcore, mapcore, axial
   implicit none
!=======================================================================
!
!  Subroutines for writing 1D and 2D maps to output
!
!  Copyright (c) 2014-2017 Core Physics, Inc.
!
!  Distributed under the MIT license.
!  See the LICENSE file in the main directory for details.
!
!=======================================================================
!
!  Subroutine to print 3D statistics
!
!  General version that considers qtr-symmetry
!
!  2017/05/18 - added 2PIN edit
!
!=======================================================================
      integer, intent(in) :: npin, kd, nassm
      real(8), intent(in) :: power(npin,npin,kd,nassm)
      real(8), intent(out):: pave   ! average
      real(8), intent(out):: pmax   ! max 3pin
      real(8), intent(out):: pmax2d ! max 2pin
      character(len=*), intent(in) :: title

!--- local

      integer  :: i, j, k
      integer  :: ia, ja, na, kpin, nexit
      real(8)  :: p2d, z2d
      real(8)  :: pp
      real(8)  :: zlen, zave      ! axial values
      real(8)  :: c3min, c3max    ! 3D values
      real(8)  :: exave           ! exit average
      integer  :: k2max(3)        ! 2D max locations
      integer  :: k3min(4)        ! 3D min locations
      integer  :: k3max(4)        ! 3D max locations

!--- calculate 3D statistics

      pave=0.0d0   ! output average
      pmax=0.0d0   ! output maximum 3pin
      pmax2d=0.0d0 ! output maximum 2pin

      kpin=0
      zave=0.0d0
      zlen=0.0d0
      c3min=1.0d20   ! core min/max
      c3max=0.0d0
      k3min(:)=0
      k3max(:)=0
      k2max(:)=0

      exave=0.0d0
      nexit=0

      do ja=1, jcore     ! loop over assemblies in full-core
        do ia=1, icore   ! loop over assemblies in full-core
        na=mapcore(ia,ja)
        if (na.eq.0) cycle

        do j=1, npin
          do i=1, npin
            p2d=0.0d0   ! reset for each rod
            z2d=0.0d0
            do k=1, kd
              pp=power(i,j,k,na)
              if (pp.gt.0.0d0) then
                p2d=p2d+axial(k)*pp
                z2d=z2d+axial(k)
                zave=zave+axial(k)*pp
                zlen=zlen+axial(k)
                kpin=kpin+1
                if (pp.lt.c3min) then
                  c3min=pp
                  k3min(1)=i
                  k3min(2)=j
                  k3min(3)=k
                  k3min(4)=na
                endif
                if (pp.gt.c3max) then
                  c3max=pp
                  k3max(1)=i
                  k3max(2)=j
                  k3max(3)=k
                  k3max(4)=na
                endif

              endif
            enddo
            if (z2d.gt.0.0d0) then
              p2d=p2d/z2d
              if (p2d.gt.pmax2d) then
                pmax2d=p2d
                k2max(1)=i
                k2max(2)=j
                k2max(3)=na
              endif
            endif
            pp=power(i,j,kd,na)   ! exit value
            exave=exave+pp
            nexit=nexit+1
          enddo
        enddo
      enddo    ! ia
      enddo    ! ja
      if (zlen.gt.0.0d0) then
        zave=zave/zlen
      else
        c3min=0.0d0    ! avoid overflow
      endif

      write (*,'(/,1x,2a)') '3D Statistics - ', trim(title)

      if (c3max.gt.0.0d0 .and. c3max.lt.0.001d0) then ! print exponent
        write (*,190) 'Min  ', c3min
        write (*,190) 'Max  ', c3max
        write (*,192) 'Max2d', pmax2d
        write (*,192) 'Ave  ', zave
      else
        write (*,180) 'Min  ', c3min, k3min(:)
        write (*,180) 'Max  ', c3max, k3max(:)
        write (*,181) 'Max2D', pmax2d, k2max(:)
        write (*,182) 'Ave  ', zave
      endif
      write (*,'(a,i10)  ')  '  Num  ', kpin

      write (*,*) 'exit average ', exave/dble(nexit), nexit

  180 format (2x, a,f10.4,'  at (i,j,k,na)', 4i4)
  181 format (2x, a,f10.4,'  at (i,j,  na)', 2i4,4x,i4)
  182 format (2x, a,f10.4)
  190 format (2x, a,1p,e14.5,'  at (i,j,k,na)', 0p, 4i4)
  192 format (2x, a,1p,e14.5)

      pave=zave    ! return max value
      pmax=c3max   ! return max value

!--- check for very small values

!     do k=1, kd
!       do na=1, nassm
!         do j=1, npin
!           do i=1, npin
!             pp=power(i,j,k,na)
!             if (pp.gt.0.0d0) then
!               if (pp.lt.zave*1.0d-4) then
!                  write (*,*) 'small value detected ', i, j, k, na, pp
!               endif
!             endif
!           enddo
!         enddo
!       enddo
!     enddo

      return
      end subroutine stat3d
!=======================================================================
!
!  Subroutine to print 3D statistics
!
!  Simple version - does straight average without considering qtr-symmetry
!
!=======================================================================
      subroutine stat3d_simple (title, npin, kd, nassm, axial, power, pmax)
      implicit none
      integer, intent(in) :: npin, kd, nassm
      real(8), intent(in) :: axial(kd)
      real(8), intent(in) :: power(npin,npin,kd,nassm)
      real(8), intent(out):: pmax
      character(len=*), intent(in) :: title

!--- local

      integer  :: i, j, k
      integer  :: ia, kpin
      real(8)  :: pp
      real(8)  :: zlen, zave      ! axial values
      real(8)  :: c3min, c3max    ! 3D values
      integer  :: k3min(4)        ! 3D min locations
      integer  :: k3max(4)        ! 3D max locations

!--- calculate 3D statistics

      pmax=0.0d0   ! output maximum

      kpin=0
      zave=0.0d0
      zlen=0.0d0
      c3min=1.0d20
      c3max=0.0d0
      k3min(:)=0
      k3max(:)=0

      do k=1, kd
        do ia=1, nassm
          do j=1, npin
            do i=1, npin
              pp=power(i,j,k,ia)
              if (pp.gt.0.0d0) then
                zave=zave+axial(k)*pp
                zlen=zlen+axial(k)
                kpin=kpin+1
                if (pp.lt.c3min) then
                  c3min=pp
                  k3min(1)=i
                  k3min(2)=j
                  k3min(3)=k
                  k3min(4)=ia
                endif
                if (pp.gt.c3max) then
                  c3max=pp
                  k3max(1)=i
                  k3max(2)=j
                  k3max(3)=k
                  k3max(4)=ia
                endif

              endif
            enddo
          enddo
        enddo
      enddo
      if (zlen.gt.0.0d0) then
        zave=zave/zlen
      else
        c3min=0.0d0    ! avoid overflow
      endif

      write (*,'(/,1x,2a)') '3D Statistics - ', trim(title)

      if (c3max.gt.0.0d0 .and. c3max.lt.0.001d0) then ! print exponent
        write (*,190) 'Min ', c3min
        write (*,190) 'Max ', c3max
        write (*,192) 'Ave ', zave
      else
        write (*,180) 'Min ', c3min, k3min(:)
        write (*,180) 'Max ', c3max, k3max(:)
        write (*,182) 'Ave ', zave
      endif
      write (*,'(a,i10)  ')  '  Num ', kpin

  180 format (2x, a,f10.4,'  at (i,j,k,na)', 4i4)
  182 format (2x, a,f10.4)
  190 format (2x, a,1p,e14.5,'  at (i,j,k,na)', 0p, 4i4)
  192 format (2x, a,1p,e14.5)

      pmax=c3max   ! return max value

!--- check for very small values

!     do k=1, kd
!       do ia=1, nassm
!         do j=1, npin
!           do i=1, npin
!             pp=power(i,j,k,ia)
!             if (pp.gt.0.0d0) then
!               if (pp.lt.zave*1.0d-4) then
!                  write (*,*) 'small value detected ', i, j, k, ia, pp
!               endif
!             endif
!           enddo
!         enddo
!       enddo
!     enddo

      return
      end subroutine stat3d_simple

!=======================================================================
!
!  Subroutine to print 3D exit values by assembly using pretty output
!  (useful for looking at exit temperatures/densities)
!
!=======================================================================
      subroutine print_exit_map(title, npin, kd, nassm, power)
      implicit none
      integer, intent(in) :: npin, kd, nassm
      real(8), intent(in) :: power(npin,npin,kd,nassm)
      character(len=*), intent(in) :: title

!--- local

      integer            :: ia, klev

!--- print exit values (level kd)

      klev=kd

      do ia=1, nassm        ! loop over assemblies
        write (*,*)
        write (*,*) 'Exit values level ', kd
        call print_single_pin_map(title, npin, kd, nassm, ia, klev, power)
      enddo      ! ia

      return
      end subroutine print_exit_map
!=======================================================================
!
!  Subroutine to print 3D pin powers by assembly using pretty output
!
!=======================================================================
      subroutine print_3D_pin_map(title, npin, kd, nassm, power)
      implicit none
      integer, intent(in) :: npin, kd, nassm
      real(8), intent(in) :: power(npin,npin,kd,nassm)
      character(len=*), intent(in) :: title

!--- local

      integer            :: ia, klev

!--- start

      do ia=1, nassm        ! loop over assemblies
        do klev=kd, 1, -1   ! loop over axial levels
          call print_single_pin_map(title, npin, kd, nassm, ia, klev, power)
        enddo    ! klev
      enddo      ! ia

      return
      end subroutine print_3D_pin_map
!=======================================================================
!
!  Subroutine to print a single pin power map of an assembly using pretty output
!
!  Pass in entire 3D array (nassm,kd), but only print one assembly (ia,klev)
!
!=======================================================================
      subroutine print_single_pin_map(title, npin, kd, nassm, ia, klev, power)
      implicit none
      integer, intent(in) :: npin, kd, nassm, klev, ia
      real(8), intent(in) :: power(npin,npin,kd,nassm)
      character(len=*), intent(in) :: title

!--- local

      integer            :: i, j, m
      integer            :: kpin
      real(8)            :: pp, pave, pmin, pmax
      real(8)            :: zmax
      character(len=200) :: line
      character(len=8)   :: fmt

      if (npin.gt.30) then
        write (*,*) '**** too many pins across to print nice maps ****'
        return
      endif

!--- calculate maximum and set the format correctly

      pp=0.0d0
      zmax=0.0d0
      do j=1, npin
        do i=1, npin
          pp=power(i,j,klev,ia)
          zmax=max(zmax,pp)
        enddo
      enddo

      fmt='(f7.4)'
      if (zmax.gt.  10.0d0) fmt='(f7.3)'
      if (zmax.gt. 100.0d0) fmt='(f7.2)'
      if (zmax.gt.1000.0d0) fmt='(f7.1)'

!d    write (*,*) 'debug: zmax = ', zmax
!d    write (*,*) 'debug: fmt  = ', fmt

!--- write map

      write (*,'(/,1x,a)') trim(title)
!!    write (*,'(1x,a,i3,2a)') 'Assembly ', ia,'    type ', trim(assmmap(nassm))   ! no types to mpact
      if (kd.eq.1) then
        write (*,30)  ia
      else
        write (*,35)  ia, klev
      endif
      pave=0.0
      kpin=0
      pmin=2.0d20
      pmax=0.0d0
      do j=1, npin
        line=' '
        m=0
        do i=1, npin
          pp=power(i,j,klev,ia)
          if (pp.eq.0.0) then
            line(m+1:m+7)='   --- '
          else
            write (line(m+1:m+7),fmt) pp
            pave=pave+pp
            kpin=kpin+1
            pmin=min(pmin,pp)
            pmax=max(pmax,pp)
          endif
          m=m+7
        enddo
        write (*,'(1x,a)') line(1:m)
      enddo
      if (kpin.gt.0) pave=pave/real(kpin)
      if (pmin.gt.1.0d20) pmin=0.0d0   ! protect if no data
      write (*,210) kpin, pave, pmin, pmax
!d    write (*,*) 'debug average = ', pave

! 210 format (5x,'number of hot pins',i4,'    average=',f7.4,'   min=', f7.4,'   max=', f7.4)
  210 format (5x,'number of hot pins',i4,'    average=',f12.4,'   min=', f12.4,'   max=', f12.4)

   30 format (1x,'Assembly ',i3,2a)
   35 format (1x,'Assembly ',i3,'  Level',i3)

      return
      end subroutine print_single_pin_map
!=======================================================================
!
!  Subroutine to print axial edits
!
!=======================================================================
      subroutine print1d(title, npin, kd, nassm, power)
      use mod_coregeom, only : icore, jcore, mapcore, axial
      implicit none
      character(len=*), intent(in) :: title
      integer, intent(in) :: npin, kd, nassm
      real(8), intent(in) :: power(npin,npin,kd,nassm)

      integer          :: i, j, k, nn
      integer          :: ia, ja
      integer          :: na
      real(8)          :: z1, z2
      real(8)          :: pp
      real(8)          :: pave
      real(8)          :: zave
      real(8)          :: ztot, zheat
      real(8)          :: axpow(kd)    ! automatic
      real(8)          :: axmid(kd)    ! automatic
      real(8)          :: axmin(kd)    ! automatic
      real(8)          :: axmax(kd)    ! automatic
      integer          :: numax(kd)    ! automatic

      if (kd.le.1) return   ! skip 1D edits

!--- calculate axial averages

      axpow(:)=0.0d0
      axmin(:)=1.0d20
      axmax(:)=0.0d0
      numax(:)=0

      do ja=1, jcore     ! loop over assemblies in full-core
        do ia=1, icore   ! loop over assemblies in full-core
          na=mapcore(ia,ja)
          if (na.eq.0) cycle

          do k=1, kd       ! loop over axial levels
            do j=1, npin
              do i=1, npin
                pp=power(i,j,k,na)
                if (pp.gt.0.0d0) then
                  axpow(k)=axpow(k)+pp
                  numax(k)=numax(k)+1
                  axmin(k)=min(pp,axmin(k))  ! don't count zero power locations
                  axmax(k)=max(pp,axmax(k))
                endif
              enddo
            enddo
          enddo    ! k
        enddo
      enddo

!  Normalize each axial level
!  Note: this will give "odd" results if different axial levels have different number of pins

      pave=0.0d0   ! core average
      zave=0.0d0   ! height of fuel
      nn=0
      do k=1, kd
        if (numax(k).eq.0) then
           axmin(k)=0    ! protect overflow
        else
           axpow(k)=axpow(k)/dble(numax(k))
           pave=pave+axial(k)*axpow(k)
           zave=zave+axial(k)
           if (nn.eq.0) nn=numax(k)       ! save first non-zero level
           if (nn.ne.numax(k)) nn=-1000   ! set error
        endif
      enddo
      if (zave.gt.0.0d0) pave=pave/zave
      if (nn.lt.0) then
        write (*,*) '*** 1D axial distributions not calculated because some ***'
        write (*,*) '***    axial levels have different numbers of pins     ***'
        return
      endif

!--- calculate heated length

      ztot=0.0d0
      zheat=0.0d0
      do k=1, kd
        ztot=ztot+axial(k)
        if (axpow(k).gt.0.0d0) zheat=zheat+axial(k)
      enddo

!--- find axial midpoints for plotting

      z1=0.0d0
      do k=1, kd
        z2=z1+axial(k)
        axmid(k)=0.5d0*(z1+z2)
        z1=z2
      enddo

      write (*,'(/,1x,a,/,1x,a)') trim(title), '1D Axial Distribution'
      write (*,'(a)')      '    K   Elev         Delta      Average     Minimum     Maximum'
      if (pave.lt.0.001d0) then  ! print exponent
        do k=kd, 1, -1
          write (*,210) k, axmid(k), axial(k), axpow(k), axmin(k), axmax(k)
        enddo
      else
        do k=kd, 1, -1
          write (*,200) k, axmid(k), axial(k), axpow(k), axmin(k), axmax(k)
        enddo
      endif
      write (*,205) pave

      write (*,140) 'total ', ztot
      write (*,140) 'heated', zheat
  140 format (2x,a,' axial length ', f12.4,' cm')

  200 format (i5,f10.4,6f12.4)
  210 format (i5,f10.4,f12.4,1p,6e14.5)
  205 format (2x,'ave',22x,f12.4)

      return
      end subroutine print1d

!=======================================================================
!
!  Collapse 3D edits to 2D
!
! 05/19/2017 - took out 3D max/min edits, these are calculated elsewhere
!            - added option to collapse with loadings
!
!=======================================================================
      subroutine collapse2d(label, npin, kd, nassm, axial, power, pow2d, ifload)
      use mod_coregeom, only : pinload
      implicit none
      character(len=*), intent(in) :: label
      integer, intent(in) :: npin, kd, nassm
      real(8), intent(in) :: power(npin,npin,kd,nassm)
      real(8), intent(in) :: axial(kd)
      real(8), intent(out):: pow2d(npin,npin,nassm)
      logical, intent(in) :: ifload   ! flag to weight with loadings

!--- local

      integer :: i, j, k
      integer :: na
      integer :: nhot
      integer :: k2min(3)
      integer :: k2max(3)
      real(8) :: pp, pl
      real(8) :: zave, zrod
      real(8) :: c2min, c2max

      pow2d(:,:,:)=0.0d0

      write (*,'(/,1x,a,a)') 'Collapsing 3D distributions to 2D - ', trim(label)

!--- collapse to 2D  (loading)

      if (ifload) then

        write (*,*) 'performing axial collapse with pin loadings'
        do na=1, nassm     ! loop over assemblies
          do j=1, npin
            do i=1, npin
              zrod=0.0d0
              zave=0.0d0   ! 2D axial height with power
              do k=1, kd   ! loop over axial levels
                pp=power(i,j,k,na)
                if (pp.gt.0.0d0) then
                  pl=pinload(i,j,k,na)
                  zrod=zrod+axial(k)*pl*pp
                  zave=zave+axial(k)*pl
                  if (pl.eq.0.0d0) write (*,*) 'WARNING: zero loading used in collapse'
                endif
              enddo
              if (zave.gt.0.0d0) pow2d(i,j,na)=zrod/zave
            enddo
          enddo
        enddo      ! na

      else

!--- collapse to 2D  (no loading)

        do na=1, nassm     ! loop over assemblies
          do j=1, npin
            do i=1, npin
              zrod=0.0d0
              zave=0.0d0   ! 2D axial height with power
              do k=1, kd   ! loop over axial levels
                pp=power(i,j,k,na)
                if (pp.gt.0.0d0) then
                  zrod=zrod+axial(k)*pp
                  zave=zave+axial(k)
                endif
              enddo
              if (zave.gt.0.0d0) pow2d(i,j,na)=zrod/zave
            enddo
          enddo
        enddo      ! na

      endif   ! pin loadings

!--- find max/min

      c2max=0.0d0        ! 2D core max
      c2min=1.0d20       ! 2D core min
      k2max(:)=0         ! 2D core max i, j, n
      k2min(:)=0         ! 2D core min i, j, n
      nhot=0

      do na=1, nassm     ! loop over assemblies
        do j=1, npin
          do i=1, npin
            pp=pow2d(i,j,na)
            if (pp.gt.0.0d0) then
              nhot=nhot+1
              if (pp.lt.c2min) then
                c2min=pp
                k2min(1)=i
                k2min(2)=j
                k2min(3)=na
              endif
              if (pp.gt.c2max) then
                c2max=pp
                k2max(1)=i
                k2max(2)=j
                k2max(3)=na
              endif
            endif
          enddo
        enddo
      enddo      ! na

      if (nhot.eq.0) then    ! protect edits
        c2min=0.0d0
      endif

      write (*,190) 'max', c2max, k2max(:)
      write (*,190) 'min', c2min, k2min(:)
  190 format ('   2D ',a,' in core =', f10.4,' at (i,j,na)  ', 2i4,4x,i4)

      return
      end subroutine collapse2d

!=======================================================================
!
!  Subroutine to print 2D Assembly Maps (assembly average values)
!
!=======================================================================
      subroutine print2d_assm_map(title, npin, nassm, pow2)
      use mod_coregeom, only : icore, jcore, mapcore, xlabel, ylabel
      implicit none
      integer, intent(in) :: npin, nassm
      real(8), intent(in) :: pow2(npin,npin,nassm)
      character(len=*), intent(in) :: title

      character(len=8) :: fmt

      integer          :: i, j, nn
      integer          :: nnsave
      integer          :: ia, ja, na
      real(8)          :: pmin, pmax
      real(8)          :: pp
      real(8)          :: passm(nassm)  ! automatic

      logical  :: ifbw  ! flag if assemblies have different number of pins

      ifbw=.false.

      write (*,'(/,1x,2a,/)') trim(title), ' - 2D assembly average (2RPF)'

!--- calculate 2D assembly averages

      passm(:)=0.0d0

      nnsave=0
      do na=1, nassm
        nn=0
        do j=1, npin
          do i=1, npin
            pp=pow2(i,j,na)
            if (pp.gt.0.0d0) then
              nn=nn+1
              passm(na)=passm(na)+pp
            endif
          enddo
        enddo
        if (nn.gt.0) passm(na)=passm(na)/dble(nn)
        if (nnsave.eq.0) nnsave=nn    ! save first time
        if (nn.ne.nnsave) then
          ifbw=.true.
!         write (0,*) 'debug: assembly ', na, nn, nnsave
        endif
      enddo

!--- calculate average - use mapcore in case this is qtr-core

      nn=0
      pp=0.0d0
      pmin=passm(1)
      pmax=passm(1)

      do ja=1, jcore
        do ia=1, icore
          na=mapcore(ia,ja)
          if (na.gt.0) then
            pp=pp+passm(na)
            nn=nn+1
            pmax=max(pmax,passm(na))
            pmin=min(pmin,passm(na))
          endif
        enddo    ! ia
      enddo      ! ja

      if (nn.gt.0) pp=pp/dble(nn)
      write (*,125) nn
      if (ifbw) then
        write (*,*) 'WARNING: average of assemblies cannot be calculated because assemblies have different number of rods'
      else
        write (*,130) 'average', pp
!!      if (abs(pp-1.0d0).gt.0.0001) write (0,*) '***** check normalization *****'   ! msg only valid for pin powers
      endif
      write (*,130) 'maximum', pmax
      write (*,130) 'minimum', pmin
  125 format (4x,'number of assemblies in full-core ', i0)
  130 format (4x,a,' assembly power  ', f8.4)

      write (*,*)

      fmt='(f8.4)'
      if (pp.gt.100.0d0) fmt='(f8.2)'

!--- print map

      write (*,'("  **  ",50(4x,a2,2x))') (xlabel(i),i=1,icore)   ! labels are 2

      do j=1, jcore
        write (*,'(2x,a2,"- ")',advance='no') ylabel(j)
        do i=1, icore
           if (mapcore(i,j).eq.0) then
             write (*,'(8x)',advance='no')
           else
             write (*,fmt,advance='no') passm(mapcore(i,j))
           endif
        enddo
        write (*,*)
      enddo

      return
      end subroutine print2d_assm_map
!=======================================================================
!
!  Subroutine to print 2D max pin in assembly maps (2PIN)
!
!  2017/05/18 - added for Jim
!
!=======================================================================
      subroutine print2d_2pin_map(title, npin, nassm, power, kd)
      use mod_coregeom, only : icore, jcore, mapcore, axial, xlabel, ylabel
      implicit none
      integer, intent(in) :: npin, nassm
      integer, intent(in) :: kd
      real(8), intent(in) :: power(npin,npin,kd,nassm)
      character(len=*), intent(in) :: title

      character(len=8) :: fmt

      integer          :: i, j, k
      integer          :: ia, ja, na
      real(8)          :: pp
      real(8)          :: zave, zlen
      real(8)          :: pin2d(icore,jcore)  ! automatic
      real(8)          :: core2d      ! core max

      write (*,'(/,1x,2a,/)') trim(title), ' - max 2D rod in each assembly (2PIN)'

!--- calculate max pin in each assembly - use mapcore in case this is qtr-core

      pin2d(:,:)=0.0d0
      core2d=0.0d0

      do ja=1, jcore
        do ia=1, icore
          na=mapcore(ia,ja)
          if (na.gt.0) then

            do j=1, npin
              do i=1, npin
                zave=0.0d0
                zlen=0.0d0
                do k=1, kd   ! loop over axial levels
                  pp=power(i,j,k,na)
                  if (pp.gt.0.0d0) then
                    zave=zave+axial(k)*pp
                    zlen=zlen+axial(k)
                  endif
                enddo
                if (zlen.gt.0.0d0) then
                  zave=zave/zlen
                  pin2d(ia,ja)=max(pin2d(ia,ja),zave)
                endif
              enddo
            enddo

            core2d=max(core2d,pin2d(ia,ja))   ! core max

          endif  ! na
        enddo    ! ia
      enddo      ! ja

      fmt='(f8.4)'
      if (core2d.gt.100.0d0) fmt='(f8.2)'

!--- print map

      write (*,'("  **  ",50(4x,a2,2x))') (xlabel(i),i=1,icore)   ! labels are 2

      do ja=1, jcore
        write (*,'(2x,a2,"- ")',advance='no') ylabel(ja)
        do ia=1, icore
           if (mapcore(ia,ja).eq.0) then
             write (*,'(8x)',advance='no')
           else
             write (*,fmt,advance='no') pin2d(ia,ja)
           endif
        enddo
        write (*,*)
      enddo

      write (*,*) 'WARNING: 2PIN maps do not use loadings if collapsing exposures'

      return
      end subroutine print2d_2pin_map
!=======================================================================
!
!  Subroutine to calculate core-wide axial offset
!
!  General version that considers qtr-symmetry
!
!=======================================================================
      subroutine calc_axoff (title, npin, kd, nassm, power, axoff)
      use mod_coregeom, only : icore, jcore, mapcore, axial
      implicit none
      integer, intent(in) :: npin, kd, nassm
      real(8), intent(in) :: power(npin,npin,kd,nassm)
      real(8), intent(out):: axoff  ! Axial offset
      character(len=*), intent(in) :: title

!--- local

      real(8)  :: zmesh(0:kd)    ! automatic

      integer  :: i, j, k
      integer  :: ia, ja, na
      integer  :: kmid
      real(8)  :: chi
      real(8)  :: pp
      real(8)  :: ptop, pbot
      real(8)  :: ztop, zbot
      real(8)  :: zmid
      real(8)  :: zold

      axoff=0.0d0
      if (kd.eq.1) return   ! exit for 2D problems

!--- calculate midpoint

      zold=0.0d0
      do k=1, kd
        zmesh(k)=zold+axial(k)
        zold=zmesh(k)
      enddo
      zmid=zmesh(kd)*0.5d0

!d    write (*,*) 'debug: core height ', zmesh(kd)
!d    write (*,*) 'debug: core mid    ', zmid

      chi=1.0d0    ! set if boundaries align
      kmid=0
      zold=0.0d0
      do k=1, kd
        if (zmid.ge.zold .and. zmid.lt.zmesh(k)) then    ! fix 2019/02/21
          kmid=k
          chi=(zmid-zmesh(k-1))/axial(k)
        endif
        zold=zmesh(k)  ! protect bounds when k=1
      enddo

!d    write (*,*) 'debug: kmid ', kmid
!d    write (*,*) 'debug: chi  ', chi

!--- calculate 3D statistics

      ptop=0.0d0     ! axial offset edits
      pbot=0.0d0
      ztop=0.0d0
      zbot=0.0d0

      do ja=1, jcore     ! loop over assemblies in full-core
        do ia=1, icore   ! loop over assemblies in full-core
        na=mapcore(ia,ja)
        if (na.eq.0) cycle

        do k=1, kd
          do j=1, npin
            do i=1, npin
              pp=power(i,j,k,na)
              if (pp.gt.0.0d0) then
                if (k.lt.kmid) then
                  pbot=pbot+axial(k)*pp
                  zbot=zbot+axial(k)
                elseif (k.gt.kmid) then
                  ptop=ptop+axial(k)*pp
                  ztop=ztop+axial(k)
                else    ! core midplane in node
                  pbot=pbot+axial(k)*chi*pp
                  zbot=zbot+axial(k)*chi
                  ptop=ptop+axial(k)*(1.0d0-chi)*pp
                  ztop=ztop+axial(k)*(1.0d0-chi)
                endif
              endif
            enddo  ! i
          enddo    ! j
        enddo      ! k

      enddo    ! ia
      enddo    ! ja

!d    write (*,*) 'debug: zbot ', zbot
!d    write (*,*) 'debug: ztop ', ztop

      ptop=ptop/ztop
      pbot=pbot/zbot
      axoff=(ptop-pbot)/(ptop+pbot)*100.0d0

      write (*,210) trim(title), axoff
  210 format (2x,a,' A/O  ',f10.4,' %')

      return
      end subroutine calc_axoff

max_stars_repo_name: palmtag/VERAout-tools

max_issues_repo_name: palmtag/VERAout-tools

max_forks_repo_name: palmtag/VERAout-tools