     module soil_module

     implicit none

       type soilayer
        real :: ec = 0.
        real :: cal = 0.
        real :: ph = 0.
        real :: alb = 0.         !! none          albedo when soil is moist
        real :: usle_k = 0.      !!               USLE equation soil erodibility (K) factor
        real ::conk = 0.         !! mm/hr          lateral saturated hydraulic conductivity for each profile layer in a give HRU.
        real ::flat = 0.         !! mm H2O         lateral flow storage array
        real :: prk = 0.         !! mm H2O         percolation from soil layer on current day
        real :: volcr = 0.       !! mm             crack volume for soil layer
        real :: tillagef = 0.
        real :: rtfr = 0.        !! none           root fraction
        real :: watp = 0.
        integer :: a_days = 0
        integer :: b_days = 0
        real :: psp_store = 0.
        real :: ssp_store = 0.
        real :: percc = 0.       !!
        real :: latc = 0.        !!
        real :: vwt = 0.         !!
      end type soilayer
      type (soilayer), dimension(:), allocatable :: layer1

      type soil_physical_properties
        real :: d = 0.            !! mm            depth to bottom of soil layer
        real :: thick = 0.        !! mm            thichness of soil layer
        real :: bd = 0.           !! Mg/m**3       bulk density of the soil
        real :: k = 0.            !! mm/hr         saturated hydraulic conductivity of soil layer. Index:(layer,HRU)
        real :: clay = 0.         !! none          fraction clay content in soil material (UNIT CHANGE!)
        real :: silt = 0.         !! %             percent silt content in soil material
        real :: sand = 0.         !! none          fraction of sand in soil material
        real :: rock = 0.         !! %             percent of rock fragments in soil layer
        real :: conv_wt = 0.       !! none          factor which converts kg/kg to kg/ha
        real :: crdep = 0.         !! mm            maximum or potential crack volume
        real :: awc = 0.           !! mm H20/mm     soil available water capacity of soil layer
        real :: fc = 0.           !! mm H2O         amount of water available to plants in soil layer at field capacity (fc - wp),Index:(layer,HRU)
        real :: hk = 0.           !! none           beta coefficent to calculate hydraulic conductivity
        real :: por = 0.         !! none           total porosity of soil layer expressed as a fraction of the total volume, Index:(layer,HRU)
        real :: st = 0.          !! mm H2O         amount of water stored in the soil layer on any given day (less wp water)
        real :: tmp = 0.         !! deg C          daily average temperature of second soil layer
        real :: ul = 0.          !! mm H2O         amount of water held in the soil layer at saturation (sat - wp water)
        real :: up = 0.          !! mm H2O/mm      soil water content of soil at -0.033 MPa (field capacity)
        real :: wp = 0.          !! mm H20/mm      soil water content of soil at -1.5 MPa (wilting point)
        real :: wpmm = 0.        !! mm H20         water content of soil at -1.5 MPa (wilting point)
      end type soil_physical_properties

      type soil_profile
        character(len=16) :: snam = ""     !! NA            soil series name
        character(len=16) :: hydgrp = ""    !! NA            hydrologic soil group
        character(len=16) :: texture = ""
        integer ::  nly  = 0               !! none          number of soil layers
        type (soil_physical_properties),dimension (:), allocatable::phys
        type (soilayer), dimension (:), allocatable :: ly
        real, dimension(:),allocatable :: pest              !! kg/ha    total pesticide in the soil profile
        real :: zmx = 0.                   !! mm            maximum rooting depth in soil
        real :: anion_excl = 0.            !! none          fraction of porosity from which anions are excluded
        real :: crk = 0.                   !! none          crack volume potential of soil
        real :: alb = 0.                   !! none          albedo when soil is moist
        real :: usle_k = 0.                !!               USLE equation soil erodibility (K) factor
        real :: det_san = 0.
        real :: det_sil = 0.
        real :: det_cla = 0.
        real :: det_sag = 0.
        real :: det_lag = 0.
        real :: sumul = 0.                 !! mm H2O         amount of water held in soil profile at saturation
        real :: sumfc = 0.                 !! mm H2O         amount of water held in the soil profile at field capacity
        real :: sw = 0.                    !! mm H2O         amount of water stored in soil profile
        real :: sw_300 = 0.                !! mm H2O         amount of water stored to 300 mm
        real :: sumwp = 0.                 !!
        real :: swpwt = 0.                 !!
        real :: ffc = 0.                   !! none           initial HRU soil water content expressed as fraction of field capacity
        real :: wat_tbl = 0.               !!
        real :: avpor = 0.                 !! none           average porosity for entire soil profile
        real :: avbd = 0.                  !! Mg/m^3         average bulk density for soil profile
      end type soil_profile
      type (soil_profile), dimension(:), allocatable :: soil
      type (soil_profile), dimension(:), allocatable :: soil_init

      type soil_hru_database
         character(len=16) :: snam = ""     !! NA            soil series name
         character(len=16) :: hydgrp = ""    !! NA            hydrologic soil group
         character(len=16) :: texture = ""
         type (soil_profile) :: s
         type (soil_physical_properties),dimension(:), allocatable::phys
         type (soilayer), dimension(:), allocatable :: ly
      end type soil_hru_database
      type (soil_hru_database), dimension(:), allocatable :: sol

     end module soil_module


text:
!< Wenoof objects factory.
module wenoof_objects_factory
!< Wenoof factory.

use wenoof_alpha_coefficients
use wenoof_alpha_coefficients_js
use wenoof_alpha_coefficients_m
use wenoof_alpha_coefficients_z
use wenoof_base_object
use wenoof_optimal_weights
use wenoof_optimal_weights_js
use wenoof_polynomials
use wenoof_polynomials_js
use wenoof_smoothness_indicators
use wenoof_smoothness_indicators_js

implicit none
private
public :: objects_factory

type :: objects_factory
  !< Factory, create an instance of concrete extension of [[base_object]] given its constructor.
  contains
    ! public methods
    generic :: create => create_alpha_coefficients, &
                         create_optimal_weights,    &
                         create_polynomials,        &
                         create_smoothness_indicators !< Create a concrete instance of [[alpha_coefficients]] or
                                                      !< [[optimal_weights]] or [[polynomials]] or [[smoothness_indicators]].
    procedure, nopass :: create_base_object           !< Create a concrete instance of [[base_object]].
    ! private methods
    procedure, nopass, private :: create_alpha_coefficients    !< Create a concrete instance of [[alpha_coefficients]].
    procedure, nopass, private :: create_optimal_weights       !< Create a concrete instance of [[optimal_weights]].
    procedure, nopass, private :: create_polynomials           !< Create a concrete instance of [[polynomials]].
    procedure, nopass, private :: create_smoothness_indicators !< Create a concrete instance of [[smoothness_indicators]].
endtype objects_factory

contains
  subroutine create_base_object(constructor, object)
  !< Create an instance of concrete extension of [[base_object]] given its constructor.
  class(base_object_constructor),  intent(in)  :: constructor !< Constructor.
  class(base_object), allocatable, intent(out) :: object      !< Object.

  select type(constructor)
  type is(alpha_coefficients_js_constructor)
    allocate(alpha_coefficients_js :: object)
  type is(alpha_coefficients_m_constructor)
    allocate(alpha_coefficients_m :: object)
  type is(alpha_coefficients_z_constructor)
    allocate(alpha_coefficients_z :: object)
  type is(optimal_weights_js_constructor)
    allocate(optimal_weights_js :: object)
  type is(polynomials_js_constructor)
    allocate(polynomials_js :: object)
  type is(smoothness_indicators_js_constructor)
    allocate(smoothness_indicators_js :: object)
  class default
    error stop 'error: WenOOF object factory do NOT know the constructor given'
  endselect
  call object%create(constructor=constructor)
  endsubroutine create_base_object

  subroutine create_alpha_coefficients(constructor, object)
  !< Create an instance of concrete extension of [[alpha_coefficients]] given its constructor.
  class(alpha_coefficients_constructor),  intent(in)  :: constructor !< Constructor.
  class(alpha_coefficients), allocatable, intent(out) :: object      !< Object.

  select type(constructor)
  type is(alpha_coefficients_js_constructor)
    allocate(alpha_coefficients_js :: object)
  type is(alpha_coefficients_m_constructor)
    allocate(alpha_coefficients_m :: object)
  type is(alpha_coefficients_z_constructor)
    allocate(alpha_coefficients_z :: object)
  class default
    error stop 'error: WenOOF object factory do NOT know the constructor given'
  endselect
  call object%create(constructor=constructor)
  endsubroutine create_alpha_coefficients

  subroutine create_optimal_weights(constructor, object)
  !< Create an instance of concrete extension of [[optimal_weights]] given its constructor.
  class(optimal_weights_constructor),  intent(in)  :: constructor !< Constructor.
  class(optimal_weights), allocatable, intent(out) :: object      !< Object.

  select type(constructor)
  type is(optimal_weights_js_constructor)
    allocate(optimal_weights_js :: object)
  class default
    error stop 'error: WenOOF object factory do NOT know the constructor given'
  endselect
  call object%create(constructor=constructor)
  endsubroutine create_optimal_weights

  subroutine create_polynomials(constructor, object)
  !< Create an instance of concrete extension of [[polynomials]] given its constructor.
  class(polynomials_constructor),  intent(in)  :: constructor !< Constructor.
  class(polynomials), allocatable, intent(out) :: object      !< Object.

  select type(constructor)
  type is(polynomials_js_constructor)
    allocate(polynomials_js :: object)
  class default
    error stop 'error: WenOOF object factory do NOT know the constructor given'
  endselect
  call object%create(constructor=constructor)
  endsubroutine create_polynomials

  subroutine create_smoothness_indicators(constructor, object)
  !< Create an instance of concrete extension of [[smoothness_indicators]] given its constructor.
  class(smoothness_indicators_constructor),  intent(in)  :: constructor !< Constructor.
  class(smoothness_indicators), allocatable, intent(out) :: object      !< Object.

  select type(constructor)
  type is(smoothness_indicators_js_constructor)
    allocate(smoothness_indicators_js :: object)
  class default
    error stop 'error: WenOOF object factory do NOT know the constructor given'
  endselect
  call object%create(constructor=constructor)
  endsubroutine create_smoothness_indicators
endmodule wenoof_objects_factory

max_stars_repo_name: giacrossi/WenOOF

max_issues_repo_name: giacrossi/WenOOF

max_forks_repo_name: giacrossi/WenOOF


text:
      PROGRAM TEZMAP
C
C Define the error file, the Fortran unit number, the workstation type,
C and the workstation ID to be used in calls to GKS routines.
C
C     PARAMETER (IERRF=6, LUNIT=2, IWTYPE=1,  IWKID=1)   ! NCGM
C     PARAMETER (IERRF=6, LUNIT=2, IWTYPE=8,  IWKID=1)   ! X Windows
C     PARAMETER (IERRF=6, LUNIT=2, IWTYPE=11, IWKID=1)   ! PDF
C     PARAMETER (IERRF=6, LUNIT=2, IWTYPE=20, IWKID=1)   ! PostScript
C
      PARAMETER (IERRF=6, LUNIT=2, IWTYPE=1,  IWKID=1)
C
C OPEN GKS, OPEN WORKSTATION OF TYPE 1, ACTIVATE WORKSTATION
C
      CALL GOPKS (IERRF, ISZDM)
      CALL GOPWK (IWKID, LUNIT, IWTYPE)
      CALL GACWK (IWKID)
C
C INVOKE DEMO DRIVER
C
      CALL EZMAP(IERR)
C
C DEACTIVATE AND CLOSE WORKSTATION, CLOSE GKS.
C
      CALL GDAWK (IWKID)
      CALL GCLWK (IWKID)
      CALL GCLKS
C
      STOP
      END
C
      SUBROUTINE EZMAP (IERROR)
C
C PURPOSE                To provide a simple demonstration of
C                        the mapping utility, EZMAP.
C
C USAGE                  CALL EZMAP (IERROR)
C
C ARGUMENTS
C
C ON OUTPUT              IERROR
C                          An integer variable
C                          = 0, if the test was successful,
C                          = 1, the test was not successful.
C
C I/O                    The map dataset must be connected to a
C                        FORTRAN unit.  See the open statement in
C                        this test deck.
C
C                        If the test is successful, the message
C
C              EZMAP TEST EXECUTED--SEE PLOTS TO CERTIFY
C
C                        is printed on unit 6.  In addition, 10
C                        frames are produced on the machine graphics
C                        device.  In order to determine if the test
C                        was successful, it is necessary to examine
C                        the plots.
C
C PRECISION              Single
C
C LANGUAGE               FORTRAN 77
C
C REQUIRED ROUTINES      EZMAP, GRIDAL
C
C REQUIRED GKS LEVEL     0A
C
C ALGORITHM              TEZMAP calls routine MAPROJ once for each of
C                        the 10 types of projections:
C                          stereographic,
C                          orthographic,
C                          Lambert conformal conic with 2 standard
C                            parallels,
C                          Lambert equal area,
C                          gnomonic,
C                          azimuthal equidistant,
C                          cylindrical equidistant,
C                          Mercator,
C                          Mollweide type, and
C                          satellite view.
C
C HISTORY                Written October, 1976, main entry SUPMAP.
C                        Updated September, 1985, main entry EZMAP.
C
      DIMENSION PLM1(2),PLM2(2),PLM3(2),PLM4(2)
      COMMON /ERROR/ IFRAME, IERRR
C
      DATA PLM1,PLM2,PLM3,PLM4 /8*0.0/
C
C Define the center of a plot title string on a square grid of size
C 0. to 1.
C
      DATA TX/0.5/, TY/0.9765/
C
C Initialize the error parameter.
C
      IERRR = 0
C
C Turn on the error recovery mode.
C
      CALL ENTSR(IDUM,1)
C
C
C     Frame 1 -- The stereographic projection.
C
       IFRAME = 1
C
C Set the projection-type parameter.
C
      CALL MAPROJ('ST',80.0,-160.0,0.0)
C
C Set the limit parameters.
C
      CALL MAPSET('MA',PLM1,PLM2,PLM3,PLM4)
C
C Set the outline-dataset parameter.
C
      CALL MAPSTC('OU','PS')
C
C Draw the map.
C
      CALL MAPDRW
C
C Report any errors encountered.
C
       IF( NERRO(IERR).NE.0) CALL RPTERR
C
C Select normalization transformation 0.
C
      CALL GSELNT(0)
C
C Call PLCHLQ to write the plot title.
C
      CALL PLCHLQ(TX,TY,'EZMAP DEMONSTRATION:  STEREOGRAPHIC PROJECTION'
     +           ,16.,0.,0.)
      CALL FRAME
 10   CONTINUE
C
C
C     Frame 2 -- The orthographic projection.
C
      IFRAME = 2
C
C
C Set the projection-type parameter.
C
      CALL MAPROJ('OR',60.0,-120.0,0.0)
C
C Draw the map.
C
      CALL MAPDRW
C Report any errors encountered.
C
      IF( NERRO(IERR).NE.0) CALL RPTERR
C
C Write the title.
C
C
C Select normalization transformation 0.
C
      CALL GSELNT(0)
C
C Call PLCHLQ to write the plot title.
C
      CALL PLCHLQ(TX,TY,
     +     'EZMAP DEMONSTRATION:  ORTHOGRAPHIC PROJECTION',16.,0.,0.)
      CALL FRAME
 20   CONTINUE
C
C
C     Frame 3 -- The Lambert conformal conic projection.
C
      IFRAME = 3
C
C
C Set the projection-type, limits, and outline-dataset parameters.
C
      CALL MAPROJ('LC',45.0,-100.0,45.0)
      CALL MAPSET('CO',50.0,-130.0,20.0,-75.0)
      CALL MAPSTC('OU','US')
C
C Draw the map.
C
      CALL MAPDRW
C
C Report any errors encountered.
C
      IF( NERRO(IERR).NE.0) CALL RPTERR
C
C Select normalization transformation 0.
C
      CALL GSELNT(0)
C
C Call PLCHLQ to write the plot title.
C
      CALL PLCHLQ(TX,TY,
     +     'EZMAP DEMONSTRATION: LAMBERT CONFORMAL CONIC PROJECTION',
     +           16.,0.,0.)
      CALL FRAME
 30   CONTINUE
C
C
C     Frame 4 -- The Lambert equal area projection.
C
      IFRAME = 4
C
C Set the projection-type, limits, and outline-dataset parameters.
C
      CALL MAPROJ('LE',20.0,-40.0,0.0)
      CALL MAPSET('MA',PLM1,PLM2,PLM3,PLM4)
      CALL MAPSTC('OU','CO')
C
C Draw the map.
C
      CALL MAPDRW
C
C Report any errors encountered.
C
      IF( NERRO(IERR).NE.0) CALL RPTERR
C
C Select normalization transformation 0.
C
      CALL GSELNT(0)
C
C Call PLCHLQ to write the plot title.
C
      CALL PLCHLQ(TX,TY,
     +     'EZMAP DEMONSTRATION:  LAMBERT EQUAL AREA PROJECTION',
     +     16.,0.,0.)
      CALL FRAME
 40   CONTINUE
C
C
C     Frame 5 -- The gnomonic projection.
C
      IFRAME = 5
C
C Set the projection-type parameter.
C
      CALL MAPROJ('GN',0.0,0.0,0.0)
C
C Draw the map.
C
      CALL MAPDRW
C Report any errors encountered.
C
      IF( NERRO(IERR).NE.0) CALL RPTERR
C
C Select normalization transformation 0.
C
      CALL GSELNT(0)
C
C Call PLCHLQ to write the plot title.
C
      CALL PLCHLQ(TX,TY,'EZMAP DEMONSTRATION:  GNOMONIC PROJECTION',
     +           16.,0.,0.)
      CALL FRAME
 50   CONTINUE
C
C
C     Frame 6 -- The azimuthal equidistant projection.
C
      IFRAME = 6
C
C Set the projection-type parameter.
C
      CALL MAPROJ('AE',-20.0,40.0,0.0)
C
C Set the grid spacing.
C
      CALL MAPSTR('GR',5.0)
C
C Draw the map.
C
      CALL MAPDRW
C
C Report any errors encountered.
C
      IF( NERRO(IERR).NE.0) CALL RPTERR
C
C Select normalization transformation 0.
C
      CALL GSELNT(0)
C
C Call PLCHLQ to write the plot title.
C
      CALL PLCHLQ(TX,TY,
     +     'EZMAP DEMONSTRATION:  AZIMUTHAL EQUIDISTANT PROJECTION',
     +           16.,0.,0.)
      CALL FRAME
 60   CONTINUE
C
C
C     Frame 7 -- The cylindrical equidistant projection.
C
      IFRAME = 7
C
C Set the map projection type parameter.
C
      CALL MAPROJ('CE',-40.0,80.0,0.0)
C
C Draw the map.
C
      CALL MAPDRW
C
C Report any errors encountered.
C
      IF( NERRO(IERR).NE.0) CALL RPTERR
C
C Select normalization transformation 0.
C
      CALL GSELNT(0)
C
C Call PLCHLQ to write the plot title.
C
      CALL PLCHLQ(TX,TY,
     +     'EZMAP DEMONSTRATION:  CYLINDRICAL EQUIDISTANT PROJECTION',
     +           16.,0.,0.)
      CALL FRAME
 70   CONTINUE
C
C
C     Frame 8 -- The mercator projection.
C
      IFRAME = 8
C
C Set the map projection type parameter.
C
      CALL MAPROJ('ME',-60.0,120.0,0.0)
C
C Draw the map.
C
      CALL MAPDRW
C
C Report any errors encountered.
C
      IF( NERRO(IERR).NE.0) CALL RPTERR
C
C Select normalization transformation 0.
C
      CALL GSELNT(0)
C
C Call PLCHLQ to write the plot title.
C
      CALL PLCHLQ(TX,TY,'EZMAP DEMONSTRATION: MERCATOR PROJECTION',
     +           16.,0.,0.)
      CALL FRAME
 80   CONTINUE
C
C
C     Frame 9 -- The Mollyweide-type projection.
C
      IFRAME = 9
C
C Set the map projection type parameter.
C
      CALL MAPROJ('MO',-80.0,160.0,0.0)
C
C Draw the map.
C
      CALL MAPDRW
C
C Report any errors encountered.
C
      IF( NERRO(IERR).NE.0) CALL RPTERR
C
C Select normalization transformation 0.
C
      CALL GSELNT(0)
C
C Call PLCHLQ to write the plot title.
C
      CALL PLCHLQ(TX,TY,
     +     'EZMAP DEMONSTRATION:  MOLLWEIDE-TYPE PROJECTION',
     +     16.,0.,0.)
      CALL FRAME
 90   CONTINUE
C
C
C     Frame 10 -- The satellite view projection.
C
      IFRAME = 10
C
C Set the map projection type parameter.
C
      CALL MAPROJ('SV',0.0,-135.0,0.0)
C
C Set the satellite distance and supress grid lines.
C
      CALL MAPSTR('SA',6.631)
      CALL MAPSTI('GR',0)
C
C Draw the map.
C
      CALL MAPDRW
C
C Report any errors encountered.
C
      IF( NERRO(IERR).NE.0) CALL RPTERR
C
C Select normalization transformation 0.
C
      CALL GSELNT(0)
C
C Call PLCHLQ to write the plot title.
C
      CALL PLCHLQ(TX,TY,
     +     'EZMAP DEMONSTRATION: SATELLITE VIEW PROJECTION',16.,0.,0.)
      CALL FRAME
 100  CONTINUE
C
      IF(IERRR.EQ.0) WRITE(6,1000)
      IF(IERRR.EQ.1) WRITE(6,1001)
      IERROR = IERR
      RETURN
C
 1000 FORMAT(' EZMAP TEST EXECUTED--SEE PLOTS TO CERTIFY')
 1001 FORMAT(' EZMAP TEST UNSUCCESSFUL')
      END
      SUBROUTINE RPTERR
C
C ROUTINE TO REPORT ERROR MESSEGES
C
      COMMON /ERROR/ IFRAME,IERRR
C
      WRITE(6,1000) IFRAME
      CALL EPRIN
      WRITE(6,1001)
      CALL ERROF
      IERRR = 1
      RETURN
C
 1000 FORMAT(' ERROR IN FRAME ',I2,5X,'ERROR MESSAGE FOLLOWS:')
 1001 FORMAT(' ******',//)
      END

ext: f

max_stars_repo_name: tenomoto/ncl

max_issues_repo_name: tenomoto/ncl

max_forks_repo_name: tenomoto/ncl


text:
c
c ----------------------------------------------------------------------
c
      subroutine iosgld(fi,fo,sc,gc,smass,mNLSP,idecay , vertex)
c
c This subroutine computes an amplitude of the fermion-fermion-scalar
c SUSY Goldstino coupling. In this routine, the NLSP is decaying to a
c boson and a Goldstino. The h.c. of the NLSP decay is handled via the
c input parameter idecay.
c
c input:
c       complex fi(6)          : flow-in  fermion                   |fi>
c       complex fo(6)          : flow-out fermion                   <fo|
c       complex sc(3)          : input    scalar                      s
c       complex gc(2)          : coupling constants                  gsf
c       real    mNLSP          : mass of the NLSP
c       real    smass          : mass of the scalar boson
c       integer idecay         :  1 for NLSP decay to Goldstino
c                              : -1 for Goldstino to NLSP (h.c. of above)
c
c output:
c       complex vertex         : amplitude                     <fo|s|fi>
c
      implicit none
      double complex  fi(6), fo(6), gc(2), sc(3), vertex
      double complex  p14p, p14m, p23p, p23m
      double precision  mNLSP, smass
      integer idecay

      double complex ci, cZero
      parameter( ci = ( 0.0d0, 1.0d0 ), cZero = ( 0.0d0, 0.0d0 ) )
c
      if ( idecay.ne.1 .or. idecay.ne.-1 ) then
         write(6,*) 'error in idecay of IOSGLD'
         stop
      end if

      p14p = dble(sc(2)) + dimag(sc(2))
      p14m = dble(sc(2)) - dimag(sc(2))
      p23p = dble(sc(3)) + dimag(sc(3))*ci
      p23m = dble(sc(3)) - dimag(sc(3))*ci

      vertex = gc(1) *
     &         ( ( ( fo(3)*p14p + fo(4)*p23p )*fi(1)
     &            +( fo(3)*p23m + fo(4)*p14m )*fi(2) )*mNLSP
     &          -( fo(1)*fi(1) + fo(2)*fi(2) )*idecay*smass**2 )

      if ( gc(2).ne.cZero ) then
         vertex = vertex + gc(2) *
     &            ( ( ( fo(1)*p14m - fo(2)*p23p )*fi(3)
     &               -( fo(1)*p23m - fo(2)*p14p )*fi(4) )*mNLSP
     &             -( fo(3)*fi(3) + fo(4)*fi(4) )*idecay*smass**2 )
      end if

      vertex = vertex * sc(1)
c
      return
      end


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

! Program to compute nodal factors and equilibrium arguements.
! Based on "tide_fac", implemented additional feature to change the time
! inside "bctides.in" for SCHISM easily.
! Many thanks to the author of "tide_fac".
! Feel free to email kai.li.jx#gmail.com for advice, issues and criticism.
!

integer,parameter :: ncnst=37
character :: cname(ncnst)*8, tmp, arg*20
common /cnsnam/ cname
real :: nodfac,month
integer :: ncon(200)
common /cnst/ nodfac(ncnst),grterm(ncnst),speed(ncnst),p(ncnst)
logical :: lselfe=.false.

open(unit=11,file='tide_fac.out',status='unknown')
if (iargc().eq.0) then
!  write(*,*) 'Pass "-h" option for command line arguments instruction'
  write(*,*) 'SCHISM trunk/src/Utility/Tide/tf_main.f90 -- calculate nodal factors and re-write bctides.in'
  write(*,*) 'enter length of run time (days)'
  read(*,*) xdays
  write(*,*) 'enter start time - bhr,iday,imo,iyr (iyr e.g. 1992)'
  read(*,*) bhr,iday,imo,iyr
  write(*,*) 'generate "bctides.in.out" file for SCHISM? ( "bctides.in" is required for input )'
  read(*,*)tmp
  if (tmp=='y'.or.tmp=='Y') lselfe=.true.
else
  write(*,*) 'command line parsing not supported yet'
  stop
endif
rhrs=xdays*24.
yr=iyr
month=imo
day=iday
hrm=bhr+rhrs/2.
write(11,10) bhr,iday,imo,iyr
write(*,10) bhr,iday,imo,iyr
10 format(' tidal factors starting: ',' hr-',f5.2,',  day-',i3,',  month-',i3,'  year-',i5,/)
write(*,11) xdays
write(11,11) xdays
11 format(' for a run lasting ',f8.2,' days',//)
!-- determine the julian time at beginning and middle of record
dayj=dayjul(yr,month,day)
!-- determine node factors at middle of record
call nfacs(yr,dayj,hrm)
!-- determine greenwich equil. terms at beginning of record
call gterms(yr,dayj,bhr,dayj,hrm)
if (lselfe) then
  call selecttides(numcon,ncon)
else
  numcon=8
  ncon(1)=4
  ncon(2)=6
  ncon(3)=30
  ncon(4)=26
  ncon(5)=3
  ncon(6)=1
  ncon(7)=2
  ncon(8)=35
endif
write(11,*) 'const   node     eq arg (ref gm)'
write(11,1300)
1300 format(' name   factor    (deg) ',//)
do 20 nc=1,numcon
  ic=ncon(nc)
! equilibrium arguement is referenced to the greenwich meridian
  if (ic==0) then
    write(11,2001) 'z0  ',1.0,0.0
  else
    write(11,2001) cname(ic),nodfac(ic),grterm(ic)
  endif
  2001 format(1x,a4,2x,f8.5,4x,f7.2,2x,f7.4)
20 continue
if (lselfe) call writebctides(ncon,bhr,iday,imo,iyr)
stop
end

subroutine nfacs(yr,dayj,hr)

!-- calculates node factors for constituent tidal signal

!-- the equations used in this routine come from:
!         "manual of harmonic analysis and prediction of tides"
!         by paul schureman, special publication #98, us coast
!         and geodetic survey, department of commerce (1958).

!-- if daym and hrm correspond to midyear, then this routine
!-- returns the same values as found in table 14 of schureman.
!---------------------------------------------------------------------

  character*8   cst(37)
  real          i,n,nu

  common/orbitf/ds,dp,dh,dp1,dn,di,dnu,dxi,dnup,dnup2,dpc
  common/ cnst /fndcst(37),eqcst(37),acst(37),pcst(37)
  common/cnsnam/cst

!-- constituent names:
  data cst     /'M2      ','S2      ','N2      ','K1      ',&
                    'M4      ','O1      ','M6      ','MK3     ',&
                    'S4      ','MN4     ','nu2     ','S6      ',&
                    'mu2     ','2N2     ','OO1     ','lambda2 ',&
                    'S1      ','M1      ','J1      ','Mm      ',&
                    'Ssa     ','Sa      ','MSf     ','Mf      ',&
                    'rho1    ','Q1      ','T2      ','R2      ',&
                    '2Q1     ','P1      ','2SM2    ','M3      ',&
                    'L2      ','2MK3    ','K2      ','M8      ',&
                    'MS4     '/  ! to be added: 2MS6 2SM6 SK3 2MN6 eta2 ups1 MNS2 SN4

!-- orbital speeds (degrees/hour):
  data acst/28.984104252,30.0000000,28.439729568,15.041068632,&
            57.968208468,13.943035584,86.952312720,44.025172884,&
            60.0,57.423833820,28.5125831,90.0,&
            27.9682084,27.8953548,16.139101680,29.4556253,&
            15.0,14.496693984,15.5854433,0.5443747,&
            0.0821373,0.0410686,1.0158957720,1.0980331,&
            13.4715145,13.3986609,29.9589333,30.0410667,&
            12.854286252,14.9589314,31.01589576,43.476156360,&
            29.5284789,42.927139836,30.0821373,115.936416972,&
            58.984104240/

!-- number of tide cycles per day per constituent:
  data pcst/2.,2.,2.,1.,4.,1.,6.,3.,4.,4.,2.,6.,2.,2.,1.,2.,1.,1.,&
  1.,0.,0.,0.,0.,0.,1.,1.,2.,2.,1.,1.,2.,3.,2.,3.,2.,8.,4./

  pi180=3.14159265/180.
  call orbit(yr,dayj,hr)
  n=dn*pi180
  i=di*pi180
  nu=dnu*pi180
  xi=dxi*pi180
  p=dp*pi180
  pc=dpc*pi180
  sini=sin(i)
  sini2=sin(i/2.)
  sin2i=sin(2.*i)
  cosi2=cos(i/2.)
  tani2=tan(i/2.)
!-- equation 197, schureman
  qainv=sqrt(2.310+1.435*cos(2.*pc))
!-- equation 213, schureman
  rainv=sqrt(1.-12.*tani2**2*cos(2.*pc)+36.*tani2**4)
!-- variable names refer to equation numbers in schureman
  eq73=(2./3.-sini**2)/.5021
  eq74=sini**2/.1578
  eq75=sini*cosi2**2/.37988
  eq76=sin(2*i)/.7214
  eq77=sini*sini2**2/.0164
  eq78=(cosi2**4)/.91544
  eq149=cosi2**6/.8758
  eq207=eq75*qainv
  eq215=eq78*rainv
  eq227=sqrt(.8965*sin2i**2+.6001*sin2i*cos(nu)+.1006)
  eq235=.001+sqrt(19.0444*sini**4+2.7702*sini**2*cos(2.*nu)+.0981)
!-- node factors for 37 constituents:
  fndcst(1)=eq78
  fndcst(2)=1.0
  fndcst(3)=eq78
  fndcst(4)=eq227
  fndcst(5)=fndcst(1)**2
  fndcst(6)=eq75
  fndcst(7)=fndcst(1)**3
  fndcst(8)=fndcst(1)*fndcst(4)
  fndcst(9)=1.0
  fndcst(10)=fndcst(1)**2
  fndcst(11)=eq78
  fndcst(12)=1.0
  fndcst(13)=eq78
  fndcst(14)=eq78
  fndcst(15)=eq77
  fndcst(16)=eq78
  fndcst(17)=1.0
!** equation 207 not producing correct answer for m1
!**set node factor for m1 = 0 until can further research
  fndcst(18)=0.
!  fndcst(18)=eq207
  fndcst(19)=eq76
  fndcst(20)=eq73
  fndcst(21)=1.0
  fndcst(22)=1.0
  fndcst(23)=eq78
  fndcst(24)=eq74
  fndcst(25)=eq75
  fndcst(26)=eq75
  fndcst(27)=1.0
  fndcst(28)=1.0
  fndcst(29)=eq75
  fndcst(30)=1.0
  fndcst(31)=eq78
  fndcst(32)=eq149
!** equation 215 not producing correct answer for l2
!** set node factor for l2 = 0 until can further research
  fndcst(33)=0.
!  fndcst(33)=eq215
  fndcst(34)=fndcst(1)**2*fndcst(4)
  fndcst(35)=eq235
  fndcst(36)=fndcst(1)**4
  fndcst(37)=eq78
end subroutine nfacs

subroutine gterms(yr,dayj,hr,daym,hrm)
!-- calculates equilibrium arguments v0+u for constituent tide

!-- the equations used in this routine come from:
!         "manual of harmonic analysis and prediction of tides"
!         by paul schureman, special publication #98, us coast
!         and geodetic survey, department of commerce (1958).

!-- if daym and hrm correspond to midyear, then this routine
!-- returns the same values as found in table 15 of schureman.
!---------------------------------------------------------------------
  real nu,nup,nup2,i
  common /orbitf/ds,dp,dh,dp1,dn,di,dnu,dxi,dnup,dnup2,dpc
  common /cnst/ fndcst(37),eqcst(37),acst(37),pcst(37)
  pi180=3.14159265/180.
!* obtaining orbital values at beginning of series for v0
  call orbit(yr,dayj,hr)
  s=ds
  p=dp
  h=dh
  p1=dp1
  t=angle(180.+hr*(360./24.))
!** obtaining orbital values at middle of series for u
  call orbit(yr,daym,hrm)
  nu=dnu
  xi=dxi
  nup=dnup
  nup2=dnup2
!* summing terms to obtain equilibrium arguments
  eqcst(1)=2.*(t-s+h)+2.*(xi-nu)
  eqcst(2)=2.*t
  eqcst(3)=2.*(t+h)-3.*s+p+2.*(xi-nu)
  eqcst(4)=t+h-90.-nup
  eqcst(5)=4.*(t-s+h)+4.*(xi-nu)
  eqcst(6)=t-2.*s+h+90.+2.*xi-nu
  eqcst(7)=6.*(t-s+h)+6.*(xi-nu)
  eqcst(8)=3.*(t+h)-2.*s-90.+2.*(xi-nu)-nup
  eqcst(9)=4.*t
  eqcst(10)=4.*(t+h)-5.*s+p+4.*(xi-nu)
  eqcst(11)=2.*t-3.*s+4.*h-p+2.*(xi-nu)
  eqcst(12)=6.*t
  eqcst(13)=2.*(t+2.*(h-s))+2.*(xi-nu)
  eqcst(14)=2.*(t-2.*s+h+p)+2.*(xi-nu)
  eqcst(15)=t+2.*s+h-90.-2.*xi-nu
  eqcst(16)=2.*t-s+p+180.+2.*(xi-nu)
  eqcst(17)=t
  i=di*pi180
  pc=dpc*pi180
  top=(5.*cos(i)-1.)*sin(pc)
  bottom=(7.*cos(i)+1.)*cos(pc)
  q=arctan(top,bottom,1)
  eqcst(18)=t-s+h-90.+xi-nu+q
  eqcst(19)=t+s+h-p-90.-nu
  eqcst(20)=s-p
  eqcst(21)=2.*h
  eqcst(22)=h
  eqcst(23)=2.*(s-h)
  eqcst(24)=2.*s-2.*xi
  eqcst(25)=t+3.*(h-s)-p+90.+2.*xi-nu
  eqcst(26)=t-3.*s+h+p+90.+2.*xi-nu
  eqcst(27)=2.*t-h+p1
  eqcst(28)=2.*t+h-p1+180.
  eqcst(29)=t-4.*s+h+2.*p+90.+2.*xi-nu
  eqcst(30)=t-h+90.
  eqcst(31)=2.*(t+s-h)+2.*(nu-xi)
  eqcst(32)=3.*(t-s+h)+3.*(xi-nu)
  r=sin(2.*pc)/((1./6.)*(1./tan(.5*i))**2-cos(2.*pc))
  r=atan(r)/pi180
  eqcst(33)=2.*(t+h)-s-p+180.+2.*(xi-nu)-r
  eqcst(34)=3.*(t+h)-4.*s+90.+4.*(xi-nu)+nup
  eqcst(35)=2.*(t+h)-2.*nup2
  eqcst(36)=8.*(t-s+h)+8.*(xi-nu)
  eqcst(37)=2.*(2.*t-s+h)+2.*(xi-nu)
  do 1 ih=1,37
    1 eqcst(ih)=angle(eqcst(ih))
end subroutine gterms

subroutine orbit(yr,dayj,hr)

!-- determination of primary and secondary orbital functions

!-- the equations programmed here are not represented by equations in
!   schureman.  the coding in this routine derives from a program by
!   the national oceanic and atmospheric administration (noaa).
!   however, tabular values of the orbital functions can be found in
!   table 1 of schureman.

!---------------------------------------------------------------------
  real i,n,nu,nup,nup2
  common /orbitf/ds,dp,dh,dp1,dn,di,dnu,dxi,dnup,dnup2,dpc

  pi180=3.14159265/180.
  x=aint((yr-1901.)/4.)
  dyr=yr-1900.
  dday=dayj+x-1.
!-- dn is the moon's node (capital n, table 1, schureman)
  dn=259.1560564-19.328185764*dyr-.0529539336*dday-.0022064139*hr
  dn=angle(dn)
  n=dn*pi180
!-- dp is the lunar perigee (small p, table 1)
  dp=334.3837214+40.66246584*dyr+.111404016*dday+.004641834*hr
  dp=angle(dp)
  p=dp*pi180
  i=acos(.9136949-.0356926*cos(n))
  di=angle(i/pi180)
  nu=asin(.0897056*sin(n)/sin(i))
  dnu=nu/pi180
  xi=n-2.*atan(.64412*tan(n/2.))-nu
  dxi=xi/pi180
  dpc=angle(dp-dxi)
!-- dh is the mean longitude of the sun (small h, table 1)
  dh=280.1895014-.238724988*dyr+.9856473288*dday+.0410686387*hr
  dh=angle(dh)
!-- dp1 is the solar perigee (small p1, table 1)
  dp1=281.2208569+.01717836*dyr+.000047064*dday+.000001961*hr
  dp1=angle(dp1)
!-- ds is the mean longitude of the moon (small s, table 1)
  ds=277.0256206+129.38482032*dyr+13.176396768*dday+.549016532*hr
  ds=angle(ds)
  nup=atan(sin(nu)/(cos(nu)+.334766/sin(2.*i)))
  dnup=nup/pi180
  nup2=atan(sin(2.*nu)/(cos(2.*nu)+.0726184/sin(i)**2))/2.
  dnup2=nup2/pi180
end subroutine orbit

function angle(arg)
!
!*** this routine places an angle in 0-360 (+) format
!
  m=-ifix(arg/360.)
  angle=arg+float(m)*360.
  if(angle .lt. 0.) angle=angle+360.
end function angle

function arctan(top,bottom,key)
!** determine arctangent and place in correct quadrant
!   if key eq 0  no quadrant selection made
!   if key .ne. 0 proper quadrant is selected

  if(bottom .ne. 0.0) go to 4
  if(top) 2,9,3
2 arctan=270.
  return
3 arctan=90.
  return
4 arctan=atan(top/bottom)*57.2957795
  if(key.eq.0) return
  if(top) 5,5,7
5 if(bottom) 6,9,8
6 arctan=arctan+180.
  return
7 if(bottom) 6,3,10
8 arctan=arctan+360.
  return
9 arctan=0.
10 return
end function arctan

function dayjul(yr,xmonth,day)
!
!*** this routine computes the julian day (as a real variable)
!
  dimension dayt(12),days(12)
  data dayt/0.,31.,59.,90.,120.,151.,181.,212.,243.,273.,304.,334./
  data days(1),days(2) /0.,31./
  dinc=0.
  yrlp=mod((yr-1900.),4.)
  if(yrlp .eq. 0.) dinc=1.
  do 1 i=3,12
1 days(i)=dayt(i)+dinc
  dayjul=days(ifix(xmonth))+day
end function dayjul

max_stars_repo_name: rustychris/schism

max_issues_repo_name: rustychris/schism

max_forks_repo_name: rustychris/schism