text:
#ifdef PICO

    module Calculator_PICO
    use CosmologyTypes
    use CosmoTheory
    use Calculator_CAMB
    use CAMB, only : CAMBParams_Set, highL_unlensed_cl_template, ThermoDerivedParams, nthermo_derived, &
        BackgroundOutputs, lmin, CAMBParams
    use Errors !CAMB
    use settings
    use likelihood
    use Calculator_Cosmology
    use GeneralTypes
    use fpico
    implicit none
    private

    Type, extends(CAMB_Calculator) :: PICO_Calculator
    contains
    procedure :: CMBToCAMB => PICO_CMBToCAMB
    procedure :: GetNewTransferData => PICO_GetNewTransferData
    procedure :: GetNewPowerData => PICO_GetNewPowerData
    procedure :: GetCosmoTheoryForImportance => PICO_GetTheoryForImportance
    procedure :: VersionTraceOutput => PICO_VersionTraceOutput
    procedure :: ReadParams => PICO_ReadParams
    end type PICO_Calculator

    public PICO_Calculator
    contains

    subroutine PICO_CMBToCAMB(this,CMB,P)
    class(PICO_Calculator) :: this
    class(CMBParams) CMB
    type(CAMBParams)  P

    call this%CAMB_Calculator%CMBToCamb(CMB,P)
    P%Reion%use_optical_depth = .true.
    P%OnlyTransfers = .true.
    if (.not. CosmoSettings%CMB_Lensing) call MpiStop('PICO assumes lensed CMB')

    end subroutine PICO_CMBToCAMB

    subroutine PICO_GetNewTransferData(this, CMB,Info,Theory,error)
    class(PICO_Calculator) :: this
    class(CMBParams) CMB
    class(TTheoryIntermediateCache), pointer :: Info
    class(TCosmoTheoryPredictions) :: Theory
    integer error
    !Do nothing
    end subroutine PICO_GetNewTransferData

    subroutine PICO_GetOutputArray(Theory,i,j, name)
    class(TCosmoTheoryPredictions) :: Theory
    integer i,j
    character(LEN=*) :: name
    integer(fpint) lmax, lmn

    if (CosmoSettings%cl_lmax(i,j)>0) then
        lmax = min(CosmoSettings%lmax_computed_cl,CosmoSettings%cl_lmax(i,j))
        lmn = lmin
        call fpico_read_output(name,Theory%Cls(i,j)%CL(lmin:),lmn,lmax)
    end if

    end subroutine PICO_GetOutputArray

    subroutine PICO_GetNewPowerData(this, CMB, Info, Theory, error)
    use ModelParams
    use ModelData
    class(PICO_Calculator) :: this
    class(CMBParams) :: CMB
    class(TTheoryIntermediateCache), pointer :: Info
    class(TCosmoTheoryPredictions) :: Theory
    integer error,i,j, lmaxCL,lmx
    type(CAMBParams) P
    integer(fpint) :: success
    real(mcp) :: highL_norm = 0
    integer, save, allocatable :: indicesT(:,:)

    if (.not. allocated(indicesT)) then
        allocate(indicesT(3,3))
        indicesT=0
        indicesT(1,1) =  C_Temp
        indicesT(2,2) =  C_E
        indicesT(2,1) =  CT_Cross
        indicesT(3,3) =  CT_B
    end if

    call this%CMBToCAMB(CMB, P)
    call CAMBParams_Set(P)
    call this%SetBackgroundTheoryData(CMB,Theory,error)

    if (p%InitPower%n_runrun(1)/=0 .or. p%InitPower%nt_run(1)/=0 .or. p%InitPower%ant(1)/=0) &
        & call MpiStop('PICO: currently unsupported initial power parameter')

    if (CosmoSettings%Use_LSS .or. CosmoSettings%get_sigma8) then
        call MpiStop('PICO: currently no MPK or sigma8')
    end if

    call fpico_reset_params()
    call fpico_set_param("ombh2", CMB%ombh2)
    call fpico_set_param("omch2", CMB%omch2)
    call fpico_set_param("omnuh2", CMB%omnuh2)
    call fpico_set_param("omvh2", CMB%omv*CMB%h**2)
    call fpico_set_param("omk", CMB%omk)
    call fpico_set_param("hubble", CMB%H0)
    call fpico_set_param("Alens", CMB%Alens)
    call fpico_set_param("w", CMB%w)
    call fpico_set_param("theta", CosmomcTheta())
    call fpico_set_param("helium_fraction", p%yhe)

    call fpico_set_param("massive_neutrinos", CMB%nnu)
    call fpico_set_param("scalar_spectral_index(1)",p%InitPower%an(1))
    call fpico_set_param("scalar_nrun(1)",p%InitPower%n_run(1))
    call fpico_set_param("initial_ratio(1)",p%InitPower%rat(1))
    call fpico_set_param("scalar_amp(1)",p%InitPower%ScalarPowerAmp(1))
    call fpico_set_param("pivot_scalar",p%InitPower%k_0_scalar)
    call fpico_set_param("pivot_tensor",p%InitPower%k_0_tensor)
    call fpico_set_param("re_optical_depth",CMB%tau)
    !    call fpico_set_param("force",1.d0)

    call fpico_reset_requested_outputs()
    if (P%WantCls) then
        call fpico_request_output("cl_TT")
        if (CosmoSettings%num_cls>1) then
            call fpico_request_output("cl_TE")
            call fpico_request_output("cl_EE")
            if (CosmoSettings%num_cls>2) then
                if (CosmoSettings%cl_lmax(3,3)>0) call fpico_request_output("cl_BB")
                if (CosmoSettings%num_cls>3) then
                    !                    call MpiStop('PICO: lensing output currently innaccurate')
                    !                    if (CosmoSettings%cl_lmax(4,4)>0) call fpico_request_output("cl_pp")
                end if
            end if
        end if
    end if
   ! if (P%WantTransfer) then
   !     call fpico_request_output("k")
   !     call fpico_request_output("pk")
   ! end if

    call fpico_compute_result(success)
    if (success == 0) then
        ! for now just reject if out of pico bounds
        !call mpiStop('PICO failed to get result')
        if (Feedback>1) write(*,*) 'PICO out of bounds'
        error = 1
        return
    end if
    error = 0
    if (P%WantCls) then
        call PICO_GetOutputArray(Theory,1,1, "cl_TT")
        if (CosmoSettings%num_cls>1) then
            call PICO_GetOutputArray(Theory,2,1, "cl_TE")
            call PICO_GetOutputArray(Theory,2,2, "cl_EE")
            if (CosmoSettings%num_cls>2) call PICO_GetOutputArray(Theory,3,3, "cl_BB")
            !            if (CosmoSettings%num_cls>3) call PICO_GetOutputArray(Theory,4,4, "cl_pp") need to change units, as in scalCls
        end if
    end if

    do i=1, min(3,CosmoSettings%num_cls)
        do j= i, 1, -1
            lmaxCL = CosmoSettings%cl_lmax(i,j)
            lmx = min(CosmoSettings%lmax_computed_cl, lmaxCL)
            if (lmx/=0) then
                associate( CL => Theory%Cls(i,j)%CL)
                    if (CosmoSettings%lmax_computed_cl < lmaxCL) then
                        if (highL_norm ==0) & !normally normalize off TT
                        & highL_norm = CL(lmx)/this%highL_lensedCL_template(lmx,indicesT(i,j))
                        CL(lmx+1:lmaxCL) =  highL_norm*this%highL_lensedCL_template(lmx+1:lmaxCL,indicesT(i,j))
                    end if
                end associate
            end if
        end do
    end do

    end subroutine PICO_GetNewPowerData

    subroutine PICO_GetTheoryForImportance(this, CMB, Theory, error)
    class(PICO_Calculator) :: this
    class(CMBParams) CMB
    class(TCosmoTheoryPredictions) Theory
    integer :: error

    call this%GetNewPowerData(CMB, null(), Theory, error)
    call this%GetNewBackgroundData(CMB,Theory,error)
    if (error==0) call this%SetDerived(Theory)

    end subroutine PICO_GetTheoryForImportance

    subroutine PICO_VersionTraceOutput(this, ReadValues)
    use GaugeInterface, only : Eqns_name
    class(PICO_Calculator) :: this
    class(TNameValueList) :: ReadValues

    call ReadValues%Add('Compiled_PICO_version', '3.3.0') !dynamic version variable?
    call this%CAMB_Calculator%VersionTraceOutput(ReadValues)

    end subroutine PICO_VersionTraceOutput

    subroutine PICO_ReadParams(this,Ini)
    class(PICO_Calculator) :: this
    class(TSettingIni) :: Ini

    call this%CAMB_Calculator%ReadParams(Ini)
    this%calcName ='PICO'

    !$ write(*,*) '**WARNING**: pico may not work when CosmoMC compiled with -openmp (why??)'

    call fpico_init(1_fpint)
    call fpico_load(Ini%Read_String_Default("pico_datafile", EnvDefault=.true.))
    call fpico_set_verbose(int(IfThenElse(Ini%Read_Logical("pico_verbose",.false.),1,0),fpint))

    end subroutine PICO_ReadParams

    end module Calculator_PICO
#endif

max_stars_repo_name: sjoudaki/CosmoJBD

max_issues_repo_name: sjoudaki/CosmoJBD

max_forks_repo_name: sjoudaki/CosmoJBD


text:
c===========================================================================
c
c   This file is part of TISEAN
c
c   Copyright (c) 1998-2007 Rainer Hegger, Holger Kantz, Thomas Schreiber
c
c   TISEAN is free software; you can redistribute it and/or modify
c   it under the terms of the GNU General Public License as published by
c   the Free Software Foundation; either version 2 of the License, or
c   (at your option) any later version.
c
c   TISEAN is distributed in the hope that it will be useful,
c   but WITHOUT ANY WARRANTY; without even the implied warranty of
c   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
c   GNU General Public License for more details.
c
c   You should have received a copy of the GNU General Public License
c   along with TISEAN; if not, write to the Free Software
c   Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
c
c===========================================================================
c   help.f
c   Utilities for usage message
c   author T. Schreiber (1998)
c===========================================================================
      subroutine whatido(text,iverb)
      character*72 progname
      character*(*) text

      call getarg(0,progname)
      call argdel(0)
      iverb=igetv(iverb)
      if(iv_io(iverb).eq.1) then
         write(istderr(),'()')
         write(istderr(),'(a)')
     .              "TISEAN 3.0.1 (C) R. Hegger, H. Kantz, T. Schreiber
     .(1998-2007)"
         write(istderr(),'()')
         write(istderr(),'(a,a,a)')
     .      progname(1:index(progname," ")-1), ": ", text
      endif
      if(lopt("h",1).eq.1) call usage()
      end

      subroutine whatineed(text)
      character*72 progname
      character*(*) text

      call getarg(0,progname)
      write(istderr(),'()')
      write(istderr(),'(a,a,x,a)')
     .   "Usage: ", progname(1:index(progname," ")-1),  text
      end

      subroutine popt(c,text)
      character*(*) c,text

      write(istderr(),'(5h    -,a,x,1h<,a,1h>)') c, text
      end

      subroutine ptext(text)
      character*(*) text

      write(istderr(),'(3x,a)') text
      end

      subroutine pout(text)
      character*(*) text

      write(istderr(),'(8h    -o <,a,a,1h>)')
     .   "output file name, just -o means ", text
      end

      subroutine pall()

      call popt("V","verbosity level (0 = only fatal errors)")
      call popt("h","show this message")
      write(istderr(),'()')
      end

ext: f

max_stars_repo_name: homiak/jabaws

max_issues_repo_name: homiak/jabaws

max_forks_repo_name: homiak/jabaws


text:
!   Copyright 2014 College of William and Mary
!
!   Licensed under the Apache License, Version 2.0 (the "License");
!   you may not use this file except in compliance with the License.
!   You may obtain a copy of the License at
!
!     http://www.apache.org/licenses/LICENSE-2.0
!
!   Unless required by applicable law or agreed to in writing, software
!   distributed under the License is distributed on an "AS IS" BASIS,
!   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
!   See the License for the specific language governing permissions and
!   limitations under the License.

!
!****************************************************************************************
!											*
!	Read binary format v5.0 (hybrid S-Z models) of *hvel.67 (sidecenter vel.) and
!       compute fluxes at the borders between multiple regions (excluding open/land boundaries).
!       WARNING: does not work for lat/lon!
!       Works with mixed grids and LSC2
!
!       Inputs: vgrid.in, *_elev.61, *hvel.67, fluxflag.prop (depth=-1,0,1,...);
!               screen inputs. If itsflux/=0, *temp.70, *salt.70 (not completed)
!       Outputs: fluxes.out (if the difference between region #s=1, and neith=-1)
!       History:
!****************************************************************************************
!     ifort -cpp -O2 -mcmodel=medium -CB -Bstatic -o compute_fluxes_ns.WW ../UtilLib/schism_geometry.f90 ../UtilLib/compute_zcor.f90 compute_fluxes_ns.f90
      program read_out
      use compute_zcor
      use schism_geometry_mod
      parameter(nbyte=4)
      character*30 file63
      character*12 it_char
      character*48 start_time,version,variable_nm,variable_dim
      character*48 data_format
      integer,allocatable :: elnode(:,:),elside(:,:),i34(:)
      allocatable :: sigma(:),cs(:),ztot(:),x(:),y(:),dp(:),kbp(:),kfp(:)
      allocatable :: icum(:,:,:),eta2(:),ztmp(:,:)
      allocatable :: xcj(:),ycj(:),dps(:),kbs(:),xctr(:),yctr(:),dpe(:),kbe(:)
      allocatable :: zs(:,:),ze(:,:),idry(:),idry_s(:),idry_e(:)
      allocatable :: ic3(:,:),isdel(:,:),isidenode(:,:),nwild(:),iflux_e(:)
      allocatable :: fluxes(:),suv2(:,:,:),tsel(:,:,:),kbp2(:),kbs2(:),kbe2(:),sigma_lcl(:,:)

!      print*, 'Outputs in Cartesian (1) or lat/lon (2)?'
!      read(*,*) ics
      ics=1
      print*, 'Input start and end stack #:'
      read(*,*)iday1,iday2

      print*, 'Do you want to compute T,S fluxes? 0: no; 1: yes'
      read(*,*)itsflux

!...  Header
      open(63,file='1_elev.61',status='old',access='direct',recl=nbyte)
      open(64,file='1_hvel.67',status='old',access='direct',recl=nbyte)
      if(itsflux==1) then
        open(65,file='1_temp.70',status='old',access='direct',recl=nbyte)
        open(66,file='1_salt.70',status='old',access='direct',recl=nbyte)
      endif
      irec=0
      do m=1,48/nbyte
        read(63,rec=irec+m) data_format(nbyte*(m-1)+1:nbyte*m)
      enddo
      if(data_format.ne.'DataFormat v5.0') then
        print*, 'This code reads only v5.0:  ',data_format
        stop
      endif
      irec=irec+48/nbyte
      do m=1,48/nbyte
        read(63,rec=irec+m) version(nbyte*(m-1)+1:nbyte*m)
      enddo
      irec=irec+48/nbyte
      do m=1,48/nbyte
        read(63,rec=irec+m) start_time(nbyte*(m-1)+1:nbyte*m)
      enddo
      irec=irec+48/nbyte
      do m=1,48/nbyte
        read(63,rec=irec+m) variable_nm(nbyte*(m-1)+1:nbyte*m)
      enddo
      irec=irec+48/nbyte
      do m=1,48/nbyte
        read(63,rec=irec+m) variable_dim(nbyte*(m-1)+1:nbyte*m)
      enddo
      irec=irec+48/nbyte

      write(*,'(a48)')data_format
      write(*,'(a48)')version
      write(*,'(a48)')start_time
      write(*,'(a48)')variable_nm
      write(*,'(a48)')variable_dim

      read(63,rec=irec+1) nrec
      read(63,rec=irec+2) dtout
      read(63,rec=irec+3) nspool
      read(63,rec=irec+4) ivs
      read(63,rec=irec+5) i23d
      irec=irec+5

      print*, 'ivs=',ivs,'; i23d=',i23d,' ;nrec= ',nrec

!     Vertical grid (obsolete; will be overwritten later)
      read(63,rec=irec+1) nvrt
      read(63,rec=irec+2) kz
      read(63,rec=irec+3) h0
      read(63,rec=irec+4) h_s
      read(63,rec=irec+5) h_c
      read(63,rec=irec+6) theta_b
      read(63,rec=irec+7) theta_f
      irec=irec+7

!      do k=1,kz-1
!        read(63,rec=irec+k) ztot(k)
!      enddo
!      do k=kz,nvrt
!        kin=k-kz+1
!        read(63,rec=irec+k) sigma(kin)
!        cs(kin)=(1-theta_b)*sinh(theta_f*sigma(kin))/sinh(theta_f)+ &
!     &theta_b*(tanh(theta_f*(sigma(kin)+0.5))-tanh(theta_f*0.5))/2/tanh(theta_f*0.5)
!      enddo
      irec=irec+nvrt

      irec_tmp=irec !save for non-standard outputs

!     Horizontal grid
      read(63,rec=irec+1) np
      read(63,rec=irec+2) ne
      irec=irec+2

      allocate(x(np),y(np),dp(np),kbp(np),kfp(np),elnode(4,ne),idry(np),idry_e(ne), &
     &eta2(np),xctr(ne),yctr(ne),dpe(ne),kbe(ne),ztmp(nvrt,np),ze(nvrt,ne), &
     &tsel(2,nvrt,ne),kbp2(np),kbe2(ne),i34(ne),stat=istat)
      if(istat/=0) stop 'Falied to allocate (2)'

      do m=1,np
        read(63,rec=irec+1)x(m)
        read(63,rec=irec+2)y(m)
        read(63,rec=irec+3)dp(m)
        read(63,rec=irec+4)kbp(m)
        irec=irec+4
      enddo !m=1,np

      do m=1,ne
        read(63,rec=irec+1)i34(m)
        irec=irec+1
        do mm=1,i34(m)
          read(63,rec=irec+1)elnode(mm,m)
          irec=irec+1
        enddo !mm
      enddo !m
      irec0=irec

!     hvel.67
      read(64,rec=irec_tmp+1) ns
      read(64,rec=irec_tmp+2) ns_e
      print*, 'side grid:',ns,ns_e
      allocate(dps(ns),kbs(ns),kbs2(ns),zs(nvrt,ns),idry_s(ns),suv2(2,nvrt,ns),stat=istat)
      if(istat/=0) stop 'Allocation error: side(8)'
      irec0_u=irec_tmp+2
      do m=1,ns
        read(64,rec=irec0_u+3)dps(m)
        read(64,rec=irec0_u+4)kbs(m)
        irec0_u=irec0_u+4
      enddo !m

!      do m=1,ns_e
!        read(64,rec=irec0_u+1)i34
!        irec0_u=irec0_u+1
!        do mm=1,i34
!          read(64,rec=irec0_u+1)itmp
!          irec0_u=irec0_u+1
!          if(m==ns_e) print*, 'last elem in side grid:',itmp
!        enddo !mm
!      enddo !m

      irec0_u=irec0_u+4*ns_e !pure tri for sidecenters.gr3

!     *.70
      if(itsflux==1) then
        read(65,rec=irec_tmp+1) ne1
        read(65,rec=irec_tmp+2) ne_e
        if(ne1/=ne) stop 'mismatch in elem.'
        print*, 'Elem. grid:',ne1,ne_e
        irec0_TS=irec_tmp+2
        do m=1,ne
          read(65,rec=irec0_TS+3)dpe(m)
          read(65,rec=irec0_TS+4)kbe(m)
          irec0_TS=irec0_TS+4
        enddo !m
        irec0_TS=irec0_TS+4*ne_e
      endif !itsflux

!     Compute geometry
      call compute_nside(np,ne,i34,elnode,ns2)
      if(ns/=ns2) stop 'mismatch in sides'
      print*, 'done compute_nside'
      allocate(ic3(4,ne),elside(4,ne),isdel(2,ns),isidenode(2,ns),xcj(ns),ycj(ns),stat=istat)
      if(istat/=0) stop 'Allocation error: side(7)'
      call schism_geometry_single(np,ne,ns2,x,y,i34,elnode,ic3,elside,isdel,isidenode,xcj,ycj)
      print*, 'done schism_geometry'

      print*, 'last element',(elnode(j,ne),j=1,i34(ne))

!     Read vgrid.in
      allocate(ztot(nvrt),sigma(nvrt),sigma_lcl(nvrt,np),kbp(np),stat=istat)
      if(istat/=0) stop 'Falied to allocate (1)'
      call get_vgrid_single('vgrid.in',np,nvrt,ivcor,kz,h_s,h_c,theta_b,theta_f,ztot,sigma,sigma_lcl,kbp)

!     Read in fluxflag.prop
      allocate(iflux_e(ne),stat=istat)
      open(32,file='fluxflag.prop',status='old')
      do i=1,ne
        read(32,*)j,tmp1
        itmp=tmp1
        if(itmp<-1) stop 'fluxflag.prop wrong'
        iflux_e(i)=itmp
      enddo
      close(32)

!     Only compute if diff(region #)=1; fluxes(i) is from i to i-1
      max_reg=maxval(iflux_e)
      allocate(fluxes(max_reg),stat=istat)
      if(istat/=0) stop 'Falied to allocate (6)'
      open(60,file='fluxes.out')

!...  Compute relative record # for a node and level for 3D outputs
!...
!      icount=0
!      do i=1,np
!        do k=max0(1,kbp(i)),nvrt
!          do m=1,ivs
!            icount=icount+1
!            icum(i,k,m)=icount
!          enddo !m
!        enddo !k
!      enddo !i=1,np

!...  Time iteration
!...
      it_tot=0
      do iday=iday1,iday2
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      write(it_char,'(i12)')iday
      open(63,file=it_char//'_elev.61',status='old',access='direct',recl=nbyte)
      open(64,file=it_char//'_hvel.67',status='old',access='direct',recl=nbyte)
      if(itsflux==1) then
        open(65,file=it_char//'_temp.70',status='old',access='direct',recl=nbyte)
        open(66,file=it_char//'_salt.70',status='old',access='direct',recl=nbyte)
      endif

      irec=irec0
      irec_u=irec0_u
      irec_TS=irec0_TS
      do it1=1,nrec
        read(63,rec=irec+1) time
        read(63,rec=irec+2) it
        irec=irec+2
        it_tot=it_tot+1
!        time=it_tot*dtout

        print*, 'time in days =',time/86400

        do i=1,np
          read(63,rec=irec+i) eta2(i)
        enddo !i
        irec=irec+np+np !done with (63 for this step

        !Compute z coordinates
        do i=1,np
          if(ivcor==1) then !localized
            if(dp(i)+eta2(i)<=h0) then
              idry(i)=1
            else !wet
              idry(i)=0
              do k=kbp(i),nvrt
                ztmp(k,i)=(eta2(i)+dp(i))*sigma_lcl(k,i)+eta2(i)
              enddo !k
            endif !wet/dry
          else if(ivcor==2) then !SZ
            call zcor_SZ_single(dp(i),eta2(i),h0,h_s,h_c,theta_b,theta_f,kz,nvrt,ztot, &
     &sigma,ztmp(:,i),idry(i),kbpl)
            if(kbpl/=kbp(i)) stop 'mismatch (2)'
          endif
        enddo !i=1,np

!       idry_e, idry_s
        do i=1,ne
          idry_e(i)=maxval(idry(elnode(1:i34(i),i)))
          kbe2(i)=minval(kbp(elnode(1:i34(i),i)))
        enddo !i
        idry_s=1 !init.
        do i=1,ns
          kbs2(i)=minval(kbp(isidenode(1:2,i)))
          do j=1,2
            ie=isdel(j,i)
            if(ie/=0.and.idry_e(max(1,ie))==0) idry_s(i)=0
          enddo !j
          !Check
          if(idry_s(i)==0) then
            itmp=maxval(idry(isidenode(1:2,i)))
            if(itmp/=0) then
              write(*,*)'Wet side has dry nodes:',i,isidenode(1:2,i)
              stop
            endif
          endif
        enddo !i=1,ns

!       hvel.67
        irec_u=irec_u+ns+2
        suv2=0 !for below bottom
        do i=1,ns
          do k=max0(1,kbs(i)),nvrt
            do m=1,2
              read(64,rec=irec_u+1) suv2(m,k,i)
              irec_u=irec_u+1
            enddo !m
          enddo !k

          !Debug
          !if(iday==ndays.and.it1==nrec) write(97,*)xcj(i),ycj(i),suv2(1:2,nvrt,i)
        enddo !i

        !Compute zcor for wet sides (deal with some inconsistency)
        do i=1,ns
          if(idry_s(i)==1) cycle

          n1=isidenode(1,i); n2=isidenode(2,i)
          do k=kbs2(i),nvrt
            k1=max(k,kbp(n1))
            k2=max(k,kbp(n2))
            zs(k,i)=(ztmp(k1,n1)+ztmp(k2,n2))/2
            if(abs(zs(k,i))>1.e18) then
              write(*,*)'Dry side (1):',i,k,n1,n2,ztmp(k,n1),ztmp(k,n2)
              stop
            endif
          enddo !k
        enddo !i=1,ns

!       T,S
        if(itsflux==1) then
          irec_TS=irec_TS+ne+2
          tsel=-99
          do i=1,ne
            do k=max0(1,kbe(i)),nvrt
              read(65,rec=irec_TS+1) tsel(1,k,i)
              read(66,rec=irec_TS+1) tsel(2,k,i)
              irec_TS=irec_TS+1
            enddo !k

            !Debug
            !if(iday==ndays.and.it1==nrec) write(99,*)i,tsel(2,nvrt,i)
          enddo !i

          do i=1,ne
            if(idry_e(i)==1) cycle

            !Compute zcor for wet elem.
            do k=kbe2(i),nvrt
              ze(k,i)=0
              do j=1,i34(i)
                nd=elnode(j,i)
                kl=max(k,kbp(nd))
                ze(k,i)=ze(k,i)+ztmp(kl,nd)/i34(i)
              enddo !j
              if(abs(ze(k,i))>1.e18) then
                write(*,*)'Dry elem:',i
                stop
              endif
            enddo !k

            !Extend valid values to bottom for wet elem.
            !After this all wet elem. have valid values from 1:nvrt
            kk0=-1
            do k=1,nvrt
              if(tsel(1,k,i)>=0) then
                kk0=k; exit
              endif
            enddo !k
            if(kk0<=0) then
              write(*,*)'Weird T:',i,kk0,tsel(1,:,i)
              stop
            else if(tsel(2,kk0,i)<0) then
              write(*,*)'Weird S:',i,kk0,tsel(2,:,i)
              stop
            endif
            do k=1,kk0-1
              tsel(1:2,k,i)=tsel(1:2,kk0,i)
            enddo !k
          enddo !i=1,ne
        endif !itsflux

!       Fluxes
        fluxes=0
        do i=1,ns
          if(isdel(2,i)==0.or.idry_s(i)==1) cycle

          !Internal wet sides
          ie1=isdel(1,i); ie2=isdel(2,i)
          if(iflux_e(ie1)==-1.or.iflux_e(ie2)==-1) cycle
          if(iabs(iflux_e(ie1)-iflux_e(ie2))/=1) cycle

          n1=isidenode(1,i); n2=isidenode(2,i)
          distj=sqrt((x(n1)-x(n2))**2+(y(n1)-y(n2))**2)
          thetan=atan2(x(n1)-x(n2),y(n2)-y(n1))
          snx=cos(thetan); sny=sin(thetan)
          ireg_hi=max(iflux_e(ie1),iflux_e(ie2))
          ireg_lo=min(iflux_e(ie1),iflux_e(ie2))
          if(ireg_hi==iflux_e(ie1)) then
            isgn=1 !'hi' to 'lo' is same as elem. '1' to '2'
          else
            isgn=-1
          endif

          do k=kbs2(i),nvrt-1
            !Normal vel. from 1 to 2
            vnn=(suv2(1,k+1,i)+suv2(1,k,i))/2*snx+(suv2(2,k+1,i)+suv2(2,k,i))/2*sny
            fluxes(ireg_hi)=fluxes(ireg_hi)+isgn*vnn*distj*(zs(k+1,i)-zs(k,i)) !m^3/s
          enddo !k
        enddo !i=1,ns

!       Output. fluxes(i) is from region i to i-1
        write(60,'(f16.6,60000(1x,e14.4))')time/86400,fluxes(:)
      enddo !it1=1,nrec

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      enddo !iday

      stop
      end

max_stars_repo_name: rustychris/schism

max_issues_repo_name: rustychris/schism

max_forks_repo_name: rustychris/schism


text:
! Contain the functions related to brute force method
! and the test functionas
! use brute force method to solve
    module BruteForce
    use constpara
    use mysolclass
    implicit none

    integer::totalnumfeasible
    integer, allocatable, dimension(:,:)::pool
    real*8, allocatable, dimension(:,:)::fre_pool

    contains

    subroutine bf_enumerate_fleet(basenwk)
    use GraphLib
    use mysolclass
    implicit none
    integer::p,l
    real*8::inifre(nline)
    type(solclass)::sol
    real*8,ALLOCATABLE::tmp_fre(:,:)
    type(graphclass),intent(in)::basenwk

    call get_pool
    if (nline.ne.4) then
        write(*,*) " Warnning: the bruteforce method only applies when the number of line equals 4"
    end if
    allocate(tmp_fre(totalnumfeasible+1,4))
    caseindex = 0
    call sol%inisol(basenwk)

    do l = 1, nline
        call sol%mylines(l)%copylines(basenwk%mylines(l))
    enddo

    open(1, file="C:\GitCodes\BTNDP\Input\TestNetwork\inifre.txt")
    do l = 1, nline
        read(1, *) inifre(l)
    enddo
    close(1)

    ! write cases
    open(1, file = "C:\GitCodes\BTNDP\Input\TestNetwork\numcases.txt")

    write(1,*) totalnumfeasible+1
    ! evaluate base case
    do l =1, nline
        tmp_fre(1,l) =  inifre(l)
    end do

    call sol%evaluate(basenwk)
    call sol%dp%outputod(sol%dp%xfa,sol%dp%fx,nl,ndest)
    call sol%dp%outputx
    caseindex= caseindex + 1

    do p = 1, totalnumfeasible
        write(*,*) "Test Fre = ",pool(p,:)
        call sol%update_fleet_and_fre(pool(p,:))
        do l =1, nline
            tmp_fre(p+1,l) =  sol%mylines(l)%fre
        end do
        write(*,*) "Output Fre = ",sol%mylines(1)%fre,sol%mylines(2)%fre,sol%mylines(3)%fre,sol%mylines(4)%fre
        call sol%evaluate(basenwk)
        call sol%dp%outputod(sol%dp%xfa,sol%dp%fx,nl,ndest)
        call sol%dp%outputx
        caseindex= caseindex + 1
    enddo

    open(2, file = "C:\GitCodes\BTNDP\Input\TestNetwork\setfre.txt")
    open(3, file = "C:\GitCodes\BTNDP\Results\enumeratefleet.txt")

    do p = 1, totalnumfeasible+1
        write(2,"(f8.3,a,f8.3,a,f8.3,a,f8.3)") tmp_fre(p,1),",",tmp_fre(p,2),",",tmp_fre(p,3),",",tmp_fre(p,4)
    enddo
    do p = 1, totalnumfeasible
       write(3,"(i4,a,i4,a,i4,a,i4)") pool(p,1),",",pool(p,2),",",pool(p,3),",",pool(p,4)
    enddo
    close(1)
    close(2)
    close(3)

    deallocate(tmp_fre)
    end subroutine

    subroutine bf_given_fre(basenwk)
    use GraphLib
    use mysolclass
    implicit none
    integer::p,l
    type(solclass)::sol
    type(graphclass),intent(in)::basenwk
    call read_pool
    caseindex = 0
    call sol%inisol(basenwk)
    do l = 1, nline
        call sol%mylines(l)%copylines(basenwk%mylines(l))
    enddo
    !do p = totalfea, 1,-1
    do p = 1, totalnumfeasible
        ! write(*,"(4(f6.2,2X))") fre_pool(p,:)
        write(*,*) " Test Fre = ",fre_pool(p,:)
        do l =1, nline
            sol%mylines(l)%fre = fre_pool(p,l)
        end do
        call sol%evaluate(basenwk)
        call sol%dp%outputod(sol%dp%xfa,sol%dp%fx,nl,ndest)
        call sol%dp%outputx
        caseindex= caseindex + 1
    enddo
    end subroutine

    subroutine read_pool
    implicit none

    real*8::tff(4)
    integer::i
    if (nline.ne.4) then
        write(*,*) "Warnning: read pool, the size of line should equal 4"
    end if
    open(1, file='C:\GitCodes\BTNDP\Input\TestNetwork\numcases.txt')
    read(1,*) totalnumfeasible
    close(1)
    allocate(fre_pool(totalnumfeasible,4))
    open(1, file='C:\GitCodes\BTNDP\Input\TestNetwork\setfre.txt')
    do i = 1, totalnumfeasible
        read(1,*) tff(:)
        fre_pool(i,:) =  tff(:)
    enddo

    end subroutine

    ! get the solution pool
    subroutine get_pool
    implicit none
    integer::l1,l2,l3,l4
    totalnumfeasible = 0
    do l1 = fleet_lb(1), fleet_ub(1)
        do l2 =  fleet_lb(2),fleet_ub(2)
            do l3 = fleet_lb(3), fleet_ub(3)
                do l4 =  fleet_lb(4),fleet_ub(4)
                    if (l1+l2+l3+l4.le.fleetsize) then
                        totalnumfeasible =  totalnumfeasible + 1
                    end if
                end do
            enddo
        end do
    end do

    write(*,*) "total feasible pool = ",totalnumfeasible

    allocate(pool(totalnumfeasible,4))
    totalnumfeasible = 0
    do l1 = fleet_lb(1), fleet_ub(1)
        do l2 =  fleet_lb(2),fleet_ub(2)
            do l3 = fleet_lb(3), fleet_ub(3)
                do l4 =  fleet_lb(4),fleet_ub(4)
                    !if (l1+l2+l3+l4.eq.fleetsize) then
                    if (l1+l2+l3+l4.le.fleetsize) then
                        totalnumfeasible = totalnumfeasible + 1
                        pool(totalnumfeasible,1) = l1
                        pool(totalnumfeasible,2) = l2
                        pool(totalnumfeasible,3) = l3
                        pool(totalnumfeasible,4) = l4
                    end if
                end do
            enddo
        end do
    end do

    end subroutine

    end module

max_stars_repo_name: hkujy/BTNDP

max_issues_repo_name: hkujy/BTNDP

max_forks_repo_name: hkujy/BTNDP