text:
D Street House, Davis:

For 12 years, the Yolo Center for Families (http://thecenterforfamilies.org/) has been strengthening and supporting families and individuals in Yolo County by connecting them to community resources and services. Recently, in partnership with Davis Community Church, the Center opened up a resource center in central Davis, “D Street House” (441 D St.; 5307571261).  While many of the services offered at D Street House provide significant support for lowerincome families and individuals, the Center for Families also provides programs and connections to services that support individuals and families from all socioeconomic backgrounds.

The services and programs currently offered out of D St House cover four main areas: Family Support, Children’s Services, Health Services, and Access to Benefits. Examples of services/programs include case management, infant and toddler playgroups, parent education, child development screenings, Zumba classes, and food distribution. Please see https://www.facebook.com/dstreethouse for uptodate offerings.

Individuals and families can both give and receive at D Street House, so please consider serving as a volunteer.  Go to http://bit.ly/DstreetHouse for more information.

The Center for Families Sites:

Woodland Center (Main Office)
828 Court St., Woodland, CA 95695
5304067221

Woodland High  La Semilla
21 N. West St., Woodland, CA 95695
5306624678 x5422

Knights Landing Center
9586 Mill St., Knights Landing, CA 95645
5307351776 x102

West Sacramento Center
2455 W. Capitol Ave., West Sacramento, CA 95691
9163724252

D Street House
441 D St., Davis, CA 95616
5307571261

Marguerite Montgomery Elementary Center
1441 Danbury St., Davis, CA 95618
5307592159

ABOUT US
Our mission is to provide prevention and intervention services that empower individuals and families to achieve their goals and be active stewards of our community.

The Center for Families core cultural values:
Value and nurture family
Honor and celebrate culture
Invest in community
Be creative with your resources
Be passionate and have fun
Build open and honest relationships
Communication, collaboration, and partnership are key
Be welcoming and friendly
Be willing to teach and learn from one another
Dedicated to individuals and families in need

Our center provides a holistic, onestop access points for family support, child, youth development, and health resources.
We serve individuals of all ages.
Every child and youth has a safe, nurturing and stimulating environment in which to grow and develop to his or her potential.
Our goal is to improve the health, social, educational, and economic outcomes for Yolo County residents.
Every family is guided to be stable and selfreliant.
The Center for Families honors clients cultural practices and traditions.
We recognize our clients as equal partners whose strengths are appreciated, enhanced, and utilized.

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
      SUBROUTINE LLIKE(NX,X,INTS,NINT,DOUBLS,NDOUB,TMAT,NTMAT,IMAT,&
     &                 NIMAT,OMAT,NOMAT,NOBS,NSET,NMISST,FPEN,FPRIOR,&
     &                 EPSM,JOB,FUNK,INFO,VX0)

      INTEGER          NX,NINT,INTS(NINT),NDOUB,NTMAT,NIMAT,NOMAT,NSET,&
     &                 NOBS(NSET),NMISST,JOB,INFO,VX0
      DOUBLE PRECISION X(NX),DOUBLS(NDOUB),TMAT(NTMAT),IMAT(NIMAT),&
     &                 OMAT(NOMAT),EPSM,FUNK

      END

max_stars_repo_name: perNyfelt/renjinSamplesAndTests

max_issues_repo_name: perNyfelt/renjinSamplesAndTests

max_forks_repo_name: perNyfelt/renjinSamplesAndTests


text:
      PROGRAM CCPCIR
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)

      PARAMETER (M=30,N=30,LRWK=3500,LIWK=3500,LZDT=2000)
      PARAMETER (THETMN=25.,THETMX=125.,RHOMN=0.5,RHOMX=5.)

      REAL Z(M,N), RWRK(LRWK)
      INTEGER IWRK(LIWK)
C
C Generate some data to contour
C
      CALL MKDAT (2, Z, M, N)
C
C  Open GKS, open and activate a workstation.
C
      CALL GOPKS (IERRF, ISZDM)
      CALL GOPWK (IWKID, LUNIT, IWTYPE)
      CALL GACWK (IWKID)
C
C Turn off clipping.
C
      CALL GSCLIP (0)
C
C
C Tell CONPACK that the SET call has been done, force it to generate X
C coordinates that are longitudes and Y coordinates that are latitudes,
C turn on mapping to an EZMAP background. Define the out-of-range value
C (returned by MAPTRN for an unprojectable point).
C
      CALL SET (0.05,0.95,0.05,0.95,
     1     -RHOMX,RHOMX,-RHOMX,RHOMX,1)
      CALL CPSETI ('SET - DO SET-CALL FLAG',0)
      CALL CPSETR ('XC1 - X COORDINATE AT INDEX 1',RHOMN)
      CALL CPSETR ('XCM - X COORDINATE AT INDEX M',RHOMX)
      CALL CPSETR ('YC1 - Y COORDINATE AT INDEX 1',THETMN)
      CALL CPSETR ('YCN - Y COORDINATE AT INDEX N',THETMX)
      CALL CPSETI ('MAP - MAPPING FLAG',2)
C
C Initialize the drawing of the first contour plot.
C
      CALL CPRECT (Z,M,M,N,RWRK,LRWK,IWRK,LIWK)
C
C Draw Contours
C
      CALL CPCLDR(Z,RWRK,IWRK)
C
C Close frame
C
      CALL FRAME
C
C Deactivate and close workstation, close GKS.
C
      CALL GDAWK (IWKID)
      CALL GCLWK (IWKID)
      CALL GCLKS
      STOP
      END

      SUBROUTINE MKDAT(NCTFR,ZDAT,M,N)
C
C Create some data to contour
C
      REAL ZDAT(M,N)
C
C NCTFR is used to generate a random number, ZDAT is the data array
C
      CALL GGDINI (0.,1.,NCTFR,.9)
      ZMIN= 1.E36
      ZMAX=-1.E36
      DO 10 I=1,M
         RLON=.017453292519943*(-180.+360.*REAL(I-1)/REAL(M-1))
         DO 20 J=1,N
            RLAT=.017453292519943*(-90.+180.*REAL(J-1)/REAL(N-1))
            ZDAT(I,J)=GGDPNT(RLAT,RLON)+.5*COS(4.*RLAT)
            ZMIN=MIN(ZMIN,ZDAT(I,J))
            ZMAX=MAX(ZMAX,ZDAT(I,J))
 20      CONTINUE
 10   CONTINUE
      DO 30 I=1,M
         DO 40 J=1,N
            ZDAT(I,J)=((ZDAT(I,J)-ZMIN)/(ZMAX-ZMIN))*130.-50.
 40      CONTINUE
 30   CONTINUE

      RETURN
      END

ext: f

max_stars_repo_name: tenomoto/ncl

max_issues_repo_name: tenomoto/ncl

max_forks_repo_name: tenomoto/ncl


text:
      SUBROUTINE MTRIRG
C
C
C*****
C THIS ROUTINE COMPUTES THE   MASS    MATRIX FOR A AXI-SYMMETRIC RING
C WITH A TRIANGULAR CROSS SECTION
C*****
C
C
C                        ECPT FOR THE TRIANGULAR RING
C
C
C                                                      TYPE
C ECPT( 1) ELEMENT IDENTIFICATION                        I
C ECPT( 2) SCALAR INDEX NO. FOR GRID POINT A             I
C ECPT( 3) SCALAR INDEX NO. FOR GRID POINT B             I
C ECPT( 4) SCALAR INDEX NO. FOR GRID POINT C             I
C ECPT( 5) MATERIAL ORIENTATION ANGLE(DEGREES)           R
C ECPT( 6) MATERIAL IDENTIFICATION                       I
C ECPT( 7) COOR. SYS. ID. FOR GRID POINT A               I
C ECPT( 8) X-COOR. OF GRID POINT A (IN BASIC COOR.)      R
C ECPT( 9) Y-COOR. OF GRID POINT A (IN BASIC COOR.)      R
C ECPT(10) Z-COOR. OF GRID POINT A (IN BASIC COOR.)      R
C ECPT(11) COOR. SYS. ID. FOR GRID POINT B               I
C ECPT(12) X-COOR. OF GRID POINT B (IN BASIC COOR.)      R
C ECPT(13) Y-COOR. OF GRID POINT B (IN BASIC COOR.)      R
C ECPT(14) Z-COOR. OF GRID POINT B (IN BASIC COOR.)      R
C ECPT(15) COOR. SYS. ID. FOR GRID POINT C               I
C ECPT(16) X-COOR. OF GRID POINT C (IN BASIC COOR.)      R
C ECPT(17) Y-COOR. OF GRID POINT C (IN BASIC COOR.)      R
C ECPT(18) Z-COOR. OF GRID POINT C (IN BASIC COOR.)      R
C ECPT(19) EL. TEMPERATURE FOR MATERIAL PROPERTIES       R
C
C
      DOUBLE PRECISION CONSTD, D2PI
      DOUBLE PRECISION
     1                   D ,  GAMBQ,    R,   Z
     2,                            DELINT,   AK,  AKI
     3,                  AKT,  AM
      DOUBLE PRECISION   R1,  R2,  R3,  Z1,  Z2,  Z3,  ZMIN, DGAMA
     1,                  DR,  RH,  DZ,  ZH,  RA,  ZA,  AREA
     2,                  TWOPI,    RHOD,     DKI
C
      DIMENSION          IECPT(19)
      DIMENSION  AM(36)
C
      COMMON /CONDAD/ CONSTD(5)
      COMMON   /SMA2IO/
     1                   DUM1(10)
     2,                  IFMGG
     3,                  DUM2(25)
      COMMON   /SMA2CL/
     1                   DUM3(2)
     2,                  NPVT
     3,                  DUM4(7)
     4,                  LINK(10)           ,NOGO
      COMMON   /SMA2ET/
     1                   ECPT(19)
     2,                  DUM5(81)
      COMMON   /MATIN/
     1                   MATIDC             ,MATFLG
     2,                  ELTEMP             ,STRESS
     3,                  SINTH              ,COSTH
      COMMON   /MATOUT/
     1                   E(3)               ,ANU(3)
     2,                  RHO                ,G(3)
     3,                  ALF(3)             ,TZERO
      COMMON   /SMA2DP/
     1                   D(36) ,   GAMBQ(36),     R(3) ,    Z(3)
     2,                                      DELINT(8),     AK(36)
     3,                  AKI(36),  AKT(9)
     4,                  DGAMA,    ZMIN,     RHOD,     TWOPI
     5,                  DR,  RH,  DZ,  ZH,  RA,  ZA,  AREA
     8,                  IGP(3) , ICS(3) , SP(18)
     9,                  TEMPE
C
      EQUIVALENCE        (IECPT(1) , ECPT(1))
      EQUIVALENCE   (R(1),R1),     (R(2),R2),     (R(3),R3)
     1,             (Z(1),Z1),     (Z(2),Z2),     (Z(3),Z3)
      EQUIVALENCE  (AM(1) , AK(1))
      EQUIVALENCE ( CONSTD(2) , D2PI   )
C
C ----------------------------------------------------------------------
C
C STORE ECPT PARAMETERS IN LOCAL VARIABLES
C
      IDEL  = IECPT(1)
      IGP(1)= IECPT(2)
      IGP(2)= IECPT(3)
      IGP(3)= IECPT(4)
      MATID = IECPT(6)
      ICS(1)= IECPT(7)
      ICS(2)= IECPT(11)
      ICS(3)= IECPT(15)
      R(1)  = ECPT(8)
      D(1)  = ECPT(9)
      Z(1)  = ECPT(10)
      R(2)  = ECPT(12)
      D(2)  = ECPT(13)
      Z(2)  = ECPT(14)
      R(3)  = ECPT(16)
      D(3)  = ECPT(17)
      Z(3)  = ECPT(18)
      TEMPE = ECPT(19)
      DGAMA = ECPT(5)
C
C
C CHECK INTERNAL GRID POINTS FOR PIVOT POINT
C
      IPP = 0
      DO 100 I = 1,3
      IF (NPVT .EQ. IGP(I)) IPP = I
  100 CONTINUE
      IF (IPP .EQ. 0) CALL MESAGE (-30,34,IDEL)
C
C
C TEST THE VALIDITY OF THE GRID POINT COORDINATES
C
      IEROR1 = 0
      DO 200 I = 1,3
      IF (R(I).GT.0.0D0) GO TO 200
      IF (IEROR1.NE.0) GO TO 200
      CALL MESAGE (30, 211, IDEL)
      IEROR1 = 1
  200 CONTINUE
      IEROR2 = 0
      DO 210 I = 1, 3
      IF (D(I).EQ.0.0D0) GO TO 210
      IF (IEROR2.NE.0) GO TO 210
      CALL MESAGE (30, 212, IDEL)
      IEROR2 = 1
  210 CONTINUE
C
C
C COMPUTE THE ELEMENT COORDINATES
C
      IF (IEROR1.EQ.0.AND.IEROR2.EQ.0) GO TO 220
      NOGO = 2
      RETURN
  220 ZMIN = DMIN1 (Z1, Z2, Z3)
      Z1 = Z1 - ZMIN
      Z2 = Z2 - ZMIN
      Z3 = Z3 - ZMIN
C
C
C
C FORM THE TRANSFORMATION MATRIX (6X6) FROM FIELD COORDINATES TO GRID
C POINT DEGREES OF FREEDOM
C
      DO 300 I = 1,36
      GAMBQ(I) = 0.0D0
  300 CONTINUE
      GAMBQ( 1) = 1.0D0
      GAMBQ( 2) = R1
      GAMBQ( 3) = Z1
      GAMBQ(10) = 1.0D0
      GAMBQ(11) = R1
      GAMBQ(12) = Z1
      GAMBQ(13) = 1.0D0
      GAMBQ(14) = R2
      GAMBQ(15) = Z2
      GAMBQ(22) = 1.0D0
      GAMBQ(23) = R2
      GAMBQ(24) = Z2
      GAMBQ(25) = 1.0D0
      GAMBQ(26) = R3
      GAMBQ(27) = Z3
      GAMBQ(34) = 1.0D0
      GAMBQ(35) = R3
      GAMBQ(36) = Z3
C
C
C     NONEED TO COMPUTE DETERMINANT SINCE IT IS NOT USED SUBSEQUENTLY.
      ISING = -1
      CALL INVERD (6, GAMBQ(1),6 , D(10), 0, D(11) , ISING , SP)
C
      IF (ISING.EQ.2) GO TO 920
C
C
C
C CALCULATE THE INTEGRAL VALUES IN ARRAY DELINT WHERE THE ORDER IS
C INDICATED BY THE FOLLOWING TABLE
C
C              DELINT( 1) - ( 1,0)
C              DELINT( 2) - ( 1,1)
C              DELINT( 3) - ( 1,2)
C              DELINT( 4) - ( 2,0)
C              DELINT( 5) - ( 2,1)
C              DELINT( 6) - ( 0,2)
C              DELINT( 7) - ( 3,0)
C              DELINT( 8) - (-1,2)
C
C
C TEST FOR RELATIVE SMALL AREA OF INTEGRATION
C AND IF AREA IS SMALL THEN APPROXIMATE INTEGRALS
C
      DR = DMAX1 ( DABS(R1-R2) , DABS(R2-R3) , DABS(R3-R1) )
      RH = DMIN1 ( R1 , R2 , R3 ) / 10.0D0
      DZ = DMAX1 ( DABS(Z1-Z2) , DABS(Z2-Z3) , DABS(Z3-Z1) )
      ZH = DMIN1 ( Z1 , Z2 , Z3 ) / 10.0D0
      RA = (R1 + R2 + R3) / 3.0D0
      ZA = (Z1 + Z2 + Z3) / 3.0D0
      AREA =(R1*(Z2-Z3) + R2*(Z3-Z1) + R3*(Z1-Z2)) / 2.0D0
      KODE = 0
      IF (DABS( (R2-R1)/R2 ) .LT. 1.0D-5) KODE = 1
      IF ( DR .LE. RH  .OR.  DZ .LE. ZH ) KODE = -1
C
C
  310 CONTINUE
      I1 = 0
      DO 400 I = 1,3
      IP = I
      DO 350 J = 1,3
      IQ = J - 1
      IF (IP.EQ.2 .AND. IQ.EQ.2) IP = 0
      IF (IP.EQ.3 .AND. IQ.EQ.2) IP = -1
      IF (IP.EQ.3 .AND. IQ.EQ.1) GO TO 350
      I1 = I1 + 1
      IF (KODE) 320,330,340
  320 DELINT(I1) =((RA) ** IP)*((ZA) ** IQ) * AREA
      GO TO 350
  330 DELINT(I1) = DKI(1,3,1,2,1,3,IP,IQ,R,Z)
     1           + DKI(3,2,1,2,3,2,IP,IQ,R,Z)
      GO TO 350
  340 CONTINUE
      DELINT(I1) = DKI(1,3,3,2,1,3,IP,IQ,R,Z)
  350 CONTINUE
  400 CONTINUE
C
C
C TEST FOR EXCESSIVE ROUND-OFF ERROR IN INTEGRAL CALCULATIONS
C AND IF IT EXIST APPROXIMATE INTEGRALS
C
      IF (KODE .LT. 0) GO TO 500
      DO 450 I = 1,8
      IF (DELINT(I) .LT. 0.0D0) GO TO 475
  450 CONTINUE
      IF (DELINT(3) .LE. DELINT(6)) GO TO 475
      IF (DELINT(8) .GE. DELINT(3)) GO TO 475
      IF (DELINT(8) .GT. DELINT(6)) GO TO 475
      GO TO 500
  475 CONTINUE
      KODE = -1
      GO TO 310
  500 CONTINUE
C
C
C
C LOCATE THE MATERIAL PROPERTIES IN THE MAT1 OR MAT3 TABLE
C
      MATIDC = MATID
      MATFLG = 7
      ELTEMP = TEMPE
      CALL  MAT (IDEL)
C
C
C SET MATERIAL PROPERTIES IN DOUBLE PRECISION VARIABLES
C
      RHOD = RHO
C
C GENERATE THE CONSISTENT MASS MATRIX IN FIELD COORDINATES
C
      DO 600 I = 1,36
      AM(I) = 0.0D0
  600 CONTINUE
      TWOPI = D2PI * RHOD
      AM( 1)= TWOPI *       DELINT(1)
      AM( 2)= TWOPI *       DELINT(4)
      AM( 3)= TWOPI *       DELINT(2)
      AM( 7)= AM( 2)
      AM( 8)= TWOPI *       DELINT(7)
      AM( 9)= TWOPI *       DELINT(5)
      AM(13)= AM( 3)
      AM(14)= AM( 9)
      AM(15)= TWOPI *       DELINT(3)
      AM(22)= AM( 1)
      AM(23)= AM( 2)
      AM(24)= AM( 3)
      AM(28)= AM(23)
      AM(29)= AM( 8)
      AM(30)= AM( 9)
      AM(34)= AM(24)
      AM(35)= AM(30)
      AM(36)= AM(15)
C
C TRANSFORM THE ELEMENT   MASS    MATRIX FROM FIELD COORDINATES
C TO GRID POINT DEGREES OF FREEDOM
C
      CALL GMMATD (GAMBQ , 6, 6, 1, AK , 6, 6, 0, D )
      CALL GMMATD (D  , 6, 6, 0, GAMBQ , 6, 6, 0, AK)
C
C
C
C ZERO OUT THE (6X6) MATRIX USED AS INPUT TO THE INSERTION ROUTINE
C
      DO 700 I = 1,36
      AKI(I) = 0.0D0
  700 CONTINUE
C
C
C LOCATE THE TRANSFORMATION MATRICES FOR THE THREE GRID POINTS
C
      DO 800 I = 1,3
      IF (ICS(I) .EQ. 0) GO TO 800
      K = 9 * (I-1) + 1
      CALL TRANSD (ICS(I) , D(K))
  800 CONTINUE
C
C
C
C START THE LOOP FOR INSERTION OF THE THREE (6X6) MATRICES
C INTO THE MASTER MASS MATRIX
C
      IR1  = 2 * IPP - 1
      IAPP = 9 * (IPP-1) + 1
      DO 900 I = 1,3
C
C PLACE THE APPROIATE (2X2) SUBMATRIX OF THE MASS MATRIX
C IN A (3X3) MATRIX FOR TRANSFORMATION
C
      IC1 = 2 * I - 1
      IRC = (IR1 - 1) * 6 + IC1
      AKT(1) = AK(IRC)
      AKT(2) = 0.0D0
      AKT(3) = AK(IRC+1)
      AKT(4) = 0.0D0
      AKT(5) = 0.0D0
      AKT(6) = 0.0D0
      AKT(7) = AK(IRC+6)
      AKT(8) = 0.0D0
      AKT(9) = AK(IRC+7)
C
C TRANSFORM THE (3X3) MASS MATRIX
C
      IF (ICS(IPP) .EQ. 0) GO TO 820
      CALL GMMATD (D(IAPP) , 3, 3, 1, AKT(1) , 3, 3, 0, D(28) )
      DO 810 J = 1,9
      AKT(J) = D(J+27)
  810 CONTINUE
  820 CONTINUE
      IF (ICS(I) .EQ. 0) GO TO 840
      IAI = 9 * (I - 1) + 1
      CALL GMMATD (AKT(1) , 3, 3, 0, D(IAI) , 3, 3, 0, D(28) )
      DO 830 J = 1,9
      AKT(J) = D(J+27)
  830 CONTINUE
  840 CONTINUE
C
C PLACE THE TRANSFORMED (3X3) MATRIX INTO A (6X6) MATRIX FOR
C THE INSERTION ROUTINE
C
      J = 0
      DO 850 J1 = 1,18,6
      DO 850 J2 = 1,3
      J = J + 1
      K = J1 + J2 - 1
      AKI(K) = AKT(J)
  850 CONTINUE
C
C CALL THE INSERTION ROUTINE
C
      CALL SMA2B (AKI(1) , IGP(I), -1, IFMGG, 0.0D0)
  900 CONTINUE
      RETURN
C
C  SET FLAG FOR FATAL ERROR WHILE ALLOWING ERROR MESSAGES TO ACCUMULATE
C
  915 NOGO=1
      RETURN
  920 CALL MESAGE(30,37,IDEL)
      GO TO 915
C
      END


text:
MODULE USR_SMlib_Matrix_Arithmetic
  !
  !!  This module is part of SMLIB v. 1.1.  It collects all the matrices
  !!  modules.  It also contains a routine for conversion of an MSR matrix
  !!  to a CSR matrix by an overloaded assignment(=) operator.
  !
  !! ------------------------------------------------------------------------
  !
  !!  Copyright (C) 1996, 1997, 1998,  Ernst A. Meese
  !!    Refer to the file copyright.doc for details and important disclaimer.
  !
  use USR_SMlib_CSR        !!((17-B-USR_SMlib_CSR.f90))
  use USR_SMlib_MSR        !!((17-B-USR_SMlib_MSR.f90))
  use USR_SMlib_LU_CSR_MSR !!((18-B-USR_SMlib_LU_CSR_MSR.f90))
  use USR_SMlib_Band_LU    !!((17-B-USR_SMlib_Band_LU.f90))

  INTERFACE ASSIGNMENT(=)
     MODULE PROCEDURE CSR_from_MSR
  END INTERFACE

CONTAINS

  SUBROUTINE CSR_from_MSR( C, M )
    TYPE(CSR), INTENT(INOUT) :: C
    TYPE(MSR), INTENT(IN)    :: M
    INTEGER :: I

    IF ( ALLOCATED_MATRIX(C) ) CALL DEALLOCATE_MATRIX(C)
    CALL ALLOCATE_MATRIX( C, M % N, M % JA ( M % N + 1 ) )

    DO I = 1, M % N
       Call setRow ( C, I, (/ I,         M % JA (M % JA (I) : M % JA(I+1)-1)  /), &
                           (/ M % A (I), M % A  (M % JA (I) : M % JA(I+1)-1)  /) )
    END DO

  END SUBROUTINE CSR_from_MSR

END MODULE USR_SMlib_Matrix_Arithmetic

max_stars_repo_name: wawiesel/tapack

max_issues_repo_name: wawiesel/tapack

max_forks_repo_name: wawiesel/tapack


text:
c read a few eigenvectors, a spectrum, and try to fit z

      program fitz

      include 'specarray.h'

      parameter(NMAXVECTS=12)
      parameter(NMAXSTEP=2*NLOGWMAX)
c scale sky by estimating read and photon noise ?
      parameter (IFSCALENOISE=1)
c plot intermediate steps?
      parameter (IFPLOTALL=0)
c read spectrum names from a file?
      parameter (IFFILE=1)
c overplot the actual z, when known?
      parameter (IFREALZ=1)
c max number of minima to find and zestimates to return
      parameter (NMAXEST=12)

c      real evects(NMAXVECTS,NLOGWMAX)
      integer koffset,nvmax,nwlmax
      real z
      common /vectorfit/ koffset,z,nvmax,nwlmax,
     $     evects(NMAXVECTS,NLOGWMAX)

      real spec1(NWAVEMAX),espec1(NWAVEMAX)
      real spec2(NLOGWMAX),espec2(NLOGWMAX)
      real spec3(NLOGWMAX),espec3(NLOGWMAX)
      real tmpspec(NLOGWMAX),wavelog(NLOGWMAX),wavelin(NWAVEMAX)
      integer ifit(NLOGWMAX)
      real wfit(NLOGWMAX),sfit(NLOGWMAX),efit(NLOGWMAX)
      real rifit(NLOGWMAX),fitvect(NLOGWMAX)
      real waverest1(NLOGWMAX)
      real wrestlin(NWAVEMAX)
      real zp1log(NMAXSTEP),zarr(NMAXSTEP)
      real rnparr(NMAXSTEP)
      integer indkarr(NMAXEST)
      real diffarr(NMAXEST),zestarr(NMAXEST)
c      real xpl(2),ypl(2)
c      real zout(NLOGWMAX),zchisq(NLOGWMAX)
      character sname*60,esname*60,vname*60
      character answ*3,tlabel*60,xlabel*60
      character fsname*60,fename*60,fzname*60
      integer isize(7)

c  stuff for svdfit
      real uu(NLOGWMAX,NMAXVECTS)
      real vv(NMAXVECTS,NMAXVECTS),ww(NMAXVECTS)
      real acoeff(NMAXVECTS),evals(NMAXVECTS)
      real chifit(NMAXSTEP),rchifit(NMAXSTEP)

      external funcs

      include 'pcredshift.h'

c      call pgbeg(0,'?',2,2)
      call pgbeg(0,'?',1,1)
      call pgscf(2)
      call pgsch(1.3)

      write(*,'("Full continuum fit subtraction [1,y-yes]? ",$)')
      read(*,'(a1)') answ
      if (answ(1:1) .eq. '0' .or. answ(1:1) .eq. 'n' .or.
     $     answ(1:1) .eq. 'N') then
         ifcontsub = 0
      else
         ifcontsub = 1
      end if
      write(*,'("Do mean subtraction [1,y-yes]? ",$)')
      read(*,'(a1)') answ
      if (answ(1:1) .eq. '0' .or. answ(1:1) .eq. 'n' .or.
     $     answ(1:1) .eq. 'N') then
         ifmeansub = 0
      else
         ifmeansub = 1
      end if
 100  write(*,
     $     '("Diameter for median smoothing, pixels [0,1=none]: ",$)')
      read(*,'(a3)') answ
      if (answ(1:3) .eq. '   ') then
         ndiamed = 0
      else
         read(answ,*,err=100) ndiamed
      end if

c      write(*,
c     $     '("Restwave spectrum region to use in PCA, min, max: ",$)')
c      read(*,*) pcawmin,pcawmax
c      pwminlog = log10(pcawmin)
c      pwmaxlog = log10(pcawmax)

      nvmax=NMAXVECTS
      nwlmax=NLOGWMAX

c open file and read eigenvectors. nvects is returned as
c the # of vectors to use, nw is the actual # of wavelength pixels
      call readevects(evects,nvmax,nwlmax,waverest1,nv,nw)
      nwlog = nw

c figure out w0 and dw for the eigenvectors.  This is in log w
c There is no check yet to make sure it's evenly spaced.
      dwrlog = (waverest1(nw) - waverest1(1)) / (nwlog-1.)
      w0rlog = waverest1(1) - dwrlog
      do i=1,nwlog
         wrestlin(i) = 10**waverest1(i)
      end do
c      write(*,*) "w0rlog, dwrlog = ", w0rlog,dwrlog

c plot the eigenvectors
      call pgsubp(2,2)
      do i=1,nv
         do j=1,nwlog
            tmpspec(j) = evects(i,j)
         end do
         call showspec(nwlog,waverest1,tmpspec)
         write(tlabel,'(a,1x,i3)') "eigenvector",i-1
         call pglabel("log wavelength"," ",tlabel)
      end do
c      call pgsubp(1,1)
      call pgsubp(2,2)

 200  continue
      if (IFFILE .eq. 1) then
         write(*,'("File with spectrum names: ",$)')
         read(*,'(a)') fsname
         open(10,file=fsname,status='old',err=200)
         write(*,'("File with sky/rms names: ",$)')
         read(*,'(a)') fename
         open(11,file=fename,status='old',err=200)
         if (IFREALZ .eq. 1) then
            write(*,'("File with actual zs for plot [none]: ",$)')
            read(*,'(a)') fzname
            if (fzname(1:3) .ne. '   ') then
               open(12,file=fzname,status='old',err=202)
               ifzfile=1
            else
c  no actual z file
               ifzfile=0
            end if
 202        continue
         end if
      end if

      open(2,file='fitz.out',status='unknown')
      open(3,file='fitz.out1',status='unknown')
      open(4,file='fitz.out2',status='unknown')

 220  continue
c open files and read spec and error, with wavelength calib.
      if (IFFILE .eq. 1) then
         read(10,'(a)',err=666,end=666) sname
         read(11,'(a)',err=666,end=666) esname
         zreal=1.0e6
         if (IFREALZ .eq. 1 .and. ifzfile .eq. 1) then
            read(12,*,err=222,end=222) zreal
         end if
 222     continue
      else
         write(*,'("fits file with spectrum [quit]: ",$)')
         read(*,'(a)') sname
         if (sname(1:3) .eq. '   ') go to 666
         write(*,'("fits file with sky/rms [quit]: ",$)')
         read(*,'(a)') esname
         if (esname(1:3) .eq. '   ') go to 666
         if (IFREALZ .eq. 1) then
            write(*,'("actual z for plotting [none]: ",$)')
            read(*,'(a)') fzname
            zreal=1.0e6
            if (fzname(1:3) .ne. '   ') then
               read(fzname,*,err=232) zreal
            end if
 232        continue
         end if
      end if

      call imopen(sname,1,imgs,ier)
      if (ier .ne. 0) go to 220
      call imgsiz(imgs,isize,idim,itype,ier)
c trim a buffer at end to get rid of the region where Marc encoded
c the slit number
      ibuff2 = 25
      nwspec = isize(1) - ibuff2
      call imopen(esname,1,imge,ier)
c      call imgsiz(imge,isize,idim,itype,ier)
      if (ier .ne. 0) go to 220

      call imgl1r(imgs,spec1,ier)
      call imgl1r(imge,espec1,ier)

c  find wavelength of pixel 0, and dw/dpix
      call imgkwr(imgs,'CRPIX1',refpix,ier)
c      if (ier .ne. 0) refpix = badset
      if (ier .ne. 0) then
         ier=0
         refpix = 0.
      end if
      call imgkwr(imgs,'CRVAL1',refw,ier)
      if (ier .ne. 0) refw = badset
      call imgkwr(imgs,'CDELT1',dwsp,ier)
      if (ier .ne. 0 .or. abs(dwsp) .lt. 1.e-3) then
         ier=0
         call imgkwr(imgs,'CD1_1',dwsp,ier)
         if (ier .ne. 0) dwsp = badset
      end if
      if (refpix .gt. bad .and. refw .gt. bad .and. dwsp .gt. bad)
     $     then
         w0sp = refw - refpix*dwsp
c         dws = dwsp
      else
c  we should really do something here
c         w0s = badset
c         dws = badset
         write(*,
     $     '("Couldnt get w0 and dw for ",a," enter w0,dw: ",$)') sname
         read(*,*) w0sp,dwsp
      end if
      write(*,*) "w0sp, dwsp = ", w0sp,dwsp

      call imclos(imgs,ier)
      call imclos(imge,ier)

c zero out anything beyond the actual spectrum
      do i=nwspec+1,nw
         spec1(i) = 0.0
         espec1(i) = 0.0
      end do

c try to clean out bad values
      badmax = 5000.
      do i=1,nw
         if(spec1(i) .lt. bad .or. spec1(i) .gt. badmax)
     $        spec1(i) = badset
         if(espec1(i) .lt. bad .or. espec1(i) .gt. badmax)
     $        espec1(i) = badset
      end do

c Figure out the ref wavelength for the log rebinning so that
c it falls on the wl scale of the eigenvectors.
      tmp = log10(w0sp)
      itmp = int( (tmp-w0rlog)/dwrlog )
      w0log = w0rlog + dwrlog*itmp
c k0offset is the offset in integer index of log wavelength
c from the beginning of the eigenvectors to the spectrum to be fit.
      k0offset = itmp
      dwlog = dwrlog
c      write(*,*) "w0log, dwlog = ", w0log,dwlog

c convert sky spectrum into std.dev spectrum.  This might be better
c done after smoothing and log rebinning?
c Placeholder assumptions: 2x30 min exposures, 5 pixel diam
c extraction window.
      if (IFSCALENOISE .eq. 1) then
         nexp = 2
         exptime = 30.*60.
c         dpix = 5.
         dpix = 1.
         call scalenoise(espec1,nwspec,nexp,exptime,dpix,espec2)
      else
c or if the sky spectrum is actually a rms/per pixel spectrum
c we could do nothing, or scale down by a factor of
c sqrt(#pixels).  #pixels is most likely 7.
         do i=1,nwspec
            espec2(i) = espec1(i) / sqrt(7.0)
         end do
      end if

      do i=1,nw
         wavelin(i) = w0sp+i*dwsp
         wavelog(i) = log10(wavelin(i))
      end do
      if (IFPLOTALL .ne. 0) then
         call showspec(nwspec,wavelin,spec1)
         call pgqci(indexc)
         call pgsci(3)
         call pgline(nwspec,wavelin,espec2)
         call pgsci(indexc)
         call pglabel("wavelength","counts","spectrum and error")
      end if

c blank out?  blankoutsky is supposed to run on a 2-d array

      call blankoutsky(spec1,1,1,nwspec,wavelin)

cc find region of good data
c      call findends(nwspec,spec1,imin,imax)
c      imingood = imin
c      nwspecgood = imax-imin+1

c median smooth
      if (ndiamed .gt. 1) then
         call medsmooth(nwspec,spec1,ndiamed,spec2)
c         call medsmooth(nwspecgood,spec1(imingood),ndiamed,
c     $        spec2(imingood))
c attempt to compensate the errors.  I divided ndiamed by 2 as
c a hack since adjacent pixels are not independent (assume the #
C of indep measurements is about pixels/2 if well sampled)
         tmp = sqrt(max(ndiamed/2.,1.))
         do i=1,nwspec
            espec2(i) = espec2(i) / tmp
         end do
      else
         do i=1,nwspec
            spec2(i) = spec1(i)
            espec2(i) = espec2(i)
         end do
      end if

      if (IFPLOTALL .ne. 0) then
         call showspec(nwspec,wavelin,spec2)
         call pglabel("wavelength","counts","median smoothed")
         call pgqci(indexc)
         call pgsci(3)
         call pgline(nwspec,wavelin,espec2)
         call pgsci(indexc)
      end if

c log rebin both.  log rebinning the std dev the same way is a
c total hack, thus it is better done on the variance

      call logrebin(spec2,nwspec,w0sp,dwsp,nwlog,w0log,dwlog,spec3)
c      call logrebin(spec2(imingood),nwspecgood,w0sp,dwsp,
c     $     nwlog,w0log,dwlog,spec3)
c convert std dev to variance
      do i=1,nwlog
         espec3(i) = espec2(i)**2
      end do
      call logrebin(espec3,nwspec,w0sp,dwsp,nwlog,w0log,dwlog,espec2)
c      call logrebin(espec3(imingood),nwspecgood,w0sp,dwsp,
c     $     nwlog,w0log,dwlog,espec2)
c convert back
      do i=1,nwlog
         espec3(i) = sqrt(max(espec2(i),1.e-2))
      end do

      if (IFPLOTALL .ne. 0) then
         call showspec(nwlog,wavelog,spec3)
         call pglabel("log wavelength","counts","log rebinned")
         call pgqci(indexc)
         call pgsci(3)
         call pgline(nwlog,wavelog,espec3)
         call pgsci(indexc)
      end if

c find&clean ends
      call findends(nwlog,spec3,imin,imax)
      call cleanends(nwlog,spec3,imin,imax)
      write(*,*) "imin, imax = ", imin,imax

c continuum subtract
      if (ifcontsub .ne. 0) then
         contwrad = 100.*dwlog
         call contsubmed(spec3,nwlog,dwlog,contwrad)
      else
         call contsubconst(spec3,nwlog)
      end if

      if (IFPLOTALL .ne. 0) then
         call showspec(nwlog,wavelog,spec3)
         call pgqci(indexc)
         call pgsci(3)
         call pgline(nwlog,wavelog,espec3)
         call pgsci(indexc)
         call pglabel("log wavelength","counts","continuum subtracted")
      end if

c copy wave and spec to temp array, cleaning out bad points
      npfit=0
      do i=1,nwlog
         if(spec3(i) .gt. bad .and. espec3(i) .gt. bad) then
            npfit=npfit+1
            ifit(npfit) = i
            rifit(npfit) = real(i)
            wfit(npfit) = w0log + i*dwlog
            sfit(npfit) = spec3(i)
            efit(npfit) = espec3(i)
         end if
      end do
      write (*,*) npfit, " points to fit"

c      if (IFPLOTALL .ne. 0) then
         call showspec(npfit,wfit,sfit)
         call pgqci(indexc)
         call pgsci(3)
         call pgline(npfit,wfit,efit)
         call pgsci(indexc)
         call pglabel("log wavelength","counts",
     $        "spectrum and error to fit, "//sname)
c      end if

c What are the good data regions?  The spectrum was
c good from imin to imax (before cleanends).  Eigenvector 1
c (actually the mean) is good from jmin to jmax.
      do j=1,nwlog
         spec2(j) = evects(1,j)
      end do
      call findends(nwlog,spec2,jmin,jmax)
c      write(*,*) 'imin,imax,jmin,jmax, k0off: ',imin,imax,jmin,jmax,
c     $     k0offset

c loop over steps in log w.  the spectrum, indexed by i,
c is offset to the red of the eigenvectors, indexed by j,
c by k steps: j+k = i, k = i-j
c  and note that the two scales are offset by k0offset, so when
c  z=0, i=j-k0offset.  Typically the spectrum starts redder
c  than the eigenv. so k0offset>0.  If we started at
c  z=0, k would start at -koffset.
c log w = log wrest + log(1+z), so  k*dwlog = log(1+z)
c Start with an offset of z=-0.01
      kmin = int(log10(1-0.01)/dwlog) - k0offset
c For pos. redshift, i+k0offset>j (i is position in spectrum, j in evect)
c Go to the point where only 10 pts overlap, which doesn't
c depend on k0offset.
      kmax = imax-10-jmin
      nk=kmax-kmin+1
c      write(*,*) "kmin, kmax: ",kmin,kmax
      tmp1 = (kmin+k0offset)*dwlog
      tmp2 = (kmax+k0offset)*dwlog
c      write(*,*) "log(1+z): ",tmp1,tmp2,"  z: ",
c     $     10**tmp1-1.0,10**tmp2-1.0

c      open(2,file='fitz.out1',status='unknown')

      do k=kmin,kmax
c  kindex is 1 when k=kmin
         kindex = k-kmin+1
c  for passing k to the funcs subroutine
         koffset=k
         zp1log(kindex) = (k+k0offset)*dwlog
         zarr(kindex) = 10**zp1log(kindex) - 1.0
         z = zarr(kindex)

c  we should only use the data that overlaps the eigenvectors,
c  the eigenvectors extend from jmin to jmax, and i=j+k.
c  But, since we cleaned out bad points, the spectrum to fit
c  is no longer evenly spaced in dwlog, so ifitmax is not
c  trivial to calculate.  Bummer!
c  that is, what I've been calling "i" is the index of spec3(i),
c  and the contents of ifit() and rifit(), but it is not the
c  index of ifit
         iminorig = jmin+k
         imaxorig = jmax+k
c now find the corresponding index of ifit - i.e.
c we want iminfit where ifit(iminfit) >= iminorig.  Confused yet?
         call findindex(npfit,ifit,iminorig,imaxorig,iminfit,imaxfit)

         npfittmp = imaxfit-iminfit+1
         rnparr(kindex) = real(npfittmp)
c         write(*,*) "iminfit,imaxfit: ",iminfit,imaxfit,npfittmp

c now we gotta only use the matching points in the eigenvectors,
c which is a PITA.  Actually we've solved this through the way
c that funcs() works.

c do fit using svdfit or something like it.
c Each of the functions will be an eigenvector

c         call bigsvdfit(wfit,sfit,efit,npfit,acoeff,nv,uu,vv,ww,
c     $        nwlmax,nvmax,chisq1,funcs)
c         call bigsvdfit(rifit,sfit,efit,npfit,acoeff,nv,uu,vv,ww,
c     $        nwlmax,nvmax,chisq1,funcs)
         call bigsvdfit(rifit(iminfit),sfit(iminfit),efit(iminfit),
     $        npfittmp,acoeff,nv,uu,vv,ww,nwlmax,nvmax,chisq1,funcs)
         chifit(kindex) = chisq1
c to get rid of division by zero when there were no points
c to fit.
         rchifit(kindex) = chisq1/max(real(npfittmp),1.e-4)
c         write(2,*) k,kindex,npfittmp,zp1log(kindex),zarr(kindex),
c     $        chifit(kindex),rchifit(kindex)

      end do
c      close(2)

c find the absolute minimum in chi-squared.  there are probably
c better ways to do this in the long run
c      call findabsmin(nk,chifit,indkmin,chimin)
c      call findabsmin(nk,rchifit,indkrmin,rchimin)
c  try finding the point with the biggest drop below "continuum"
c  i.e. deepest local minimum
      call findlocmin(nk,chifit,indkmin,chimin,chiminavg)
      call findlocmin(nk,rchifit,indkrmin,rchimin,rchiminavg)

c find the n deepest local minima in chi squared
      nfind=5
      call findmultmin(nk,chifit,nfind,indkarr,diffarr)

      zmin=zarr(indkmin)
      zrmin=zarr(indkrmin)
      do i=1,nfind
         zestarr(i) = zarr(indkarr(i))
      end do

      write(*,*) "chisq min at ",indkmin,zmin,chimin,chiminavg
      write(*,*) "reduced chisq min at ",indkrmin,zrmin,rchimin,
     $     rchiminavg

c write stuff to output files.
c 2 - fitz.out gets # of z-estimates, z-estimates, spec name
c 3 - fitz.out1 gets specname(truncated), z(deepest chisq min) ,
c          z(deepest rchisq min), chi-min, chi-"continuum",
c          rchi-min, rchi-"continuum"
c 4 - fitz.out2 gets for each z-est: k-index, z-est, depth in chisq

c      write(3,'(a,f7.4,3x,f7.4)') sname,zmin,zrmin
      write(3,'(a25,2(3x,f7.4),3x,4(1x,1pe10.3))')
     $     sname,zmin,zrmin,chimin,chiminavg,rchimin,rchiminavg
      write(2,'(i2,$)') nfind
      do i=1,nfind
         write(2,'(2x,f7.4$)') zestarr(i)
         write(4,'(i5,2x,f7.4,2x,f9.1,$)') indkarr(i),zestarr(i),
     $        diffarr(i)
      end do
      write(2,'(2x,a25)') sname
      write(4,'(2x,a25)') sname

c plot
c      call showspec(nk,zp1log,chifit)
c      call plotz(log10(1.+zreal),log10(1.+zmin),log10(1+zrmin))
c      write(xlabel,1010) "log(1+z)",zreal,zmin,zrmin
c      call pglabel(xlabel,"chi-squared",sname)
      call showspec(nk,zarr,chifit)
      call plotz(zreal,zmin,zrmin)
      write(xlabel,1010) "z",zreal,zmin,zrmin
      call pglabel(xlabel,"chi-squared",sname)
      call showspec(nk,zarr,rchifit)
      call plotz(zreal,zmin,zrmin)
      write(xlabel,1010) "z",zreal,zmin,zrmin
      call pglabel(xlabel,"reduced chi-squared",sname)
 1010 format(a,",  z=",f7.4,"  zest1=",f7.4,"  zest2=",f7.4)
c      call showspec(nk,zarr,rnparr)
c      call pglabel("z","number of points in fit",sname)

c  calculate the fit at chimin.  Redo the fit to get the coeffs
c  for the best-fit spectrum
      koffset=indkmin+kmin-1
      call findindex(npfit,ifit,jmin+koffset,jmax+koffset,
     $     iminfit,imaxfit)
      npfittmp=imaxfit-iminfit+1
c      write(*,*) "iminfit,imaxfit: ",iminfit,imaxfit,npfittmp
c      call bigsvdfit(rifit,sfit,efit,npfit,acoeff,nv,uu,vv,ww,
c     $     nwlmax,nvmax,chisq1,funcs)
      call bigsvdfit(rifit(iminfit),sfit(iminfit),efit(iminfit),
     $     npfittmp,acoeff,nv,uu,vv,ww,nwlmax,nvmax,chisq1,funcs)
      write(*,*) "coeffs ",(acoeff(i),i=1,nv)
      do i=1,npfit
         call funcs(rifit(i),evals,nv)
         fitvect(i) = 0.0
         do j=1,nv
            fitvect(i) = fitvect(i) + acoeff(j)*evals(j)
         end do
      end do

      call showspec(npfit,wfit,sfit)
      call pglabel("log wavelength","counts"," ")
      call pgmtxt('T',2.5,0.0,0.0,sname)
      call pgmtxt('T',1.5,0.0,0.0,"spectrum and best fits")
      call pgqci(indexc)
      call pgsci(3)
      call pgline(npfit,wfit,fitvect)
      write(tlabel,'(a,f7.4)') "chi-sq, z=",zmin
      call pgmtxt('T',2.5,1.0,1.0,tlabel)
      call pgsci(indexc)

c  calculate the fit at rchimin
      koffset=indkrmin+kmin-1
      call findindex(npfit,ifit,jmin+koffset,jmax+koffset,
     $     iminfit,imaxfit)
      npfittmp=imaxfit-iminfit+1
c      call bigsvdfit(rifit,sfit,efit,npfit,acoeff,nv,uu,vv,ww,
c     $     nwlmax,nvmax,chisq1,funcs)
      call bigsvdfit(rifit(iminfit),sfit(iminfit),efit(iminfit),
     $     npfittmp,acoeff,nv,uu,vv,ww,nwlmax,nvmax,chisq1,funcs)
c      write(*,*) "coeffs ",(acoeff(i),i=1,nv)
      do i=1,npfit
         call funcs(rifit(i),evals,nv)
         fitvect(i) = 0.0
         do j=1,nv
            fitvect(i) = fitvect(i) + acoeff(j)*evals(j)
         end do
      end do
      call pgqci(indexc)
      call pgsci(2)
      call pgline(npfit,wfit,fitvect)
      write(tlabel,'(a,f7.4)') "red.chisq, z=",zrmin
      call pgmtxt('T',1.5,1.0,1.0,tlabel)
      call pgsci(indexc)

      go to 220

 666  continue

      close(2)
      close(3)
      close(4)
      close(10)
      close(11)
      if (ifzfile .ne. 0) close(12)

      call pgend()

      end

cccccccccccccccccccc
c  funcs returns the values of the nv eigenvectors
c  evaluated at position xfit, in the array evals

c  if svdfit is called with argument wfit, xfit is the
c  log wavelength.
c  if svdfit is called with argument rifit, xfit is the
c  index in the log wavelength array, as a real.
c

      subroutine funcs(xfit,evals,nv)

      include 'specarray.h'
c having NMAXVECTS in here as well as fitz is a bad kludge
      parameter(NMAXVECTS=12)

      real evals(nv)

      common /vectorfit/ koffset,z,nvmax,nwlmax,
     $     evects(NMAXVECTS,NLOGWMAX)

c find the index corresponding to the position xfit
c and retrieve the eigenvector values

c  this is for use with rifit
      ispec = nint(xfit)
      jvect = ispec-koffset
      if(jvect .ge. 1 .and. jvect .le. nwlmax) then
         do i=1,nv
            evals(i) = evects(i,jvect)
         end do
      else
         do i=1,nv
            evals(i) = 0.0
         end do
      end if

      return
      end

cccccccccccccccccccccccccccccccccccccccc
c  Draw vertical lines for the actual z and z-estimates
c  on the chi-squared plots.

      subroutine plotz(zreal,zest1,zest2)

      real xpl(2),ypl(2)
      integer indexc,indexls

      ypl(1) = -100.
      ypl(2) = 1.0e10
      call pgqci(indexc)
      call pgqls(indexls)
c      call pgsls(1)
      xpl(1) = zreal
      xpl(2) = zreal
      call pgsci(4)
      call pgsls(5)
      call pgline(2,xpl,ypl)
      xpl(1) = zest1
      xpl(2) = zest1
      call pgsci(3)
      call pgsls(4)
      call pgline(2,xpl,ypl)
      xpl(1) = zest2
      xpl(2) = zest2
      call pgsci(2)
      call pgsls(3)
      call pgline(2,xpl,ypl)
      call pgsci(indexc)
      call pgsls(indexls)

      return
      end


text:
      program info1
c----------------------------------------------------------------------
c usage    info1.ex < HB_file
c
c where info1 is the executable generated by makefile, HB_file is a
c file containing a matrix stored in Harwell-Boeing matrices.
c Info1 will then dump the information into the standard output.
c
c To use with larger matrices, increase nmax and nzmax.
c----------------------------------------------------------------------
      implicit none
      integer nmax, nzmax
      parameter (nmax = 30000, nzmax = 800000)
      integer ia(nmax+1),ia1(nmax+1),ja(nzmax),ja1(nzmax)
      real*8  a(nzmax),a1(nzmax),rhs(1)
      character title*72, type*3, key*8, guesol*2
      logical valued
c
      integer job, iin, nrow,ncol,nnz,ierr, nrhs, iout
c--------------
      data iin /5/, iout/6/
c--------------
      job = 2
      nrhs = 0
      call readmt (nmax,nzmax,job,iin,a,ja,ia, rhs, nrhs,
     *     guesol,nrow,ncol,nnz,title,key,type,ierr)
c---- if not readable return
      if (ierr .ne. 0) then
         write (iout,100) ierr
 100     format(' **ERROR: Unable to read matrix',/,
     *        ' Message returned fom readmt was ierr =',i3)
         stop
      endif
      valued = (job .ge. 2)
c-------
      call dinfo1(ncol,iout,a,ja,ia,valued,title,key,type,a1,ja1,ia1)
c--------------------end------------------------------------------------
c-----------------------------------------------------------------------
      end

max_stars_repo_name: xywei/Mirana

max_issues_repo_name: xywei/Mirana

max_forks_repo_name: xywei/Mirana


text:
module da_geoamv

   use module_domain, only : domain

   use da_control, only : obs_qc_pointer,max_ob_levels,missing_r, &
      v_interp_p, v_interp_h, check_max_iv_print, &
      missing, max_error_uv, max_error_t, rootproc, kms,kme,kts,kte, &
      max_error_p,max_error_q, check_max_iv_unit,check_max_iv, trace_use_dull, &
      max_stheight_diff,missing_data,max_error_bq,max_error_slp,fails_error_max, &
      max_error_bt, max_error_buv, geoamv, anal_type_verify,qcstat_conv_unit,ob_vars, &
      convert_fd2uv, convert_uv2fd, max_error_spd, max_error_dir, max_omb_spd, max_omb_dir, pi, qc_rej_both,&
      wind_sd_geoamv, wind_stats_sd, write_rej_obs_conv
   use da_grid_definitions, only : da_ffdduv, da_ffdduv_model,da_ffdduv_diagnose
   use da_physics, only : da_uv_to_sd_lin, da_uv_to_sd_adj
   use da_define_structures, only : maxmin_type, iv_type, y_type, jo_type, &
      bad_data_type, x_type, number_type, bad_data_type, &
      maxmin_type
  use da_interpolation, only : da_interp_lin_3d, da_to_zk, &
      da_interp_lin_3d_adj
   use da_par_util, only : da_proc_stats_combine
   use da_par_util1, only : da_proc_sum_int
   use da_statistics, only : da_stats_calculate
   use da_tools, only : da_max_error_qc, da_residual, da_convert_zk, da_get_print_lvl
   use da_tracing, only : da_trace_entry, da_trace_exit

   ! The "stats_geoamv_type" is ONLY used locally in da_geoamv:

   type residual_geoamv1_type
      real          :: u                        ! u-wind.
      real          :: v                        ! v-wind.
   end type residual_geoamv1_type

   type maxmin_geoamv_stats_type
      type (maxmin_type)         :: u, v, t, q
   end type maxmin_geoamv_stats_type

   type stats_geoamv_type
      type (maxmin_geoamv_stats_type)  :: maximum, minimum
      type (residual_geoamv1_type)     :: average, rms_err
   end type stats_geoamv_type

contains

#include "da_ao_stats_geoamv.inc"
#include "da_jo_and_grady_geoamv.inc"
#include "da_residual_geoamv.inc"
#include "da_oi_stats_geoamv.inc"
#include "da_print_stats_geoamv.inc"
#include "da_transform_xtoy_geoamv.inc"
#include "da_transform_xtoy_geoamv_adj.inc"
#include "da_check_max_iv_geoamv.inc"
#include "da_get_innov_vector_geoamv.inc"
#include "da_calculate_grady_geoamv.inc"

end module da_geoamv

max_stars_repo_name: matzegoebel/WRF-fluxavg

max_issues_repo_name: teb-model/wrf-teb

max_forks_repo_name: teb-model/wrf-teb


text:
subroutine trtrou(lundia    ,nmax      ,mmax      ,nmaxus    ,kmax      , &
                & cfrou     ,rouflo    ,linit     ,gdis_dp   ,gdis_zet  , &
                & huv       ,kcuv      ,uvdir     ,uvperp    ,sig       , &
                & z0rou     ,deltadir  ,idir      ,gdp       )
!----- GPL ---------------------------------------------------------------------
!
!  Copyright (C)  Stichting Deltares, 2011-2016.
!
!  This program is free software: you can redistribute it and/or modify
!  it under the terms of the GNU General Public License as published by
!  the Free Software Foundation version 3.
!
!  This program is distributed in the hope that it will be useful,
!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!  GNU General Public License for more details.
!
!  You should have received a copy of the GNU General Public License
!  along with this program.  If not, see <http://www.gnu.org/licenses/>.
!
!  contact: delft3d.support@deltares.nl
!  Stichting Deltares
!  P.O. Box 177
!  2600 MH Delft, The Netherlands
!
!  All indications and logos of, and references to, "Delft3D" and "Deltares"
!  are registered trademarks of Stichting Deltares, and remain the property of
!  Stichting Deltares. All rights reserved.
!
!-------------------------------------------------------------------------------
!  $Id: trtrou.f90 6319 2016-07-18 12:19:10Z ye $
!  $HeadURL: https://svn.oss.deltares.nl/repos/delft3d/tags/6686/src/engines_gpl/flow2d3d/packages/kernel/src/compute/trtrou.f90 $
!!--description-----------------------------------------------------------------
!
! Calculate rougness due to trachytopes.
! Routine is called for U/V-direction
! respectively.
!
!!--pseudo code and references--------------------------------------------------
! NONE
!!--declarations----------------------------------------------------------------
    use precision
    use mathconsts
    use globaldata
    use string_module
    !
    implicit none
    !
    type(globdat),target :: gdp
    !
    ! The following list of pointer parameters is used to point inside the gdp structure
    !
    integer                    , pointer :: i50
    integer                    , pointer :: i90
    integer                    , pointer :: iarea_avg
    integer                    , pointer :: nttaru
    integer                    , pointer :: nttarv
    integer                    , pointer :: ntrt
    integer , dimension(:,:)   , pointer :: ittaru
    integer , dimension(:,:)   , pointer :: ittarv
    integer , dimension(:,:)   , pointer :: ittdef
    real(fp)                   , pointer :: eps
    real(fp)                   , pointer :: rhow
    real(fp)                   , pointer :: ag
    real(fp)                   , pointer :: z0
    real(fp)                   , pointer :: vonkar
    real(fp)                   , pointer :: vicmol
    real(fp)                   , pointer :: dryflc
    real(fp)                   , pointer :: alf_area_ser
    real(fp)                   , pointer :: trtminh
    real(fp), dimension(:,:)   , pointer :: rgcalu
    real(fp), dimension(:,:)   , pointer :: rgcalv
    real(fp), dimension(:)     , pointer :: rttaru
    real(fp), dimension(:)     , pointer :: rttarv
    real(fp), dimension(:,:)   , pointer :: rttdef
    real(fp), dimension(:,:,:) , pointer :: rttfu
    real(fp), dimension(:,:,:) , pointer :: rttfv
    real(fp), dimension(:,:)   , pointer :: vegh2d
    real(fp), dimension(:,:)   , pointer :: vden2d
    logical                    , pointer :: flsedprop_rqrd
    logical                    , pointer :: spatial_bedform
    logical                    , pointer :: v2dwbl
    logical                    , pointer :: waqol
    !
    real(fp), dimension(:,:)   , pointer :: dxx
    real(fp), dimension(:)     , pointer :: rhosol
    real(fp), dimension(:)     , pointer :: bedformD50
    real(fp), dimension(:)     , pointer :: bedformD90
    real(fp), dimension(:)     , pointer :: rksr
    real(fp), dimension(:)     , pointer :: rksmr
    real(fp), dimension(:)     , pointer :: rksd
!
! Local parameters
!
    integer , parameter :: max_cl   = 8
    integer , parameter :: ch_type  = 1
    integer , parameter :: kn_type  = 0
    integer , parameter :: area_rgh = 1
    integer , parameter :: line_rgh = 2
    integer , parameter :: pnt_rgh  = 3
    integer , parameter :: spec_rgh = 0
    integer , parameter :: skip_rgh = -999
!
! Global variables
!
    integer                                                            , intent(in)  :: idir   !  Flag for direction, 1=U, 2=V
    integer                                                            , intent(in)  :: kmax
    integer                                                                          :: lundia
    integer                                                            , intent(in)  :: mmax
    integer                                                            , intent(in)  :: nmax
    integer                                                            , intent(in)  :: nmaxus
    integer, dimension(gdp%d%nlb:gdp%d%nub, gdp%d%mlb:gdp%d%mub)       , intent(in)  :: kcuv
    logical                                                            , intent(in)  :: linit
    real(fp), dimension(kmax)                                          , intent(in)  :: sig
    real(fp), dimension(gdp%d%nlb:gdp%d%nub, gdp%d%mlb:gdp%d%mub)      , intent(in)  :: deltadir
    real(fp), dimension(gdp%d%nlb:gdp%d%nub, gdp%d%mlb:gdp%d%mub)      , intent(in)  :: gdis_dp
    real(fp), dimension(gdp%d%nlb:gdp%d%nub, gdp%d%mlb:gdp%d%mub)      , intent(in)  :: gdis_zet
    real(fp), dimension(gdp%d%nlb:gdp%d%nub, gdp%d%mlb:gdp%d%mub)      , intent(in)  :: huv
    real(fp), dimension(gdp%d%nlb:gdp%d%nub, gdp%d%mlb:gdp%d%mub)      , intent(in)  :: z0rou
    real(fp), dimension(gdp%d%nlb:gdp%d%nub, gdp%d%mlb:gdp%d%mub, 3)                 :: cfrou
    real(fp), dimension(gdp%d%nlb:gdp%d%nub, gdp%d%mlb:gdp%d%mub, kmax), intent(in)  :: uvdir
    real(fp), dimension(gdp%d%nlb:gdp%d%nub, gdp%d%mlb:gdp%d%mub, kmax), intent(in)  :: uvperp
    character(4)                                                                     :: rouflo
!
! Local variables
!
    integer                     :: ifrom
    integer                     :: ilist
    integer                     :: ircod
    integer                     :: ita
    integer                     :: ito
    integer                     :: itrt
    integer, dimension(max_cl)  :: itrt_list
    integer                     :: itrt_user
    integer                     :: k
    integer                     :: kmaxx
    integer                     :: m
    integer                     :: mc
    integer                     :: md
    integer                     :: ml
    integer                     :: mu
    integer                     :: n
    integer                     :: nc
    integer                     :: nd
    integer                     :: nm
    integer                     :: nmu
    integer                     :: nl
    integer                     :: nlist
    integer                     :: nu
    integer                     :: num
    integer, dimension(2)       :: numlen
    integer                     :: rgh_geom
    integer                     :: rgh_type
    logical                     :: lfound
    logical                     :: lnew
    real(fp)                    :: a0
    real(fp)                    :: a1
    real(fp)                    :: alfa
    real(fp)                    :: alfam
    real(fp)                    :: betam
    real(fp)                    :: bfh
    real(fp)                    :: bfl
    real(fp)                    :: calca1
    real(fp)                    :: calca2
    real(fp)                    :: cbed
    real(fp)                    :: ch_icode
    real(fp)                    :: ch_lin_ser
    real(fp)                    :: ch_pnt_ser
    real(fp)                    :: ch_sum_par
    real(fp)                    :: ch_sum_ser
    real(fp)                    :: d50
    real(fp)                    :: d90
    real(fp)                    :: densit
    real(fp)                    :: depth
    real(fp)                    :: drag
    real(fp)                    :: dstar
    real(fp)                    :: e1
    real(fp)                    :: f
    real(fp)                    :: fracbu
    real(fp)                    :: fraccu
    real(fp), dimension(max_cl) :: fraccu_list
    real(fp)                    :: fracto
    real(fp)                    :: hk
    real(fp)                    :: kbed
    real(fp)                    :: kn_icode
    real(fp)                    :: kn_sum
    real(fp)                    :: rc0
    real(fp)                    :: rc3
    real(fp)                    :: rcgrn
    real(fp)                    :: rchmin
    real(fp)                    :: rcmain
    real(fp)                    :: relden
    real(fp)                    :: relh2
    real(fp)                    :: rkgrn
    real(fp)                    :: rkmain
    real(fp)                    :: rktra
    real(fp)                    :: rksdu
    real(fp)                    :: rksmru
    real(fp)                    :: rksru
    real(fp)                    :: rleng
    real(fp)                    :: rm0
    real(fp)                    :: rmndrf
    real(fp)                    :: rmup
    real(fp)                    :: t1
    real(fp)                    :: t2
    real(fp)                    :: t3
    real(fp)                    :: t4
    real(fp)                    :: t5
    real(fp)                    :: thetac
    real(fp)                    :: thetag
    real(fp)                    :: thetam
    real(fp)                    :: tsp
    real(fp)                    :: u2dh
    real(fp)                    :: ubsvg2
    real(fp)                    :: ucbsv2
    real(fp)                    :: umod
    real(fp)                    :: uuu
    real(fp)                    :: uv0
    real(fp)                    :: vheigh
    real(fp)                    :: vd2d
    real(fp)                    :: vh2d
    real(fp)                    :: vvv
    real(fp)                    :: vz0
    real(fp)                    :: zcc
    character(12), dimension(2) :: cnum
    character(132)              :: cmsg
!
    integer                   , pointer :: nttar
    integer , dimension(:,:)  , pointer :: ittar
    real(fp), dimension(:,:)  , pointer :: rgcal
    real(fp), dimension(:,:,:), pointer :: rttf
    real(fp), dimension(:)    , pointer :: rttar
!
!! executable statements -------------------------------------------------------
!
    rhow            => gdp%gdphysco%rhow
    ag              => gdp%gdphysco%ag
    z0              => gdp%gdphysco%z0
    vonkar          => gdp%gdphysco%vonkar
    vicmol          => gdp%gdphysco%vicmol
    eps             => gdp%gdconst%eps
    dryflc          => gdp%gdnumeco%dryflc
    v2dwbl          => gdp%gdnumeco%v2dwbl
    !
    alf_area_ser    => gdp%gdtrachy%alf_area_ser
    trtminh         => gdp%gdtrachy%trtminh
    iarea_avg       => gdp%gdtrachy%iarea_avg
    nttaru          => gdp%gdtrachy%nttaru
    nttarv          => gdp%gdtrachy%nttarv
    ntrt            => gdp%gdtrachy%ntrt
    ittaru          => gdp%gdtrachy%ittaru
    ittarv          => gdp%gdtrachy%ittarv
    ittdef          => gdp%gdtrachy%ittdef
    rgcalu          => gdp%gdtrachy%rgcalu
    rgcalv          => gdp%gdtrachy%rgcalv
    rttaru          => gdp%gdtrachy%rttaru
    rttarv          => gdp%gdtrachy%rttarv
    rttdef          => gdp%gdtrachy%rttdef
    rttfu           => gdp%gdtrachy%rttfu
    rttfv           => gdp%gdtrachy%rttfv
    vegh2d          => gdp%gdtrachy%vegh2d
    vden2d          => gdp%gdtrachy%vden2d
    flsedprop_rqrd  => gdp%gdtrachy%flsedprop_rqrd
    waqol           => gdp%gdwaqpar%waqol
    !
    dxx             => gdp%gderosed%dxx
    rhosol          => gdp%gdsedpar%rhosol
    i50             => gdp%gdmorpar%i50
    i90             => gdp%gdmorpar%i90
    bedformD50      => gdp%gdbedformpar%bedformD50
    bedformD90      => gdp%gdbedformpar%bedformD90
    spatial_bedform => gdp%gdbedformpar%spatial_bedform
    rksr            => gdp%gdbedformpar%rksr
    rksmr           => gdp%gdbedformpar%rksmr
    rksd            => gdp%gdbedformpar%rksd
    !
    ! refer to appropriate arrays depending on direction
    !
    if (idir==1) then
       !
       ! U-direction
       !
       nttar => nttaru
       ittar => ittaru
       rgcal => rgcalu
       rttf  => rttfu
       rttar => rttaru
    else
       !
       ! V-direction
       !
       nttar => nttarv
       ittar => ittarv
       rgcal => rgcalv
       rttf  => rttfv
       rttar => rttarv
    endif
    !
    ! initialize ...
    !
    ! Reset RTTF to zero
    !
    do k = 1, kmax
       do m = 1, mmax
          do n = 1, nmax
             rttf(n, m, k) = 0.0_fp
          enddo
       enddo
    enddo
    !
    ! Copy CFROU. When in init-mode, make backup of defaults,
    ! (CFROU(..,..,3) = CFROU(..,..,2))
    ! otherwise retrieve backup (CFROU(..,..,2) = CFROU(..,..,3)).
    !
    if (linit) then
       ifrom = 2
       ito   = 3
    else
       ifrom = 3
       ito   = 2
    endif
    do m = 1, mmax
       do n = 1, nmax
          cfrou(n, m, ito) = cfrou(n, m, ifrom)
       enddo
    enddo
    !
    ! Initialize locals
    !
    fracbu     = 0.0_fp
    fracto     = 0.0_fp
    ch_sum_par = 0.0_fp
    ch_sum_ser = 0.0_fp
    kn_sum     = 0.0_fp
    depth      = 0.0_fp
    ch_lin_ser = 0.0_fp
    ch_pnt_ser = 0.0_fp
    ml         = -99999
    nl         = -99999
    lnew       = .false.
    !
    ! Main loop over area data
    !
    do ita = 1, nttar
       !
       ! Get the next point
       !
       nc = ittar(ita, 1)
       mc = ittar(ita, 2)
       !
       ! Check if it is a new cell
       !
       lnew = nc/=nl .or. mc/=ml
       !
       if (lnew) then
          !
          ! Calculate the roughness for the old grid cell
          !
          if (nl>0 .and. ml>0) then
             !
             ! Calculate the roughness for the old grid cell
             !
             call calrou(kn_sum          ,fracto          ,fracbu    ,depth     ,ch_lin_ser, &
                       & cfrou(nl, ml, 2),cfrou(nl, ml, 3),rouflo    , &
                       & iarea_avg       ,ch_sum_par      ,ch_sum_ser,ch_pnt_ser, &
                       & alf_area_ser    )
          endif
          !
          ! Check if the new point was a block separator
          ! or an inactive point.
          ! If so reset last point and go to the next point
          !
          if (nc== - 1 .and. mc== - 1) then
             !
             ! Reset and get next point
             !
             nl = -99999
             ml = -99999
             cycle
          elseif (kcuv(nc, mc)==0) then
             !
             ! Reset and get next point
             ! This must be a separate elseif branch: kcuv(-1,-1) is not allowed to be read
             !
             nl = -99999
             ml = -99999
             cycle
          endif
          !
          ! Initialise new cell
          !
          nl         = nc
          ml         = mc
          nu         = min(nc + 1, nmax)
          mu         = min(mc + 1, mmax)
          nd         = max(nc - 1, 1)
          md         = max(mc - 1, 1)
          fracbu     = 0.0_fp
          fracto     = 0.0_fp
          ch_sum_par = 0.0_fp
          ch_sum_ser = 0.0_fp
          kn_sum     = 0.0_fp
          ch_lin_ser = 0.0_fp
          ch_pnt_ser = 0.0_fp
          depth      = max(trtminh , huv(nc, mc))
          !
          kmaxx = kmax
          !
          if (v2dwbl) then
             do k = kmax, 1, -1
                zcc = (1.0 + sig(k))*depth
                if (zcc>deltadir(nc, mc) .or. zcc>0.5*depth) then
                   kmaxx = k
                   exit
                endif
             enddo
          endif
          !
          if (idir==1) then
             call n_and_m_to_nm(nc, mc, nm, gdp)
             call n_and_m_to_nm(nc, mu, nmu, gdp)
             !
             ! U-direction
             !
             if (flsedprop_rqrd) then
                !
                ! get D50, D90, relative density from morphology module
                !
                if (associated(gdp%gderosed%dxx))  then
                   d50 = 0.5_fp*(dxx(nm, i50) + dxx(nmu, i50))
                   d90 = 0.5_fp*(dxx(nm, i90) + dxx(nmu, i90))
                   relden = (rhosol(1)-rhow)/rhow ! assume density equal for all sediment fractions
                else
                   if (spatial_bedform) then
                      d50 = 0.5_fp*(bedformD50(nm) + bedformD50(nmu))
                      d90 = 0.5_fp*(bedformD90(nm) + bedformD90(nmu))
                   else
                      d50 = bedformD50(1)
                      d90 = bedformD90(1)
                   endif
                   relden = 1.65_fp
                endif
             endif
             !
             ! Average perpendicular velocity
             !
             vvv = 0.25_fp*(  uvperp(nc, mc, kmaxx) + uvperp(nc, mu, kmaxx)  &
                 &          + uvperp(nd, mc, kmaxx) + uvperp(nd, mu, kmaxx))
             !
             ! Riples, Megaripples, Dunes
             !
             rksru  = (rksr (nm) + rksr (nmu))*0.5_fp
             rksmru = (rksmr(nm) + rksmr(nmu))*0.5_fp
             rksdu  = (rksd (nm) + rksd (nmu))*0.5_fp
             if (waqol) then
                !
                ! 2D Vegetation characterics (coming from WAQ)
                !
                vd2d = vden2d(nc,mc) + vden2d(nc,mu)
                vh2d = (vegh2d(nc,mc)*vden2d(nc,mc) + vegh2d(nc,mu)*vden2d(nc,mu)) / max(vd2d,eps)
                vd2d = vd2d * 0.5_fp
             endif
          else
             call n_and_m_to_nm(nc, mc, nm, gdp)
             call n_and_m_to_nm(nu, mc, num, gdp)
             !
             ! V-direction
             !
             if (flsedprop_rqrd) then
                !
                ! get D50, D90, relative density from morphology module
                !
                if (associated(gdp%gderosed%dxx))  then
                   d50 = 0.5_fp*(dxx(nm, i50) + dxx(num, i50))
                   d90 = 0.5_fp*(dxx(nm, i90) + dxx(num, i90))
                   relden = (rhosol(1)-rhow)/rhow ! assume density equal for all sediment fractions
                else
                   if (spatial_bedform) then
                      d50 = 0.5_fp*(bedformD50(nm) + bedformD50(num))
                      d90 = 0.5_fp*(bedformD90(nm) + bedformD90(num))
                   else
                      d50 = bedformD50(1)
                      d90 = bedformD90(1)
                   endif
                   relden = 1.65_fp
                endif
             endif
             !
             ! Average perpendicular velocity
             !
             vvv = 0.25_fp*(  uvperp(nc, mc, kmaxx) + uvperp(nc, md, kmaxx)  &
                 &          + uvperp(nu, mc, kmaxx) + uvperp(nu, md, kmaxx))
             !
             ! Riples, Megaripples, Dunes
             !
             rksru  = (rksr (nm) + rksr (num))*0.5_fp
             rksmru = (rksmr(nm) + rksmr(num))*0.5_fp
             rksdu  = (rksd (nm) + rksd (num))*0.5_fp
             if (waqol) then
                !
                ! 2D Vegetation characterics (coming from WAQ)
                !
                vd2d = vden2d(nc,mc) + vden2d(nu,mc)
                vh2d = (vegh2d(nc,mc)*vden2d(nc,mc) + vegh2d(nu,mc)*vden2d(nu,mc)) / max(vd2d,eps)
                vd2d = vd2d * 0.5_fp
             endif
          endif
          !
          ! Depth-average velocity (similar as in TAUBOT)
          !
          uuu  = uvdir(nc, mc, kmaxx)
          umod = sqrt(uuu**2 + vvv**2)
          if (kmax==1) then
             u2dh = umod
          else
             u2dh = (umod/depth*((depth + z0rou(nc, mc))         &
                  &              *log(1.0_fp + depth/max(z0rou(nc, mc),1.0e-20_fp)) &
                  &              - depth)                         ) &
                  & /log(1.0_fp + (1.0_fp + sig(kmaxx))*depth/max(z0rou(nc, mc),1.0e-20_fp))
          endif
       endif
       !
       ! Get the trachytope and roughness characteristic
       !
       itrt_list(1)   = ittar(ita, 3)
       fraccu_list(1) = rttar(ita)
       nlist          = 1
       !
       ilist=1
       do while (ilist<=nlist)
          itrt_user = itrt_list(ilist)
          fraccu    = fraccu_list(ilist)
          !
          lfound = .false.
          do itrt = 1, ntrt
             if (ittdef(itrt, 1)==itrt_user) then
                lfound = .true.
                exit
             endif
          enddo
          if (.not.lfound) then
             call prterr(lundia    ,'J001'    ,'TRTROU: Trachytope not found.' )
             call d3stop(1, gdp)
          endif
          !
          ircod = ittdef(itrt, 2)
          !
          if (ircod==2) then
             !
             ! Combination of area roughnesses: expand
             !
             nlist = nlist+1
             if (nlist>max_cl) then
                call prterr(lundia    ,'J001'    ,'TRTROU: Maximum recursion depth.' )
                call d3stop(1, gdp)
             endif
             itrt_list(ilist)   = nint(rttdef(itrt, 1))
             itrt_list(nlist)   = nint(rttdef(itrt, 2))
             fraccu_list(ilist) = fraccu*rttdef(itrt, 3)
             fraccu_list(nlist) = fraccu*rttdef(itrt, 4)
          else
             !
             ! Other roughness formulation: fill in itrt
             ! this is the effective internal trachytope
             ! number where the formula known to the user
             ! as trachytope type itrt_user is stored.
             !
             itrt_list(ilist) = itrt
             ilist            = ilist+1
          endif
       enddo
       !
       ! Calculate the roughnes according to the given types
       !
       do ilist = 1, nlist
          itrt   = itrt_list(ilist)
          fraccu = fraccu_list(ilist)
          ircod  = ittdef(itrt, 2)
          !
          rgh_type = kn_type
          rgh_geom = spec_rgh
          if (ircod==1) then
             !
             ! Water free area
             !
             fracbu = fracbu + fraccu
          elseif (ircod==2) then
             !
             ! Combination of area roughnesses: parsed above
             ! So, the program should never come here.
             !
          elseif (ircod==51) then
             !
             ! Constant White-Colebrook / Nikuradse value
             !
             kn_icode = rttdef(itrt, 1)
             rgh_geom = area_rgh
          elseif (ircod==52) then
             !
             ! Constant Chezy value
             !
             ch_icode = rttdef(itrt, 1)
             rgh_type = ch_type
             rgh_geom = area_rgh
          elseif (ircod==53) then
             !
             ! Constant Manning value
             !
             ch_icode = (depth**(1.0_fp/6.0_fp))/rttdef(itrt, 1)
             rgh_type = ch_type
             rgh_geom = area_rgh
          elseif (ircod==54) then
             !
             ! Constant Z0 value
             !
             kn_icode = 30.0_fp*rttdef(itrt, 1)
             rgh_geom = area_rgh
          elseif (ircod==101) then
             !
             ! WAQUA roughness predictor
             !         (Van Rijn simplified)
             !
             alfam    = rttdef(itrt, 1)
             betam    = rttdef(itrt, 2)
             rkmain   = alfam*depth**0.7_fp*(1.0_fp - exp( - betam*depth**( -0.3_fp)))
             kn_icode = max(rkmain, 0.00001_fp)
             rgh_geom = area_rgh
          elseif (ircod==102) then
             !
             ! power law: C = A * H**B
             !
             ch_icode = rttdef(itrt, 1)*depth**rttdef(itrt, 2)
             rgh_type = ch_type
             rgh_geom = area_rgh
          elseif (ircod==103) then
             !
             ! Van Rijn roughness predictor
             !
             dstar = d50*(ag*relden/vicmol**2)**(1.0_fp/3.0_fp)
             !
             if (dstar<=4.0_fp) then
                thetac = 0.24_fp/dstar
             elseif (dstar<=10.0_fp) then
                thetac = 0.14_fp*dstar**( - 0.64_fp)
             elseif (dstar<=20.0_fp) then
                thetac = 0.04_fp*dstar**( - 0.10_fp)
             elseif (dstar<=150.0_fp) then
                thetac = 0.013_fp*dstar**(0.29_fp)
             else
                thetac = 0.055_fp
             endif
             !
             ucbsv2 = ag*relden*d50*thetac
             !
             rkgrn  = 3.0_fp*d90
             rcgrn  = 18.0_fp*log10(12.0_fp*depth/rkgrn)
             ubsvg2 = ag*u2dh**2/rcgrn**2
             tsp    = (ubsvg2 - ucbsv2)/ucbsv2
             !
             if (tsp<=0.0_fp .or. tsp>=25.0_fp) then
                bfh = 0.0_fp
             else
                t1  = 0.11_fp*depth
                t2  = (d50/depth)**(0.3_fp)
                t3  = 1.0_fp - exp( - 0.5_fp*tsp)
                t4  = 25.0_fp - tsp
                bfh = t1*t2*t3*t4
             endif
             !
             bfl      = max(7.3_fp*depth,1e-6_fp)
             !
             t1       = 1.0_fp - exp(-25.0_fp*bfh/bfl)
             kn_icode = rkgrn + 1.1_fp*bfh*t1
             rgh_geom = area_rgh
          elseif (ircod==104) then
             !
             ! Struiksma roughness predictor
             !
             calca1 = rttdef(itrt, 1)
             calca2 = rttdef(itrt, 2)
             thetac = rttdef(itrt, 3)
             thetam = rttdef(itrt, 4)
             rchmin = rttdef(itrt, 5)
             !
             rcgrn = calca1*log10(calca2*depth/d90)
             !
             thetag = u2dh**2/(rcgrn**2*relden*d50)
             if (u2dh<=0.0_fp) then
                thetag = 0.001_fp
             elseif (thetag<=0.0_fp) then
                write (cnum(1), '(i12)') nc
                call remove_leading_spaces(cnum(1)   ,numlen(1) )
                write (cnum(2), '(i12)') mc
                call remove_leading_spaces(cnum(2)   ,numlen(2) )
                cmsg = cnum(1)(1:numlen(1)) // ', ' // cnum(2)(1:numlen(2))
                call prterr(lundia    ,'J015'    ,cmsg      )
                call d3stop(1, gdp)
                !
                thetag = 0.001_fp
             endif
             !
             if (thetag>=1.0_fp) then
                write (cnum(1), '(i12)') nc
                call remove_leading_spaces(cnum(1)   ,numlen(1) )
                write (cnum(2), '(i12)') mc
                call remove_leading_spaces(cnum(2)   ,numlen(2) )
                cmsg = cnum(1)(1:numlen(1)) // ', ' // cnum(2)(1:numlen(2))
                call prterr(lundia    ,'J016'    ,cmsg      )
                call d3stop(1, gdp)
                !
                thetag = 0.999_fp
             endif
             !
             if (thetag<=thetac) then
                rcmain = rcgrn
             else
                t1     = thetac*(thetag - thetac)
                t2     = thetag - (thetam**2/thetac)
                t3     = ((rcgrn/rchmin)**2) - 1.0_fp
                t4     = thetag*(thetam - thetac)**2
                rmndrf = 1.0_fp - (t1*t2*t3/t4)
                rcmain = rcgrn/sqrt(rmndrf)
             endif
             !
             ch_icode = rcmain
             rgh_type = ch_type
             rgh_geom = area_rgh
          elseif (ircod==105) then
             !
             ! Quadratic combination of bedform heights (Van Rijn 2004)
             !
             rkmain   = sqrt(rksru**2 + rksmru**2 + rksdu**2)
             kn_icode = min(rkmain,0.5_fp*depth)
             rgh_geom = area_rgh
          elseif (ircod==106) then
             !
             ! Linear addition of bedform heights
             !
             rkmain   = rksru + rksmru + rksdu
             kn_icode = min(rkmain,0.5_fp*depth)
             rgh_geom = area_rgh
          elseif (ircod==151 .or. ircod==152) then
             !
             ! Vegetation roughness WAQUA / Van Barneveld
             !
             vheigh = rttdef(itrt, 1)
             densit = rttdef(itrt, 2)
             if (ircod==152) then
                drag = rttdef(itrt, 3)
                kbed = rttdef(itrt, 4)
                cbed = 18.0_fp*log10(12.0_fp*depth/kbed)
             else
                drag = 1.65_fp
                cbed = 1.0e5_fp
             endif
             !
             if (vheigh>=depth) then
                !
                ! For flow through vegetation only
                !
                if (ircod==152) then
                   rc0  = sqrt(1.0_fp/((drag*densit*vheigh)/(2.0_fp*ag)+1.0_fp/(cbed*cbed)))
                else
                   rc0  = sqrt((2.0_fp*ag)/(drag*densit*depth))
                endif
             else
                !
                ! For flow through and over vegetation
                !
                if (ircod==152) then
                   uv0  = sqrt(vheigh/((drag*densit*vheigh)/(2.0_fp*ag)+1.0_fp/(cbed*cbed)))
                   !
                   alfa = 0.0227_fp*vheigh**0.7_fp
                else
                   uv0  = sqrt((2.0_fp*ag)/(drag*densit))
                   !
                   alfa = 0.01_fp*sqrt(depth*vheigh)
                endif
                alfa = max(alfa, 0.001_fp)
                a0   = (densit*drag)/(2.0_fp*alfa)
                !
                t1   = sqrt(2.0_fp*a0)
                rc3  = (2.0_fp*ag*(depth - vheigh))                   &
                     & /(alfa*t1*(exp(vheigh*t1) + exp( - vheigh*t1)))
                t2   = rc3*exp(vheigh*t1)
                t3   = sqrt(t2 + (uv0*uv0))
                t4   = sqrt(rc3 + (uv0*uv0))
                e1   = (t1*t2)/(2.0_fp*t3)
                a1   = (1.0_fp + sqrt(1.0_fp + ((4.0_fp*(e1**2)*(vonkar**2)*(depth - vheigh))/ag)))&
                     & /((2.0_fp*(e1**2)*(vonkar**2))/ag)
                a1   = min(a1,vheigh)
                t5   = depth - (vheigh - a1)
                f    = (vonkar*t3)/sqrt(ag*t5)
                vz0  = a1*exp( - f)
                rc0  = 1.0_fp/(depth**1.5_fp)                                          &
                     & *((2.0_fp/t1)*(t3 - t4) + (uv0/t1)*log(((t3 - uv0)*(t4 + uv0))  &
                     & /((t3 + uv0)*max((t4 - uv0), eps_fp)))                          &
                     & + (sqrt(ag*t5)/vonkar)                                          &
                     & *((t5*log(t5/vz0)) - (a1*log(a1/vz0)) - (depth - vheigh)))
             endif
             !
             if (ircod==152) then
                ch_icode = rc0
                rgh_type = ch_type
             else
                rktra = (12.0_fp*depth)/10.0_fp**(rc0/18.0_fp)
                !
                ! 0.50 is a kind of bed rougness under trachytopes
                !
                kn_icode = max(rktra, 0.50_fp)
             endif
             rgh_geom = area_rgh
          elseif (ircod==153 .or. ircod==154) then
             !
             ! Baptist vegetation formulation
             !
             vheigh = rttdef(itrt, 1)
             densit = rttdef(itrt, 2)
             if (vheigh < 0.0_fp) then
                vheigh = vh2d
                densit = vd2d
             endif
             drag   = rttdef(itrt, 3)
             cbed   = rttdef(itrt, 4)
             !
             rgh_type = ch_type
             rgh_geom = area_rgh
             !
             if (vheigh < eps) then
                rgh_geom = skip_rgh
             elseif (ircod==153) then
                if (depth>vheigh) then
                   ch_icode = 1.0_fp/sqrt(1.0_fp/(cbed*cbed) + &
                            &          (drag*densit*vheigh)/(2.0_fp*ag)) &
                            & + sqrt(ag)*log(depth/vheigh)/vonkar
                else
                   ch_icode = 1.0_fp/sqrt(1.0_fp/(cbed*cbed) + &
                            &          (drag*densit*depth)/(2.0_fp*ag))
                endif
             else
                hk     = max(1.0_fp,depth/vheigh)
                ch_icode = cbed + sqrt(ag)/vonkar*log(hk)* &
                         & sqrt(1.0_fp+(drag*densit*vheigh*cbed**2)/(2.0_fp*ag))
                rttf(nc, mc, 1) = rttf(nc, mc, 1) + fraccu * &
                         & drag*densit/hk*(cbed*cbed)/(ch_icode*ch_icode)
             endif
             !
             !rgh_type = ch_type
             !rgh_geom = area_rgh
          elseif (ircod==201) then
             !
             ! Get coefficients for hedges
             !
             vheigh = rttdef(itrt, 1)
             densit = rttdef(itrt, 2)
             !
             ! Calculate hedges
             !
             rleng = gdis_zet(nc, mc)
             if (vheigh>depth) then
                !
                ! For flow through vegetation only
                !
                rmup = 1.0_fp + 0.175_fp*densit*(depth/vheigh - 2.0_fp)
             else
                !
                ! For flow through and over vegetation
                !
                rmup = 1.0_fp - 0.175_fp*densit*vheigh/depth
             endif
             rmup   = min(0.999_fp, max(0.001_fp, rmup))
             ch_icode = sqrt((2.0_fp*ag*rleng*rmup**2)/(fraccu*depth*(1.0_fp - rmup**2)))
             !
             rgh_type = ch_type
             rgh_geom = line_rgh
          elseif (ircod==202) then
             !
             ! Get coefficients for hedges
             !
             vheigh = rttdef(itrt, 1)
             densit = rttdef(itrt, 2)
             drag   = rttdef(itrt, 3)
             rm0    = rttdef(itrt, 4)
             !
             ! Calculate hedges
             !
             rleng = gdis_zet(nc, mc)
             if (vheigh>depth) then
                !
                ! For flow through vegetation only
                !
                ch_icode = sqrt((2.0_fp*ag*rleng)/(fraccu*drag*densit*depth))
             else
                !
                ! For flow through and over vegetation
                !
                relh2 = ((depth-vheigh)/depth)**2
                ch_icode = sqrt((2.0_fp*ag*rleng)/(fraccu*depth))* &
                         & (vheigh/(depth*sqrt(drag*densit)) + &
                         &  rm0*sqrt(relh2/(1.0_fp-relh2)))
             endif
             !
             rgh_type = ch_type
             rgh_geom = line_rgh
          elseif (ircod==251) then
             !
             ! Get coefficients for trees
             !
             vheigh = rttdef(itrt, 1)
             drag   = rttdef(itrt, 2)
             densit = fraccu
             !
             ! Calculate trees
             !
             vheigh   = min(vheigh,depth)
             ch_icode = sqrt((2.0_fp*ag)/(drag*densit*vheigh))
             !
             rgh_type = ch_type
             rgh_geom = pnt_rgh
          else
             !
             ! Specified roughness type not implemented
             !
             call prterr(lundia    ,'J001'    ,'TRTROU: Specified roughness type not implemented.'   )
             call d3stop(1, gdp)
          endif
          !
          if (rgh_geom == area_rgh) then
             if (rgh_type == ch_type) then
                kn_icode = (12.0_fp*depth)/10.0_fp**(ch_icode/18.0_fp)
             elseif (rgh_type == kn_type) then
                ch_icode = 18.0_fp*log10(12.0_fp*depth/kn_icode)
             endif
             fracto     = fracto + fraccu
             kn_sum     = kn_sum + fraccu*kn_icode
             ch_sum_par = ch_sum_par + fraccu*ch_icode
             ch_sum_ser = ch_sum_ser + fraccu/(ch_icode**2)
          elseif (rgh_geom == line_rgh) then
             ch_lin_ser = ch_lin_ser + 1.0_fp/(ch_icode**2)
          elseif (rgh_geom == pnt_rgh) then
             ch_pnt_ser = ch_pnt_ser + 1.0_fp/(ch_icode**2)
          endif
       enddo
    enddo
    !
    ! Process the last pending cell
    !
    if (nttar>0 .and. nc>0 .and. mc>0) then
       depth = max(dryflc , huv(nc, mc))
       call calrou(kn_sum    ,fracto    ,fracbu    ,depth     ,ch_lin_ser, &
                 & cfrou(nc, mc, 2)     ,cfrou(nc, mc, 3)     ,rouflo    , &
                 & iarea_avg ,ch_sum_par,ch_sum_ser,ch_pnt_ser, &
                 & alf_area_ser         )
    endif
    !
    ! Multiply roughness with calibration factor
    !
    do m = 1, mmax
       do n = 1, nmaxus
          if (cfrou(n, m, 2)>0.0) then
             cfrou(n, m, 2) = cfrou(n, m, 2)*rgcal(n, m)
          endif
       enddo
    enddo
end subroutine trtrou

max_stars_repo_name: liujiamingustc/phd

max_issues_repo_name: liujiamingustc/phd

max_forks_repo_name: liujiamingustc/phd