text:
PROGRAM trim_copy
IMPLICIT NONE
CHARACTER(len=32) :: in_name, out_name, close_status
CHARACTER(len=160) :: line
CHARACTER(len=80) :: msg
INTEGER :: istat
LOGICAL :: lexist = .FALSE.
CHARACTER :: choice

WRITE(*, *) "Enter input filename: "
READ(*, *) in_name

OPEN(UNIT=7, FILE=in_name, STATUS='OLD', ACTION='READ', IOSTAT=istat, IOMSG=msg)
!INQUIRE(UNIT=7, EXIST=lexist)
input_open_ok: IF (istat /= 0) THEN
  WRITE(*, *) "Error: Could not open file ", TRIM(in_name)
  WRITE(*, *) msg
ELSE
  WRITE(*, *) "Enter output filename: "
  READ(*, *) out_name
  INQUIRE(FILE=out_name, EXIST=lexist)
  exists: IF (lexist) THEN
    WRITE(*, *) "File ", TRIM(out_name), " already exists.", &
                " Would you like to overwrite it? "
    READ(*, *) choice
    overwrite: IF (choice == 'Y') THEN
      OPEN(UNIT=8, FILE=out_name, STATUS='REPLACE', ACTION='WRITE', &
      IOSTAT=istat, IOMSG=msg)
    ELSE
      WRITE(*, *) "Exiting program."
    END IF overwrite
  ELSE
    OPEN(UNIT=8, FILE=out_name, STATUS='NEW', ACTION='WRITE', &
    IOSTAT=istat, IOMSG=msg)
  END IF exists
  output_open_ok: IF (istat /= 0) THEN
    WRITE(*, *) "Error: Could not open file ", TRIM(in_name)
    WRITE(*, *) msg
  ELSE
    DO
      IF (istat /= 0) EXIT
      READ(7, '(A160)', IOSTAT=istat) line
      IF (istat /= 0) EXIT
      WRITE(8, *, IOSTAT=istat) TRIM(line)
    END DO
    WRITE(*, *) "Copy complete. Would you like to delete ", TRIM(in_name), "?"
    READ(*, *) choice
    IF (choice == 'Y') THEN
      close_status = 'DELETE'
    ELSE
      close_status = 'KEEP'
    END IF
    CLOSE(7, STATUS=close_status)
    CLOSE(8)
  END IF output_open_ok
END IF input_open_ok
END PROGRAM trim_copy

max_stars_repo_name: evanmacbride/fortran-practice

max_issues_repo_name: evanmacbride/fortran-practice

max_forks_repo_name: evanmacbride/fortran-practice


SUBROUTINE formke(km,kp,c,ke,theta)
!
! This subroutine creates the ke matrix for incremental Biot.
!
 IMPLICIT NONE
 INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15)
 REAL(iwp),INTENT(IN)::km(:,:),kp(:,:),c(:,:),theta
 REAL(iwp),INTENT(OUT)::ke(:,:)
 INTEGER::ndof
 ndof=UBOUND(km,1)
 ke(:ndof,:ndof)=km
 ke(:ndof,ndof+1:)=c
 ke(ndof+1:,:ndof)=TRANSPOSE(c)
 ke(ndof+1:,ndof+1:)=-theta*kp
RETURN
END SUBROUTINE formke


text:
program t
  implicit none
! open stmt replace status
  character*22::a,b

  open (unit=95, file='tmpfile', status='new')
  write (95, '(a)') 'this is not the answer'
  rewind (95)
  read (95, '(a)') a
  close (95)

  open(unit=85, file='tmpfile', status='replace')
  write (85, '(a)') 'but this is the answer'
  close(85)

  open(85,file='tmpfile',status='old')
  read (85, '(a)') b
  print *,a
  print *,b

  close(85,status='delete')
endprogram t

max_stars_repo_name: maddenp/ppp

max_issues_repo_name: maddenp/ppp

max_forks_repo_name: maddenp/ppp


text:
Cafe California closed on December 17, 2005. It is now Aioli Bodega Espanola.

The modest white doors concealed a surprisingly open and roomy dining area, decked with Art Galleries local artwork. Filled with Upscale Restaurants fancy meals, the Cafe Californias menu contained things like Cajun Style Pork Ribs, Lamb Stew, and New York Steak, with prices ranging from $8.95 to $23ish for the main course alone. The entrees came with soup of the day or house salad, which was excellent.  Most of their salads and desserts contain at least one edible flower. Their music was the Big Band and Swing channel from megachannel tv service. On certain Fridays the Dick Livingston Trio performed classic and original jazz, blues and beebop music. They had a full bar in the same room as the dining area.

In December 2005, Cafe California was sold to Sacramento restaurant owners Reda Bellarbi and Jennifer Sparks, who will change it to an upscale Spanish restaurant in January 2006. They currently own http://www.sacramentomidtown.com/aioli/aioli.shtml Aïoli Bodega Española in Midtown Sacramento.

 The only thing fancy and upscale about this restaurant is the prices.  while the atmosphere is nice, the service is terible.  the waitress was cold and the waiter was rude and brash.  the food was down right bad.  the chicken was undercooked.  the menu did not list important ingredients that changed the entire flavor of the dish.  their pasta dish was bland, their assortive olives were worse than the ones I would buy at a grocery store, and the snapper tasted worse than microwaved frozen fish.  also the chicken breast and the fish had little bones in them.  the only thing good was the dressing to the side salad that came with the meal, and perhalps their mash patatoes.  this is perhalps the worst food ive had in davis (including all the places that serve food for under 5 dollars).  if you want to go to an upscale restaurant try Seasons or Mustard Seed; they are excellant.  however if you want low quality food and bad service for upscale prices try california cafe. Users/MattHh

 When I ate here, the food and service was pretty damn good.  Just too expensive though. Users/PatrickSing

 When we went there the food was indeed unimpressive, and the service was definitely lacking.  I would agree that the mashed potatoes were good, though.  Users/MattJurach

 I think this place tries to be fancy but cant back it up.  Atmosphere is not that nice, food is too expensive for how it tastes.  I havent been impressed.  Users/StephanieLiese

 I was asked 6 times (yes, that many) if I wanted wine, which I didnt, by service personnel who sneered when I didnt order any. Even the dishware is crappy old Corning Ware. Poseurs, all the way around! Users/NoelBruening

 Last night was the second time that I went there and although the food tends to be pricey, I have never been let down. The baked goat cheese starter is oh so yummy but I find it weird that you only get three crostini w/ it. However, since you also get a basket of bread, it doesnt really matter. The baslamic vinegrette that comes with the salad is amazing. Last night the soup de jour was this tex mex chowder which was so complex plus a nice spicy kick in the background. The chipotle prawn salad was good (they gave you a whole HALF of a Haas avocado!!!) however I thought the dressing would be more of a creamy aioli. Its presentation was excellent. The spicy peanut chicken was pretty good too (the sauce was so good!) and it came with some couscous and greens. However the most amazing dish that Ive tried there has to be the cowboy steak. The flavors were so rich and yummy that you didnt want it to end. I like the atmosphere b/c it feels very intimate and I like the kind of music they play.

My biggest issue with Cafe California has to be the fact that their menus are on pieces of paper, the dinnerware could be improved to go w/ the atmosphere of the restaurant, and the menu is pretty much the same. Users/KarlaFung

20051123 14:59:37 nbsp Recently went here for supper because they advertized clams. We arrived, looked at the menu  no clams, and asked the waitthing if they had clams. No, she snapped, Were sold out.  This just pisses me off that they advertize something and then dont have it. And then no apology and get treated like we are ignorant or stupid. Grrrrrrrrrrrrrrrr. Despite all that the food here is quite good  no real complaint there. An aside: Im 58 and when I go here I feel young  I wouldnt be surprised if they started carding people and refusing entry to anyone under 70. Users/SorryAboutThat

20051203 09:54:36 nbsp A recent article in the Davis Enterprise explained that Cafe California has been sold.  The new owners will be converting the space into a Spanish restaurant that will serve tapas. Users/SamanthaGrant

20051203 10:42:00 nbsp Any relation to Hotel California? Users/SteveDavison

20051206 00:47:45 nbsp We went there tonight.  The menu we were handed was four sheets of loose paper.  Dessert menu was recited without pricesa chocolate mousse, bowl of raspberries/strawberries with some sauce, or HaagenDaaz.  My lamb stew seemed to be even thirds lamb, mushrooms, and onions, in a heavy brandycream broth.  It was tasty, though.  Veggiestwo green beans and 1/2 carrot, sliced lengthwise. Users/JudithTruman

20110617 10:45:52 nbsp I hear some locals are opening a new restaurant in this location, since its predecessor Aioli has closed. Its great to see Davis welcome a new eating establishment NOT based on Thai cuisine! I do hope the new owners create a friendly hospitable atmosphere. . .I went to Aioli once and never went back after the bartender was TOTALLY rude, and speaking of Cafe California. . .I once thought I would treat my mom to an upscale dinner and took her there at the recommendation of a friend. The waitress (who I was later told was the owners sister) didnt even try to be polite to us, in fact she was a snide sourpuss who made no effort to even pretend she liked her job, in fact, she acted like we were a huge inconvenience and was really snide to my mom, which I think is despicable. Had to get that off my chest. Whenever I pay through the nose to eat at a restaurant where the staff puts on airs or cops attitude about simple questions, I am reminded that its healthier, more fun, and more rewarding not to mention easier on the budget, to cook and eat home. Users/MT
 Fair enough.  Colin OHara, bartender at this former establishment and the new present one, Our house is great. Users/BruceHansen

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
! Rice test
! check that a 'stop' statement without a stop-code is handled correctly

program stop_no_stop_code
  stop
end program

max_stars_repo_name: maurizioabba/rose

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


text: This Memorials memorial Address(benches bench, 38.564853, 121.762647) is located on the Senda Nueva Greenbelt just east of the entrance to Antelope Avenue.


text:
      PROGRAM lsl2rdf

c     ***********************************************************************
c     * Convert a trajectory file into a lat/lon/p RDF-netCDF file          *
c     * RDF = Reverse Domain Filling                                        *
c     * Michael Sprenger / Spring, summer 2010                              *
c     ***********************************************************************

      use netcdf

      implicit none

c     ----------------------------------------------------------------------
c     Declaration of parameters
c     ----------------------------------------------------------------------

c     Real comparison
      real       eps
      parameter (eps=0.001)

c     ----------------------------------------------------------------------
c     Declaration of variables
c     ----------------------------------------------------------------------

c     Input and output format for trajectories (see iotra.f)
      character*80                           inpfile     ! Input filename
      character*80                           outfile     ! Output filename

c     Trajectories
      integer                                ntra        ! Number of trajectories
      integer                                ntim        ! Number of times
      integer                                ncol        ! Number of columns
      real,allocatable, dimension (:,:,:) :: tra         ! Trajectories (ntra,ntim,ncol)
      character*80                           vars(100)   ! Variable names
      integer                                refdate(6)  ! Reference date

c     Output RDF-netCDF file
      character*80                           cdfname
      character*80                           varname
      character*80                           longname
      character*80                           unit
      integer                                nx,ny,nz
      real,allocatable, dimension (:) ::     lon,lat,lev
      real,allocatable, dimension (:,:,:) :: arr
      real                                   time
      integer                                crefile,crevar,cretime
      integer,allocatable, dimension (:,:,:) :: ind
      integer                                indx,indy,indz

c     Auxiliary variables
      integer                                inpmode
      integer                                stat
      integer                                fid
      integer                                i,j,k,ivar,itim
      real                                   tmp(10000)
      integer                                ok
      integer                                count

c     ----------------------------------------------------------------------
c     Read input data
c     ----------------------------------------------------------------------

c     Write start message
      print*,'========================================================='
      print*,'              *** START OF PROGRAM LSL2RDF ***'
      print*

c     Read parameters
      open(10,file='lsl2rdf.param')
       read(10,*) inpfile
       read(10,*) outfile
       read(10,*) ntra,ntim,ncol
      close(10)

c     Determine the formats
      call mode_tra(inpmode,inpfile)
      if (inpmode.eq.-1) inpmode=1

c     Allocate memory for input trajectory
      allocate(tra(ntra,ntim,ncol),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array tra      ***'

c     Read inpufile
      call ropen_tra(fid,inpfile,ntra,ntim,ncol,refdate,vars,inpmode)
      call read_tra (fid,tra,ntra,ntim,ncol,inpmode)
      call close_tra(fid,inpmode)

c     ----------------------------------------------------------------------
c     Determine output grid and get mapping {output array} <-> {trajectory}
c     ----------------------------------------------------------------------

c     Get a list of longitudes at time 0
      nx = 0
      do i=1,ntra
        ok = 1
        do j=1,nx
           if ( abs(tmp(j)-tra(i,1,2)).lt.eps ) then
              ok = 0
              goto 100
           endif
        enddo
 100    if ( ok.eq.1 ) then
          nx      = nx + 1
          tmp(nx) = tra(i,1,2)
        endif
      enddo

      if ( nx.gt.1 ) then
        call sort(nx,tmp)
      endif

      allocate(lon(nx),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array lon      ***'
      do i=1,nx
        lon(i) = tmp(i)
      enddo

c     Get a list of latitudes at time 0
      ny = 0
      do i=1,ntra
        ok = 1
        do j=1,ny
           if ( abs(tmp(j)-tra(i,1,3)).lt.eps ) then
              ok = 0
           endif
        enddo
 101    if ( ok.eq.1 ) then
          ny      = ny + 1
          tmp(ny) = tra(i,1,3)
        endif
      enddo

      if ( ny.gt.1 ) then
        call sort(ny,tmp)
      endif

      allocate(lat(ny),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array lat      ***'
      do i=1,ny
        lat(i) = tmp(i)
      enddo

c     Get a list of pressure levels at time 0
      nz = 0
      do i=1,ntra
        ok = 1
        do j=1,nz
           if ( abs(tmp(j)-tra(i,1,4)).lt.eps ) then
              ok = 0
              goto 102
           endif
        enddo
 102    if ( ok.eq.1 ) then
          nz      = nz + 1
          tmp(nz) = tra(i,1,4)
        endif
      enddo

      if ( nz.gt.1 ) then
         call sort(nz,tmp)
      endif

      allocate(lev(nz),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array lev      ***'
      do i=1,nz
        lev(i) = tmp(i)
      enddo

c     Write output grid
      print*,'---- OUTPUT GRID ----------------------------------------'
      print*
      write(*,'(a10,1000f5.0)') 'lon:',(lon(i),i=1,nx)
      write(*,'(a10,1000f5.0)') 'lat:',(lat(i),i=1,ny)
      write(*,'(a10,1000f5.0)') 'lev:',(lev(i),i=1,nz)
      print*

c     Decide whether several levels are really needed
      if ( (nx*ny).eq.ntra ) then
        if ( nz.ne.1 ) then
            print*
            print*,'All levels will be mapped to index level 1'
            nz = 1
        endif
      endif

c     Allocate memory for output RDF-netCDF
      allocate(arr(nx,ny,nz),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array arr      ***'
      allocate(ind(nx,ny,nz),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array ind      ***'

c     Init the mapping array <-> trajectory
      do i=1,nx
        do j=1,ny
            do k=1,nz
                ind(i,j,k) = 0
            enddo
        enddo
      enddo

c     Now get for each grid point the right trajectory index
      do i=1,ntra

         indx = 0
         do j=1,nx
            if ( abs(lon(j)-tra(i,1,2)).lt.eps  ) then
              indx = j
              goto 200
            endif
         enddo
 200     continue

         indy = 0
         do j=1,ny
            if ( abs(lat(j)-tra(i,1,3)).lt.eps  ) then
              indy = j
              goto 201
            endif
         enddo
 201     continue

         if ( nz.gt.1 ) then
           indz = 0
           do j=1,nz
            if ( abs(lev(j)-tra(i,1,4)).lt.eps  ) then
              indz = j
              goto 202
            endif
           enddo
 202       continue
         else
           indz = 1
         endif

         if ( (indx.ne.0).and.(indy.ne.0).and.(indz.ne.0) ) then
            ind(indx,indy,indz) = i
         endif

      enddo

c     Check whether all grid points are linked to a traejctory
      count = 0
      do i=1,nx
        do j=1,ny
            do k=1,nz
                if ( ind(i,j,k).eq.0 ) count = count + 1
            enddo
        enddo
      enddo
      if ( count.ne.0 ) then
        print*,' WARNING: not all grid points in RDF domain linked'
        print*,'          to a trajectory... proceed anyway.'
      endif

c     ----------------------------------------------------------------------
c     Fill the output array
c     ----------------------------------------------------------------------

      print*,'---- MAPPING TO RDF GRID --------------------------------'
      print*

c     Loop over all variables
      crefile = 1
      do ivar=2,ncol

c        Loop over all times
         cretime = 1
         crevar  = 1
         do itim=1,ntim

c            Get the time from the first trajectory
             time = tra(1,itim,1)
             write(*,'(a20,a10,f8.2)')
     >            ' variable,time : ',trim( vars(ivar) ),time

c            Do the remapping
             do i=1,nx
                do j=1,ny
                    do k=1,nz
                       if ( ind(i,j,k).ne.0 ) then
                        arr(i,j,k) = tra( ind(i,j,k), itim, ivar )
                       endif
                    enddo
                enddo
             enddo

c            Save the array to RDF-netCDF
             call writecdf2D_cf
     >          (outfile,vars(ivar),nx,lon,ny,lat,nz,lev,
     >           arr,time,crefile,crevar,cretime)

             crefile = 0
             crevar  = 0

         enddo

      enddo

c     Write end message
      print*
      print*,'              *** END OF PROGRAM LSL2RDF ***'
      print*,'========================================================='

      end

c     ********************************************************************
c     * SUBROUTINES                                                      *
c     ********************************************************************

c     --------------------------------------------------------------------
c     Subroutines to write the CF netcdf output file
c     --------------------------------------------------------------------

      subroutine writecdf2D_cf
     >          (cdfname,varname,nx,lon,ny,lat,nz,lev,
     <           arr,time,crefile,crevar,cretime)

      USE netcdf

      IMPLICIT NONE

c     Declaration of input parameters
      character*80 cdfname
      character*80 varname,longname,unit
      integer      nx,ny,nz
      real         lon(nx),lat(ny),lev(nz)
      real         arr(nx,ny,nz)
      real         time
      integer      crefile,crevar,cretime

c     Local variables
      integer      ierr
      integer      ncID
      integer      LonDimId,    varLonID
      integer      LatDimID,    varLatID
      integer      LevDimID,    varLevID
      integer      TimeDimID,   varTimeID
      real         timeindex
      integer      i
      integer      nvars,varids(100)
      integer      ndims,dimids(100)
      real         timelist(1000)
      integer      ntimes
      integer      ind
      integer      varID
      real         xmin,xmax,ymin,ymax,zmin,zmax

c     Quick an dirty solution for fieldname conflict
      if ( varname.eq.'time' ) varname = 'TIME'

c     Initially set error to indicate no errors.
      ierr = 0

c     ---- Create the netCDF - skip if <crefile=0> ----------------------
      if ( crefile.ne.1 ) goto 100

c     Create the file
      ierr = nf90_create(trim(cdfname), NF90_CLOBBER, ncID)

c     Define dimensions
      ierr=nf90_def_dim(ncID,'longitude' ,nx            , LonDimID )
      ierr=nf90_def_dim(ncID,'latitude'  ,ny            , LatDimID )
      ierr=nf90_def_dim(ncID,'level    ' ,nz            , LevDimID )
      ierr=nf90_def_dim(ncID,'time'      ,nf90_unlimited, TimeDimID)

c     Define coordinate Variables
      ierr = nf90_def_var(ncID,'longitude',NF90_FLOAT,
     >     (/ LonDimID /),varLonID)
      ierr = nf90_put_att(ncID, varLonID, "standard_name","longitude")
      ierr = nf90_put_att(ncID, varLonID, "units"      ,"degree_east")

      ierr = nf90_def_var(ncID,'latitude',NF90_FLOAT,
     >     (/ LatDimID /),varLatID)
      ierr = nf90_put_att(ncID, varLatID, "standard_name", "latitude")
      ierr = nf90_put_att(ncID, varLatID, "units"    ,"degree_north")

      ierr = nf90_def_var(ncID,'level',NF90_FLOAT,
     >     (/ LevDimID /),varLevID)
      ierr = nf90_put_att(ncID, varLevID, "standard_name", "level")
      ierr = nf90_put_att(ncID, varLevID, "units"    ,"nil")

      ierr = nf90_def_var(ncID,'time',NF90_FLOAT,
     >     (/ TimeDimID /), varTimeID)

c     Write general global attributes
      ierr = nf90_put_att(ncID, NF90_GLOBAL, 'Conventions', 'CF-1.0')
      ierr = nf90_put_att(ncID, NF90_GLOBAL, 'title',
     >     'RDF trajectory')
      ierr = nf90_put_att(ncID, NF90_GLOBAL, 'source',
     >     'Lagranto Trajectories')
      ierr = nf90_put_att(ncID, NF90_GLOBAL, 'institution',
     >     'ETH Zurich, IACETH')

c     Write grid extension
      xmin = lon(1)
      xmax = lon(1)
      do i=2,nx
        if ( lon(i).lt.xmin ) xmin = lon(i)
        if ( lon(i).gt.xmax ) xmax = lon(i)
      enddo

      ymin = lat(1)
      ymax = lat(1)
      do i=2,ny
        if ( lat(i).lt.ymin ) ymin = lat(i)
        if ( lat(i).gt.ymax ) ymax = lat(i)
      enddo

      zmin = lev(1)
      zmax = lev(1)
      do i=2,nz
        if ( lev(i).lt.zmin ) zmin = lev(i)
        if ( lev(i).gt.zmax ) zmax = lev(i)
      enddo

      ierr = nf90_put_att(ncID, NF90_GLOBAL, 'nx',nx )
      ierr = nf90_put_att(ncID, NF90_GLOBAL, 'ny',ny )
      ierr = nf90_put_att(ncID, NF90_GLOBAL, 'nz',nz )
      ierr = nf90_put_att(ncID, NF90_GLOBAL, 'xmin',xmin )
      ierr = nf90_put_att(ncID, NF90_GLOBAL, 'xmax',xmax )
      ierr = nf90_put_att(ncID, NF90_GLOBAL, 'ymin',ymin )
      ierr = nf90_put_att(ncID, NF90_GLOBAL, 'ymax',ymax )
      ierr = nf90_put_att(ncID, NF90_GLOBAL, 'zmin',zmin )
      ierr = nf90_put_att(ncID, NF90_GLOBAL, 'zmax',zmax )

c     Check whether the definition was successful
      ierr = nf90_enddef(ncID)
      if (ierr.gt.0) then
         print*, 'An error occurred while attempting to ',
     >        'finish definition mode.'
         stop
      endif

c     Write coordinate data
      ierr = nf90_put_var(ncID,varLonID ,lon )
      ierr = nf90_put_var(ncID,varLatID ,lat )
      ierr = nf90_put_var(ncID,varLevID ,lev )

c     Close netCDF file
      ierr = nf90_close(ncID)

 100  continue

c     ---- Define a new variable - skip if <crevar=0> -----------------------

      if ( crevar.ne.1 ) goto 110

c     Open the file for read(write access
      ierr = nf90_open  (trim(cdfname), NF90_WRITE  , ncID)

c     Get the IDs for dimensions
      ierr = nf90_inq_dimid(ncID,'longitude'     , LonDimID )
      ierr = nf90_inq_dimid(ncID,'latitude'      , LatDimID )
      ierr = nf90_inq_dimid(ncID,'level'         , LevDimID )
      ierr = nf90_inq_dimid(ncID,'time'          , TimeDimID)

c     Enter define mode
      ierr = nf90_redef(ncID)

c     Write definition and add attributes
      ierr = nf90_def_var(ncID,varname,NF90_FLOAT,
     >             (/ LonDimID, LatDimID, LevDimID, TimeDimID /),varID)
      ierr = nf90_put_att(ncID, varID, '_FillValue', -999.99  )

c     Check whether definition was successful
      ierr = nf90_enddef(ncID)
      if (ierr.gt.0) then
         print*, 'An error occurred while attempting to ',
     >           'finish definition mode.'
         stop
      endif

c     Close netCDF file
      ierr = nf90_close(ncID)

 110  continue

c     ---- Create a new time (unlimited dimension) - skip if <cretime=0> ------

      if ( cretime.ne.1 ) goto 120

c     Open the file for read/write access
      ierr = nf90_open  (trim(cdfname), NF90_WRITE, ncID)

c     Get the list of times on the netCDF file
      ierr = nf90_inq_dimid(ncID,'time', TimeDimID)
      if ( ierr.ne.0 ) then
         print*,'Time dimension is not defined on ',trim(cdfname),
     >          ' .... Stop'
         stop
      endif
      ierr = nf90_inquire_dimension(ncID, TimeDimID, len = ntimes)
      ierr = nf90_inq_varid(ncID,'time', varTimeID)
      if ( ierr.ne.0 ) then
         print*,'Variable time is not defined on ',trim(cdfname),
     >          ' ... Stop'
         stop
      endif
      ierr = nf90_get_var(ncID,varTimeID,timelist(1:ntimes))

c     Decide whether a new time must be written
      ind = 0
      do i=1,ntimes
         if ( time.eq.timelist(i) ) ind = i
      enddo

c     Extend the time list if required
      if ( ind.eq.0 ) then
         ntimes           = ntimes + 1
         timelist(ntimes) = time
         ierr = nf90_put_var(ncID,varTimeID,timelist(1:ntimes))
      endif

c     Close netCDF file
      ierr = nf90_close(ncID)

 120  continue

c     ---- Write data --------------------------------------------------

c     Open the file for read/write access
      ierr = nf90_open  (trim(cdfname), NF90_WRITE , ncID)

c     Get the varID
      ierr = nf90_inq_varid(ncID,varname, varID )
      if (ierr.ne.0) then
         print*,'Variable ',trim(varname),' is not defined on ',
     >          trim(cdfname)
         stop
      endif

c     Get the time index
      ierr = nf90_inq_dimid(ncID,'time', TimeDimID)
      if ( ierr.ne.0 ) then
         print*,'Time dimension is not defined on ',trim(cdfname),
     >          ' .... Stop'
         stop
      endif
      ierr = nf90_inquire_dimension(ncID, TimeDimID, len = ntimes)
      ierr = nf90_inq_varid(ncID,'time', varTimeID)
      if ( ierr.ne.0 ) then
         print*,'Variable time is not defined on ',trim(cdfname),
     >          ' ... Stop'
         stop
      endif
      ierr = nf90_get_var(ncID,varTimeID,timelist(1:ntimes))
      ind = 0
      do i=1,ntimes
         if ( time.eq.timelist(i) ) ind = i
      enddo
      if (ind.eq.0) then
         print*,'Time',time,' is not defined on the netCDF file',
     >          trim(cdfname),' ... Stop'
         stop
      endif

c     Write data block
      ierr = nf90_put_var(ncID,varID,arr,
     >                    start = (/  1,  1,  1, ind /),
     >                    count = (/ nx, ny, nz,   1 /) )

c     Check whether writing was successful
      ierr = nf90_close(ncID)
      if (ierr.ne.0) then
         write(*,*) trim(nf90_strerror(ierr))
         write(*,*) 'An error occurred while attempting to ',
     >              'close the netcdf file.'
         write(*,*) 'in clscdf_CF'
      endif

      end

c     --------------------------------------------------------------------
c     Sort an array (from Muerical Recipes)
c     --------------------------------------------------------------------

      SUBROUTINE SORT(N,RA)
      DIMENSION RA(N)
      L=N/2+1
      IR=N
10    CONTINUE
        IF(L.GT.1)THEN
          L=L-1
          RRA=RA(L)
        ELSE
          RRA=RA(IR)
          RA(IR)=RA(1)
          IR=IR-1
          IF(IR.EQ.1)THEN
            RA(1)=RRA
            RETURN
          ENDIF
        ENDIF
        I=L
        J=L+L
20      IF(J.LE.IR)THEN
          IF(J.LT.IR)THEN
            IF(RA(J).LT.RA(J+1))J=J+1
          ENDIF
          IF(RRA.LT.RA(J))THEN
            RA(I)=RA(J)
            I=J
            J=J+J
          ELSE
            J=IR+1
          ENDIF
        GO TO 20
        ENDIF
        RA(I)=RRA
      GO TO 10
      END

max_stars_repo_name: GEOS-ESM/AeroApps

max_issues_repo_name: GEOS-ESM/AeroApps

max_forks_repo_name: GEOS-ESM/AeroApps


text:
      subroutine epflag(matype,ipfl,ipelem,lpb,numel)
c     implicit double precision (a-h,o-z)                                    dp
c
c     set element print flags
c
      dimension ipelem(2,*),matype(*),ipfl(*)
      do 20 i=1,numel
      do 10 j=1,lpb
      k=ipelem(1,j)
      l=ipelem(2,j)
      if (i.ge.k.and.i.le.l) go to 20
   10 continue
      ipfl(i)=1
   20 continue
c
      return
c
      end

max_stars_repo_name: imohame/LabCode

max_issues_repo_name: imohame/LabCode

max_forks_repo_name: imohame/LabCode


subroutine initialise(z, y, x, state_a_m1)
    implicit none

    ! spaceshape
    integer, intent(in) :: z, y, x
    ! component states
    real(kind=8), intent(inout), dimension(z, y, x) :: state_a_m1

    state_a_m1 = 0

end subroutine initialise

subroutine run(z, y, x, &
               transfer_j, transfer_m, &
               ancillary_b, &
               parameter_c, &
               state_a_m1, state_a_0, &
               constant_c, &
               transfer_l, transfer_n, transfer_o, &
               output_x, output_y)

    implicit none

    ! spaceshape
    integer, intent(in) :: z, y, x
    ! from exchanger
    real(kind=8), intent(in), dimension(z, y, x) :: transfer_j, transfer_m
    ! component ancillary data
    real(kind=8), intent(in), dimension(12, z, y, x) :: ancillary_b
    ! component parameters
    real(kind=8), intent(in) :: parameter_c
    ! component states
    real(kind=8), intent(in), dimension(z, y, x) :: state_a_m1
    real(kind=8), intent(inout), dimension(z, y, x) :: state_a_0
    ! component constants
    real(kind=8), intent(in) :: constant_c
    ! to exchanger
    real(kind=8), intent(out), dimension(z, y, x) :: &
        transfer_l, transfer_n, transfer_o
    ! component outputs
    real(kind=8), intent(out), dimension(z, y, x) :: &
        output_x, output_y

    state_a_0 = state_a_m1 + 1

    transfer_l = (ancillary_b(12,:,:,:) * transfer_m) + state_a_0
    transfer_n = parameter_c * transfer_j
    transfer_o = parameter_c + transfer_j

    output_x = (parameter_c * transfer_j) + constant_c
    output_y = (ancillary_b(12,:,:,:) * transfer_m) - state_a_0

end subroutine run

subroutine finalise()
    implicit none

end subroutine finalise


!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
!! 3D oil reservoir Simulator project
!! by Mehrdad Gharib Shirangi From Spring 2010 until Fall 2010
!! mehrdad.ghsh@gmail.com
!! Copy Right by Mehrdad Gharib Shirangi
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

Module Derivative

! contains transmissibility functions
use globalvar
use Trans

implicit none

!real :: Press

contains
   real*8 function DTrxa_DP(i,j,k,m,wrt) ! the value of the derivative of x transmissibility between a gridblock and a next grid
   integer :: ind , N_ind ! Neighbor grid index
   integer :: i,j,k
   real*8 ::  delx_avg , Press
   real*8 :: denomenator, press_term, Dip_cor
   Character*1 :: m ! 'o' or 'w' : oil or water
   Character*1 :: wrt ! 'a' , 'c'  ! wrt to after , or before or itself

      ind = (k-1) * Nx * Ny + (j-1) * Nx + i
      N_ind = ind + 1
      delx_avg = (Deltax(i+1) + Deltax(i)) / 2.0
      denomenator = 1 / 2.0 * ( Deltax(i+1) / ( kx(ind+1) * Deltay(j) * Deltaz(ind+1) ) + Deltax(i)/( kx(ind)*Deltay(j)*Deltaz(ind) ) )
      Dip_cor = delx_avg ** 2 / ((Zblock(ind)-Zblock(N_ind)) ** 2 + delx_avg ** 2)

      if ( D_Poten( ind , N_ind , m ) > 0 ) then !( Potential(ind,m) > Potential(N_ind,m) ) then ! upstream block
         if (wrt=='a') then
            press_term = 0.0
         Else
            if ( m == 'w') then
			    press_term = - k_rw(ind) / ( mu_w (ind) * B_w(ind ) ) * ( B_der(ind ,'w') / B_w(ind ) + mu_der(ind ,'w') / mu_w(ind) )
            Else
				press_term = - k_ro(ind) / ( mu_o (ind) * B_o(ind ) ) * ( B_der(ind ,'o') / B_o(ind ) + mu_der(ind , 'o') / mu_o(ind) )
            EndIf
         EndIf
      else
         if (wrt=='c') then
            press_term = 0.0
         Else
             if ( m == 'w') then
			    press_term = - k_rw( N_ind ) / ( mu_w (N_ind) * B_w(N_ind ) ) * ( B_der(N_ind ,'w') / B_w(N_ind ) + mu_der(N_ind ,'w') / mu_w(N_ind) )
             Else
				press_term = - k_ro( N_ind) / ( mu_o ( N_ind) * B_o( N_ind ) ) * ( B_der( N_ind ,'o') / B_o( N_ind ) + mu_der( N_ind , 'o') / mu_o( N_ind) )
             EndIf
         Endif
      EndIf

      DTrxa_DP = betac * Dip_cor * press_term /denomenator

   Endfunction DTrxa_DP
!*******************************************************************************************
   real*8 function DTrxb_DP(i,j,k,m,wrt) ! the value of the derivative of x transmissibility between a gridblock and a next grid
   integer :: ind , N_ind ! Neighbor grid index
   integer :: i,j,k
   real*8 ::  delx_avg , Press
   real*8 :: denomenator, press_term, Dip_cor
   Character*1 :: m ! 'o' or 'w' : oil or water
   Character*1 :: wrt ! 'b' , 'c'  ! wrt to after , or before or itself

      ind = (k-1) * Nx * Ny + (j-1) * Nx + i
      N_ind = ind - 1
      delx_avg = (Deltax(i-1) + Deltax(i))/2.0
      denomenator = 1 / 2.0 * ( Deltax(i-1) / ( kx(ind-1) * Deltay(j) * Deltaz(ind-1) ) + Deltax(i)/( kx(ind)*Deltay(j)*Deltaz(ind) ) )
      Dip_cor = delx_avg ** 2 / ((Zblock(ind)-Zblock(N_ind)) ** 2 + delx_avg ** 2)

      if ( D_Poten( ind , N_ind , m ) > 0 ) then ! upstream block
         if (wrt=='b') then
            press_term = 0.0
         Else
            if ( m == 'w') then
			    press_term = - k_rw(ind) / ( mu_w (ind) * B_w(ind ) ) * ( B_der(ind ,'w') / B_w(ind ) )
            Else
				press_term = - k_ro(ind) / ( mu_o (ind) * B_o(ind ) ) * ( B_der(ind ,'o') / B_o(ind ) + mu_der(ind , 'o') / mu_o(ind) )
            EndIf
         EndIf
      else
         if (wrt=='c') then
            press_term = 0.0
         Else
             if ( m == 'w') then
			    press_term = - k_rw( N_ind ) / ( mu_w (N_ind) * B_w(N_ind ) ) * ( B_der(N_ind ,'w') / B_w(N_ind ) )
             Else
				press_term = - k_ro( N_ind) / ( mu_o ( N_ind) * B_o( N_ind ) ) * ( B_der( N_ind ,'o') / B_o( N_ind ) + mu_der( N_ind , 'o') / mu_o( N_ind) )
             EndIf
         Endif
      EndIf

      DTrxb_DP = betac * Dip_cor * press_term /denomenator

   Endfunction DTrxb_DP
!*******************************************************************************************
   real*8 function DTryb_DP(i,j,k,m,wrt) ! the value of the derivative of x transmissibility between a gridblock and a next grid
   integer :: ind , N_ind ! Neighbor grid index
   integer :: i,j,k
   real*8 :: dely_avg , Press
   real*8 :: denomenator, press_term, Dip_cor
   Character*1 :: m ! 'o' or 'w' : oil or water
   Character*1 :: wrt ! 'b' , 'c'  ! wrt to after , or before or itself

      ind = (k-1) * Nx * Ny + (j-1) * Nx + i
      N_ind = ind - Nx
      dely_avg = (Deltay(j-1) + Deltay(j))/2.0
      denomenator = 1 / 2.0 * ( Deltay(j-1) / ( ky(N_ind) * Deltax(i) * Deltaz(N_ind) ) + Deltay(j)/( ky(ind) * Deltax(i) * Deltaz(ind) ) )
      Dip_cor = dely_avg ** 2 / ((Zblock(ind)-Zblock(N_ind)) ** 2 + dely_avg ** 2)

      if ( D_Poten( ind , N_ind , m ) > 0 ) then ! upstream block
         if (wrt=='b') then
            press_term = 0.0
         Else
            if ( m == 'w') then
			    press_term = - k_rw(ind) / ( mu_w (ind) * B_w(ind ) ) * ( B_der(ind ,'w') / B_w(ind ) )
            Else
				press_term = - k_ro(ind) / ( mu_o (ind) * B_o(ind ) ) * ( B_der(ind ,'o') / B_o(ind ) + mu_der(ind , 'o') / mu_o(ind) )
            EndIf
         EndIf
      else
         if (wrt=='c') then
            press_term = 0.0
         Else
             if ( m == 'w') then
			    press_term = - k_rw( N_ind ) / ( mu_w (N_ind) * B_w(N_ind ) ) * ( B_der(N_ind ,'w') / B_w(N_ind ) )
             Else
				press_term = - k_ro( N_ind) / ( mu_o ( N_ind) * B_o( N_ind ) ) * ( B_der( N_ind ,'o') / B_o( N_ind ) + mu_der( N_ind , 'o') / mu_o( N_ind) )
             EndIf
         Endif
      EndIf

      DTryb_DP = betac * Dip_cor * press_term /denomenator

   Endfunction DTryb_DP
!*******************************************************************************************
   real*8 function DTrya_DP(i,j,k,m,wrt) ! the value of the derivative of x transmissibility between a gridblock and a next grid
   integer :: ind , N_ind ! Neighbor grid index
   integer :: i,j,k
   real*8 :: dely_avg, Press
   real*8 :: denomenator, press_term, Dip_cor
   Character*1 :: m ! 'o' or 'w' : oil or water
   Character*1 :: wrt ! 'a' , 'c'  ! wrt to after , or before or itself

      ind = (k-1) * Nx * Ny + (j-1) * Nx + i
      N_ind = ind + Nx
      dely_avg = (Deltay(j+1) + Deltay(j))/2.0
      denomenator = 1 / 2.0 * ( Deltay(j+1) / ( ky(N_ind) * Deltax(i) * Deltaz(N_ind) ) + Deltay(j)/( ky(ind) * Deltax(i) * Deltaz(ind) ) )
      Dip_cor = dely_avg ** 2 / ((Zblock(ind)-Zblock(N_ind)) ** 2 + dely_avg ** 2)

      if ( D_Poten( ind , N_ind , m ) > 0.0 ) then ! upstream block
         if (wrt=='a') then
            press_term = 0.0
         Else
            if ( m == 'w') then
			    press_term = - k_rw(ind) / ( mu_w (ind) * B_w(ind ) ) * ( B_der(ind ,'w') / B_w(ind ) )
            Else
				press_term = - k_ro(ind) / ( mu_o (ind) * B_o(ind ) ) * ( B_der(ind ,'o') / B_o(ind ) + mu_der(ind , 'o') / mu_o(ind) )
            EndIf
         EndIf
      else
         if (wrt=='c') then
            press_term = 0.0
         Else
             if ( m == 'w') then
			    press_term = - k_rw( N_ind ) / ( mu_w (N_ind) * B_w(N_ind ) ) * ( B_der(N_ind ,'w') / B_w(N_ind ) )
             Else
				press_term = - k_ro( N_ind) / ( mu_o ( N_ind) * B_o( N_ind ) ) * ( B_der( N_ind ,'o') / B_o( N_ind ) + mu_der( N_ind , 'o') / mu_o( N_ind) )
             EndIf
         Endif
      EndIf

      DTrya_DP = betac * Dip_cor * press_term /denomenator

   Endfunction DTrya_DP

!*******************************************************************************************

   real*8 function DTrzb_DP(i,j,k,m,wrt) ! the value of the derivative of x transmissibility between a gridblock and a next grid
   integer :: ind , N_ind ! Neighbor grid index
   integer :: i,j,k
   real*8 :: denomenator, press_term , Press
   Character*1 :: m ! 'o' or 'w' : oil or water
   Character*1 :: wrt ! 'a' , 'c'  ! wrt to after , or before or itself

      ind = (k-1) * Nx * Ny + (j-1) * Nx + i
      N_ind = ind - Nx * Ny
      denomenator = 1 / 2.0 * ( Deltaz(N_ind) / ( kz(N_ind) * Deltax(i) * Deltay(j) ) + Deltaz(ind)/( kz(ind) * Deltax(i) * Deltay(j) ) )

      if ( D_Poten( ind , N_ind , m ) > 0 ) then ! upstream block
         if (wrt=='b') then
            press_term = 0.0
         Else
            if ( m == 'w') then
			    press_term = - k_rw(ind) / ( mu_w (ind) * B_w(ind ) ) * ( B_der(ind ,'w') / B_w(ind ) )
            Else
				press_term = - k_ro(ind) / ( mu_o (ind) * B_o(ind ) ) * ( B_der(ind ,'o') / B_o(ind ) + mu_der(ind , 'o') / mu_o(ind) )
            EndIf
         EndIf
      else
         if (wrt=='c') then
            press_term = 0.0
         Else
             if ( m == 'w') then
			    press_term = - k_rw( N_ind ) / ( mu_w (N_ind) * B_w(N_ind ) ) * ( B_der(N_ind ,'w') / B_w(N_ind ) )
             Else
				press_term = - k_ro( N_ind) / ( mu_o ( N_ind) * B_o( N_ind ) ) * ( B_der( N_ind ,'o') / B_o( N_ind ) + mu_der( N_ind , 'o') / mu_o( N_ind) )
             EndIf
         Endif
      EndIf

      DTrzb_DP = betac * press_term /denomenator

   Endfunction DTrzb_DP

!*******************************************************************************************

   real*8 function DTrza_DP(i,j,k,m,wrt) ! the value of the derivative of x transmissibility between a gridblock and a next grid
   integer :: ind , N_ind ! Neighbor grid index
   integer :: i,j,k
   real*8 ::delz_avg , Press
   real*8 :: denomenator, press_term
   Character*1 :: m ! 'o' or 'w' : oil or water
   Character*1 :: wrt ! 'a' , 'c'  ! wrt to after , or before or itself

      ind = (k-1) * Nx * Ny + (j-1) * Nx + i
      N_ind = ind + Nx * Ny
      denomenator = 1 / 2.0 * ( Deltaz(N_ind) / ( kz(N_ind) * Deltax(i) * Deltay(j) ) + Deltaz(ind)/( kz(ind) * Deltax(i) * Deltay(j) ) )

      if ( D_Poten( ind , N_ind , m ) > 0 ) then ! upstream block
         if (wrt=='a') then
            press_term = 0.0
         Else
            if ( m == 'w') then
			    press_term = - k_rw(ind) / ( mu_w (ind) * B_w(ind ) ) * ( B_der(ind ,'w') / B_w(ind ) )
            Else
				press_term = - k_ro(ind) / ( mu_o (ind) * B_o(ind ) ) * ( B_der(ind ,'o') / B_o(ind ) + mu_der(ind , 'o') / mu_o(ind) )
            EndIf
         EndIf
      else
         if (wrt=='c') then
            press_term = 0.0
         Else
             if ( m == 'w') then
			    press_term = - k_rw( N_ind ) / ( mu_w (N_ind) * B_w(N_ind ) ) * ( B_der(N_ind ,'w') / B_w(N_ind ) )
             Else
				press_term = - k_ro( N_ind) / ( mu_o ( N_ind) * B_o( N_ind ) ) * ( B_der( N_ind ,'o') / B_o( N_ind ) + mu_der( N_ind , 'o') / mu_o( N_ind) )
             EndIf
         Endif
      EndIf

      DTrza_DP = betac * press_term /denomenator

   Endfunction DTrza_DP

!*******************************************************************************************
     real*8 function DTrxa_DS(i,j,k,m,wrt) ! the value of the derivative of x transmissibility between a gridblock and a next grid
   integer :: ind , N_ind ! Neighbor grid index
   integer :: i,j,k
   real*8 :: delx_avg , Press
   real*8 :: denomenator, press_term, Dip_cor
   Character*1 :: m ! 'o' or 'w' : oil or water
   Character*1 :: wrt ! 'a' , 'c'  ! wrt to after , or before or itself

      ind = (k-1) * Nx * Ny + (j-1) * Nx + i
      N_ind = ind + 1
      delx_avg = (Deltax(i+1) + Deltax(i))/2.0
      denomenator = 1 / 2.0 * ( Deltax(i+1) / ( kx(ind+1) * Deltay(j) * Deltaz(ind+1) ) + Deltax(i)/( kx(ind)*Deltay(j)*Deltaz(ind) ) )
      Dip_cor = delx_avg ** 2 / ((Zblock(ind)-Zblock(N_ind)) ** 2 + delx_avg ** 2)

      if ( D_Poten( ind , N_ind , m ) > 0 ) then ! upstream block
         if (wrt=='a') then
            press_term = 0
         Else
             if ( m == 'w') then
				press_term = DkrwdS( ind ) / ( mu_w ( ind ) * B_w( ind ) ) + k_rw(ind) / ( mu_w ( ind) * B_w(ind) ) * ( B_der(ind ,'w') / B_w(ind ) + mu_der(ind ,'w') / mu_w(ind) ) * PCO_der(ind)
             Else
				press_term = DkrodS( ind ) / ( mu_o ( ind ) * B_o( ind ) )
             EndIf
         EndIf
      else
         if (wrt=='c') then
            press_term = 0
         Else
            if ( m == 'w') then
				press_term = DkrwdS( N_ind ) / ( mu_w ( N_ind ) * B_w( N_ind ) ) + k_rw(N_ind) / ( mu_w ( N_ind) * B_w(N_ind) ) * ( B_der(N_ind ,'w') / B_w(N_ind ) + mu_der(N_ind ,'w') / mu_w(N_ind) ) * PCO_der(N_ind)
            Else
				press_term = DkrodS( N_ind ) / ( mu_o ( N_ind ) * B_o( N_ind ) )
            EndIf
         Endif
      EndIf

      DTrxa_DS = betac * Dip_cor * press_term /denomenator

   Endfunction DTrxa_DS
!*******************************************************************************************
       real*8 function DTrxb_DS(i,j,k,m,wrt) ! the value of the derivative of x transmissibility between a gridblock and a next grid
   integer :: ind , N_ind ! Neighbor grid index
   integer :: i,j,k
   real*8 :: delx_avg , Press
   real*8 :: denomenator, press_term, Dip_cor
   Character*1 :: m ! 'o' or 'w' : oil or water
   Character*1 :: wrt ! 'b' , 'c'  ! wrt to after , or before or itself

      ind = (k-1) * Nx * Ny + (j-1) * Nx + i
      N_ind = ind - 1
      delx_avg = (Deltax(i-1) + Deltax(i))/2.0
      denomenator = 1 / 2.0 * ( Deltax(i-1) / ( kx(N_ind) * Deltay(j) * Deltaz(N_ind) ) + Deltax(i)/( kx(ind)*Deltay(j)*Deltaz(ind) ) )
      Dip_cor = delx_avg ** 2 / ((Zblock(ind)-Zblock(N_ind)) ** 2 + delx_avg ** 2)

      if ( D_Poten( ind , N_ind , m ) > 0 ) then ! upstream block
         if (wrt=='b') then
            press_term = 0
         Else
             if ( m == 'w') then
				press_term = DkrwdS( ind ) / ( mu_w ( ind ) * B_w( ind ) ) + k_rw(ind) / ( mu_w ( ind) * B_w(ind) ) * ( B_der(ind ,'w') / B_w(ind ) ) * PCO_der(ind)
             Else
				press_term = DkrodS( ind ) / ( mu_o ( ind ) * B_o( ind ) )
             EndIf
         EndIf
      else
         if (wrt=='c') then
            press_term = 0
         Else
            if ( m == 'w') then
				press_term = DkrwdS( N_ind ) / ( mu_w ( N_ind ) * B_w( N_ind ) ) + k_rw(N_ind) / ( mu_w ( N_ind) * B_w(N_ind) ) * ( B_der(N_ind ,'w') / B_w(N_ind ) ) * PCO_der(N_ind)
            Else
				press_term = DkrodS( N_ind ) / ( mu_o ( N_ind ) * B_o( N_ind ) )
            EndIf
         Endif
      EndIf

      DTrxb_DS = betac * Dip_cor * press_term /denomenator

   Endfunction DTrxb_DS
!*******************************************************************************************

   real*8 function DTryb_DS(i,j,k,m,wrt) ! the value of the derivative of x transmissibility between a gridblock and a next grid
   integer :: ind , N_ind ! Neighbor grid index
   integer :: i,j,k
   real*8 ::dely_avg , Press
   real*8 :: denomenator, press_term, Dip_cor
   Character*1 :: m ! 'o' or 'w' : oil or water
   Character*1 :: wrt ! 'b' , 'c'  ! wrt to after , or before or itself

      ind = (k-1) * Nx * Ny + (j-1) * Nx + i
      N_ind = ind - Nx
      dely_avg = (Deltay(j-1) + Deltay(j))/2.0
      denomenator = 1 / 2.0 * ( Deltay(j-1) / ( ky(N_ind) * Deltax(i) * Deltaz(N_ind) ) + Deltay(j)/( ky(ind) * Deltax(i) * Deltaz(ind) ) )
      Dip_cor = dely_avg ** 2 / ((Zblock(ind)-Zblock(N_ind)) ** 2 + dely_avg ** 2)

      if ( D_Poten( ind , N_ind , m ) > 0 ) then ! upstream block
         if (wrt=='b') then
            press_term = 0
         Else
             if ( m == 'w') then
				press_term = DkrwdS( ind ) / ( mu_w ( ind ) * B_w( ind ) )+ k_rw(ind) / ( mu_w ( ind) * B_w(ind) ) * ( B_der(ind ,'w') / B_w(ind ) ) * PCO_der(ind)
             Else
				press_term = DkrodS( ind ) / ( mu_o ( ind ) * B_o( ind ) )
             EndIf
         EndIf
      else
         if (wrt=='c') then
            press_term = 0
         Else
            if ( m == 'w') then
				press_term = DkrwdS( N_ind ) / ( mu_w ( N_ind ) * B_w( N_ind ) ) + k_rw(N_ind) / ( mu_w ( N_ind) * B_w(N_ind) ) * ( B_der(N_ind ,'w') / B_w(N_ind ) ) * PCO_der(N_ind)
            Else
				press_term = DkrodS( N_ind ) / ( mu_o ( N_ind ) * B_o( N_ind ) )
            EndIf
         Endif
      EndIf

      DTryb_DS = betac * Dip_cor * press_term /denomenator

   Endfunction DTryb_DS

!*******************************************************************************************

   real*8 function DTrya_DS(i,j,k,m,wrt) ! the value of the derivative of x transmissibility between a gridblock and a next grid
   integer :: ind , N_ind ! Neighbor grid index
   integer :: i,j,k
   real*8 :: dely_avg , Press
   real*8 :: denomenator, press_term, Dip_cor
   Character*1 :: m ! 'o' or 'w' : oil or water
   Character*1 :: wrt ! 'a' , 'c'  ! wrt to after , or before or itself

      ind = (k-1) * Nx * Ny + (j-1) * Nx + i
      N_ind = ind + Nx
      dely_avg = (Deltay(j+1) + Deltay(j) )/2.0
      denomenator = 1 / 2.0 * ( Deltay(j+1) / ( ky(N_ind) * Deltax(i) * Deltaz(N_ind) ) + Deltay(j)/( ky(ind) * Deltax(i) * Deltaz(ind) ) )
      Dip_cor = dely_avg ** 2 / ((Zblock(ind)-Zblock(N_ind)) ** 2 + dely_avg ** 2)

      if ( D_Poten( ind , N_ind , m ) > 0 ) then ! upstream block
         if (wrt=='a') then
            press_term = 0
         Else
             if ( m == 'w') then
				press_term = DkrwdS( ind ) / ( mu_w ( ind ) * B_w( ind ) ) + k_rw(ind) / ( mu_w ( ind) * B_w(ind) ) * ( B_der(ind ,'w') / B_w(ind ) ) * PCO_der(ind)
             Else
				press_term = DkrodS( ind ) / ( mu_o ( ind ) * B_o( ind ) )
             EndIf
         EndIf
      else
         if (wrt=='c') then
            press_term = 0
         Else
            if ( m == 'w') then
				press_term = DkrwdS( N_ind ) / ( mu_w ( N_ind ) * B_w( N_ind ) ) + k_rw(N_ind) / ( mu_w ( N_ind) * B_w(N_ind) ) * ( B_der(N_ind ,'w') / B_w(N_ind ) ) * PCO_der(N_ind)
            Else
				press_term = DkrodS( N_ind ) / ( mu_o ( N_ind ) * B_o( N_ind ) )
            EndIf
         Endif
      EndIf

      DTrya_DS = betac * Dip_cor * press_term /denomenator

   Endfunction DTrya_DS
!*******************************************************************************************

   real*8 function DTrzb_DS(i,j,k,m,wrt) ! the value of the derivative of x transmissibility between a gridblock and a next grid
   integer :: ind , N_ind ! Neighbor grid index
   integer :: i,j,k
   real*8 :: denomenator, press_term , Press
   Character*1 :: m ! 'o' or 'w' : oil or water
   Character*1 :: wrt ! 'b' , 'c'  ! wrt to after , or before or itself

      ind = (k-1) * Nx * Ny + (j-1) * Nx + i
      N_ind = ind - Nx * Ny
      denomenator = 1 / 2.0 * ( Deltaz(N_ind) / ( kz(N_ind) * Deltax(i) * Deltay(j) ) + Deltaz(ind)/( kz(ind) * Deltax(i) * Deltay(j) ) )

      if ( D_Poten( ind , N_ind , m ) > 0 ) then ! upstream block
         if (wrt=='b') then
            press_term = 0
         Else
             if ( m == 'w') then
				press_term = DkrwdS( ind ) / ( mu_w ( ind ) * B_w( ind ) ) + k_rw(ind) / ( mu_w ( ind) * B_w(ind) ) * ( B_der(ind ,'w') / B_w(ind ) ) * PCO_der(ind)
             Else
				press_term = DkrodS( ind ) / ( mu_o ( ind ) * B_o( ind ) )
             EndIf
         EndIf
      else
         if (wrt=='c') then
            press_term = 0
         Else
            if ( m == 'w') then
				press_term = DkrwdS( N_ind ) / ( mu_w ( N_ind ) * B_w( N_ind ) ) + k_rw(N_ind) / ( mu_w ( N_ind) * B_w(N_ind) ) * ( B_der(N_ind ,'w') / B_w(N_ind ) ) * PCO_der(N_ind)
            Else
				press_term = DkrodS( N_ind ) / ( mu_o ( N_ind ) * B_o( N_ind ) )
            EndIf
         Endif
      EndIf

      DTrzb_DS = betac * press_term /denomenator

   Endfunction DTrzb_DS
!*******************************************************************************************

   real*8 function DTrza_DS(i,j,k,m,wrt) ! the value of the derivative of x transmissibility between a gridblock and a next grid
   integer :: ind , N_ind ! Neighbor grid index
   integer :: i,j,k
   real*8 :: denomenator, press_term , Press
   Character*1 :: m ! 'o' or 'w' : oil or water
   Character*1 :: wrt ! 'a' , 'c'  ! wrt to after , or before or itself

      ind = (k-1) * Nx * Ny + (j-1) * Nx + i
      N_ind = ind + Nx * Ny
      denomenator = 1 / 2.0 * ( Deltaz(N_ind) / ( kz(N_ind) * Deltax(i) * Deltay(j) ) + Deltaz(ind)/( kz(ind) * Deltax(i) * Deltay(j) ) )

      if ( D_Poten( ind , N_ind , m ) > 0 ) then ! upstream block
         if (wrt=='a') then
            press_term = 0
         Else
             if ( m == 'w') then
				press_term = DkrwdS( ind ) / ( mu_w ( ind ) * B_w( ind ) ) + k_rw(ind) / ( mu_w ( ind) * B_w(ind) ) * ( B_der(ind ,'w') / B_w(ind ) ) * PCO_der(ind)
             Else
				press_term = DkrodS( ind ) / ( mu_o ( ind ) * B_o( ind ) )
             EndIf
         EndIf
      else
         if (wrt=='c') then
            press_term = 0
         Else
            if ( m == 'w') then
				press_term = DkrwdS( N_ind ) / ( mu_w ( N_ind ) * B_w( N_ind ) ) + k_rw(N_ind) / ( mu_w ( N_ind) * B_w(N_ind) ) * ( B_der(N_ind ,'w') / B_w(N_ind ) ) * PCO_der(N_ind)
            Else
				press_term = DkrodS( N_ind ) / ( mu_o ( N_ind ) * B_o( N_ind ) )
            EndIf
         Endif
      EndIf

      DTrza_DS = betac * press_term /denomenator

   Endfunction DTrza_DS
!*******************************************************************************************
!*******************************************************************************************

  real*8 function B_der(ind , m )
    integer :: ind , n2 , m2 , i4   ! Neighbor grid index
    Character*1 :: m ! 'o' or 'w' : oil or water
    real*8 :: p2 , X

           if (m=='o') then
			  p2 = pressure(ind)
              n2 = N_Bo
              m2 = 1
              do while (abs(n2-m2) .NE. 1)
                i4 = FLOOR(m2 + (n2-m2+1)/2.0)
                if (p2 < Bo_p(i4,1) )  then
                    n2 = FLOOR(m2+(n2-m2+1)/2.0);
                else
                    m2 = FLOOR(m2+(n2-m2+1)/2.0);
                end if
              end do
              B_der = (Bo_p(n2,2)-Bo_p(m2,2))/(Bo_p(n2,1)-Bo_p(m2,1))
!			   B_der = -5.0d-5
           Else
              p2 = W_press(ind)
              !B_der = - PVTW(2) / ( ( 1 + PVTW(3) * (p2 - PVTW(1) ) )**2 ) * PVTW(3)
              X =  PVTW(3) * (p2 - PVTW(1) )
              B_der = - PVTW(2) / ( ( 1 + X + ( X ** 2 )/2.0 )**2 ) * ( PVTW(3) + ( PVTW(3) ** 2 ) * (p2 - PVTW(1) ) )

           EndIf
  Endfunction B_der

  !*******************************************************************************************

  real*8 function mu_der(ind , m )
    integer :: ind ! Neighbor grid index
    Character*1 :: m ! 'o' or 'w' : oil or water
    real*8 :: p1
    Integer :: m1,n1 , i5

    p1 = pressure(ind)
    if (m=='w') then
        mu_der = 0.0
    else
        n1 = N_mup
        m1 = 1
        do while (abs(n1-m1) .NE. 1)
            i5 = floor(m1 + (n1-m1+1)/2.0)
            if (p1 < mu_p(i5,1) )  then
                n1 = floor(m1+(n1-m1+1)/2.0);
            else
                m1 = floor(m1+(n1-m1+1)/2.0);
            end if
        end do
        mu_der = (mu_p(n1,2)-mu_p(m1,2))/(mu_p(n1,1)-mu_p(m1,1))
    Endif

  Endfunction mu_der

  !*******************************************************************************************

  real*8 function gama_der(ind , m)
    integer :: ind ! Ind is the block number which we want to have the derivative wrt the pressure of that block
    Character*1 :: m ! 'o' or 'w' : oil or water
    !Character*1 :: wrt ! 'a' , 'c'  ! wrt to after , or before or itself

    If ( m == 'o' ) then
       gama_der = - den_o / ( 144.0 * 2.0 * ( B_o(ind) ) **2 ) *  B_der(ind , 'o' )
    Else !( m == 'w' )
       gama_der = - den_w / ( 144.0 * 2.0 * ( B_w(ind) ) **2 ) *  B_der(ind , 'w' )
    EndIf

  Endfunction gama_der
  !*******************************************************************************************

  real*8 function Dkr_dS(index , m1)   ! this function calculates the Derivative of Kr wrt the saturation of the block

    Implicit None
    !
    integer :: index ! Ind is the block number which we want to have the derivative wrt the saturation of that block
    Character*1 :: m1 ! 'o' or 'w' : oil or water
    real*8 :: swater
    Integer :: m ,n

        swater = sw(index)
        n = n_Kr
        m = 1
        do while (abs(n-m) .NE. 1)
            if ( swater < Sw_mat(floor(m+(n-m+1)/2.0)))  then
                n = floor(m+(n-m+1)/2.0);
            else
                m = floor(m+(n-m+1)/2.0);
            end if
        end do
        if  ( m1 == 'o' ) then
			if ((swater > Sw_mat(n_Kr - 1) ) .OR. (swater .LT. Sw_mat(1) ) ) then
				Dkr_dS = 0.0
			Else
				Dkr_dS = (kro(n)-kro(m))/(Sw_mat(n)-Sw_mat(m))
			EndIf
        else
			if ( (swater > 1.0 ) .OR. (swater .LT. Sw_mat(1) ) ) then
				Dkr_dS = 0.0
			Else
				Dkr_dS = (krw(n)-krw(m))/(Sw_mat(n)-Sw_mat(m))
			EndIf
        EndIf

  Endfunction Dkr_dS

  !*******************************************************************************************
  real*8 function PCO_der(ind)
    integer :: ind
    real*8 :: Swater
    Integer :: m,n

        Swater = sw(ind)

		if ( Swater .LE. Sw_mat(1) ) then
			n = 2
			m = 1
		elseif ( Swater .GE. 1.0   ) then
			n = n_Kr
			m = n_Kr - 1
		else
			n = n_Kr
			m=1
			do while (abs(n-m) .NE. 1)
				if (Swater < Sw_mat( floor(m+(n-m+1)/2.0) ) ) then
					n = floor(m+(n-m+1)/2.0);
				else
					m = floor(m+(n-m+1)/2.0);
				end if
			end do
		EndIf
        PCO_der = (Pcow(n) - Pcow(m))/(Sw_mat(n) - Sw_mat(m))
	   ! PCO_der = - 3 * 34.98 * Swater ** 2 + 2 * 83.96 * Swater - 67.17

  Endfunction PCO_der
  !*******************************************************************************************
!*******************************************************************************************

  real*8 function B_der_P(p , m )
    integer :: ind , n2 , m2 , i4   ! Neighbor grid index
    Character*1 :: m ! 'o' or 'w' : oil or water
    real*8 :: p , p2 , X

           p2 = p
           if (m=='o') then
              n2 = N_Bo
              m2 = 1
              do while (abs(n2-m2) .NE. 1)
                i4 = FLOOR(m2 + (n2-m2+1)/2.0)
                if (p2 < Bo_p(i4,1) )  then
                    n2 = FLOOR(m2+(n2-m2+1)/2.0);
                else
                    m2 = FLOOR(m2+(n2-m2+1)/2.0);
                end if
              end do
              B_der_P = (Bo_p(n2,2)-Bo_p(m2,2))/(Bo_p(n2,1)-Bo_p(m2,1))
!			   B_der = -5.0d-5
           Else
              !B_der = - PVTW(2) / ( ( 1 + PVTW(3) * (p2 - PVTW(1) ) )**2 ) * PVTW(3)
              X =  PVTW(3) * (p2 - PVTW(1) )
              B_der_P = - PVTW(2) / ( ( 1 + X + ( X ** 2 )/2.0 )**2 ) * ( PVTW(3) + ( PVTW(3) ** 2 ) * (p2 - PVTW(1) ) )
           EndIf

  Endfunction B_der_P

  !*******************************************************************************************
  !*******************************************************************************************
  real*8 function PCO_der_Sw(Swater)
    integer :: ind
    real*8 :: Swater
    Integer :: m,n

		if ( Swater .LE. Sw_mat(1) ) then
			n = 2
			m = 1
		elseif ( Swater .GE. 1.0   ) then
			n = n_Kr
			m = n_Kr - 1
		else
			n = n_Kr
			m=1
			do while (abs(n-m) .NE. 1)
				if (Swater < Sw_mat( floor(m+(n-m+1)/2.0) ) ) then
					n = floor(m+(n-m+1)/2.0);
				else
					m = floor(m+(n-m+1)/2.0);
				end if
			end do
		EndIf
        PCO_der_Sw = (Pcow(n) - Pcow(m))/(Sw_mat(n) - Sw_mat(m))
	   ! PCO_der = - 3 * 34.98 * Swater ** 2 + 2 * 83.96 * Swater - 67.17

  Endfunction PCO_der_Sw
  !*******************************************************************************************

  EndModule Derivative