<html><head><meta name="color-scheme" content="light dark"></head><body><pre style="word-wrap: break-word; white-space: pre-wrap;">package gem;

//   wlf_GEM6.f
import java.io.File;
import java.io.FileNotFoundException;
import java.io.PrintWriter;
import java.util.Arrays;
import java.util.Scanner;
import java.util.logging.Level;
import java.util.logging.Logger;

//     ##################################################################
//     ##################################################################
//     ######                                                      ######
//     ######                                                      ######
//     ######            PROGRAM FOR DAILY WEATHER SIMULATION      ######
//     ######               IN THE CONTIGOUS UNITED STATES         ######
//     ######                                                      ######
//     ######                     ARS - 114                        ######
//     ######                   (c) July  1994                     ######
//     #####                        Feburay 2000                   ######
//     ######                                                      ######
//     ######        UNITED STATES DEPARTMENT OF AGRICULTUE        ######
//     ######            AGRICULTURAL RESEARCH SERVICE.            ######
//     ######                  All rights reserved.                ######
//     ######                                                      ######
//     ##################################################################
//     ##################################################################
public class GEM6 {

//#######################################################################
//
//     PURPOSE:
//
//    This program provides precipitation probabilities and simulates
//    data for daily precipitation, maximum and minimum temperature, and 
//    solar radiation for an n-year period at a given location within the 
//    contiguous United States.  Wind Speed and Dewpoint were added in
//    Feb. 2000.  The model is designed to preserve the  dependence in
//    time, the internal correlation, and the seasonal characteristic
//    that exist in actual weather data. Daily maximum and minimum
//    temperture, dew point, wind speed, and solar radiation data are
//    simulated using a weakly stationary generating process conditioned
//    on the precipitation process described by a Markov chain-mixed
//    exponential model.
//
//    Paratmeters for special stations within a region can be accessed
//    directly.  An estimation routine(s) will be added for estimation of
//    data (and / or parameters) for points between existing stations.
//
//    The seasonal variations of parameters are described by the Fourier
//    series. Information on the type of equipment needed to run the model
//    and an example of running the model are provided. 
//
//#######################################################################
//
//     AUTHOR: C.L.Hanson,K.A.Cumming,D.A.Woolhiser and C.W.Richardson
//
//     Modify history:
//     Yunyun Lu (November,1994)
//     transfer program from QBASIC to FORTRAN 77
//
//     Y. Lu (Dec. 1994)
//     added the minimal documentation
//
//######################################################################
//     Will Frymire (Dec 1999)
//     transefered to PC version:
//     No graphics.
//     Input Name of Station location file (i.e.  "Boise" for Boise Idaho)
//     Added comments
//     Output file 'filename.out' starts in January.(i.e. Boise.out)
//#######################################################################
//
//     Variable Declarations.
//
//#######################################################################
//
//     OUTPUT:
//
//     p00        Transition probability of a dry day following a dry day
//     p10        Transition probability of a dry day following a wet day
//     beta       The mean of the smaller exponential distributions
//     delta      The mean of the larger exponential distributions
//     albar      The weighting function
//
//     anp        Expected annual precipitation (inches)
//     sum1       number of wet days
//
//     p          Simulated M years of daily rainfall
//     tmax       Simulated M years of daily max. temperature
//     tmin       Simulated M years of daily min. temperature
//     rad        Simulated M years of daily total radiation
//
//#######################################################################
//
//     Variable Declarations.
//
//#######################################################################
//
    public static final int idyr = 365;
    public static final int MAX_NUM_STA = 360;
    // common -- const1
    static double pi = Math.PI;
    // common -- sited2
    static double[] delta = new double[idyr];
    static double[] p00 = new double[idyr];
    static double[] p10 = new double[idyr];
    static double[] beta = new double[idyr];
    static double albar;
    // common -- rainf1
    static double sum1, anp;
    // common -- tem11
    static double[] tmad = new double[idyr];
    static double[] txs = new double[idyr];
    static double[] txm1 = new double[idyr];
    static double[] txs1 = new double[idyr];
    static double[] dmad = new double[idyr];
    static double[] dxs = new double[idyr];
    static double[] dxm1 = new double[idyr];
    static double[] dxs1 = new double[idyr];
    static double[] wmad = new double[idyr];
    static double[] wxs = new double[idyr];
    static double[] wxm1 = new double[idyr];
    static double[] wxs1 = new double[idyr];
    static double[] tnm = new double[idyr];
    static double[] tns = new double[idyr];
    static double[] tnm1 = new double[idyr];
    static double[] tns1 = new double[idyr];
    static double[] rm0 = new double[idyr];
    static double[] rs0 = new double[idyr];
    static double[] rm1 = new double[idyr];
    static double[] rs1 = new double[idyr];
    static double[] rm = new double[idyr];
    static double[] chipre = new double[6];
    // common -- temp2
    static double p, tmax, tmin, dew, wind, rad;
    // common -- temp3
    static double[] fran = new double[12];
    // common -- temp4
    static double[][][] a = new double[6][6][12];
    static double[][][] b = new double[6][6][12];

////
////%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
////
////     Beginning of executable code...
////
////%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
////
    public static void main(String[] args) {
        System.out.println(
                "###############################################################\n" +
                "###############################################################\n" +
                "#####                                                     #####\n" +
                "#####                  Welcome to                         #####\n" +
                "#####                                                     #####\n" +
                "#####        UNITED STATES WEATHER SIMULATOR GEM6         #####\n" +
                "#####                                                     #####\n" +
                "#####                      ARS - 114                      #####\n" +
                "#####                    (c) July  1994                   #####\n" +
                "#####                        FEB.  2000                   #####\n" +
                "#####                                                     #####\n" +
                "#####   This program provides easy access to informate    #####\n" +
                "#####       daily precipitation, maximum and minimum      #####\n" +
                "#####            temperture, and solar radiation.         #####\n" +
                "#####                                                     #####\n" +
                "#####         UNITED STATES DEPARTMENT OF AGRICULTUE      #####\n" +
                "#####             AGRICULTURAL RESEARCH SERVICE.          #####\n" +
                "#####                   All rights reserved.              #####\n" +
                "#####                                                     #####\n" +
                "###############################################################\n" +
                "###############################################################");

        input_parameters();
//      CALL INPUT_PARAMETERS   ! input fourier coefficients for
//  		              ! each station
    }
////     ##################################################################
////     ##################################################################
////     ######                                                      ######
////     ######            SUBROUTINE INPUT_PARAMETERS               ######
////     ######                                                      ######
////     ######                     (c) 1994                         ######
////     ######                                                      ######
////     ##################################################################
////     ##################################################################
//
//        SUBROUTINE INPUT_PARAMETERS
////
////#######################################################################
////
////     PURPOSE:
////
////     Input Fourier coefficients for each station.
////
////     Parameter interpolation for user chosen station. One option used.
////     (Choose Closest station)
////
////     Other options to be implemented at a later date: (after Feb. 2000)
////     (1) to use arithmetic average of parameters from stations
////         within a radius of 100 miles.
////     (2) to use the parameters of the nearest neighbor.
////
////#######################################################################
////
////     AUTHOR: C.L.Hanson,K.A.Cumming,D.A.Woolhiser and C.W.Richardson
////
////     Modify history:
////     Yunyun Lu (November,1994)
////     transfer program from QBASIC to FORTRAN 77
////     set up the algorism for distance oon curved earth surface
////
////     Y. Lu (Dec. 1994)
////     added the full documentation
////
////#######################################################################
////
////     INPUT:   (Inputs from AGUA46 ) This version of GEM uses only 3 harmonics
////
////      Total 360 stations in US region
////
////       x        station's longitude.
////       y        station's latitude.
////      stname    station's name.
////       z        parameter (p00,p10,beta,mu) matrix for
////                MCME model (upto 6 harmonics).
////                1. Annual mean
////                2. Amplitude of first harmonic
////                3. Phase angle of first harmonic
////                4. Amplitude of second harmonic
////                5. Phase angle of second harmonic
////                6. Amplitude of third harmonic
////                7. Phase angle of third harmonic
////                8. Amplitude of fourth harmonic
////                9. Phase angle of fourth harmonic
////               10. Amplitude of fifth harmonic
////               11. Phase angle of fifth harmonic
////               12. Amplitude of sixth harmonic
////               13. Phase angle of sixth harmonic
////       zz      weighting function
////
////     OUTPUT:
////
////       average parameter for the station
////
////     p00        Transition probability of a dry day following a dry day
////     p10        Transition probability of a dry day following a wet day
////     beta       The mean of the smaller exponential distributions
////     delta      The mean of the larger exponential distributions
////     albar      The weighting function
////
////     WORKING ARRAY:
////
////       th
////       za
////        r        distance between center and stations
////
////#######################################################################
////
////     Variable Declarations.
////
////      n nh   (Aqual coefficients array size (n by,nh) )
////
////#######################################################################
////
    private static final int n = 4;
    private static final int nh = 7;
    // common -- filname
    static String site;
    static String temp;
    static String matrix;
    static String gem6out;
    static String filename;
    static Scanner inBR = new Scanner(System.in);

    static void input_parameters() {        // start method after variable declarations
        // that are in common blocks

        double[] th = new double[idyr];
        double[][] za = new double[n][nh];
        String[] staname = new String[MAX_NUM_STA];
        String[] stna = new String[3];
        String char1;
        String AguaHeader;

        ////#######################################################################
////
//// VarName filename.extension PROGRAM   Comments
//// site     Location.site - AGUA generated precipitation file
//// temp        "    .temp - GENPAR generated annual means (coefficients)
//// matrix      "    .max  - GENPAR generated A &amp; B matrix (March 1)

//// gem6out     "    .out  - GEM6  X num years of generated weather (Jan1)
////             "    .txt  - SAMSON (CD) --&gt; SAMFORMT--&gt; *.txt file
//// filename - location name of file (i.e. Boise for Boise.txt)
////  *.site  contains coeficients for stations
////#######################################################################
////       Read data from file "sites.dat" by station name,latitude,
////       longitude and parameter matrix. Calculate the distance between
////       center (provided by user) and each station, store it in array
////       r(m).
////#######################################################################
//
        for (int ddx = 0; ddx &lt; th.length; ddx++) {
            th[ddx] = 0;
        }

        albar = 0.;

////ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

        System.out.println("Enter Location file name (no extension-no blanks) ");
//        filename = inBR.nextLine();
//        if (filename.trim().length() == 0) {
//            filename = "OR_Astoria";
//            filename = "KS_Dodge";
//        }
        filename = "C:/weps.svn/GEM/data/states/Id/ID_Boise/ID_Boise";

        System.out.println("input &gt;&gt; " + filename);
//         write (6,*)''
//         read(5,*) filename
////
//// Fortran pads the character string 'filename' with blanks.  Check for
//// the 1st blank character. This is the end of the file name string.  Look
//// for 'filename.site' etc.
////
//         ilen = index(filename,' ')              ! string position of 1st ' '
////
//         site=filename(1:ilen-1)// '.site'       ! from AGUA46 - ppt output
        site = filename + ".site";
////
//         temp = filename(1:ilen-1) // '.temp'    ! from GENPAR - annual means
        temp = filename + ".temp";
////
//         matrix = filename(1:ilen-1) // '.max'   ! GENPAR - A&amp;B matrix
        matrix = filename + ".max";
////
//         gem6out = filename(1:ilen-1) // '.out'  ! GEN6 output file
        gem6out = filename + ".out2";
//
////
////#######################################################################
////       Begining the reading process
////#######################################################################
//
////cc	open(2,file='boise61a.site',status='old') ! Input SITES (Agua)
////
//	open(2,file=site,status='old')            ! Input SITES (Agua)
        Scanner siteS = null;
        try {
            siteS = new Scanner(new File(site));
        } catch (FileNotFoundException fnfe) {
            System.err.println("" + fnfe);
            System.exit(1);
        }
//
////cc	plat1=43.+43./60.                         !HARD CODED for BOISE
////
        AguaHeader = siteS.nextLine();
        double plat1 = siteS.nextDouble();      // get the station latitude
        siteS.nextLine();         // skip rest of line
//        read(2,*)AguaHeader
//        read(2,*)plat1
//
        for (int jdx = 0; jdx &lt; nh; jdx++) {
            for (int idx = 0; idx &lt; n; idx++) {
                za[idx][jdx] = siteS.nextDouble();
            }
        }
        albar = siteS.nextDouble();
        siteS.close();
        System.out.println("" + albar);
//	read(2,*) ((za(i1,i2),i1=1,n),i2=1,nh)    ! Agua coefficients
//	read(2,*) albar                 ! read weighting coefficient
//	close(2)
////
////#######################################################################
////       Calculate the seasonal variations (see orange book Eq. (3) Pg 2
////#######################################################################
////
        double t = 2 * pi / idyr;
//	t=2*pi/float(idyr)
//
        for (int i1 = 0; i1 &lt; n; i1++) {
//        DO 23 i1=1,n
            for (int i = 0; i &lt; idyr; i += 4) {
//	DO 24 i=1,idyr,4
                th[i] = za[i1][0];
//          th(i)=za(i1,1)
                if (za[i1][1] &gt; 0.0001) {
                    th[i] += za[i1][1] * Math.sin(t * (i + 1) + za[i1][2]);
                }
//	  IF(za(i1,2).gt.0.0001)
//     !      th(i)=th(i)+za(i1,2)*sin(t*i+za(i1,3))
                if (za[i1][3] &gt; 0.0001) {
                    th[i] += za[i1][3] * Math.sin(t * (i + 1) * 2 + za[i1][4]);
                }
//	  IF(za(i1,4).gt.0.0001)
//     !      th(i)=th(i)+za(i1,4)*sin(t*i*2+za(i1,5))
                if (za[i1][5] &gt; 0.001) {
                    th[i] += za[i1][5] * Math.sin(t * (i + 1) * 3 + za[i1][6]);
                }
//	  IF(za(i1,6).gt.0.0001)
//     !      th(i)=th(i)+za(i1,6)*sin(t*i*3+za(i1,7))
////
////cccccccccccccccccccccccccccccccccccccccccccc
//
                if (i &gt; 0) {
                    for (int k = 0; k &lt; 3; k++) {
                        th[i - k - 1] = th[i - 4] + (th[i] - th[i - 4]) * (3 - k) / 4;
                    }
                }
//	  IF(i.gt.1) THEN
//	    DO 27 k=1,3
//	      th(i-k)=th(i-4)+(th(i)-th(i-4))*(4-k)/4.
//27          CONTINUE
//          END IF
            }
//24      CONTINUE
            for (int i = 0; i &lt; idyr; i++) {
                switch (i1) {
                    case 0:
                        p00[i] = th[i];
                        break;
                    case 1:
                        p10[i] = th[i];
                        break;
                    case 2:
                        beta[i] = th[i];
                        break;
                    case 3:
                        delta[i] = (th[i] - albar * beta[i]) / (1 - albar);
                }

            }
//        DO 26 i=1,idyr
//	  IF(i1.eq.1) p00(i)=th(i)
//	  IF(i1.eq.2) p10(i)=th(i)
//	  IF(i1.eq.3) beta(i)=th(i)
//	  IF(i1.eq.4) delta(i)=(th(i)-albar*beta(i))/(1.-albar)
//                           ! the calculation of delta is according to orange
//                           ! book Eq. (6), page 3
//26      CONTINUE
        }
//23      CONTINUE
//
        System.out.println("Do you want to calculate average annual precipitation (Y/N)? ");
        //todo: String res = inBR.nextLine().trim();
        String res = "n";
        double eer;
        double dedp1;
//        write(6,'(/ 22(/5x,a)//)')
//     :'                                                              ',
//     :'                                                              ',
//     :'Do you want to calculate average annual precipitation (Y/N)?  '
//        read(5,'(a1)') char1
        if (res.equalsIgnoreCase("n")) {
//            System.out.println("hit 'return'/'enter' key please.");
//            inBR.nextLine();
        } else {
//	IF(char1.eq.'N'.or.char1.eq.'n') THEN
//          write(6,'(5x,a)') 'hit "return"/"enter" key please.'
//	  read(5,*)
//	  GOTO 30
//        END IF
            annual_rainfall();
//        CALL ANNUAL_RAINFALL      ! average annual rainfall
            System.out.println("If you know average annual precipitation please enter: ");
            System.out.println("(inches) in form of xxx.xx.                            ");
            double avp = inBR.nextDouble();
//        write(6,'(/ 22(/5x,a)//)')
//     :'If you know average annual precipitation please enter:        ',
//     :'(inches) in form of xxx.xx.                                   '
//        read(5,'(f10.5)') avp
////       read(5,*) avp
            if (avp &lt;= 0.0) {
//        IF(avp.eq.0.0) GOTO 30
//
////#######################################################################
////       A NEWTON-RAPHSON iterative procedure is used to adjust
////       alfarbar and P10 bar, so that the theoretical mean (from
////       interpolation isohyetal map) within +/- 0.1 percent of
////       the known average annual precipitation
////#######################################################################
                while (true) {
                    eer = avp - anp;
//2       eer=avp-anp
                    if (Math.abs(eer) &lt; 0.001) {
                        break;
                    }
//        IF(abs(eer).lt.0.001) GOTO 30
                    dedp1 = (1 - za[0][0]) / Math.pow(1 + za[1][0] - za[0][0], 2) * (albar * za[2][0] + (1 - albar) * za[3][0]) * idyr;
//	dedp1=(1-za(1,1))/(1+za(2,1)-za(1,1))**2
//     ! *(albar*za(3,1)+(1-albar)*za(4,1))*idyr
//		                   ! orange book Eq. (25), page 11
                    double deda = -(1 - za[0][0]) * (za[2][0] - za[3][0]) * idyr / (1 + za[1][0] - za[0][0]);
//	deda=-(1-za(1,1))*(za(3,1)-za(4,1))*idyr/(1+za(2,1)
//     ! -za(1,1))
//		                   ! orange book Eq. (26), page 12
                    double da = 0.33 * eer / deda;
//	da=-0.33*eer/deda          ! orange book Eq. (23), page 11
                    double dp1 = -0.33 * eer / dedp1;
//	dp1=-0.33*eer/dedp1        ! orange book Eq. (24), page 11
//
                    za[1][0] = za[1][0] + dp1;
//	za(2,1)=za(2,1)+dp1
                    albar -= da;
//	albar=albar+da
                    for (int i = 0; i &lt; idyr; i++) {
                        p10[i] += dp1;
                    }
//	DO i=1,idyr
//	  p10(i)=p10(i)+dp1
//        END DO
                    annual_rainfall();
//
//	CALL ANNUAL_RAINFALL
//	GOTO 2
                }
            }
        }
//30      CONTINUE
        System.out.println("");
        System.out.println("");
        System.out.println("Do you want to estimate the probability of x or less inches of");
        System.out.println("precipitation over N days ? (Y/N)                             ");
//        write(6,'(/ 22(/5x,a)//)')
//     :'                                                              ',
//     :'                                                              ',
//     :'Do you want to estimate the probability of x or less inches of',
//     :'precipitation over N days ? (Y/N)                             '
        //todo: res = inBR.nextLine().trim();
        res = "n";
//        read(5,'(a1)') char1
        if (res.equalsIgnoreCase("y")) {
//	IF(char1.eq.'Y'.or.char1.eq.'y') THEN
            total_rainfall();
//          CALL TOTAL_RAINFALL
//			     ! for n-day total precipitation
        } else {
//	ELSE
//            System.out.println("hit 'return'/'enter' key please.");
////          write(6,'(5x,a)') 'hit "return"/"enter" key please.'
//            inBR.nextLine();
//	  read(5,*)
//
        }
//        END IF
        System.out.println("Do you want to start the simulation process? (Y/N)");
//        write(6,'(5x,a)')
        //todo: res = inBR.nextLine().trim();
        res = "y";
//     :'Do you want to start the simulation process? (Y/N)            '
//        read(5,'(a1)') char1
        if (res.equalsIgnoreCase("Y")) {
//	IF(char1.eq.'Y'.or.char1.eq.'y') THEN
////
////C        write(6,'(5x,a)')
////C     :'Simulated data will be saved in file "Yourfilename.out" '
////
            simulation(plat1);
//	  CALL SIMULATION(plat1)
//                      ! for simulation of p, tmax, tmin and radiation
        }
//        END IF
//31      RETURN
//	END
    }
//
////     ##################################################################
////     ##################################################################
////     ######                                                      ######
////     ######             SUBROUTINE ANNUAL_RAINFALL               ######
////     ######                                                      ######
////     ######                     (c) 1994                         ######
////     ######                                                      ######
////     ##################################################################
////     ##################################################################
//

    static void annual_rainfall() {
//        SUBROUTINE ANNUAL_RAINFALL
////
////#######################################################################
////
////     PURPOSE:
////
////     Calculate average annual (365 days) rainfall
////
////#######################################################################
////
////     AUTHOR: C.L.Hanson,K.A.Cumming,D.A.Woolhiser and C.W.Richardson
////
////     Modify history:
////     Yunyun Lu (November,1994)
////     transfer program from QBASIC to FORTRAN 77
////
////     Y. Lu (Dec. 1994)
////     added the full documentation
////
////#######################################################################
////
////     INPUT:
////
////     p00        Transition probability of a dry day following a dry day
////     p10        Transition probability of a dry day following a wet day
////     beta       The mean of the smaller exponential distributions
////     delta      The mean of the larger exponential distributions
////     albar      The weighting function
////
////    OUTPUT:
////
////     anp        Expected annual precipitation (inches)
////     sum1       Number of wet days
////
////    WORKING ARRAY:
////
////     pwet       unconditional probability of day n being wet.
////
////#######################################################################
////
////     Variable Declarations.
////
////#######################################################################
////
//
        int n = 4;
        int nh = 13;
        int idyr = 365;
//	parameter(n=4,nh=13,idyr=365)
//        common /const1/pi
//	common /sited2/delta(idyr),p00(idyr),p10(idyr),beta(idyr),albar
//	common /rainf1/sum1,anp
//	character*25 site,temp,matrix,gem6out,filename
//        common/filname/ site,temp,matrix,gem6out,filename
////%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
////
////     Beginning of executable code...
////
////%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
//
        double pwet = (1 - p00[idyr - 1]) / (1 + p10[idyr - 1] - p00[idyr - 1]);
//        pwet=(1-p00(idyr))/(1+p10(idyr)-p00(idyr))
//			       ! orange book Eq. (4), page 2
        double sum1 = 0;
        double sum2 = 0;
//	sum1=0.
//	sum2=0.
        for (int i = 0; i &lt; idyr; i++) {
//	DO 30 i=1,idyr
            pwet *= (1 - p10[i]) + (1 - pwet) * (1 - p00[i]);
            sum1 += pwet;
            sum2 += pwet * (albar * beta[i] + (1 - albar) * delta[i]);
//	  pwet=pwet*(1-p10(i))+(1-pwet)*(1-p00(i))
//	  sum1=sum1+pwet
//	  sum2=sum2+pwet*(albar*beta(i)+(1-albar)*delta(i))
        }
//30      CONTINUE
        anp = sum2 + sum1 * 0.01;
//	anp=sum2+sum1*0.01
        System.out.printf("Expected annual precipitation %6.2f inches\n", anp);
//	print*,'Expected annual precipitation', anp,'inches'
        System.out.printf("Number of wet days %f\n", sum1);
//	print*,'Number of wet days=',sum1
//
//        RETURN
    }
//        END
//
////     ##################################################################
////     ##################################################################
////     ######                                                      ######
////     ######              SUBROUTINE SIMULATION                   ######
////     ######                                                      ######
////     ######                     (c) 1994                         ######
////     ######                                                      ######
////     ##################################################################
////     ##################################################################
//
    static PrintWriter unit21 = null;
    static PrintWriter unit22 = null;
    static PrintWriter unit23 = null;

    static void simulation(double plat1) {
//        SUBROUTINE SIMULATION(PLAT1)
//
////#######################################################################
////
////     PURPOSE:
////
////     Simulate M years of daily rainfall data beginning Mar.1
////
////#######################################################################
////
////     AUTHOR: C.L.Hanson,K.A.Cumming,D.A.Woolhiser and C.W.Richardson
////
////     Modify history:
////     Yunyun Lu (November,1994)
////     transfer program from QBASIC to FORTRAN 77
////
////     Y. Lu (Dec. 1994)
////     added the full documentation
////
////     Greg Johnson (Dec. 12, 1994)
////     added the option to correct the simulated temperature in these
////     two cases: 1) tmax(day i) &lt; tmin(day i-1);
////                2) tmin(day i) &gt; tmax(day i-1);
////
////     to         1) tmax(day i) = tmin(day i-1);
////                2) tmin(day i) = tmax(day i-1).
////#######################################################################
////
////     INPUT:
////
////     p00        Transition probability of a dry day following a dry day
////     p10        Transition probability of a dry day following a wet day
////     beta       The mean of the smaller exponential distributions
////     delta      The mean of the larger exponential distributions
////     albar      The weighting function
////
////     OUTPUT:
////
////     p          Simulated M years of daily rainfall
////     tmax       Simulated M years of daily max. temperature
////     tmin       Simulated M years of daily min. temperature
////     rad        Simulated M years of daily total radiation
////
////   WORKING ARRAY
////
////     tmad       mean of tmax -- dry days
////     txs        standard deviation of tmax -- dry days
////     txm1       mean of tmax -- wet days
////     txs1       standard deviation of tmax -- wet days
////     tnm        mean of tmin -- wet and dry days
////     tns        standard deviation of tmin -- wet and dry days
////     rm0        mean of radiation -- dry days
////     rs0        standard deviation of radiation -- dry days
////     rm1        mean of radiation -- wet days
////     rs1        standard deviation of radiation -- wet days
////     rm         max. radiation along a given latitude
////    chipre      last step of chi(n), = chi(n-1)
////    chi         a vector having elements that are the standardized residuals
////                tmax, tmin, and radiation.
////
////#######################################################################
////
////     Variable Declarations.
////
////#######################################################################
////
//
// 	parameter(idyr=365,ii=360)        ! ii=number of stations
//        common /const1/pi
//	common /sited2/delta(idyr),p00(idyr),p10(idyr),beta(idyr),albar
//	common /tem11/tmad(idyr),txs(idyr),txm1(idyr),txs1(idyr),
//     !dmad(idyr),dxs(idyr),dxm1(idyr),dxs1(idyr),
//     !wmad(idyr),wxs(idyr),wxm1(idyr),wxs1(idyr),
//     !tnm(idyr),tns(idyr),tnm1(idyr),tns1(idyr),
//     !rm0(idyr),rs0(idyr),rm1(idyr),rs1(idyr),
//     !rm(idyr),chipre(6)
//        common /temp2/p,tmax,tmin,dew,wind,rad
//        common /temp3/fran(12)
//        common /temp4/a(6,6,12),b(6,6,12)
//        integer notdone,ycount,sw,mday(12)
//        character stname(ii)*20,char1*4,char2*4
        int[] mday = {1, 32, 62, 93, 123, 154, 185, 215, 246, 276, 307, 338};
//	data mday/1,32,62,93,123,154,185,215,246,276,307,338/
        String[] mname = {"MAR", "APR", "MAY", "JUN", "JUL", "AUG", "SEP",
            "OCT", "NOV", "DEC", "JAN", "FEB"};
//        character*4 mname(12)
//	data mname/'MAR','APR','MAY','JUN','JUL','AUG',
//     !'SEP','OCT','NOV','DEC','JAN','FEB'/
//                                                     ! start from March 1 !!!
////
//	character*25 site,temp,matrix,gem6out,filename
//        common/filname/site,temp,matrix,gem6out,filename
////%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
////
////     Beginning of executable code...
////
////%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
//
        initDaysToDate();
        System.out.println("How many years do you want to simulate");
//        write(6,'(/ (/5x,a)//)')
//     :'How many years do you want to simulate'
        //todo: int numYrs = inBR.nextInt();
        int numYrs = 30;
//        read(5,'(i5)') m
        numYrs++;
//        m=m+1                          !simulate additional year
        System.out.println("Please give the seed for random process (-32768 -- +32768).");
//        write(6,'(/ (/5x,a)//)')
//     :'Please give the seed for random process (-32768 -- +32768).'
        //todo: int iseed = inBR.nextInt();
        int iseed = 1;
//        read(5,'(i8)') iseed
        int leap = 0;
//	leap=0
        int ycount = 0;
//	ycount=0
        System.out.println("Do you want leap years (Y/N) ?");
//        write(6,'(/ (/5x,a)//)')
//     :'Do you want leap years (Y/N) ?'
        //todo: String res = inBR.nextLine().trim();
        String res = "n";
//	read(5,'(a1)') char1
        if (res.equalsIgnoreCase("y")) {
//	IF(char1.eq.'Y'.or.char1.eq.'y') THEN
            System.out.println("Which of the first 4 years of simulation will be a leap year?");
            leap = inBR.nextInt();
//          write(6,'(/ (/5x,a)//)')
//     :'Which of the first 4 years of simulation will be a leap year?'
//	  read(5,'(I2)') leap
        } else {
//        ELSE
            System.out.println("hit 'return'/'enter' key please.");
            //todo: inBR.nextLine();
//          write(6,'(5x,a)') 'hit "return"/"enter" key please.'
//	  read(5,*)
        }
        leap = 4;
//	END IF
//
        double ut = (1 - p00[idyr - 1]) / (1 + p10[idyr - 1] - p00[idyr - 1]);
//	ut=(1-p00(idyr))/(1+p10(idyr)-p00(idyr))
//			  ! unconditional probability of day n being
//			  ! wet. (orange book,Eq.(4), Page 4)
        double u = ran3(iseed);
//	u=ran3(iseed)     ! random number
        boolean wetDayFlg = u &lt; ut;
//	nw=0              ! preceding day is dry
//	IF(u.lt.ut) nw=1  ! proceding day is wet, if u is less than ut
//
////
        try {
            unit21 = new PrintWriter(new File(filename + "-x1.out1"));
        } catch (FileNotFoundException ex) {
            Logger.getLogger(GEM6.class.getName()).log(Level.SEVERE, null, ex);
            System.exit(1);
        }
//	  open(unit=21,file='x1.out',status='replace')
//
        read_123imation(plat1);
        try {
//	  CALL READ_123IMATION(plat1)
////
////cc	  open(unit=22,file='gem6wx.out',status='replace')  !GEM6 output
////                           (example:  gem6out ='Boise.out')
            unit22 = new PrintWriter(new File(gem6out));
            unit22.println();
        } catch (FileNotFoundException ex) {
            Logger.getLogger(GEM6.class.getName()).log(Level.SEVERE, null, ex);
        }
//	  open(unit=22,file=gem6out,status='replace')       !GEM6 output
////CCCC
////CCCC  need to add header information
////CCCC
        unit22.println("YR MO DAY prec      tmax       tmin       dewpoint    wind      rad");
        int mmd = 1;
//        int ns = 0;
//          write(22,'(/ (1x,a14,6x,a15,7x,a16,6x,a3))')
//     :'YR MO DAY prec','tmax       tmin','dewpoint    wind','rad'
////
//
        int sw = 0;
        for (int i = 0; i &lt; numYrs; i++) {
//	DO 10 i=1,m
//	  sw=0
            if (leap != 0) {
//	  IF(leap.eq.0) GOTO 30
                ycount++;
//	  ycount=ycount+1
                if (ycount == leap) {
                    sw = 1;
                    ycount = 0;
                    leap = 4;
                }
//	  IF(ycount.ne.leap) GOTO 30
//	  sw=1
//	  ycount=0
//	  leap=4
//			  ! from now on every 4th year will be a leap year
////#######################################################################
////         begin the simulation process
////#######################################################################
//
            }
//30        continue
            double pretmin = -100.;
            double pretmax = 200.;
//	  pretmin=-100.
//	  pretmax=200.
            int notdone = 0;
//	  notdone=0

            for (int j = 0; j &lt; idyr; j++) {
                do {
                    u = ran3(iseed);
                    // determine if today is wet based upon if yesterday was wet
                    if (wetDayFlg) {
                        wetDayFlg = u &gt;= p10[j];     // yesterday was wet day
                    } else {
                        wetDayFlg = u &gt;= p00[j];     // yesterday was dry day
                    }
                    if (wetDayFlg) {
                        u = ran3(iseed);
                        if (u &lt; albar) {
                            u = ran3(iseed);
                            p = -beta[j] * Math.log(u) + 0.01;
                        } else {
                            u = ran3(iseed);
                            p = -delta[j] * Math.log(u) + 0.01;
                        }
                    } else {
                        p = 0;
                    }



//	  DO 20 j=1,idyr
//40          u=ran3(iseed)
//	    IF(nw.eq.0) GOTO 60
////#######################################################################
////           if proceding day is dry, use a new generated random number
////           and do following:
////#######################################################################
//	    IF(u.lt.p10(j)) THEN
//	      GOTO 90     ! u &lt; p10(n), day n is dry
//            ELSE
//	      GOTO 70     ! u &gt;= p10(n), day n is wet
//            END IF
////#######################################################################
////           if proceding day is wet, use a new generated random number
////           and do following:
////#######################################################################
//60          IF(u.lt.p00(j)) GOTO 90
//			  ! u &lt; p00(n) day n is dry
////#######################################################################
////           for wet day the mixed exponential distribution is:
////#######################################################################
//70          nw=1
//	    u=ran3(iseed)   ! new random number
//	    IF(u.lt.albar) GOTO 80
//			  ! u &lt; alpha(n): see orange book Eq.(18),Page 7
//	    u=ran3(iseed)
//	    p=-delta(j)*alog(u)+0.01
//	    GOTO 100
//			  ! u &gt;= alpha(n): see orange book Eq.(19),Page 7
//80          u=ran3(iseed)
//	    p=-beta(j)*alog(u)+0.01
//	    GOTO 100
////#######################################################################
////           for day is dry
////#######################################################################
//90          p=0.
//	    nw=0
//	    GOTO 100
//100         CONTINUE

                    int ji = 0;
                    for (int mo = 0; mo &lt; 11; mo++) {
                        if ((j + 1) &gt;= mday[mo] &amp;&amp; (j + 1) &lt; mday[mo + 1]) {
                            ji = mo;
                        }
                    }
//            DO 120 mo=1,11
//	      IF(j.ge.mday(mo).and.j.lt.mday(mo+1)) ji=mo
//120         CONTINUE
                    if (j &gt;= mday[11] - 1) {
                        ji = 11;
                    }
//	    IF(j.ge.mday(12)) ji=12
//            mmd=j-mday(ji)+1
                    mmd = j - mday[ji] + 1;
                    if (notdone == 1) {
                        mmd++;
                        notdone = 0;
                    }
//	    IF(notdone.eq.1) THEN
//	      mmd=mmd+1
//	      notdone=0
//            END IF
                    int ijk = i;
//	    ijk=i
                    if (ji &gt;= 10) {
                        ijk = i + 1;
                    }

//	    if(ji.ge.11) ijk=i+1
////
////cc	    IF(char1.eq.'Y'.or.char1.eq.'y') THEN
////cc	      IF(char2.eq.'Y'.or.char2.eq.'y') THEN
////
                    tmax_tmin_123(j, /*ns,*/ iseed, wetDayFlg, ji);
//	      CALL TMAX_TMIN_123(j,ns,iseed,nw,ji)
//                          ! simulate tmax,tmin,rad

                    dew -= 100.;
                    wind = (wind &gt; 0) ? wind : -wind;
//	      dew=dew-100.
//	      wind=abs(wind)
////
////cc	      ELSE
////cc              CALL TMAX_TMIN_RAD(j,ns,iseed,nw)
//                          ! simulate tmax,tmin,rad
////cc              END IF
////
                    tmax -= 100;
                    tmin -= 100;
                    tmin = (tmin &lt; pretmax) ? tmin : pretmax;
                    tmax = (tmax &gt; pretmin) ? tmax : pretmin;
                    pretmin = tmin;
                    pretmax = tmax;
//	      tmax=tmax-100.
//	      tmin=tmin-100.
//	      tmin=min(tmin,pretmax)
//	      tmax=max(tmax,pretmin)
//	      pretmin=tmin
//	      pretmax=tmax
////
////cc	    END IF
////
//110         CONTINUE
                    int nmo = (ji &lt; 10) ? ji + 3 : ji - 9;
//            nmo = ji+2
//	    IF(ji.gt.10) nmo = ji-10        ! nmo -- # of month
////
////cc	    IF(char2.eq.'Y'.or.char2.eq.'y') THEN
////
//                int[] date = DayToDate.dayInYearToDate(j, false, 0);
                    unit21.printf("%3d%3d%5s%3d%6.2f%11.2f%11.2f%11.2f%11.2f%11.2f\n",
                            ijk + 1, nmo, mname[ji], mmd + 1, p, tmax, tmin, dew, wind, rad);
//              write(21,'(i3,i3,a5,i3,f6.2,5f11.2)')ijk,nmo,mname(ji),
//     !mmd,p,tmax,tmin,dew,wind,rad
                    if ((ijk == 0) || (ijk == numYrs)) {
                    } else {
//                    if ((mname[ji].equals("Feb")) &amp;&amp; (mmd == 29)) {
//                    } else {
                        unit22.printf("%3d%3d%3d%6.2f%11.2f%11.2f%11.2f%11.2f%11.2f\n",
                                ijk, nmo, mmd + 1, p, tmax, tmin, dew, wind, rad);
//                    }
                    }
//	      if(ijk.eq.1.or.ijk.eq.(m+1))goto 159
//	      if(mname(ji).eq.' Feb'.and.mmd.eq.29) goto 159
//              write(22,'(i3,i3,i3,f6.2,5f11.2)')ijk-1,nmo,
//     !mmd,p,tmax,tmin,dew,wind,rad
//159           continue
////
                    if (j &lt; (idyr - 1) || sw == 0) {
                        break;
                    } else {
                        sw = 0;
                        notdone = 1;
                    }
                } while (true);
//            IF(j.lt.idyr) GOTO 20
//	    IF(sw.eq.0)  GOTO 20
//	    sw=0
//	    notdone=1
//	    GOTO 40
            }
//20        CONTINUE
        }
//10      CONTINUE
        unit21.close();
        unit22.close();
//      close(22)
//	RETURN
    }
//	END
//
//
//
//

    /**
    // tmax_tmin_rad is never used
    //////     ##################################################################
    //////     ##################################################################
    //////     ######                                                      ######
    //////     ######              SUBROUTINE TMAX_TMIN_RAD                ######
    //////     ######                                                      ######
    //////     ######                     (c) 1994                         ######
    //////     ######                                                      ######
    //////     ##################################################################
    //////     ##################################################################
    ////
    //
    //    static void tmax_tmin_rad(int j, int ns, int iseed, boolean wetDayFlg) {
    ////        SUBROUTINE TMAX_TMIN_RAD(j,ns,iseed,nw)
    //////
    //////#######################################################################
    //////
    //////     PURPOSE:
    //////
    //////    This program simulates max and min daily temp and radiation
    //////
    //////#######################################################################
    //////
    //////     AUTHOR: C.L.Hanson,K.A.Cumming,D.A.Woolhiser and C.W.Richardson
    //////
    //////     Modify history:
    //////     Yunyun Lu (November,1994)
    //////     transfer program from QBASIC to FORTRAN 77
    //////
    //////     Y. Lu (Dec. 1994)
    //////     added the full documentation
    //////
    //////#######################################################################
    //////
    //////     INPUT:
    //////
    //////     tmad       mean of tmax -- dry days
    //////     txs        standard deviation of tmax -- dry days
    //////     txm1       mean of tmax -- wet days
    //////     txs1       standard deviation of tmax -- wet days
    //////     tnm        mean of tmin -- wet and dry days
    //////     tns        standard deviation of tmin -- wet and dry days
    //////     rm0        mean of radiation -- dry days
    //////     rs0        standard deviation of radiation -- dry days
    //////     rm1        mean of radiation -- wet days
    //////     rs1        standard deviation of radiation -- wet days
    //////     rm         max. radiation along a given latitude
    //////    chipre      last step of chi(n), =chi(n-1)
    //////
    //////     OUTPUT:
    //////
    //////     p          Simulated M years of daily rainfall
    //////     tmax       Simulated M years of daily max. temperature
    //////     tmin       Simulated M years of daily min. temperature
    //////     rad        Simulated M years of daily total radiation
    //////
    //////#######################################################################
    //////
    //////     Variable Declarations.
    //////
    //////#######################################################################
    //////
    ////	parameter(idyr=365)
    ////        common /const1/pi
    ////	common /tem11/tmad(idyr),txs(idyr),txm1(idyr),txs1(idyr),
    ////     !dmad(idyr),dxs(idyr),dxm1(idyr),dxs1(idyr),
    ////     !wmad(idyr),wxs(idyr),wxm1(idyr),wxs1(idyr),
    ////     !tnm(idyr),tns(idyr),tnm1(idyr),tns1(idyr),
    ////     !rm0(idyr),rs0(idyr),rm1(idyr),rs1(idyr),
    ////     !rm(idyr),chipre(6)
    ////        common /temp2/p,tmax,tmin,dew,wind,rad
    ////        common /temp3/fran(12)
    ////	real    a(3,3),b(3,3)
    ////	data a/.567,.253,-.006,.086,.504,-.039,-.002,-.05,.244/
    //        double[][] a = {{.567,.253,-.006},{.086,.504,-.039},{-.002,-.05,.244}};
    ////	data b/.781,.328,.238,.0,.637,-.341,.0,.0,.873/
    //        double[][] b = {{.781,.328,.238},{.0,.637,-.341},{.0,.0,.873}};
    ////	real    rd(3),rrd(3),epsilon(3),chi(3)
    //        double[] rd = new double[3];
    //        double[] rrd = new double[3];
    //        double[] epsilon = new double[3];
    //        double[] chi = new double[3];
    ////	character*25 site,temp,matrix,gem6out,filename
    ////        common/filname/site,temp,matrix,gem6out,filename
    //////%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    //////
    //////     Beginning of executable code...
    //////
    //////%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    ////
    ////	rmm=rm0(j)
    //        double rmm = rm0[j];
    ////	rs=rs0(j)
    //        double rs = rs0[j];
    ////	txxm=tmad(j)
    //        double txxm = tmad[j];
    ////	txxs=txs(j)
    //        double txxs = txs[j];
    ////        IF (nw.eq.1) THEN
    //        if (wetDayFlg) {
    ////	  rmm=rm1(j)
    //            rmm = rm1[j];
    ////	  rs=rs1(j)
    //            rs = rs1[j];
    ////	  txxm=txm1(j)
    //            txxm = txm1[j];
    ////	  txxs=txs1(j)
    //            txxs = txs1[j];
    ////        END IF
    //        }
    ////	IF(ns.le.1) THEN
    ////	  DO k=1,5,2
    ////	    u=ran3(iseed)
    ////	    u1=ran3(iseed)
    ////	    fran(k)=sqrt(-2*alog(u))*cos(2*pi*u1)
    ////	    fran(1+k)=sqrt(-2*alog(u))*sin(2*pi*u1)
    ////          END DO
    ////	  ns=3
    ////        ELSE
    ////	  ns=0
    ////        END IF
    //        if (ns &lt;= 1) {
    //            for (int kdx = 0; kdx &lt; 5; kdx += 2) {
    //                double u = ran3(iseed);
    //                double u1 = ran3(iseed);
    //                fran[kdx] = sqrt(-2*Math.log(u) * Math.cos(2*pi*u1));
    //                fran[kdx+1] = sqrt(-2*Math.log(u) * Math.sin(2*pi*u1));
    //            }
    //        }
    //////#######################################################################
    //////       for first day using last three number in the fran() array, THEN
    //////       for next day using the first three number in fran() array.
    //////       Iterately do it back and forewards.
    //////       integer ns control which three number will be used.
    //////#######################################################################
    ////
    ////	DO k=1,3
    ////	  epsilon(k)=fran(k+ns)
    ////	  rd(k)=0.
    ////	  rrd(k)=0.
    ////        END DO
    //        for (int kdx = 0; kdx &lt; 3; kdx++) {
    //            epsilon[kdx] = fran[kdx+ns];
    //            rd[kdx] = 0;
    //            rrd[kdx] = 0;
    //        }
    ////	DO k=1,3
    ////	  DO l=1,3
    ////	    rd(k)=rd(k)+b(k,l)*epsilon(l)
    ////                                     ! see orange book Eq.(13) page 5
    ////	    rrd(k)=rrd(k)+a(k,l)*chipre(l)
    ////          END DO
    ////        END DO
    //        for (int kdx = 0; kdx &lt; 3; kdx++) {
    //            for (int ldx = 0; ldx &lt; 3; ldx++) {
    //                rd[kdx] += b[kdx][ldx] * epsilon[ldx];
    //                rrd[kdx] += a[kdx][ldx] * chipre[ldx];
    //            }
    //        }
    ////	DO l=1,3
    ////	  chi(l)=rd(l)+rrd(l)
    ////	  chipre(l)=chi(l)            ! at the END of calculation
    ////                                      ! store chi as next run's chipre.
    ////        END DO
    //        for (int ldx = 0; ldx &lt; 3; ldx++) {
    //            chi[ldx] = rd[ldx] + rrd[ldx];
    //            chipre[ldx] = chi[ldx];
    //        }
    ////	tmax=chi(1)*txxs+txxm
    //        tmax = chi[0] * txxs + txxm;
    ////	tmin=chi(2)*tns(j)+tnm(j)
    //        tmin = chi[1] * tns[j] + tnm[j];
    //////#######################################################################
    //////       Adjust for tmax, tmin, radiation
    //////#######################################################################
    ////
    ////	IF(tmin.gt.tmax) THEN   ! tmax should never be smaller than tmin
    ////	  tmaxt=tmax
    ////	  tmax=tmin
    ////	  tmin=tmaxt
    ////        END IF
    //        if (tmin &gt; tmax) {
    //            double tmaxt = tmax;
    //            tmax = tmin;
    //            tmin = tmaxt;
    //        }
    ////	rad=chi(3)*rs+rmm       ! rad will not smaller than 0.1*max.radation
    ////	rmin=.1*rm(j)           ! and not larger than max. radiation
    ////	rad=max(rad,rmin)
    ////	rad=min(rad,rm(j))
    //        rad = chi[2] * rs + rmm;
    //        double rmin = .1 * rm[j];
    //        rad = (rad &gt; rmin) ? rad : rmin;
    //        rad = (rad &lt; rm[j]) ? rad : rm[j];
    ////	RETURN
    ////	END
    ////
    //    }*/
//
//
//
////     ##################################################################
////     ##################################################################
////     ######                                                      ######
////     ######            SUBROUTINE READ_123IMATION                ######
////     ######                                                      ######
////     ######                     (c) 1994                         ######
////     ######                                                      ######
////     ##################################################################
////     ##################################################################
//
    static void read_123imation(double lat) {
//        SUBROUTINE READ_123IMATION(lat)
////
////#######################################################################
////
////     PURPOSE:
////
////     Read temperture and radiation parameters from file
////     "temp.dat".
////     Temperture and radiation procedure is from Richardson, C.W. and
////     D.A. Wright 'WGEN: A model for generating daily weather varibles'
////     USDA-ARS-8, 1984: 83pp.
////
////#######################################################################
////
////     AUTHOR: C.L.Hanson,K.A.Cumming,D.A.Woolhiser and C.W.Richardson
////
////     Modify history:
////     Yunyun Lu (November,1994)
////     transfer program from QBASIC to FORTRAN 77
////
////     Y. Lu (Dec. 1994)
////     added the full documentation
////
////#######################################################################
////
////     INPUT:
////
////    instna      station's name.
////       w        temperatures in F, radiation in Langleys
////                1. Annual mean of tmax for dry days;
////                2. Amplitude of tmax for wet or dry days;
////                3. Annual mean of the coefficient of variation of tmax
////                   for wet or dry days;
////                4. Amplitude of the coefficient of variation of tmax
////                   for wet or dry days;
////                5. Annual mean of tmax for wet days;
////                6. Annual mean of tmax for wet or dry days;
////                7. Amplitude of tmin for wet or dry days;
////                8. Annual mean of the coefficient of variation of tmin
////                   for wet or dry days;
////                9. Amplitude of the coefficient of variation of tmin
////                   for wet or dry days;
////               10. Annual mean of radiation for dry days;
////               11. Amplitude of radiation for dry days;
////               12. Annual mean of radiation for wet days;
////               13. Amplitude of radiation for wet days;
////               14. Annual mean of the coefficient of variation of !?? see 16
////                   radiation for dry days;
////               15. Amplitude of the coefficient of variation of  !?? see 17
////                   radiation for dry days;
////               16. Annual mean of the coefficient of variation of
////                   radiation for dry days;
////               17. Amplitude of the coefficient of variation of
////                   radiation for dry days;
////
////     OUTPUT:
////
////     tmad       mean of tmax -- dry days
////     txs        standard deviation of tmax -- dry days
////     txm1       mean of tmax -- wet days
////     txs1       standard deviation of tmax -- wet days
////     tnm        mean of tmin -- wet and dry days
////     tns        standard deviation of tmin -- wet and dry days
////     rm0        mean of radiation -- dry days
////     rs0        standard deviation of radiation -- dry days
////     rm1        mean of radiation -- wet days
////     rs1        standard deviation of radiation -- wet days
////     rm         max. radiation along a given latitude
////    chipre      zero array
////
////     WORKING ARRAY:
////
////       sw       Averaged w by using stations' info within 100 miles
////
////#######################################################################
////
////     Variable Declarations.
////
////#######################################################################
////
//
//	parameter(m=360,n=4,nh=40,idyr=365)
        final int nh = 40;
//        common /const1/pi
//	common /sited2/delta(idyr),p00(idyr),p10(idyr),beta(idyr),albar
//	common /tem11/tmad(idyr),txs(idyr),txm1(idyr),txs1(idyr),
//     !dmad(idyr),dxs(idyr),dxm1(idyr),dxs1(idyr),
//     !wmad(idyr),wxs(idyr),wxm1(idyr),wxs1(idyr),
//     !tnm(idyr),tns(idyr),tnm1(idyr),tns1(idyr),
//     !rm0(idyr),rs0(idyr),rm1(idyr),rs1(idyr),
//     !rm(idyr),chipre(6)
//        common /temp4/a(6,6,12),b(6,6,12)
//	real     lat
        double[] sw = new double[nh];
//	real     sw(nh)                     !average of input parameters
//	character*25 site,temp,matrix,gem6out,filename
//        common/filname/site,temp,matrix,gem6out,filename
////%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
////
////     Beginning of executable code...
////
////%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
//
//        DO 10 k=1,nh  // not needed
//	  sw(k)=0.
//10      CONTINUE
//
//        np=0          // may not be needed
//
////cc	open(2,file='boise61a.temp',status='old')
        Scanner scan2 = null;
        try {
            scan2 = new Scanner(new File(temp));
        } catch (FileNotFoundException ex) {
            Logger.getLogger(GEM6.class.getName()).log(Level.SEVERE, null, ex);
            System.exit(1);
        }
//	open(2,file=temp,status='old')
//
        for (int i1 = 0; i1 &lt; nh; i1++) {
            sw[i1] = scan2.nextDouble();
        }
//	read(2,*) (sw(i1),i1=1,nh)    ! read annual means parameters
//
////#######################################################################
////   find the stations within 100 miles (ii&gt;1) or closest station (ii=1)
////#######################################################################
//
        scan2.close();
//30      close(2)
//
////#######################################################################
////      calculate the max. radiation along lat.
////#######################################################################
//
        for (int i = 0; i &lt; idyr; i++) {
//	DO 60 i=1,idyr
            double sd = 0.4102 * sind((360 / 365.2425) * (i - 20.25));
//	  sd=0.4102*sind((360/365.2425)*(i-21.25))   !corrects to solar day
            double ch = -tand(lat) * Math.tan(sd);
//	  ch=-tand(lat)*tan(sd)
            ch = (ch &lt; 1.0) ? ch : 1.;
//	  ch=min(ch,1.)
            ch = (ch &gt; -1.0) ? ch : -1.0;
//	  ch=max(ch,-1.)
            double h = Math.acos(ch);
//	  h=acos(ch)
            double dd = 1. + 0.0335 * sind((360 / 365.2425) * (i + 89.2));
//	  dd=1.+0.0335*sind((360/365.2425)*(i+88.2))  !
////cc	  rm(i)=889.2305*dd*(h*sind(lat)*sin(sd)+cos(sd)*sin(h))
////                                             ! ^ (cosd(lat)*cos(sd)...
            rm[i] = 889.2305 * dd * (h * sind(lat) * Math.sin(sd) + (cosd(lat) * Math.cos(sd) * Math.sin(h)));
//      rm(i)=889.2305*dd*(h*sind(lat)*sin(sd)+(cosd(lat)*cos(sd)*sin(h)))
////                                             !REF WGEN &amp; WGENPAR
////         rm(i)=rm(i)*0.9
            rm[i] *= 0.8;
//          rm(i)=rm(i)*0.8                    !Ref WGEN &amp; WGENPAR
////
        }
//60        CONTINUE
////
////#######################################################################
////         add 100 degrees to temperature related parameter
////#######################################################################
////
        for (int io = 0; io &lt; 21; io += 4) {
            sw[io] += 100;
        }
//	  do io=1,21,4
//	  sw(io)=sw(io)+100.
//	  end do
//
////#######################################################################
////      change solar radiation CV parameters if necessary
////#######################################################################
//// fixed in US Climate ARS-114 July 1994,
//// Hanson, Cumming, Woolhiser, Richardson
////       sw(35)=.224
////       sw(36)=-.084
////       sw(39)=.488
////       sw(40)=-.135
//
////#######################################################################
////      calculate daily parameter values for tmax, tmin and radiation
////#######################################################################
//
        for (int i = 0; i &lt; idyr; i++) {
//        DO 70 i=1,idyr
////                                             !360/365.2425 deg / solar day
            double dt = cosd((360 / 365.2425) * (i + 1 - 141));
            double dw = cosd((360 / 365.2425) * (i + 1 - 35 - 30));
            double dr = cosd((360 / 365.2425) * (i + 1 - 113));
//	  dt=cosd((360/365.2425)*float(i-141))   ! ?? 141,35,30,113
//	  dw=cosd((360/365.2425)*float(i-35-30))
//	  dr=cosd((360/365.2425)*float(i-113))
////
            tmad[i] = sw[0] + sw[1] * dt;    // mean of tmax -- dry days
            double xcr = sw[2] + sw[3] * dt;
            txs[i] = tmad[i] * xcr;        // standard deviation of tmax -- dry days
//	  tmad(i)=sw(1)+sw(2)*dt      ! mean of tmax -- dry days
//	  xcr=sw(3)+sw(4)*dt
//	  txs(i)=tmad(i)*xcr          ! standard deviation of tmax -- dry days
////
            txm1[i] = sw[4] + sw[5] * dt;    // mean of tmax-wet days
            xcr = sw[6] + sw[7] * dt;
            txs1[i] = txm1[i] * xcr;       // standard deviation of tmax -- wet days
//	  txm1(i)=sw(5)+sw(6)*dt      ! mean of tmax-wet days
//	  xcr=sw(7)+sw(8)*dt
//	  txs1(i)=txm1(i)*xcr         ! standard deviation of tmax -- wet days
////
            tnm[i] = sw[8] + sw[9] * dt;    // mean of tmin -- dry days
            xcr = sw[10] + sw[11] * dt;
//	  tnm(i)=sw(9)+sw(10)*dt      ! mean of tmin -- dry days
//	  xcr=sw(11)+sw(12)*dt
////
            tns[i] = tnm[i] * xcr;         // standard deviation of tmin -- dry days
            tnm1[i] = sw[12] + sw[13] * dt;  // mean of tmin -- wet days
            xcr = sw[14] + sw[15] * dt;
            tns1[i] = tnm1[i] * xcr;       // standard deviation of tmin -- wet days
//	  tns(i)=tnm(i)*xcr           ! standard deviation of tmin -- dry days
//	  tnm1(i)=sw(13)+sw(14)*dt    ! mean of tmin -- wet days
//	  xcr=sw(15)+sw(16)*dt
//	  tns1(i)=tnm1(i)*xcr         ! standard deviation of tmin -- wet days
////
            dmad[i] = sw[16] + sw[17] * dt;  // mean of dew -- dry days
            xcr = sw[18] + sw[19] * dt;
            dxs[i] = dmad[i] * xcr;        // standard deviation of dew -- dry days
//	  dmad(i)=sw(17)+sw(18)*dt    ! mean of dew -- dry days
//	  xcr=sw(19)+sw(20)*dt
//	  dxs(i)=dmad(i)*xcr          ! standard deviation of dew -- dry days
////
            dxm1[i] = sw[20] + sw[21] * dt;  // mean of dew-wet days
            xcr = sw[22] + sw[23] * dt;
            dxs1[i] = dxm1[i] * xcr;       // standard deviation of dew -- wet days
//	  dxm1(i)=sw(21)+sw(22)*dt    ! mean of dew-wet days
//	  xcr=sw(23)+sw(24)*dt
//	  dxs1(i)=dxm1(i)*xcr         ! standard deviation of dew -- wet days
////
            rm0[i] = sw[24] + sw[25] * dr;   // mean radiation -- dry days
            xcr = sw[26] + sw[27] * dr;
            rs0[i] = rm0[i] * xcr;         // standard deviation of radiation -- dry days
//	  rm0(i)=sw(25)+sw(26)*dr     ! mean radiation -- dry days
//	  xcr=sw(27)+sw(28)*dr
//	  rs0(i)=rm0(i)*xcr           ! standard deviation of radiation -- dry days
////
            rm1[i] = sw[28] + sw[29] * dr;   // mean radiation -- wet days
            xcr = sw[30] + sw[31] * dr;
            rs1[i] = rm1[i] * xcr;         //standard deviation of radiation -- wet days
//	  rm1(i)=sw(29)+sw(30)*dr     ! mean radiation -- wet days
//	  xcr=sw(31)+sw(32)*dr
//	  rs1(i)=rm1(i)*xcr           !standard deviation of radiation -- wet days
////
            wmad[i] = sw[32] + sw[33] * dw;  // mean of wind -- dry days
            xcr = sw[34] + sw[35] * dw;
            wxs[i] = wmad[i] * xcr;        // standard deviation of wind -- dry days
//	  wmad(i)=sw(33)+sw(34)*dw    ! mean of wind -- dry days
//	  xcr=sw(35)+sw(36)*dw
//	  wxs(i)=wmad(i)*xcr          ! standard deviation of wind -- dry days
////
            wxm1[i] = sw[36] + sw[37] * dw;  // mean of wind-wet days
            xcr = sw[38] + sw[39] * dw;
            wxs1[i] = wxm1[i] * xcr;       // standard deviation of wind -- wet days
//	  wxm1(i)=sw(37)+sw(38)*dw    ! mean of wind-wet days
//	  xcr=sw(39)+sw(40)*dw
//	  wxs1(i)=wxm1(i)*xcr         ! standard deviation of wind -- wet days
////
        }
//70      CONTINUE
//
////#######################################################################
////         subtract 100 degrees to temperature related parameter
////#######################################################################
//
        for (int io = 0; io &lt; 20; io += 4) {
            sw[io] -= 100;
        }
//	  do io=1,21,4
//	  sw(io)=sw(io)-100.
//	  end do
//
////#######################################################################
////      for the first run chi is zero. see orange book Eq. (13), Page 5
////#######################################################################
//
        Arrays.fill(chipre, 0);
//        DO 80 i=1,6
//	  chipre(i)=0.
//80      CONTINUE
//
////cc        open(unit=2,file='boise61a.max',status='old')    !was boi_m_6.mat
        try {
            scan2 = new Scanner(new File(matrix));
        } catch (FileNotFoundException ex) {
            Logger.getLogger(GEM6.class.getName()).log(Level.SEVERE, null, ex);
            System.exit(1);
        }
//        open(unit=2,file=matrix,status='old')  !was boi_m_6.mat
        scan2.nextLine();   // A matrix header
//	read(2,*)                              !win98 did not like .mat
        for (int jj = 0; jj &lt; 12; jj++) {
//	do jj=1,12
            scan2.nextLine();       // skip month header
            scan2.nextLine();       // skip blank line
//	read(2,*)
//	read(2,*)
            for (int i = 0; i &lt; 6; i++) {
                for (int j = 0; j &lt; 6; j++) {
                    a[i][j][jj] = scan2.nextDouble();
                }
            }
        }
//	read(2,'(6f7.3)') ((a(i,j,jj),j=1,6),i=1,6)
//	end do
        scan2.nextLine();   // B matrix header
//	read(2,*)                              !win98 did not like .mat
        for (int jj = 0; jj &lt; 12; jj++) {
//	do jj=1,12
            scan2.nextLine();       // skip month header
            scan2.nextLine();       // skip blank line
//	read(2,*)
//	read(2,*)
            for (int i = 0; i &lt; 6; i++) {
                for (int j = 0; j &lt; 6; j++) {
                    b[i][j][jj] = scan2.nextDouble();
                }
            }
//	read(2,*)
//	do jj=1,12
//	read(2,*)
//	read(2,*)
//	read(2,'(6f7.3)') ((b(i,j,jj),j=1,6),i=1,6)
        }
//	end do
        scan2.close();
//        close(2)
//
    }
//        END
//
//
//
//
//
////     ##################################################################
////     ##################################################################
////     ######                                                      ######
////     ######              SUBROUTINE TMAX_TMIN_123                ######
////     ######                                                      ######
////     ######                     (c) 1994                         ######
////     ######                                                      ######
////     ##################################################################
////     ##################################################################
//
    static int ns = 0;

    static void tmax_tmin_123(int j, /*int ns,*/ int iseed, boolean wetDayFlg, int mo) {
//        SUBROUTINE TMAX_TMIN_123(j,ns,iseed,nw,mo)
////
////#######################################################################
////
////     PURPOSE:
////
////    This program simulates max and min daily temp and radiation
////
////#######################################################################
////
////     AUTHOR: C.L.Hanson,K.A.Cumming,D.A.Woolhiser and C.W.Richardson
////
////     Modify history:
////     Yunyun Lu (November,1994)
////     transfer program from QBASIC to FORTRAN 77
////
////     Y. Lu (Dec. 1994)
////     added the full documentation
////
////#######################################################################
////
////     INPUT:
////
////     tmad       mean of tmax -- dry days
////     txs        standard deviation of tmax -- dry days
////     txm1       mean of tmax -- wet days
////     txs1       standard deviation of tmax -- wet days
////     tnm        mean of tmin -- wet and dry days
////     tns        standard deviation of tmin -- wet and dry days
////     rm0        mean of radiation -- dry days
////     rs0        standard deviation of radiation -- dry days
////     rm1        mean of radiation -- wet days
////     rs1        standard deviation of radiation -- wet days
////     rm         max. radiation along a given latitude
////    chipre      last step of chi(n), =chi(n-1)
////
////     OUTPUT:
////
////     p          Simulated M years of daily rainfall
////     tmax       Simulated M years of daily max. temperature
////     tmin       Simulated M years of daily min. temperature
////     rad        Simulated M years of daily total radiation
////
////#######################################################################
////
////     Variable Declarations.
////
////#######################################################################
////
        final int nn = 6;
//	parameter(idyr=365,nn=6)
//        common /const1/pi
//	common /tem11/tmad(idyr),txs(idyr),txm1(idyr),txs1(idyr),
//     !dmad(idyr),dxs(idyr),dxm1(idyr),dxs1(idyr),
//     !wmad(idyr),wxs(idyr),wxm1(idyr),wxs1(idyr),
//     !tnm(idyr),tns(idyr),tnm1(idyr),tns1(idyr),
//     !rm0(idyr),rs0(idyr),rm1(idyr),rs1(idyr),
//     !rm(idyr),chipre(nn)
//        common /temp2/p,tmax,tmin,dew,wind,rad
//        common /temp3/fran(12)
//        common /temp4/a(nn,nn,12),b(nn,nn,12)
        double[] rd = new double[nn];
        double[] rrd = new double[nn];
        double[] epsilon = new double[nn];
        double[] chi = new double[nn];
//	real    rd(nn),rrd(nn),epsilon(nn),chi(nn)
//	character*25 site,temp,matrix,gem6out,filename
//        common/filname/site,temp,matrix,gem6out,filename
////
////%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
////
////     Beginning of executable code...
////
////%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
//
        double rmm = rm0[j];
        double rs = rs0[j];
        double txxm = tmad[j];
        double txxs = txs[j];
        double tnnm = tnm[j];
        double tnns = tns[j];
        double dxxm = dmad[j];
        double dxxs = dxs[j];
        double wxxm = wmad[j];
        double wxxs = wxs[j];
//	rmm=rm0(j)
//	rs=rs0(j)
//	txxm=tmad(j)
//	txxs=txs(j)
//	tnnm=tnm(j)
//	tnns=tns(j)
//	dxxm=dmad(j)
//	dxxs=dxs(j)
//	wxxm=wmad(j)
//	wxxs=wxs(j)
        if (wetDayFlg) {
            rmm = rm1[j];
            rs = rs1[j];
            txxm = txm1[j];
            txxs = txs1[j];
            tnnm = tnm1[j];
            tnns = tns1[j];
            dxxm = dxm1[j];
            dxxs = dxs1[j];
            wxxm = wxm1[j];
            wxxs = wxs1[j];
//        IF (nw.eq.1) THEN
//	  rmm=rm1(j)
//	  rs=rs1(j)
//	  txxm=txm1(j)
//	  txxs=txs1(j)
//	  tnnm=tnm1(j)
//	  tnns=tns1(j)
//	  dxxm=dxm1(j)
//	  dxxs=dxs1(j)
//	  wxxm=wxm1(j)
//	  wxxs=wxs1(j)
        }
//        END IF
        if (ns &lt;= 1) {
            for (int k = 0; k &lt; 12; k += 2) {
                double u = ran3(iseed);
                double u1 = ran3(iseed);
                fran[k] = sqrt(-2 * Math.log(u)) * Math.cos(2 * pi * u1);
                fran[k + 1] = sqrt(-2 * Math.log(u)) * Math.sin(2 * pi * u1);
            }
            ns = nn;
        } else {
            ns = 0;
        }
//	IF(ns.le.1) THEN
//	  DO k=1,12,2
//	    u=ran3(iseed)
//	    u1=ran3(iseed)
//	    fran(k)=sqrt(-2*alog(u))*cos(2*pi*u1)
//	    fran(1+k)=sqrt(-2*alog(u))*sin(2*pi*u1)
//          END DO
//	  ns=nn
//        ELSE
//	  ns=0
//        END IF
////#######################################################################
////       for first day using last three number in the fran() array, THEN
////       for next day using the first three number in fran() array.
////       Iterately do it back and forewards.
////       integer ns control which three number will be used.
////#######################################################################
//
        for (int k = 0; k &lt; nn; k++) {
            epsilon[k] = fran[k + ns];
            rd[k] = 0.;
            rrd[k] = 0.;
        }
//	DO k=1,nn
//	  epsilon(k)=fran(k+ns)
//	  rd(k)=0.
//	  rrd(k)=0.
//        END DO
        for (int k = 0; k &lt; nn; k++) {
            for (int l = 0; l &lt; nn; l++) {
                rd[k] += b[k][l][mo] * epsilon[l];
                //! see orange book Eq.(13] page 5
                rrd[k] += a[k][l][mo] * chipre[l];
            }
        }
        for (int k = 0; k &lt; nn; k++) {
            chi[k] = rd[k] + rrd[k];
            chipre[k] = chi[k];            //! at the END of calculation
        }
//	DO k=1,nn
//	  DO l=1,nn
//	    rd(k)=rd(k)+b(k,l,mo)*epsilon(l)
//                                     ! see orange book Eq.(13) page 5
//	    rrd(k)=rrd(k)+a(k,l,mo)*chipre(l)
//          END DO
//        END DO
//	DO l=1,nn
//	  chi(l)=rd(l)+rrd(l)
//	  chipre(l)=chi(l)            ! at the END of calculation
//                                      ! store chi as next run's chipre.
//        END DO
        tmax = chi[1] * txxs + txxm;
        tmin = chi[2] * tnns + tnnm;
        dew = chi[4] * dxxs + dxxm;
        wind = chi[5] * wxxs + wxxm;
//	tmax=chi(2)*txxs+txxm
//	tmin=chi(3)*tnns+tnnm
//	dew=chi(5)*dxxs+dxxm
//	wind=chi(6)*wxxs+wxxm
////       print*,'wind',chi(5),wxxs,chi(5)*wxxs,wxxm,wind
////#######################################################################
////       Adjust for tmax, tmin, radiation
////#######################################################################
//
        if (tmin &gt; tmax) {
            double tmaxt = tmax;
            tmax = tmin;
            tmin = tmaxt;
        }
//	IF(tmin.gt.tmax) THEN   ! tmax should never be smaller than tmin
//	  tmaxt=tmax
//	  tmax=tmin
//	  tmin=tmaxt
//        END IF
        rad = chi[3] * rs + rmm;       //! rad not smaller than 0.1*max.radation
        double rmin = .1 * rm[j];           //! and not larger than max. radiation
        rad = (rad &gt; rmin) ? rad : rmin;
        rad = (rad &lt; rm[j]) ? rad : rm[j];
//	rad=chi(4)*rs+rmm       ! rad not smaller than 0.1*max.radation
//	rmin=.1*rm(j)           ! and not larger than max. radiation
//	rad=max(rad,rmin)
//	rad=min(rad,rm(j))
//	RETURN
//	END
    }
//
//
//
//
////     ##################################################################
////     ##################################################################
////     ######                                                      ######
////     ######              SUBROUTINE TOTAL_RAINFALL               ######
////     ######                                                      ######
////     ######                     (c) 1994                         ######
////     ######                                                      ######
////     ##################################################################
////     ##################################################################
//
//

    private static void total_rainfall() {
//       SUBROUTINE TOTAL_RAINFALL
////
////#######################################################################
////
////     PURPOSE:
////
////     Calculate the cdf for total rainfall in n days.
////     Using equation from ToDOrovic, P. and D.A. Woolhiser. 1975
////     A stochastic model of n-day precipitation. J.Applied Meteor.
////     14(1): 17-24. xlam=1/(mean daily rainfall)
////
////#######################################################################
////
////     AUTHOR: C.L.Hanson,K.A.Cumming,D.A.Woolhiser and C.W.Richardson
////
////     Modify history:
////     Yunyun Lu (November,1994)
////     transfer program from QBASIC to FORTRAN 77
////
////     Y. Lu (Dec. 1994)
////     added the full documentation
////
////#######################################################################
////
////
////     INPUT:
////
////     p00        Transition probability of a dry day following a dry day
////     p10        Transition probability of a dry day following a wet day
////     beta       The mean of the smaller exponential distributions
////     delta      The mean of the larger exponential distributions
////     albar      The weighting function
////
////#######################################################################
////
////     Variable Declarations.
////
////#######################################################################
////
//	parameter (idyr=365)
//	parameter (n1=4,nh=13,myr=12)
        final int n1 = 4;
        final int nh = 13;
        final int myr = 12;
//	integer   mday(myr)
//	integer xd,xm
        int xd;
        int xm;
//        common /const1/pi
//	common /sited2/delta(idyr),p00(idyr),p10(idyr),beta(idyr),albar
//        real xa(40),fx(40),psii(40),psi0(40)
//	real exp
//	character*20 char1
        String char1 = null;
//	data mday/307,338,1,32,62,93,123,154,185,215,246,276/
//                                           ! start from March 1 !!!
        int[] mday = {306, 337, 0, 31, 61, 122, 153, 185, 214, 245, 275};
//
//
////%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
////
////     Beginning of executable code...
////
////%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
////
//        print*,'Enter the month: (1,2,...,12)'
//        read(5,'(I4)') xm
//        print*,'Enter the day: (1,2,...,31)'
//        read(5,'(I4)') xd
//        print*,
//     ! 'Enter the length,n,of the period in days (1,2,...,31)'
//        read(5,'(I4)') n
//	j=mday(xm)+xd-1+n/2
//	IF (j.gt.idyr) j=j-idyr
////#######################################################################
////       calculate counting function density functions
////#######################################################################
//
//	q0=1-p00(j)
//	qi=1-p10(j)
//	psii(1)=(1-q0)**(n-1)*(1-qi)
//	psi0(1)=(1-q0)**n
//	nu=0
//	nt=n+1
//	DO 100 i=2,n
//	nu=nu+1
//	nci=aint(n+0.5-abs(2*nu-n+0.5)+0.01)
//	nc0=aint(n+0.5-abs(2*nu-n-0.5)+0.01)
//	nc=0
//	a=0
//	bc=1
//	nsw=1
//	sumi=0
//	sum0=0
//	termi=(1-qi)/(1-q0)
//	term0=q0/qi
//10      sumi=sumi+termi
//	sum0=sum0+term0
//	nc=nc+1
//	IF(nc.ne.nc0) GOTO 30
//	psi0(i)=qi**nu*(1-q0)**(n-nu)*sum0
//	IF(nc0.gt.nci) GOTO 100
//	IF(nc.eq.nci) GOTO 40
//20      IF (nsw.eq.2) GOTO 50
//	nsw=2
//	a=a+1
//	termi=termi*(nu-a+1)/a*q0/qi
//	term0=term0*(n-nu-a+1)/a*(1-qi)/(1-q0)
//	GOTO 10
//30      IF(nc.ne.nci) GOTO 20
//40      psii(i)=qi**nu*(1-q0)**(n-nu)*sumi
//	IF(nc0.gt.nci) GOTO 20
//        GOTO 100
//50      nsw=1
//	bc=bc+1
//	termi=termi*(n-nu-bc+1)/(bc-1)*(1-qi)/(1-q0)
//	term0=term0*(nu-bc+1)/(bc-1)*q0/qi
//	GOTO 10
//100     CONTINUE
//        psii(nt)=qi**n
//	psi0(nt)=qi**n*q0/qi
////#######################################################################
////       in this section we computer the CDF for n days
////       set exponential parameter to 1/mean
////#######################################################################
//	xlam=1/(albar*beta(j)+(1-albar)*delta(j))
//	nc=0
//	xa(1)=0.
//	k=j-n/2
//	IF(k.lt.1) k=idyr
//        ra=(1-p00(k))/(1+p10(k)-p00(k))
//	sbar=ra*n/xlam
//	dx=sbar/4
//        IF(dx.le.0.05) dx=0.05
//        print*,
//     :'Was it precipitating on the day before',
//     :xm,'-',xd,'? (Y,N,E(don"t know))'
//	read(5,'(a1)') char1
//	IF(char1.eq.'Y'.or.char1.eq.'y') ra=1.
//	IF(char1.eq.'N'.or.char1.eq.'n') ra=0.
//	IF(char1.eq.'E'.or.char1.eq.'e') THEN
//        write(6,'(/ 22(/5x,a)//)')
//     :'Enter the probability of precipitation
//     :on the day before in the form 0.xx'
//	  read(5,'(f10.5)')rb
//	  IF(rb.ne.0.0) ra=rb
//        END IF
//	fx(1)=(1-q0)-(ra*(qi-q0))
//	fx(1)=fx(1)*(1-q0)**(n-1)
//	sav=fx(1)
//	DO 110 i=2,40
//	  xa(i)=xa(i-1)+dx
//	  sum=0.0
//	  coef=xlam
//	  nc=nc+1
//	  DO 120 nu=1,n
//	    co=xa(i)
//	    suma=0.0
//	    terma=xlam*co**nu/nu
//	    DO 130 j=1,nu
//	      terma=terma*(nu-j+1)/(xlam*co)
//	      suma=suma-terma
//130         CONTINUE
//	    xc=terma/xlam+suma/(exp(xlam*co)*xlam)
//	    term=(ra*psii(nu+1)+(1-ra)*psi0(nu+1))*coef*xc
//	    sum=sum+term
//	    coef=coef*xlam/nu
//	    fx(i)=sav+sum
//120       CONTINUE
//	  IF (fx(i).lt.0.99) GOTO 110
//	  in=i
//	  GOTO 140
//110     CONTINUE
//        in=40
//140     print*,
//     !'Cumulative distribution function for',n,'days.'
//        print*,'Begining',xm,'-',xd
//	IF(ra.eq.0.0.or.ra.eq.1.0) THEN
//	  IF(ra.eq.1)  char1='WET'
//	  IF(ra.eq.0)  char1='DRY'
//	  print*,'Day before  ',n,' - day period was  ',char1
//        END IF
//	IF(ra.gt.0.0.or.ra.lt.1.0) THEN
//	  print*,'     x (inches)       f(x)'
//	  DO 150 i=1,in
//	    write(6,'(2f15.8)') xa(i),fx(i)
//150       CONTINUE
//        END IF
//	RETURN
    }
//	END
//
//
////
////     ##################################################################
////     ##################################################################
////     ######                                                      ######
////     ######                  FUNCTION RAN3                       ######
////     ######                                                      ######
////     ######                Copyright (c) 1995                    ######
////     ######                                                      ######
////     ##################################################################
////     ##################################################################
////
    static int iff = 0;
    static int inext;
    static int inextp;
    static int[] ma = new int[55];

    static double ran3(int iseed) {
//      FUNCTION RAN3(iseed)
////
////#######################################################################
////
////     PURPOSE:
////
////     Generates a random number between 0 and 1 by feeding
////     a negative integer iseed.
////
////     Reference: "Seminumerical Algorithms" by Donald Knuth
////
////#######################################################################
////
////     INPUT :
////
////       iseed      an arbitrary negative integer as a seed for a
////                  sequence of random numbers
////
////     OUTPUT:
////
////       ran3       A random number between 0 and 1.
////
////#######################################################################
////
//
////
////#######################################################################
////
////     Variable Declarations:
////
////#######################################################################
////
//      implicit none        ! Force explicit declarations
//
//      integer iseed        ! The seed for random number generation
//      real ran3            ! The function to generate random number.
////
////#######################################################################
////
////     Miscellaneous local variables:
////
////#######################################################################
////
//      integer mbig,mseed,mz,k,ii,inext,inextp,i,iff,mk,mj
//      real fac
//
        int mbig = 1000000000;
        int mseed = 161803398;
        int mz = 0;
        double fac = 1.0e-9;
        int ii;
        int mj;
//      parameter (mbig=1000000000,mseed=161803398,mz=0,fac=1.e-9)
//      integer ma(55)
//      save iff, inext, inextp, ma
//      data iff /0/
////
////%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
////
////     Beginning of executable code...
////
////%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
////
//
////
////######################################################################
////
////     Initialize the sequence of random numbers between 0 and 1,
////     using iseed.
////
////######################################################################
////
        if (iseed &lt; 0 || iff == 0) {
//      IF (iseed.lt.0.or.iff.eq.0) THEN
            iff = 1;
//        iff=1
            mj = mseed - Math.abs(iseed);
//        mj=mseed-iabs(iseed)
            mj = mj % mbig;
//        mj=mod(mj,mbig)
            ma[ma.length - 1] = mj;
//        ma(55)=mj
            int mk = 1;
//        mk=1
            for (int i = 0; i &lt; ma.length - 1; i++) {
//        DO 11 i=1,54
                ii = (21 * (i + 1)) % 55;
//          ii=mod(21*i,55)
                ma[ii - 1] = mk;
//          ma(ii)=mk
                mk = mj - mk;
//          mk=mj-mk
                if (mk &lt; mz) {
                    mk += mbig;
                }
//          IF (mk.lt.mz) mk=mk+mbig
                mj = ma[ii - 1];
//          mj=ma(ii)
            }
//11      CONTINUE
//
            for (int k = 0; k &lt; 4; k++) {
//        DO 13 k=1,4
                for (int i = 0; i &lt; ma.length; i++) {
//          DO 12 i=1,55
                    ma[i] -= ma[(i + 31) % 55];
//            ma(i)=ma(i)-ma(1+mod(i+30,55))
                    if (ma[i] &lt; mz) {
                        ma[i] += mbig;
                    }
//            IF (ma(i).lt.mz) ma(i)=ma(i)+mbig
                }
//12        CONTINUE
            }
//13      CONTINUE
//
            inext = 0;
            inextp = 31;
            iseed = 1;
//        inext=0
//        inextp=31
//        iseed=1
        }
//      ENDIF
////
////######################################################################
////
////     Start to generate a random number.
////
////######################################################################
////
        inext++;
//      inext=inext+1
        if (inext == 56) {
            inext = 1;
        }
//      IF (inext.eq.56) inext=1
        inextp++;
//      inextp=inextp+1
        if (inextp == 56) {
            inextp = 1;
        }
//      IF (inextp.eq.56) inextp=1
        mj = ma[inext - 1] - ma[inextp - 1];
//      mj=ma(inext)-ma(inextp)
        if (mj &lt; mz) {
            mj += mbig;
        }
//      IF (mj.lt.mz) mj=mj+mbig
        ma[inext - 1] = mj;
//      ma(inext)=mj
        return mj * fac;
//      ran3=mj*fac
//      RETURN
    }
//      END
//

    /**
     * sind - convenience method that takes degrees instead of radians
     */
    static double sind(double degrees) {
        return Math.sin(degrees / (180 / Math.PI));
    }

    /**
     * sind - convenience method that takes degrees instead of radians
     */
    static double tand(double degrees) {
        return Math.tan(degrees / (180 / Math.PI));
    }

    /**
     * sind - convenience method that takes degrees instead of radians
     */
    static double cosd(double degrees) {
        return Math.cos(degrees / (180 / Math.PI));
    }
    static int[] days = new int[365];
    static int[] mons = new int[365];
    static int[] rtnArr = new int[2];

    static void initDaysToDate() {
        int[] daysInMonths = {31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31};
        int mon = 0;
        int day = 0;
        for (int ddx = 0; ddx &lt; days.length; ddx++) {
            if (day == daysInMonths[mon]) {
                mon++;
                day = 0;
            }
            days[ddx] = day + 1;
            mons[ddx] = mon + 1;
            day++;
        }
        System.out.println("done init DayToDate");
    }

    static int[] dayInYearToDate(int day, boolean leapFlg, int offset) {
        rtnArr[0] = mons[day];
        rtnArr[1] = days[day];
        return rtnArr;
    }

    private static double sqrt(double val) {
        return Math.pow(val, 0.5);
    }
}
</pre></body></html>