module m_LESModels
  use decomp_2d
  use m_TypeDef
  use m_LogInfo
  use m_MeshAndMetries
  use m_Parameters
  use m_Variables,only: mb1,hi1
  implicit none
  private
  real(RK),dimension(:),allocatable,save::Fcoe

  public:: InitFilterCoe,Clc_SGS_Vis
contains

  !******************************************************************
  ! InitFilterCoe
  !******************************************************************
  subroutine InitFilterCoe()
    implicit none

    if(nrank==0) then
      !call MainLog%CheckForError(ErrT_Pass,"InitFilterCoe: ","NOTE: now the LES model in the code is ONLY used for channel flows(pbc in x-dir & z-dir, npbc in y-dir)!")
    endif
    allocate(Fcoe(9))
    if(FilterType==0) then     ! 0:trapezoidal type
      Fcoe = (/one, two,  one, two, four,    two,  one, two,  one/)/sixteen
    elseif(FilterType==1) then ! 1:Simpson type
      Fcoe = (/one, four, one, four,sixteen, four, one, four, one/)/36.00_RK
    else
      call MainLog%CheckForError(ErrT_Abort,"InitFilterCoe","unrecognized FilterType")
    endif
  end subroutine InitFilterCoe

  !******************************************************************
  ! Clc_SGS_Vis
  !******************************************************************
  subroutine Clc_SGS_Vis(ux,uy,uz,nut)
    implicit none
    real(RK),dimension(mb1%xmm:mb1%xpm,mb1%ymm:mb1%ypm,mb1%zmm:mb1%zpm),intent(in) ::ux,uy,uz
    real(RK),dimension(mb1%xmm:mb1%xpm,mb1%ymm:mb1%ypm,mb1%zmm:mb1%zpm),intent(inout)::nut

    SELECT CASE(LES_type)
    CASE(1)  ! 1: Smagorinsky Model

    CASE(2)  ! 2: WALE Model
      call Clc_SGS_WALE(ux,uy,uz,nut)
    CASE(3)  ! 3: MTS Model
      call Clc_SGS_MTS(ux,uy,uz,nut)
    END SELECT

    call SetBC_and_UpdateHalo_SGS(nut)
  end subroutine Clc_SGS_Vis

  !******************************************************************
  ! Calculate the SGS kinematic viscosity using WALE model
  !  WALE: Wall-adapting local eddy-vicosity model
  ! Reference:
  !  Nicoud, F., & Ducros, F. (1999).
  !   Subgrid-scale stress modelling based on the square of the velocity
  !   gradient tensor Flow, Turbulence and Combustion, 62(3), 183-200.
  !******************************************************************
  subroutine Clc_SGS_WALE(ux,uy,uz,nut)
    implicit none
    real(RK),dimension(mb1%xmm:mb1%xpm,mb1%ymm:mb1%ypm,mb1%zmm:mb1%zpm),intent(in) ::ux,uy,uz
    real(RK),dimension(mb1%xmm:mb1%xpm,mb1%ymm:mb1%ypm,mb1%zmm:mb1%zpm),intent(out)::nut

    ! locals
    real(RK),parameter::Cw_=0.325_RK
    integer::ic,jc,kc,im,jm,km,ip,jp,kp
    real(RK)::caj,cac1,cac2,cac12,SijSij,SdijSdij,SdTrace
    real(RK)::dudx,dudy,dudz,dvdx,dvdy,dvdz,dwdx,dwdy,dwdz,CwDelta2
    real(RK)::S11,S12,S13,S22,S23,S33,Sd11,Sd12,Sd13,Sd22,Sd23,Sd33

    do kc=ystart(3),yend(3)
      km=kc-1; kp=kc+1
      do jc=ystart(2),yend(2)
        jp=jc+1; jm=jc-1
        caj  = rdyp(jc)
        cac1 = rdyc(jc)
        cac2 = rdyc(jp)
        cac12= cac1-cac2
        CwDelta2= Cw_*DeltaCell(jc)
        CwDelta2= CwDelta2*CwDelta2
        do ic=ystart(1),yend(1)
          im=ic-1; ip=ic+1
          dudx= (ux(ip,jc,kc)-ux(ic,jc,kc))*rdx
          dudy=((ux(ip,jp,kc)+ux(ic,jp,kc))*cac2+(ux(ip,jc,kc)+ux(ic,jc,kc))*cac12 &
               -(ux(ip,jm,kc)+ux(ic,jm,kc))*cac1)*0.25_RK
          dudz= (ux(ip,jc,kp)+ux(ic,jc,kp)-ux(ip,jc,km)-ux(ic,jc,km))*rdz*0.25_RK

          dvdx= (uy(ip,jp,kc)-uy(im,jp,kc)+uy(ip,jc,kc)-uy(im,jc,kc))*rdx*0.25_RK
          dvdy= (uy(ic,jp,kc)-uy(ic,jc,kc))*caj
          dvdz= (uy(ic,jp,kp)+uy(ic,jc,kp)-uy(ic,jp,km)-uy(ic,jc,km))*rdz*0.25_RK

          dwdx= (uz(ip,jc,kp)-uz(im,jc,kp)+uz(ip,jc,kc)-uz(im,jc,kc))*rdx*0.25_RK
          dwdy=((uz(ic,jp,kp)+uz(ic,jp,kc))*cac2+(uz(ic,jc,kp) +uz(ic,jc,kc))*cac12 &
               -(uz(ic,jm,kp)+uz(ic,jm,kc))*cac1)*0.25_RK
          dwdz= (uz(ic,jc,kp)-uz(ic,jc,kc))*rdz

          S11= dudx
          S12=(dudy+dvdx)*half
          S13=(dudz+dwdx)*half
          S22= dvdy
          S23=(dvdz+dwdy)*half
          S33= dwdz
          SijSij= S11*S11+ S22*S22+ S33*S33+ two*(S12*S12+ S13*S13+ S23*S23)

          Sd11= dudx*dudx+ dudy*dvdx+ dudz*dwdx
          Sd12=(dudx*dudy+ dudy*dvdy+ dudz*dwdy+ dvdx*dudx+ dvdy*dvdx+ dvdz*dwdx)*half
          Sd13=(dudx*dudz+ dudy*dvdz+ dudz*dwdz+ dwdx*dudx+ dwdy*dvdx+ dwdz*dwdx)*half
          Sd22= dvdx*dudy+ dvdy*dvdy+ dvdz*dwdy
          Sd23=(dvdx*dudz+ dvdy*dvdz+ dvdz*dwdz+ dwdx*dudy+ dwdy*dvdy+ dwdz*dwdy)*half
          Sd33= dwdx*dudz+ dwdy*dvdz+ dwdz*dwdz
          SdTrace=(Sd11+Sd22+Sd33)/Three
          Sd11=Sd11-SdTrace
          Sd22=Sd22-SdTrace
          Sd33=Sd33-SdTrace
          SdijSdij= Sd11*Sd11+ Sd22*Sd22+ Sd33*Sd33+ two*(Sd12*Sd12+ Sd13*Sd13+ Sd23*Sd23)

          nut(ic,jc,kc)=CwDelta2*(SdijSdij**1.50_RK)/(SijSij**2.50_RK+SdijSdij**1.25_RK) +xnu
        enddo
      enddo
    enddo

  end subroutine Clc_SGS_WALE

  !******************************************************************
  ! Calculate the SGS kinematic viscosity using Mixed-Time-Scale model
  !
  ! Reference:
  ! Inagaki M, Kondoh T, Nagano Y, et al. (2005).
  !   A Mixed-Time-Scale SGS Model With Fixed Model-Parameters for Practical LES.
  !     Journal of Fluids Engineering-transactions of The Asme, 127(1): 1-13.
  !******************************************************************
  subroutine Clc_SGS_MTS(ux,uy,uz,nut)
    implicit none
    real(RK),dimension(mb1%xmm:mb1%xpm,mb1%ymm:mb1%ypm,mb1%zmm:mb1%zpm),intent(in) ::ux,uy,uz
    real(RK),dimension(mb1%xmm:mb1%xpm,mb1%ymm:mb1%ypm,mb1%zmm:mb1%zpm),intent(out)::nut

    ! locals
    type(HaloInfo):: hi_K
    integer::ic,jc,kc,ip,jp,kp,im,km
    real(RK)::uxf,uyf,uzf,Ksgs,Ts,rdelta
    real(RK),dimension(ystart(1):yend(1),ystart(2):yend(2),ystart(3):yend(3))::strain_cha,Kes
    real(RK),dimension(mb1%xmm:mb1%xpm,mb1%ymm:mb1%ypm,mb1%zmm:mb1%zpm)::Ktemp

    call ClcVelStrainRate(ux,uy,uz,strain_cha)

    ! filter ux, only in homogeneous directions
    DO kc=ystart(3),yend(3)
      km=kc-1
      kp=kc+1
      do jc=ystart(2),yend(2)
        do ic=ystart(1),yend(1)
          im=ic-1
          ip=ic+1
          uxf= Fcoe(1)*ux(im,jc,km) + Fcoe(2)     *ux(ic,jc,km) +Fcoe(3)*ux(ip,jc,km) + &
               Fcoe(4)*ux(im,jc,kc) +(Fcoe(5)-one)*ux(ic,jc,kc) +Fcoe(6)*ux(ip,jc,kc) + &
               Fcoe(7)*ux(im,jc,kp) + Fcoe(8)     *ux(ic,jc,kp) +Fcoe(9)*ux(ip,jc,kp)
          Ktemp(ic,jc,kc)= uxf*uxf
        enddo
      enddo
    ENDDO
    hi_K%pencil = y_pencil
    hi_K%xmh=0;  hi_K%xph=1
    hi_K%ymh=0;  hi_K%yph=0
    hi_K%zmh=0;  hi_K%zph=0
    call myupdate_halo(Ktemp, mb1, hi_K)
    DO kc=ystart(3),yend(3)
      do jc=ystart(2),yend(2)
        do ic=ystart(1),yend(1)
          ip=ic+1
          Kes(ic,jc,kc)= half*(Ktemp(ic,jc,kc) +Ktemp(ip,jc,kc))
        enddo
      enddo
    ENDDO

    ! filter uy, only in homogeneous directions
    DO kc=ystart(3),yend(3)
      km=kc-1
      kp=kc+1
      do jc=1,nyp
        do ic=ystart(1),yend(1)
          im=ic-1
          ip=ic+1
          uyf= Fcoe(1)*uy(im,jc,km) + Fcoe(2)     *uy(ic,jc,km) +Fcoe(3)*uy(ip,jc,km) + &
               Fcoe(4)*uy(im,jc,kc) +(Fcoe(5)-one)*uy(ic,jc,kc) +Fcoe(6)*uy(ip,jc,kc) + &
               Fcoe(7)*uy(im,jc,kp) + Fcoe(8)     *uy(ic,jc,kp) +Fcoe(9)*uy(ip,jc,kp)
          Ktemp(ic,jc,kc)= uyf*uyf
        enddo
      enddo
    ENDDO
    DO kc=ystart(3),yend(3)
      do jc=1,nyc
        jp=jc+1
        do ic=ystart(1),yend(1)
          Kes(ic,jc,kc)= Kes(ic,jc,kc)+ half*(Ktemp(ic,jc,kc) +Ktemp(ic,jp,kc))
        enddo
      enddo
    ENDDO

    ! filter uz, only in homogeneous directions
    DO kc=ystart(3),yend(3)
      km=kc-1
      kp=kc+1
      do jc=ystart(2),yend(2)
        do ic=ystart(1),yend(1)
          im=ic-1
          ip=ic+1
          uzf= Fcoe(1)*uz(im,jc,km) + Fcoe(2)     *uz(ic,jc,km) +Fcoe(3)*uz(ip,jc,km) + &
               Fcoe(4)*uz(im,jc,kc) +(Fcoe(5)-one)*uz(ic,jc,kc) +Fcoe(6)*uz(ip,jc,kc) + &
               Fcoe(7)*uz(im,jc,kp) + Fcoe(8)     *uz(ic,jc,kp) +Fcoe(9)*uz(ip,jc,kp)
          Ktemp(ic,jc,kc)= uzf*uzf
        enddo
      enddo
    ENDDO
    hi_K%pencil = y_pencil
    hi_K%xmh=0;  hi_K%xph=0
    hi_K%ymh=0;  hi_K%yph=0
    hi_K%zmh=0;  hi_K%zph=1
    call myupdate_halo(Ktemp, mb1, hi_K)
    DO kc=ystart(3),yend(3)
      kp=kc+1
      do jc=ystart(2),yend(2)
        do ic=ystart(1),yend(1)
          Kes(ic,jc,kc)= Kes(ic,jc,kc)+ half*(Ktemp(ic,jc,kc) +Ktemp(ic,jc,kp))
        enddo
      enddo
    ENDDO

    ! Finally calculate the nut
    DO kc=ystart(3),yend(3)
      do jc=ystart(2),yend(2)
        rdelta = one/DeltaCell(jc)
        do ic=ystart(1),yend(1)
          Ksgs= Kes(ic,jc,kc)
          Ts = sqrt(Ksgs)*rdelta + strain_cha(ic,jc,kc)/ten
          if(Ts < 1.00E-10_RK) then
            Ts= zero
          else
            Ts= one/Ts
          endif
          nut(ic,jc,kc)=0.05_RK*Ksgs*Ts + xnu
        enddo
      enddo
    ENDDO

  end subroutine Clc_SGS_MTS

  !******************************************************************
  ! ClcVelStrainTensor
  !******************************************************************
  subroutine ClcVelStrainTensor(ux,uy,uz,st1,st2,st3,st4,st5,st6)
    implicit none
    real(RK),dimension(mb1%xmm:mb1%xpm,mb1%ymm:mb1%ypm,mb1%zmm:mb1%zpm),intent(in)::ux,uy,uz
    real(RK),dimension(ystart(1):yend(1),ystart(2):yend(2),ystart(3):yend(3)),intent(out)::st1,st2,st3
    real(RK),dimension(ystart(1):yend(1),ystart(2):yend(2),ystart(3):yend(3)),intent(out)::st4,st5,st6

    ! locals
    integer::ic,jc,kc,ip,jp,kp,im,jm,km
    real(RK)::caj,cac1,cac2,cac12

    ! st1 = duxdx                st2 = duydy                st3 = duzdz
    ! st4 =(duxdy + duydx)*0.5   st5 =(duxdz + duzdx)*0.5   st6 =(duydz + duzdy)*0.5
    DO kc=ystart(3),yend(3)
      kp=kc+1
      km=kc-1
      do jc=ystart(2),yend(2)
        jp=jc+1
        jm=jc-1
        caj  = rdyp(jc)
        cac1 = rdyc(jc)
        cac2 = rdyc(jp)
        cac12= cac1 - cac2
        do ic=ystart(1),yend(1)
          ip=ic+1
          im=ic-1
          st1(ic,jc,kc)=(ux(ip,jc,kc)-ux(ic,jc,kc))*rdx
          st2(ic,jc,kc)=(uy(ic,jp,kc)-uy(ic,jc,kc))*caj
          st3(ic,jc,kc)=(uz(ic,jc,kp)-uz(ic,jc,kc))*rdz

          !rt1=(ux(ic,jc,kc)-ux(ic,jm,kc))*cac1 +(uy(ic,jc,kc)-uy(im,jc,kc))*rdx ! 2*st4 at (i,    j,    k+1/2)
          !rt2=(ux(ic,jp,kc)-ux(ic,jc,kc))*cac2 +(uy(ic,jp,kc)-uy(im,jp,kc))*rdx ! 2*st4 at (i,    j+1,  k+1/2)
          !rt3=(ux(ip,jc,kc)-ux(ip,jm,kc))*cac1 +(uy(ip,jc,kc)-uy(ic,jc,kc))*rdx ! 2*st4 at (i+1,  j,    k+1/2)
          !rt4=(ux(ip,jp,kc)-ux(ip,jc,kc))*cac2 +(uy(ip,jp,kc)-uy(ic,jp,kc))*rdx ! 2*st4 at (i+1,  j+1,  k+1/2)
          !st4(ic,jc,kc)=oneEighth*(rt1+rt2+rt3+rt4)                             !   st4 at (i+1/2,j+1/2,k+1/2)
          st4(ic,jc,kc)=oneEighth*( (uy(ip,jp,kc) -uy(im,jp,kc) +uy(ip,jc,kc) -uy(im,jc,kc))*rdx  &
                                   +(ux(ip,jp,kc) +ux(ic,jp,kc))*cac2                             &
                                   +(ux(ip,jc,kc) +ux(ic,jc,kc))*cac12                            &
                                   -(ux(ip,jm,kc) +ux(ic,jm,kc))*cac1                             )

          !rt1=(ux(ic,jc,kc)-ux(ic,jc,km))*rdz  +(uz(ic,jc,kc)-uz(im,jc,kc))*rdx ! 2*st5 at (i,    j+1/2,k    )
          !rt2=(ux(ic,jc,kp)-ux(ic,jc,kc))*rdz  +(uz(ic,jc,kp)-uz(im,jc,kp))*rdx ! 2*st5 at (i,    j+1/2,k+1  )
          !rt3=(ux(ip,jc,kc)-ux(ip,jc,km))*rdz  +(uz(ip,jc,kc)-uz(ic,jc,kc))*rdx ! 2*st5 at (i+1,  j+1/2,k    )
          !rt4=(ux(ip,jc,kp)-ux(ip,jc,kc))*rdz  +(uz(ip,jc,kp)-uz(ic,jc,kp))*rdx ! 2*st5 at (i+1,  j+1/2,k+1  )
          !st5(ic,jc,kc)=oneEighth*(rt1+rt2+rt3+rt4)                             !   st5 at (i+1/2,j+1/2,k+1/2)
          st5(ic,jc,kc)=oneEighth*( (ux(ip,jc,kp) +ux(ic,jc,kp) -ux(ip,jc,km) -ux(ic,jc,km))*rdz  &
                                   +(uz(ip,jc,kp) -uz(im,jc,kp) +uz(ip,jc,kc) -uz(im,jc,kc))*rdx  )

          !rt1=(uz(ic,jc,kc)-uz(ic,jm,kc))*cac1 +(uy(ic,jc,kc)-uy(ic,jc,km))*rdz ! 2*st6 at (i+1/2,j,    k    )
          !rt2=(uz(ic,jp,kc)-uz(ic,jc,kc))*cac2 +(uy(ic,jp,kc)-uy(ic,jp,km))*rdz ! 2*st6 at (i+1/2,j+1,  k    )
          !rt3=(uz(ic,jc,kp)-uz(ic,jm,kp))*cac1 +(uy(ic,jc,kp)-uy(ic,jc,kc))*rdz ! 2*st6 at (i+1/2,j,    k+1  )
          !rt4=(uz(ic,jp,kp)-uz(ic,jc,kp))*cac2 +(uy(ic,jp,kp)-uy(ic,jp,kc))*rdz ! 2*st6 at (i+1/2,j+1,  k+1  )
          !st6(ic,jc,kc)=oneEighth*(rt1+rt2+rt3+rt4)                             !   st6 at (i+1/2,j+1/2,k+1/2)
          st6(ic,jc,kc)=oneEighth*( (uy(ic,jp,kp) +uy(ic,jc,kp) -uy(ic,jp,km) -uy(ic,jc,km))*rdz  &
                                   +(uz(ic,jp,kp) +uz(ic,jp,kc))*cac2                             &
                                   +(uz(ic,jc,kp) +uz(ic,jc,kc))*cac12                            &
                                   -(uz(ic,jm,kp) +uz(ic,jm,kc))*cac1                             )
        enddo
      enddo
    ENDDO

  end subroutine ClcVelStrainTensor

  !******************************************************************
  ! ClcVelStrainRate
  !******************************************************************
  subroutine ClcVelStrainRate(ux,uy,uz,strain_cha)
    implicit none
    real(RK),dimension(mb1%xmm:mb1%xpm,mb1%ymm:mb1%ypm,mb1%zmm:mb1%zpm),intent(in)::ux,uy,uz
    real(RK),dimension(ystart(1):yend(1),ystart(2):yend(2),ystart(3):yend(3)),intent(out)::strain_cha

    ! locals
    integer::ic,jc,kc,ip,jp,kp,im,jm,km
    real(RK)::st1,st2,st3,st4,st5,st6,caj,cac1,cac2,cac12

    !  characteristic filtered rate of strain: strain_cha
    !  |S|= sqrt(2*S_ij*S_ij), where S_ij = 0.5*(dui/dxj + duj/dxi)
    DO kc=ystart(3),yend(3)
      kp=kc+1
      km=kc-1
      do jc=ystart(2),yend(2)
        jp=jc+1
        jm=jc-1
        caj  = rdyp(jc)
        cac1 = rdyc(jc)
        cac2 = rdyc(jp)
        cac12= cac1 - cac2
        do ic=ystart(1),yend(1)
          ip=ic+1
          im=ic-1
          st1=(ux(ip,jc,kc)-ux(ic,jc,kc))*rdx
          st2=(uy(ic,jp,kc)-uy(ic,jc,kc))*caj
          st3=(uz(ic,jc,kp)-uz(ic,jc,kc))*rdz
          st4=oneEighth*( (uy(ip,jp,kc) -uy(im,jp,kc) +uy(ip,jc,kc) -uy(im,jc,kc))*rdx  &
                         +(ux(ip,jp,kc) +ux(ic,jp,kc))*cac2                             &
                         +(ux(ip,jc,kc) +ux(ic,jc,kc))*cac12                            &
                         -(ux(ip,jm,kc) +ux(ic,jm,kc))*cac1                             )
          st5=oneEighth*( (ux(ip,jc,kp) +ux(ic,jc,kp) -ux(ip,jc,km) -ux(ic,jc,km))*rdz  &
                         +(uz(ip,jc,kp) -uz(im,jc,kp) +uz(ip,jc,kc) -uz(im,jc,kc))*rdx  )
          st6=oneEighth*( (uy(ic,jp,kp) +uy(ic,jc,kp) -uy(ic,jp,km) -uy(ic,jc,km))*rdz  &
                         +(uz(ic,jp,kp) +uz(ic,jp,kc))*cac2                             &
                         +(uz(ic,jc,kp) +uz(ic,jc,kc))*cac12                            &
                         -(uz(ic,jm,kp) +uz(ic,jm,kc))*cac1                             )
          strain_cha(ic,jc,kc)= sqrt(two*(st1*st1+ st2*st2+ st3*st3)+ four*(st4*st4 + st5*st5 + st6*st6))
        enddo
      enddo
    ENDDO

  end subroutine ClcVelStrainRate

  !******************************************************************
  ! ClcVelStrainRate
  !******************************************************************
  subroutine SetBC_and_UpdateHalo_SGS(nut)
    implicit none
    real(RK),dimension(mb1%xmm:mb1%xpm,mb1%ymm:mb1%ypm,mb1%zmm:mb1%zpm),intent(inout)::nut

    ! yp-dir
    SELECT CASE(BcOption(yp_dir))
    CASE(BC_NSLIP)
      nut(ystart(1):yend(1), nyp, ystart(3):yend(3)) =  two*xnu-nut(ystart(1):yend(1), nyc, ystart(3):yend(3))
    CASE(BC_FSLIP)
      nut(ystart(1):yend(1), nyp, ystart(3):yend(3)) =          nut(ystart(1):yend(1), nyc, ystart(3):yend(3))
    END SELECT

    ! ym-dir
    SELECT CASE(BcOption(ym_dir))
    CASE(BC_NSLIP)
      nut(ystart(1):yend(1), 0, ystart(3):yend(3)) =  two*xnu-nut(ystart(1):yend(1), 1, ystart(3):yend(3))
    CASE(BC_FSLIP)
      nut(ystart(1):yend(1), 0, ystart(3):yend(3)) =          nut(ystart(1):yend(1), 1, ystart(3):yend(3))
    END SELECT

    ! xp-dir
    SELECT CASE(myProcNghBC(y_pencil, 3))
    CASE(BC_NSLIP)
      nut(nxp, 0:nyp,ystart(3):yend(3))= two*xnu-nut(nxc, 0:nyp,ystart(3):yend(3))
    CASE(BC_FSLIP)
      nut(nxp, 0:nyp,ystart(3):yend(3))=         nut(nxc, 0:nyp,ystart(3):yend(3))
    END SELECT

    ! xm-dir
    SELECT CASE(myProcNghBC(y_pencil, 4))
    CASE(BC_NSLIP)
      nut(0,0:nyp,ystart(3):yend(3)) = two*xnu-nut(1,0:nyp,ystart(3):yend(3))
    CASE(BC_FSLIP)
      nut(0,0:nyp,ystart(3):yend(3)) =         nut(1,0:nyp,ystart(3):yend(3))
    END SELECT

    ! zp-dir
    SELECT CASE(myProcNghBC(y_pencil, 1))
    CASE(BC_NSLIP)
      nut(ystart(1):yend(1), 0:nyp, nzp) = two*xnu-nut(ystart(1):yend(1), 0:nyp, nzc)
    CASE(BC_FSLIP)
      nut(ystart(1):yend(1), 0:nyp, nzp) =         nut(ystart(1):yend(1), 0:nyp, nzc)
    END SELECT

    ! zm-dir
    SELECT CASE(myProcNghBC(y_pencil, 2))
    CASE(BC_NSLIP)
      nut(ystart(1):yend(1), 0:nyp, 0) = two*xnu-nut(ystart(1):yend(1), 0:nyp, 1)
    CASE(BC_FSLIP)
      nut(ystart(1):yend(1), 0:nyp, 0) =         nut(ystart(1):yend(1), 0:nyp, 1)
    END SELECT

    call myupdate_halo(nut, mb1, hi1)
  end subroutine SetBC_and_UpdateHalo_SGS

end module m_LESModels


text:
module sources
  ! ############################################################################
  ! Compute the sources for the solver
  ! ############################################################################

  use control, only : dp

  implicit none

  contains

  subroutine compute_source()
    ! ##########################################################################
    ! Compute the sources into group g from gp
    ! ##########################################################################

    ! Use Statements
    use state, only : mg_source, update_fission_density, sigphi, mg_sig_s, &
                      mg_mMap, mg_phi, keff, mg_chi, mg_density, mg_constant_source, &
                      scaling
    use control, only : number_cells, allow_fission, solver_type, number_groups, &
                        scatter_leg_order, use_DGM, spatial_dimension
    use dgm, only : dgm_order, phi_m, source_m

    ! Variable definitions
    integer :: &
      g,   & ! Group index
      c,   & ! Cell index
      mat, & ! Material index
      l,   & ! Legendre index
      m,   & ! Moment index
      ll,  & ! Total moment index
      ord    ! short name for scatter_legendre_order
    logical :: &
      dgm_switch    !

    ! Reset the sources
    mg_source = 0.0_8

    ! Update the fission density if needed
    if (allow_fission .or. solver_type == 'eigen') then
      if (solver_type == 'fixed' .or. dgm_order > 0) then
        call update_fission_density()
      end if
    end if

    dgm_switch = use_DGM .and. dgm_order > 0

    sigphi = 0.0_8
    ord = scatter_leg_order

    ! Compute the source
    do c = 1, number_cells
      mat = mg_mMap(c)

      ! Add the external source
      if (use_DGM) then
        mg_source(:,c) = mg_source(:,c) + source_m(:,dgm_order)
      else
        mg_source(:,c) = mg_source(:,c) + mg_constant_source
      end if

      ! Add the fission source
      if (allow_fission .or. solver_type == 'eigen') then
        mg_source(:,c) = mg_source(:,c) + scaling * mg_chi(:,mat) * mg_density(c) / keff
      end if

      ! Compute the scattering matrix
      if (spatial_dimension == 1) then
        if (dgm_switch) then
          do g = 1, number_groups
            sigphi(:,g,c) = sum(phi_m(0, 0:ord,:,c) * mg_sig_s(:,:,g,mat), 2)
          end do  ! End g loop
        else
          do g = 1, number_groups
            sigphi(:,g,c) = sum(mg_phi(0:ord,:,c) * mg_sig_s(:,:,g,mat), 2)
          end do  ! End g loop
        end if
        do l = 0, scatter_leg_order
          sigphi(l,:,c) = sigphi(l,:,c) * (2.0_8 * l + 1.0_8) * scaling
        end do
      else
        if (dgm_switch) then
          do g = 1, number_groups
            ll = 0
            do l = 0, scatter_leg_order
              do m = -l, l
                sigphi(ll,g,c) = dot_product(phi_m(0,ll,:,c), mg_sig_s(l,:,g,mat)) * scaling
                ll = ll + 1
              end do  ! End m loop
            end do  ! End l loop
          end do  ! End g loop
        else
          do g = 1, number_groups
            ll = 0
            do l = 0, scatter_leg_order
              do m = -l, l
                sigphi(ll,g,c) = dot_product(mg_phi(ll,:,c), mg_sig_s(l,:,g,mat)) * scaling
                ll = ll + 1
              end do  ! End m loop
            end do  ! End l loop
          end do  ! End g loop
        end if
      end if
    end do  ! End c loop

  end subroutine compute_source

end module sources

max_stars_repo_name: RLReed/unotran

max_issues_repo_name: RLReed/unotran

max_forks_repo_name: RLReed/unotran


      module output_ls_salt_module

      implicit none

      type salt_balance
        real :: plant = 0.          !! |kg/ha       |salt in plant foliage
        real :: soil = 0.           !! |kg/ha       |salt in soil
        real :: surq = 0.           !! |kg/ha       |amount of pesticide type lost in surface runoff in HRU
        real :: latq = 0.           !! |kg/ha       |amount of pesticide in lateral flow in HRU
        real :: tileq = 0.          !! |kg/ha       |amount of pesticide in tile flow in HRU
        real :: perc = 0.           !! |kg/ha       |amount of pesticide leached past bottom of soil
        real :: irrig = 0.          !! |kg/ha       |amount of pesticide applied on soil
      end type salt_balance
      type (salt_balance) :: saltbz

      type object_salt_balance
        type (salt_balance), dimension (:), allocatable :: salt
      end type object_salt_balance

      type (object_salt_balance), dimension (:), allocatable :: hsaltb_d
      type (object_salt_balance), dimension (:), allocatable :: hsaltb_m
      type (object_salt_balance), dimension (:), allocatable :: hsaltb_y
      type (object_salt_balance), dimension (:), allocatable :: hsaltb_a

      type (object_salt_balance), dimension (:), allocatable :: rusaltb_d
      type (object_salt_balance), dimension (:), allocatable :: rusaltb_m
      type (object_salt_balance), dimension (:), allocatable :: rusaltb_y
      type (object_salt_balance), dimension (:), allocatable :: rusaltb_a

      type (object_salt_balance) :: bsaltb_d
      type (object_salt_balance) :: bsaltb_m
      type (object_salt_balance) :: bsaltb_y
      type (object_salt_balance) :: bsaltb_a

      type output_saltbal_header
        character (len=5) :: day =      " jday"
        character (len=6) :: mo =       "   mon"
        character (len=6) :: day_mo =   "   day"
        character (len=6) :: yrc =      "    yr"
        character (len=8) :: isd =          "    unit"
        character (len=8) :: id =           "  gis_id"
        character (len=16) :: name =        "  name          "
        character (len=16) :: salt =        "  salt_ion      "
        character (len=12) :: plant =       "  plant_kg/h "
        character (len=12) :: soil =        "    soil_kg/h"
        character (len=12) :: surq =        "   surq_kg/h "
        character (len=12) :: latq =        "   latq_kg/h "
        character (len=12) :: tileq =       "  tileq_kg/h "
        character (len=12) :: perc =        "   perc_kg/h "
        character (len=16) :: irrig =       "   irrig_kg/h"
      end type output_saltbal_header
      type (output_saltbal_header) :: saltb_hdr

      interface operator (+)
        module procedure hruout_saltbal_add
      end interface

      interface operator (/)
        module procedure hruout_saltbal_div
      end interface

      interface operator (//)
        module procedure hruout_saltbal_ave
      end interface

      contains

      function hruout_saltbal_add (hru1, hru2) result (hru3)
        type (salt_balance), intent (in) :: hru1
        type (salt_balance), intent (in) :: hru2
        type (salt_balance) :: hru3
        hru3%plant = hru1%plant + hru2%plant
        hru3%soil = hru1%soil + hru2%soil
        hru3%surq = hru1%surq + hru2%surq
        hru3%latq = hru1%latq + hru2%latq
        hru3%tileq = hru1%tileq + hru2%tileq
        hru3%perc = hru1%perc + hru2%perc
        hru3%irrig = hru1%irrig + hru2%irrig
      end function hruout_saltbal_add

      function hruout_saltbal_div (hru1,const) result (hru2)
        type (salt_balance), intent (in) :: hru1
        real, intent (in) :: const
        type (salt_balance) :: hru2
        hru2%plant = hru1%plant
        hru2%soil = hru1%soil
        hru2%surq = hru1%surq / const
        hru2%latq = hru1%latq / const
        hru2%tileq = hru1%tileq / const
        hru2%perc = hru1%perc / const
        hru2%irrig = hru1%irrig / const
      end function hruout_saltbal_div

      function hruout_saltbal_ave (hru1,const) result (hru2)
        type (salt_balance), intent (in) :: hru1
        real, intent (in) :: const
        type (salt_balance) :: hru2
        hru2%plant = hru1%plant / const
        hru2%soil = hru1%soil / const
        hru2%surq = hru1%surq
        hru2%latq = hru1%latq
        hru2%tileq = hru1%tileq
        hru2%perc = hru1%perc
        hru2%irrig = hru1%irrig
      end function hruout_saltbal_ave

      end module output_ls_salt_module


text:
Module Image

   implicit none

   type RGBimage
      integer, dimension(:,:), pointer :: Red, Green, Blue
      integer                          :: width, height
   end type RGBimage

   type RGBAimage
      integer, dimension(:,:), pointer :: Red, Green, Blue, Alpha
      integer                          :: width, height
   end type RGBAimage

   type RGB
      integer :: Red, Green, Blue
   end type RGB

   type RGBA
      integer :: Red, Green, Blue, Alpha
   end type RGBA

   interface alloc_image
      module procedure alloc_imageRGB
      module procedure alloc_imageRGBA
   end interface

   interface init_image
      module procedure init_imageRGB
      module procedure init_imageRGBA
   end interface

   interface get_pixel
      module procedure get_pixelRGB
      module procedure get_pixelRGBA
   end interface

   interface set_pixel
      module procedure set_pixelRGB
      module procedure set_pixelRGBA
   end interface

   interface fill_img
      module procedure fill_imgRGB
      module procedure fill_imgRGBA
   end interface

   interface flip
      module procedure flip_RGBA
      ! module procedure flip_RGB
   end interface

   interface read_ppm
      module procedure read_ppm_RGB
      module procedure read_ppm_RGBA
   end interface read_ppm

   interface write_ppm
      module procedure write_ppm_RGB
      module procedure write_ppm_RGBA
   end interface

   interface open_image
      module procedure open_imageRGB
      module procedure open_imageRGBA
   end interface

   interface save_image
      module procedure save_imageRGB
      module procedure save_imageRGBA
   end interface

   interface operator (==)
      module procedure RGBequal
      module procedure RGBAequal
   end interface

   interface operator (/=)
      module procedure RGBnotequal
      module procedure RGBAnotequal
   end interface

   interface operator (.dist.)
      module procedure colourdistance
   end interface

   interface operator (*)
      module procedure colourmultiplyA
      module procedure colourmultiply
      module procedure colourmultiplycolour
   end interface

   interface operator (+)
      module procedure colouradd
      module procedure colouraddA
   end interface

   interface assignment (=)
      module procedure RGBimageequal
      module procedure RGBAimageequal
   end interface

   private
   public :: operator(+), assignment(=), operator(*), operator(.dist.), operator(==), operator(/=)
   public :: save_image, open_image, write_ppm, read_ppm, flip, clampInt
   public :: alloc_image, init_image, get_pixel, set_pixel, fill_img
   public :: rgb, rgba, RGBimage, RGBAimage, alpha_comp

Contains

   subroutine init_imageRGB(img)

      implicit none

      type(RGBimage), intent(OUT) :: img

      nullify(img%Red)
      nullify(img%Green)
      nullify(img%Blue)
      img%width = 0
      img%height = 0

   end subroutine init_imageRGB

  subroutine init_imageRGBA(img)

      implicit none

      type(RGBAimage), intent(OUT) :: img

      nullify(img%Red)
      nullify(img%Green)
      nullify(img%Blue)
      nullify(img%Alpha)

      img%width = 0
      img%height = 0

   end subroutine init_imageRGBA

   subroutine alloc_imageRGB(image, w, h)

      implicit none

      type(RGBimage)      :: image
      integer, intent(IN) :: w, h

      allocate(image%Red(w, h))
      allocate(image%Green(w, h))
      allocate(image%Blue(w, h))
      image%width = w
      image%height = h

   end subroutine alloc_imageRGB

   subroutine alloc_imageRGBA(image, w, h)

      implicit none

      type(RGBAimage)     :: image
      integer, intent(IN) :: w, h

      allocate(image%Red(w, h))
      allocate(image%Green(w, h))
      allocate(image%Blue(w, h))
      allocate(image%Alpha(w, h))

      image%width = w
      image%height = h

   end subroutine alloc_imageRGBA

   subroutine set_pixelRGB(img, x, y, colour)

      implicit none

      type(RGBimage), intent(INOUT) :: img
      integer,        intent(in)    :: x, y
      type(RGB),      intent(in)    :: colour

      if(x > img%width .or. x < 1)then
         return
      elseif(y > img%height .or. y < 1)then
         return
      end if

      img%Red(x, y) = colour%Red
      img%Green(x, y) = colour%Green
      img%Blue(x, y) = colour%Blue

   end subroutine set_pixelRGB

   subroutine set_pixelRGBA(img, x, y, colour)

      implicit none

      type(RGBAimage), intent(INOUT) :: img
      integer,        intent(in)     :: x, y
      type(RGBA),      intent(in)    :: colour

      if(x > img%width .or. x < 1)then
         return
      elseif(y > img%height .or. y < 1)then
         return
      end if

      img%Red(x, y) = colour%Red
      img%Green(x, y) = colour%Green
      img%Blue(x, y) = colour%Blue
      img%Alpha(x, y) = colour%Alpha

   end subroutine set_pixelRGBA

   subroutine get_pixelRGB(img, x, y, colour)

      implicit none

      type(RGBimage), intent(IN)    :: img
      type(RGB),      intent(INOUT) :: colour
      integer,        intent(in)    :: x, y

      colour%Red = img%Red(x, y)
      colour%Green = img%Green(x, y)
      colour%Blue = img%Blue(x, y)

   end subroutine get_pixelRGB

   subroutine get_pixelRGBA(img, x, y, colour)

      implicit none

      type(RGBAimage), intent(IN)    :: img
      type(RGBA),      intent(INOUT) :: colour
      integer,        intent(in)     :: x, y

      colour%Red = img%Red(x, y)
      colour%Green = img%Green(x, y)
      colour%Blue = img%Blue(x, y)
      colour%Alpha = img%Blue(x, y)

   end subroutine get_pixelRGBA

   subroutine fill_imgRGB(img, colour)

      implicit none

      type(RGBimage), intent(INOUT) :: img
      type(RGB),      intent(IN)    :: colour
      integer                       :: i, j

      do j = 1, img%height
         do i = 1, img%width
            img%Red(i,j) = colour%Red
            img%Green(i,j) = colour%Green
            img%Blue(i,j) = colour%Blue
         end do
      end do
   end subroutine fill_imgRGB

   subroutine fill_imgRGBA(img, colour)

      implicit none

      type(RGBAimage), intent(INOUT) :: img
      type(RGBA),      intent(IN)    :: colour
      integer                        :: i, j

      do j = 1, img%height
         do i = 1, img%width
            img%Red(i,j) = colour%Red
            img%Green(i,j) = colour%Green
            img%Blue(i,j) = colour%Blue
            img%Alpha(i,j) = colour%Alpha
         end do
      end do
   end subroutine fill_imgRGBA

   subroutine RGBimageequal(a1, a2)

      implicit none

      type(RGBimage), intent(OUT) :: a1
      type(RGBimage), intent(IN)  :: a2

      a1%red = a2%red
      a1%green = a2%green
      a1%blue = a2%blue

   end subroutine RGBimageequal

   subroutine RGBAimageequal(a1, a2)

      implicit none

      type(RGBAimage), intent(OUT) :: a1
      type(RGBAimage), intent(IN)  :: a2

      a1%red = a2%red
      a1%green = a2%green
      a1%blue = a2%blue
      a1%alpha = a2%alpha

   end subroutine RGBAimageequal

   logical function RGBequal(c1, c2)

      type(RGB), intent(IN) :: c1, c2

      RGBequal = .TRUE.
      if( (c1%red == c2%red) .and. (c1%Green == c2%Green) .and. (c1%Blue == c2%Blue))then
         return
      end if
      RGBequal = .FALSE.

   end function RGBequal

   logical function RGBAequal(c1, c2)

      type(RGBA), intent(IN) :: c1, c2

      RGBAequal = .TRUE.
      if( (c1%red == c2%red) .and. (c1%Green == c2%Green) .and. (c1%Blue == c2%Blue) .and. (c1%alpha == c2%alpha))then
         return
      end if
      RGBAequal = .FALSE.

   end function RGBAequal

   integer function clampInt(val, lo, hi)

      integer, intent(IN) :: lo, hi, val

      if(val < lo)then
         clampInt = lo
      elseif(val > hi)then
         clampInt = hi
      else
         clampInt = val
      end if

   end function clampInt

   logical function RGBnotequal(c1, c2)

      type(RGB), intent(IN) :: c1, c2

      RGBnotequal = .False.
      if( (c1%red == c2%red) .and. (c1%Green == c2%Green) .and. (c1%Blue == c2%Blue))then
         return
      end if
      RGBnotequal = .True.

   end function RGBnotequal

   logical function RGBAnotequal(c1, c2)

      type(RGBA), intent(IN) :: c1, c2

      RGBAnotequal = .False.
      if( (c1%red == c2%red) .and. (c1%Green == c2%Green) .and. (c1%Blue == c2%Blue) .and. (c1%alpha == c2%alpha))then
         return
      end if
      RGBAnotequal = .True.

   end function RGBAnotequal

   real function colourdistance(c1, c2)

         type(RGB), intent(IN) :: c1, c2

         colourdistance = sqrt(real(c1%red - c2%red)**2 + real(c1%Green - c2%Green)**2 &
                               + real(c1%Blue - c2%Blue)**2 )

   end function colourdistance

   function colouraddA(c1, c2)

      type(RGBA), intent(IN) :: c1, c2
      type(RGBA)             :: colouraddA

      colouraddA = RGBA(c1%red + c2%red, c1%green + c2%green, c1%blue + c2%blue, c1%alpha + c2%alpha)

   end function colouraddA

   function colouradd(c1, c2)

      type(RGB), intent(IN) :: c1, c2
      type(RGB)             :: colouradd

      colouradd = RGB(c1%red + c2%red, c1%green + c2%green, c1%blue + c2%blue)

   end function colouradd

   function colourmultiplyA(colour, scalar)

      type(RGBA),intent(IN) :: colour
      type(RGBA)            :: colourmultiplyA
      real,      intent(IN) :: scalar

      colourmultiplyA = RGBA(int(colour%red*scalar), int(colour%green*scalar), int(colour%blue*scalar), colour%alpha)

   end function colourmultiplyA

   function colourmultiply(colour, scalar)

      type(RGB),intent(IN) :: colour
      type(RGB)            :: colourmultiply
      real,      intent(IN) :: scalar

      colourmultiply = RGB(int(colour%red*scalar), int(colour%green*scalar), int(colour%blue*scalar))

   end function colourmultiply

   function colourmultiplycolour(a, b)

      type(RGB),intent(IN) :: a, b
      type(RGB)            :: colourmultiplycolour

      colourmultiplycolour = RGB(int(a%red*b%red), int(a%green*b%green), int(a%blue*b%blue))

   end function colourmultiplycolour

   function alpha_comp(ca, cb) result(co)

      implicit none

      type(RGBA), intent(IN) :: cb, ca
      type(RGBA)             :: co
      real                   :: a_tmp, b_tmp

      a_tmp = ca%alpha/255.
      b_tmp = cb%alpha/255.

      co%red =  clampInt(int(ca%red * a_tmp) + int(cb%red * (1. - a_tmp)), 0, 255)
      co%green =  clampInt(int(ca%green * a_tmp) + int(cb%green * (1. - a_tmp)), 0, 255)
      co%blue =  clampInt(int(ca%blue * a_tmp) + int(cb%blue * (1. - a_tmp)), 0, 255)
      co%alpha = clampInt(int((b_tmp + (1. - b_tmp) * a_tmp)*255.), 0, 255)

   end function alpha_comp

   subroutine read_ppm_RGB(filename, img)

      implicit none

      type(RGBImage),     intent(OUT) :: img
      integer                         :: nmax, i, j, offset
      character(len = *), intent(IN)  :: filename
      character(2)                    :: mode
      character                       :: code

      img%width = 0
      img%height = 0
      nullify(img%Red)
      nullify(img%Green)
      nullify(img%Blue)

      open(56, file = filename, access='stream',form='formatted',status= 'old')

      read(56, '(A2)') mode
      read(56, *) img%width, img%height
      read(56, *) nmax
      inquire(56,pos=offset)
      close(56)
      open(56 ,file=filename,access='stream',status='old')
      read(56, pos=offset-1) code
      call alloc_image(img, img%width, img%height)

      if(mode == 'P6')then
         do j = 1, img%height
            do i = 1, img%width
               read(56) code
               img%Red(i,j) = iachar(code)
               read(56) code
               img%Green(i,j) = iachar(code)
               read(56) code
               img%Blue(i,j) = iachar(code)
            end do
         end do
      else
         print*,'Mode not supported!'
         stop
!         do j = 1, img%height
!            do i = 1, img%width
!               read(56) img%red(i,j)
!               read(56),img%green(i,j)
!               read(56),img%blue(i,j)
!               print*,img%red(i,j),img%green(i,j),img%blue(i,j)
!               call exit(0)
!            end do
!         end do
      end if

      close(56)

   end subroutine read_ppm_RGB

   subroutine read_ppm_RGBA(filename, img)

      implicit none

      type(RGBAImage),     intent(OUT) :: img
      integer                          :: nmax, i, j, offset, io, u
      character(len = *), intent(IN)   :: filename
      character(len=2)                     :: mode
      character(len=1)                        :: code
      character(len=128) :: line

!horrible code fix!!!!!!!!!!!!!!!

      call init_image(img)

      open(newunit=u, file = filename, access='stream',form='formatted',status= 'old', iostat=io)
      if(io /= 0)then
         print*,'Cant open '//trim(filename)
      else
         read(u, '(A2)') mode
         if(mode == 'P6')then
            read(u, *) img%width, img%height
            read(u, *) nmax
            inquire(u,pos=offset)
            close(u)
            open(newunit=u ,file=filename,access='stream',status='old')
            read(u, pos=offset-1) code
            call alloc_image(img, img%width, img%height)
         elseif(mode == 'P7')then
            read(u,*)line,img%width
            read(u,*)line,img%height
            read(u,*)line!depth
            read(u,*)line!maxval
            read(u,*)line!tuple
            read(u,*)line!endhdr
            inquire(u,pos=offset)
            close(u)
            open(newunit=u ,file=filename,access='stream',status='old')
            read(u,pos=offset-1)code
            call alloc_image(img, img%width, img%height)
         end if
         if(mode == 'P6')then
            do j = 1, img%height
               do i = 1, img%width
                  read(u) code
                  img%Red(i,j) = iachar(code)
                  read(u) code
                  img%Green(i,j) = iachar(code)
                  read(u) code
                  img%Blue(i,j) = iachar(code)
               end do
            end do
         elseif(mode == 'P7')then
            do j = 1, img%height
               do i = 1, img%width
                  read(u) code
                  img%Red(i,j) = iachar(code)
                  read(u) code
                  img%Green(i,j) = iachar(code)
                  read(u) code
                  img%Blue(i,j) = iachar(code)
                  read(u) code
                  img%alpha(i,j) = iachar(code)
               end do
            end do
         else
            print*,'Mode not supported!'
            stop
         end if
      end if
      close(u)

   end subroutine read_ppm_RGBA

   subroutine flip_RGBA(img)

      implicit none

      type(RGBAimage), intent(INOUT) :: img
      type(RGBAimage) :: tmp

      integer :: i, j,k

      call init_image(tmp)
      call alloc_image(tmp, img%width, img%height)
      k = 1
      do j = img%height,1, -1
         do i = 1, img%width
            tmp%red(i, k) = img%red(i, j)
            tmp%green(i, k) = img%green(i, j)
            tmp%blue(i, k) = img%blue(i, j)
            tmp%alpha(i, k) = img%alpha(i, j)
         end do
      k = k + 1
      end do

      img = tmp

   end subroutine flip_RGBA

   subroutine write_ppm_RGB(filename, img, mode)

      implicit none

      type(RGBimage),             intent(in) :: img
      character(len=*),           intent(in) :: filename
      character(len=*), optional, intent(in) :: mode
      integer                                :: i, j, u
      logical                                :: flag

      if(present(mode))then
         if(mode == 'P6')then
            flag = .TRUE.
         else
            flag=.FALSE.
         end if
      else
         flag = .FALSE.
      end if

      open(newunit=u, file = filename)

      if(flag)then
         write(u, '(A2)') 'P6'
      else
         write(u, '(A2)') mode
         print*,mode
      end if

      write(u, '(i0, " ",i0)') img%width, img%height
      write(u, '(i0)') 255

      if(flag)then
         do j =  img%height,1, -1
            do i = 1, img%width
               write(u, '(3A1)', advance='no') achar(img%Red(i,j)), achar(img%Green(i,j)), achar(img%Blue(i,j))
            end do
         end do
      else
         do j = 1, img%height
            do i = 1, img%width
               write(u, '(3(I3.1,1X))',advance='no') img%Red(i,j), img%Green(i,j), img%Blue(i,j)
            end do
         end do
      end if
      close(u)

   end subroutine write_ppm_RGB

   subroutine write_ppm_RGBA(filename, img, mode)

      use iso_fortran_env, only : int8
      use utils,           only : str

      implicit none

      type(RGBAImage),             intent(in) :: img
      character(len=*),           intent(in)  :: filename
      character(len=*), optional, intent(in)  :: mode
      integer                                 :: i, j, u
      logical                                 :: flag
      integer(kind=int8) :: r,g,b,a

      if(present(mode))then
         if(mode == 'P7' .or. mode == 'P6')then
            flag = .TRUE.
         else
            flag = .FALSE.
         end if
      else
         flag = .FALSE.
      end if

      open(newunit=u, file=filename, status='replace')

      if(flag)then
         write(u, '(A2)') mode
         if(mode == 'P7')then
            write(u,*)'WIDTH '//str(img%width)
            write(u,*)'HEIGHT '//str(img%height)
            write(u,*)'DEPTH 4'
            write(u,*)'MAXVAL '//str(max(maxval(img%red),maxval(img%green),maxval(img%blue)))
            write(u,*)'TUPLTYPE RGB_ALPHA'
            write(u,*)'ENDHDR'
         else
            write(u, '(i0, " ",i0)') img%width, img%height
            write(u, '(i0)') 255
         end if
      else
         write(u, '(A2)') mode
      end if

      close(u)

      if(flag)then
         open(newunit=u, file = filename, access='stream', form='unformatted', position='append')
         do j = 1, img%height
            do i = 1, img%width
               r = img%Red(i,j); g = img%Green(i,j); b = img%Blue(i,j); a=img%alpha(i,j)
               if(mode == 'P7')write(u)r,g,b,a
               if(mode == 'P6')write(u)r,g,b
            end do
         end do
      else
         open(newunit=u, file = filename, position='append')
         do j = 1, img%height
            do i = 1, img%width
               write(u, '(3(I3.1,1X))',advance='no') img%Red(i,j), img%Green(i,j), img%Blue(i,j)
            end do
         end do
      end if
      close(u)

   end subroutine write_ppm_RGBA

   subroutine open_imageRGB(img, filename, format)

      implicit none

      type(RGBimage), intent(OUT) :: img
      character(*),   intent(IN)  :: filename, format

      call execute_command_line('convert '//filename//format//' '//filename//'.ppm')
      call read_ppm(filename//'.ppm', img)
      call execute_command_line('rm '//filename//'.ppm')

   end subroutine open_imageRGB

   subroutine open_imageRGBA(img, filename, format)

      implicit none

      type(RGBAimage), intent(OUT) :: img
      character(*),   intent(IN)  :: filename, format

      call execute_command_line('convert '//filename//format//' '//filename//'.ppm')
      call read_ppm(filename//'.ppm', img)
      call execute_command_line('rm '//filename//'.ppm')

   end subroutine open_imageRGBA

   subroutine save_imageRGB(img, filename, format)

      implicit none

      type(RGBimage), intent(IN) :: img
      character(*),   intent(IN) :: filename, format

      print*,'Saved image as: ',filename//format
      call write_ppm(filename//'.ppm', img, 'P6')
      call execute_command_line('convert '//filename//'.ppm '//filename//format)
      call execute_command_line('rm '//filename//'.ppm')

   end subroutine save_imageRGB

   subroutine save_imageRGBA(img, filename, format)

      implicit none

      type(RGBAimage), intent(IN) :: img
      character(*),   intent(IN) :: filename, format

      print*,'Saved image at: ',filename//format
      call write_ppm(filename//'.ppm', img, 'P6')
      call execute_command_line('convert '//filename//'.ppm '//filename//format)
      call execute_command_line('rm '//filename//'.ppm')

   end subroutine save_imageRGBA
end module image

max_stars_repo_name: lewisfish/FortGL

max_issues_repo_name: lewisfish/FortGL

max_forks_repo_name: lewisfish/FortGL


text:
C File: arrays-sect-05.f
C Illustrates the use of array sections using an implicit do loop over
C all of the array elements.  The do loop does not explicitly specify
C the bounds of the iteration.

      program main
      implicit none

      integer, dimension (5) :: A
      integer :: i

      do i = 1, 5
          A(i) = 0
      end do

      write (*,10) A(1), A(2), A(3), A(4), A(5)

 10   format(5(I5,X))

C     Use of a more complex array section to redefine all elements of A
C     This example uses an implicit do loop.
      A(:) = 17

      write (*,10) A(1), A(2), A(3), A(4), A(5)

      stop
      end program main

ext: f


text:
	SUBROUTINE HSWIND  ( iwnd, size, iwidth, itype, hdsiz, iret )
C************************************************************************
C* HSWIND - VG								*
C*									*
C* This subroutine sets the wind attributes.				*
C*									*
C* HSWIND  ( IWND, SIZE, IWIDTH, ITYPE, HDSIZ, IRET )			*
C*									*
C* Input parameters:							*
C*	IWND		INTEGER		Wind category			*
C*					  1 = Barbs			*
C*					  2 = Arrows			*
C*	SIZE		REAL		Wind size			*
C*	IWIDTH		INTEGER		Wind line width			*
C*	ITYPE		INTEGER		Wind type			*
C*	HDSIZ		REAL		Wind arrow head size		*
C*									*
C* Output parameters:							*
C*	IRET		INTEGER		Return code			*
C**									*
C* Log:									*
C* S. Jacobs/NCEP	 3/97						*
C* S. Maxwell/GSC        6/97   Documentation changes                   *
C************************************************************************
	INCLUDE		'ERROR.PRM'
	INCLUDE		'DEVCHR.CMN'
	INCLUDE		'DVWNDW.CMN'
C------------------------------------------------------------------------
	iret = NORMAL
C
	CALL VSWIND ( iwnd, size, iwidth, itype, hdsiz, iret )
C*
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
C-------GF 28.08.1980------------------------ M E M S E T S
      INTEGER*2 FSET,SETHIB
      INTEGER*2 SET   (XXSETH),SETELE(XXSETH)
     = ,EMITSH,EMITST,EMITSB   ! STA: FOR EMITT, SOURCE
     = ,EMITTH,EMITTT,EMITTB   ! STA: FOR EMITT, TARGET
     = ,LAGARH,LAGART,LAGARB   ! STA: FOR LAGAR
     = ,LAGETH,LAGETT,LAGETB   ! STA: FOR LAGET
     = ,LAPUTH,LAPUTT,LAPUTB   ! STA: FOR LAPUT
     = ,SPLITH,SPLITT,SPLITB   ! STA: FOR SPLIT
     = ,STACLH,STACLT,STACLB   ! SYM: FOR STACLO
     = ,STAGAH,STAGAT,STAGAB   ! STA: FOR STAGAR, ILL STATES
     = ,STAGBH,STAGBT,STAGBB   ! STA: FOR STAGAR, SURVIVING STATES
     = ,STASUH,STASUT,STASUB   ! STA: FOR STASUC
     = ,SYMINH,SYMINT,SYMINB   ! SYM: FOR SYMINS
      COMMON /SETS/
     =  FSET
     = ,SETHIB
     = ,EMITSH,EMITST,EMITSB   ! STA: FOR EMITT, SOURCE
     = ,EMITTH,EMITTT,EMITTB   ! STA: FOR EMITT, TARGET
     = ,LAGARH,LAGART,LAGARB   ! STA: FOR LAGAR
     = ,LAGETH,LAGETT,LAGETB   ! STA: FOR LAGET
     = ,LAPUTH,LAPUTT,LAPUTB   ! STA: FOR LAPUT
     = ,SPLITH,SPLITT,SPLITB   ! STA: FOR SPLIT
     = ,STACLH,STACLT,STACLB   ! SYM: FOR STACLO
     = ,STAGAH,STAGAT,STAGAB   ! STA: FOR STAGAR, ILL STATES
     = ,STAGBH,STAGBT,STAGBB   ! STA: FOR STAGAR, SURVIVING STATES
     = ,STASUH,STASUT,STASUB   ! STA: FOR STASUC
     = ,SYMINH,SYMINT,SYMINB   ! SYM: FOR SYMINS
      COMMON /SETTRA/ SETELE
      EQUIVALENCE (SET(1),MEMNUC(1))

max_stars_repo_name: gfis/jextra

max_issues_repo_name: gfis/jextra

max_forks_repo_name: gfis/jextra


          program demo_assert
          use M_verify, only : assert
          implicit none
          real :: a, toobig=1024
          a=2000
          call assert('myroutine', 101, a.gt.toobig, 'The value is too large', a, '.gt.', toobig)
          end program demo_assert


text:
!!
!! The Intel compiler is stumbling over the type-is statement
!! in the following example when it shouldn't be.
!!
!! % ifort --version
!! ifort (IFORT) 13.1.2 20130514
!!
!! % ifort intel-bug-20130802.f90
!! intel-bug-20130802.f90(21): error #6402: prPromoteSym : Illegal KIND & CLASS mix   [KIND]
!!   type is (integer(kind(val)))
!! -------------------^
!! compilation aborted for intel-bug-20130802.f90 (code 1)
!!

subroutine fubar (cs, val)

  class(*) :: cs
  integer :: val

  select type (cs)
  type is (integer(kind(val)))
  end select

end subroutine

max_stars_repo_name: dhnza/fortran-compiler-tests

max_issues_repo_name: dhnza/fortran-compiler-tests

max_forks_repo_name: dhnza/fortran-compiler-tests


text:
      module rrlw_kg07

      use shr_kind_mod, only: r8 => shr_kind_r8

!      use parkind ,only : jpim, jprb

      implicit none
      save

!-----------------------------------------------------------------
! rrtmg_lw ORIGINAL abs. coefficients for interval 7
! band 7:  980-1080 cm-1 (low - h2o,o3; high - o3)
!
! Initial version:  JJMorcrette, ECMWF, jul1998
! Revised: MJIacono, AER, jun2006
!-----------------------------------------------------------------
!
!  name     type     purpose
!  ----   : ----   : ---------------------------------------------
!fracrefao: real
!fracrefbo: real
! kao     : real
! kbo     : real
! kao_mco2: real
! kbo_mco2: real
! selfrefo: real
! forrefo : real
!-----------------------------------------------------------------

      integer, parameter :: no7  = 16

      real(kind=r8) , dimension(no7) :: fracrefbo
      real(kind=r8) :: fracrefao(no7,9)
      real(kind=r8) :: kao(9,5,13,no7)
      real(kind=r8) :: kbo(5,13:59,no7)
      real(kind=r8) :: kao_mco2(9,19,no7)
      real(kind=r8) :: kbo_mco2(19,no7)
      real(kind=r8) :: selfrefo(10,no7)
      real(kind=r8) :: forrefo(4,no7)

!-----------------------------------------------------------------
! rrtmg_lw COMBINED abs. coefficients for interval 7
! band 7:  980-1080 cm-1 (low - h2o,o3; high - o3)
!
! Initial version:  JJMorcrette, ECMWF, jul1998
! Revised: MJIacono, AER, jun2006
!-----------------------------------------------------------------
!
!  name     type     purpose
!  ----   : ----   : ---------------------------------------------
!fracrefa : real
!fracrefb : real
! ka      : real
! kb      : real
! ka_mco2 : real
! kb_mco2 : real
! selfref : real
! forref  : real
!
! absa    : real
!-----------------------------------------------------------------

      integer, parameter :: ng7  = 12

      real(kind=r8) , dimension(ng7) :: fracrefb
      real(kind=r8) :: fracrefa(ng7,9)
      real(kind=r8) :: ka(9,5,13,ng7) ,absa(585,ng7)
      real(kind=r8) :: kb(5,13:59,ng7),absb(235,ng7)
      real(kind=r8) :: ka_mco2(9,19,ng7)
      real(kind=r8) :: kb_mco2(19,ng7)
      real(kind=r8) :: selfref(10,ng7)
      real(kind=r8) :: forref(4,ng7)

      equivalence (ka(1,1,1,1),absa(1,1)),(kb(1,13,1),absb(1,1))

      end module rrlw_kg07

max_stars_repo_licenses: zlib-acknowledgement FTL RSA-MD

max_issues_repo_licenses: zlib-acknowledgement FTL RSA-MD

max_forks_repo_licenses: zlib-acknowledgement FTL RSA-MD


text:
! THIS VERSION: GALAHAD 2.1 - 22/03/2007 AT 09:00 GMT.

!-*-*-*-*-*-*-*-*  G A L A H A D _ S M T   M O D U L E  *-*-*-*-*-*-*-*-*

!  Copyright reserved, Gould/Orban/Toint, for GALAHAD productions
!  Principal author: Nick Gould

!  History -
!   originally released pre GALAHAD Version 1.0. December 1st 1997
!   update released with GALAHAD Version 2.0. February 16th 2005

!  For full documentation, see
!   http://galahad.rl.ac.uk/galahad-www/specs.html

   MODULE GALAHAD_SMT_double

!  ==========================
!  Sparse matrix derived type
!  ==========================

     IMPLICIT NONE

!---------------------
!   P r e c i s i o n
!---------------------

     INTEGER, PARAMETER :: wp = KIND( 1.0D+0 )

     PRIVATE
     PUBLIC :: SMT_put, SMT_get

     TYPE, PUBLIC :: SMT_type
       INTEGER :: m, n, ne
       CHARACTER, ALLOCATABLE, DIMENSION( : ) :: id
       CHARACTER, ALLOCATABLE, DIMENSION( : ) :: type
       INTEGER, ALLOCATABLE, DIMENSION( : ) :: row
       INTEGER, ALLOCATABLE, DIMENSION( : ) :: col
       INTEGER, ALLOCATABLE, DIMENSION( : ) :: ptr
       REAL ( KIND = wp ), ALLOCATABLE, DIMENSION( : ) :: val
     END TYPE

   CONTAINS

!-*-*-*-*-*-*-*-*-   S M T  _ p u t   S U B R O U T I N E  -*-*-*-*-*-*-*-*-

     SUBROUTINE SMT_put( array, string, stat )

!  Dummy arguments

     CHARACTER, ALLOCATABLE, DIMENSION( : ) :: array
     CHARACTER ( len = * ), INTENT( IN ), OPTIONAL :: string
     INTEGER, INTENT( OUT ), OPTIONAL ::  stat

!  Local variables

     INTEGER :: i, l
     LOGICAL :: ok

     l = 0
     IF ( PRESENT( string ) ) l = LEN_TRIM( string )
     IF ( PRESENT( stat ) ) THEN
!      IF ( ALLOCATED( array ) ) DEALLOCATE( array, STAT = stat )
       ALLOCATE( array( l ), STAT = stat )
       ok = stat == 0
     ELSE
!      IF ( ALLOCATED( array ) ) DEALLOCATE( array )
       ALLOCATE( array( l ) )
       ok = .true.
     END IF
     IF ( ok ) THEN
       DO i = 1, l
         array( i ) = string( i : i )
       END do
     END IF

     RETURN

!  End of SMT_put

     END SUBROUTINE SMT_put

!-*-*-*-*-*-*-*-*-*-*-   S M T  _ g e t   F U N C T I O N  -*-*-*-*-*-*-*-*-*-

     FUNCTION SMT_get( array )

!  Dummy arguments

     CHARACTER, DIMENSION( : ) :: array
     CHARACTER( SIZE( array ) ) :: SMT_get

!  Local variables

     INTEGER :: i

     DO i = 1, SIZE( array )
        SMT_get( i : i ) = array( i )
     END DO

     RETURN

!  End of function SMT_get

     END FUNCTION SMT_get

!  End of module SMT

   END MODULE GALAHAD_SMT_double

max_forks_repo_name: dynamics-of-stellar-systems/dynamite


text:
Real Function ptdis(x)
  Save
  b = 3.78
  c = 0.47
  d = 3.60
  ptdis = 1./(2.*b)*(1.-exp(-b*x**2)) - c/d*x*exp(-d*x) - c/d**2*(exp(-d*x)-1.)
  Return
End Function ptdis

max_stars_repo_name: xiaohaijin/AMPT

max_issues_repo_name: xiaohaijin/AMPT

max_forks_repo_name: xiaohaijin/AMPT


! Module providing geometry-related routines for AGN feedback.
!
! Paul Ricker (2/2013)

module agn_geometry

    use amrex_fort_module, only : rt => amrex_real
    use bl_constants_module, only: M_PI

!==============================================================================

implicit none

public

! Constants for geometric shapes

integer, parameter :: INSIDE = 0, STRADDLES = 1, OUTSIDE = 2

contains

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

! loc_sphere: For a Cartesian zone with given left and right coordinate values,
!             determine whether it lies entirely inside, entirely outside, or
!             partially inside and outside a given sphere.
!
! arguments: xl, yl, zl             Left-face coordinates of zone relative to
!                                   sphere's center
!            xr, yr, zr             Right-face coordinates of zone relative to
!                                   sphere's center
!            params(1)              Parameter array containing radius of sphere

integer function loc_sphere(xl, xr, yl, yr, zl, zr, R)

real(rt), intent(in) :: xl, xr, yl, yr, zl, zr, R

real(rt)             :: rl2, rr2, R2

R2  = R**2
rl2 = xl**2 + yl**2 + zl**2
rr2 = xr**2 + yr**2 + zr**2

if (rr2 <= R2) then
  loc_sphere = INSIDE
else if (rl2 > R2) then
  loc_sphere = OUTSIDE
else
  loc_sphere = STRADDLES
endif

return
end function loc_sphere

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

! loc_cylinder: For a Cartesian zone with given left and right coordinate
!               values, determine whether it lies entirely inside, entirely
!               outside, or partially inside and outside a given cylinder.
!
! arguments: xl, yl, zl             Left-face coordinates of zone relative to
!                                   cylinder's origin point
!            xr, yr, zr             Right-face coordinates of zone relative to
!                                   cylinder's origin point
!            params(5)              Parameter array containing radius, height,
!                                   and orientation direction cosines of
!                                   cylinder

integer function loc_cylinder(xl, xr, yl, yr, zl, zr, R, h, nxc, nyc, nzc)

real(rt), intent(in) :: xl, xr, yl, yr, zl, zr
real(rt), intent(in) :: R, h, nxc, nyc, nzc

real(rt)             :: R2, xldotn, xrdotn, rl2, rr2
real(rt)             :: xlproj, xrproj, ylproj, yrproj, zlproj, zrproj

R2  = R**2

xldotn = xl*nxc + yl*nyc + zl*nzc
xrdotn = xr*nxc + yr*nyc + zr*nzc

xlproj = xl - xldotn*nxc
xrproj = xr - xrdotn*nxc
ylproj = yl - xldotn*nyc
yrproj = yr - xrdotn*nyc
zlproj = zl - xldotn*nzc
zrproj = zr - xrdotn*nzc

rl2 = xlproj**2 + ylproj**2 + zlproj**2
rr2 = xrproj**2 + yrproj**2 + zrproj**2

if ((xldotn >= 0.) .and. (xrdotn <= h) .and. (rr2 <= R2)) then
  loc_cylinder = INSIDE
else if ((xrdotn < 0.) .or. (xldotn > h) .or. (rl2 > R2)) then
  loc_cylinder = OUTSIDE
else
  loc_cylinder = STRADDLES
endif

return
end function loc_cylinder

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

! weight_uniform: A uniform weighting function.
!
! arguments: x, y, z            Coordinates of sample point (ignored)
!            params(:)          Weighting function parameters (ignored)

real(rt) function weight_uniform(x, y, z, params)

real(rt), intent(in) :: x, y, z, params(:)

weight_uniform = 1.

return
end function weight_uniform

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

! weight_cattaneo: Cattaneo & Teyssier weighting function.
!
! arguments: x, y, z            Coordinates of sample point
!            params(:)          Weighting function parameters: R, h, nx, ny, nz, rscale

real(rt) function weight_cattaneo(x, y, z, R, height, nx, ny, nz, rscale)

real(rt), intent(in) :: x, y, z, R, height, nx, ny, nz, rscale

real(rt)             :: xdotn, xproj, yproj, zproj, r2
real(rt)             :: vol_cylinder

xdotn = x*nx + y*ny + z*nz
xproj = x - xdotn*nx
yproj = y - xdotn*ny
zproj = z - xdotn*nz
r2    = xproj**2 + yproj**2 + zproj**2

vol_cylinder = M_PI * R**2 * height

if ((xdotn >= 0.) .and. (xdotn <= height) .and. (r2 <= R**2)) then
  weight_cattaneo = dexp(-r2/(2*rscale**2)) * &
                    xdotn/height**2 / (2*M_PI*rscale**2) /(1.d0-dexp(-0.5*(R/rscale**2))) * &
                    vol_cylinder
else
  weight_cattaneo = 0.
endif

return
end function weight_cattaneo

!==============================================================================

end module agn_geometry


! { dg-do compile }
! Test the patch for PR36374 in which the different
! symbols for 'foobar' would be incorrectly flagged as
! ambiguous in foo_mod.
!
! Contributed by Salvatore Filippone  <sfilippone@uniroma2.it>
!
module s_foo_mod
  type s_foo_type
    real(kind(1.e0)) :: v
  end type s_foo_type
  interface foobar
    subroutine s_foobar(x)
      import
      type(s_foo_type), intent (inout) :: x
    end subroutine s_foobar
  end interface
end module s_foo_mod

module d_foo_mod
  type d_foo_type
    real(kind(1.d0)) :: v
  end type d_foo_type
  interface foobar
    subroutine d_foobar(x)
      import
      type(d_foo_type), intent (inout) :: x
    end subroutine d_foobar
  end interface
end module d_foo_mod

module foo_mod
  use s_foo_mod
  use d_foo_mod
end module foo_mod

subroutine s_foobar2(x)
  use foo_mod
end subroutine s_foobar2


! RUN: %python %S/test_folding.py %s %flang_fc1
! Tests folding of SCALE()
module m
  logical, parameter :: test_1 = scale(1.0, 1) == 2.0
  logical, parameter :: test_2 = scale(0.0, 1) == 0.0
  logical, parameter :: test_3 = sign(1.0, scale(-0.0, 1)) == -1.0
  logical, parameter :: test_4 = sign(1.0, scale(0.0, 0)) == 1.0
  logical, parameter :: test_5 = scale(1.0, -1) == 0.5
  logical, parameter :: test_6 = scale(2.0, -1) == 1.0
end module


text:
      subroutine clock(t)
      real*8 t
c ibm-rs clock function
      t=mclock()/100.
      return
      end

max_stars_repo_name: HRR-Sergkien/Ipydoppler

max_issues_repo_name: HRR-Sergkien/Ipydoppler

max_forks_repo_name: HRR-Sergkien/Ipydoppler


text:
      program ireal01

C     Bug found in spice.f from the PerfectClub Benchmark

C     i=MIN(MAX(k0, k1), k2) should not generate any constraint
C     but I <= k2

      real value(100)

      i = 1.2
      i = min(1.2, k)
      LWIDTH=MIN(MAX(IFIELD,72),132)
      LWIDTH=MIN(MAX(IFIELD,72.0D0),132)

      LWIDTH=MIN(MAX(IFIELD,72),132.0D0)

      LWIDTH=MIN(MAX(IFIELD,72.0D0),132.0D0)
      LWIDTH=MIN(MAX((IFIELD+IFLD),72.0D0),132.0D0)
      LWIDTH=MIN(MAX(VALUE(IFIELD+IFLD),72.0D0),132.0D0)

      print *, i

      end


text:
      SUBROUTINE BMG2_SymStd_COPY_cg_rV_G_L(
     &                Q_SER, NGx, NGy,
     &                Q, NLx, NLy, iGs, jGs,
     &                ProcCoord, NProcI, NProcJ,
     &                Nproc, MyProc
     &                )

C ==========================================================================
C  --------------------
C   DESCRIPTION:
C  --------------------
C
C   Copy the solution of the coarse-grid solve, back into the
C   MPI-based BoxMG solution vector.
C
C =======================================================================
C $license_flag$
C =======================================================================
C  --------------------
C   INPUT:
C  --------------------
C
C
C =======================================================================
C  --------------------
C   OUTPUT:
C  --------------------
C
C
C
C =======================================================================
C  --------------------
C   LOCAL:
C  --------------------
C
C
C ==========================================================================

      IMPLICIT NONE

C ----------------------------
C     Argument Declarations
C
      !
      !  Global/Local indexing
      !
      INTEGER NLx, NLy,
     &        NGx, NGy,
     &        iGs, jGs

      !
      !  Global Solution: SER
      !
      REAL*8   Q_SER(NGx,NGy)
      !
      !  Local Solution: MPI
      !
      REAL*8   Q(NLx,NLy)

      !
      !  Processor grid
      !
      INTEGER NProcI, NProcJ, NProcK, NProc, MyProc
      INTEGER ProcCoord(2, NProc)

C ----------------------------
C     Local Declarations
C
      INTEGER  iG, iL, jG, jL, kG, kL,
     &         MyProcI, MyProcJ, MyProcK

      INTEGER  iBEG, iEND, jBEG, jEND, kBEG, kEND

C =========================================================================

      MyProcI = ProcCoord(1,MyProc)
      MyProcJ = ProcCoord(2,MyProc)

      !
      !  Setup loop boundaries in x
      !
      IF ( MyProcI.EQ.1 ) THEN
         iBEG = 2
      ELSE
         iBEG = 1
      END IF

      IF ( MyProcI.EQ.NProcI ) THEN
         iEND = NLx-1
      ELSE
         iEND = NLx
      END IF

      !
      !  Setup loop boundaries in y
      !
      IF ( MyProcJ.EQ.1 ) THEN
         jBEG = 2
      ELSE
         jBEG = 1
      END IF

      IF ( MyProcJ.EQ.NProcJ ) THEN
         jEND = NLy-1
      ELSE
         jEND = NLy
      ENDIF

      DO jL=jBEG, jEND
         jG = jGs + jL - 1
         DO iL=iBEG, iEND
            iG = iGs + iL - 1
            Q(iL,jL) = Q_SER(iG,jG)
         ENDDO
      ENDDO

C =========================================================================

      RETURN
      END

max_stars_repo_name: cedar-framework/cedar

max_issues_repo_name: cedar-framework/cedar

max_forks_repo_name: cedar-framework/cedar


text:
      write (*,*)
!    !"AA"
      end

max_stars_repo_name: LChernigovskaya/fortran-plugin

max_issues_repo_name: LChernigovskaya/fortran-plugin

max_forks_repo_name: LChernigovskaya/fortran-plugin


text:
#include "eiscor.h"
subroutine u_test_passed(TIME)

  implicit none

  ! input variables
  real(8), intent(in) :: TIME

  ! print failure
  write(STDERR,'(a,F10.5,a)') 'PASSED in ',TIME,' secs'

end subroutine u_test_passed

max_stars_repo_name: trcameron/FPML

max_issues_repo_name: trcameron/FPML

max_forks_repo_name: trcameron/FPML


text:
!-----------------------------------------------------------------------
subroutine tridiagonal(N,au,bu,cu,du,fi,lt,value)
!
! !DESCRIPTION:
! A linear equation with tridiagonal matrix structure is solved here. The main
! diagonal is stored on {\tt bu}, the upper diagonal on {\tt au}, and the
! lower diagonal on {\tt cu}, the right hand side is stored on {\tt du}.
! The method used here is the simplified Gauss elimination, also called
! \emph{Thomas algorithm}.
!
! !USES:
   IMPLICIT NONE
!
! !INPUT PARAMETERS:
   integer, intent(in)                 :: N,fi,lt
   real(8), dimension(N), intent(in)   :: au, bu, cu, du
!
! !OUTPUT PARAMETERS:
   real(8), intent(out)                :: value(N)
!
!EOP
!
! !LOCAL VARIABLES:
   integer                             :: i
   real(8), dimension(N)               :: ru, qu
!
!-----------------------------------------------------------------------
!BOC
   ru(lt)=au(lt)/bu(lt)
   qu(lt)=du(lt)/bu(lt)

   do i=lt-1,fi+1,-1
      ru(i)=au(i)/(bu(i)-cu(i)*ru(i+1))
      qu(i)=(du(i)-cu(i)*qu(i+1))/(bu(i)-cu(i)*ru(i+1))
   end do

   qu(fi)=(du(fi)-cu(fi)*qu(fi+1))/(bu(fi)-cu(fi)*ru(fi+1))

   value(fi)=qu(fi)

!   write(6,*) 'fi = ',fi
!   write(6,*) 'value(fi) = ',value(fi)

   do i=fi+1,lt
      value(i)=qu(i)-ru(i)*value(i-1)
   end do

!   write(6,*) 'value(lt) = ',value(lt)

   return
end subroutine tridiagonal

max_stars_repo_name: BingzhangChen/HOT

max_issues_repo_name: BingzhangChen/HOT

max_forks_repo_name: BingzhangChen/HOT