c-------------------------------------------------------------------------------
c-------------------------------------------------------------------------------
c     TT6_NEW_VERTICAL_MSU1_ZM_VWEIGHT1.F            
c
c     Date:    10-FEB-2025.
c
c     Authors: Ben Santer, Klamath Falls.
c
c              Wolfgang Brueggemann, University of Hamburg, Germany.
c
c              Karl Taylor, LLNL, Livermore.
c
c     Version: # 1.0 / Single precision.
c
c
c
c     Computes detection times (as described in Santer et al., 2009, 2013) for 
c     raw and optimized fingerprints. Note that only the spatial properties of 
c     the fingerprint are optimized.
c
c     The user selects four different data sets:
c
c     SET1     Data set for computing noise EOFs used in optimization
c     SET2     Data set in which climate-change "fingerprint" is searched for
c     SET3     Data set for computing climate-change "fingerprint"
c     SET4     Data set for estimating natural variability statistics
c
c     In the terminology of Santer et al., 2004 ("Part II"):
c
c     SET1     is equivalent to c(x,t)  
c     SET2     is equivalent to o(x,t)
c     SET3     is equivalent to s(x,t); is used for estimating f(x) and f*m(x)
c     SET4     is equivalent to c1(x,t) 
c
c     The analysis is either in state space (nyh < nt - 1) or in sample space
c     (nyh > nt - 1), where nyh is the total number of spatial points (the
c     number of latitudes, ny, times the number of atmospheric layers, nl) and 
c     nt is the total number of time samples. Eigenvectors are always computed 
c     for the smaller scatter matrix.
c
c     Input data are zonal-mean annual-mean temperatures are either actual 
c     satellite measurements from the Stratospheric Sounding Unit (SSU) and
c     the Microwave Sounding Unit/Advanced Microwave Sounding Unit (MSU/AMSU)
c     or are "synthetic" SSU and MSU temperatures calculated from climate model 
c     simulation output. Results are for six different atmospheric layers:
c
c     KM          MB
c     45.00000    1.43 (SSU channel 3 weighting function)
c     38.00000    3.65 (SSU channel 2 weighting function)
c     30.00000   11.72 (SSU channel 1 weighting function)
c     19.02420   64.45 (MSU tls weighting function)
c      5.60216  498.01 (MSU ttt weighting function)
c      3.07098  694.78 (MSU tlt weighting function)
c
c     Note that:
c
c     tls  = Temperature Lower Stratosphere (MSU channel 4)
c     ttt  = Corrected MSU channel 2 (a.k.a. tf1)
c     tlt  = Lower Tropospheric Temperature retrieval 
c
c     The correcton to MSU channel 2 is for the influence of lower 
c     stratospheric cooling.
c
c     The input pressure coordinates for these layers are the estimated peaks 
c     of the SSU and MSU global-mean weighting functions. The user has the
c     option of weighting these six layers in different ways.
c
c     This program is capable of dealing with missing data in one or more of
c     the four input datasets.
c
c     The user supplies weights for each of the six atmospheric layers. By 
c     default, these are unit weights - i.e., each layer is weighted equally. 
c     If a layer has zero weight, the layer is not used in the fingerprint 
c     analysis.
c
c     27-JUL-22.
c     First working version of program.
c
c     30-JUL-22.
c     Implemented calculation of fractional detection time in Section 11.6.
c
c     04-AUG-22.
c     Processing of input data from CESM1 and CanESM2 LEs. Also allowed user
c     to select stipulated S/N threshold for estimating detection time.
c
c     05-AUG-22.
c     Processing of input data from CESM2 and MIROC6 LEs. 
c
c     06-AUG-22.
c     Corrected error setting path to MIROC_ES2L data.    
c
c     22-AUG-22.
c     Implemented analysis of tt6 input data - i.e., input of results for
c     SSU channels 1, 2, and 3 and MSU tls, ttt, and tlt. This version of
c     the program is no longer suitable for the fingerprint analysis of MSU
c     tt3 data!
c
c     26-AUG-22.
c     Corrected error in setting path to cesm2_waccm realizations 2 and 3
c     of the hst_ssp585 extended historical run.
c
c     25-SEP-22.
c     Implemented the imean=3 case, in which the overall global mean of the
c     four stratospheric layers (ssu3, ssu2, ssu1, and tls) is subtracted
c     from each latitude band of these four layers. For the troposphere, the
c     overall global mean of ttt and tlt is subtracted from each latitude
c     band of these two layers.
c
c     08-OCT-22.
c     Ensured that npc=216 (the combined ny * nh dimension).
c                
c     17-OCT-22.
c     Set the initial analysis period to 5 years (previously 10 years). This
c     means that the first analysis period for signal trends is 1986 to 1990.
c     The noise trends are then of length 5, 6, 7, ... 36 years.
c                
c     25-OCT-22.
c     Enabled a mass weighting option. Mass weights were calculated by Cheng-
c     Zhi Zou. The option is selected via ivweight. Weights for each of the six
c     layers are normalized by the sum of the weights.
c
c     02-NOV-22.
c     Revisited the mass weighting option. Mass weights were calculated by Cheng-
c     Zhi Zou on Oct. 31, 2022. These are different from the mass weights that 
c     were initially provided by Cheng-Zhi.            
c
c     04-NOV-22.
c     Enabled the ivweight=5 option. In this option, three layers with 
c     relatively small overlap in their weighting functions are considered: 
c     ssu3, tls, and tlt. This is the MIN_olap case.
c
c     07-NOV-22.
c     Allowed selection of update STAR MSU TLS and TTT data. These data sets
c     have been merged with STAR SSU ddata and with RSS TLT (since STAR does
c     not currently provide a TLT product).
c
c     10-NOV-22.
c     Modified program to allow selection of sm6 data sets. These data sets
c     use ttt data that have been calculated with latitudinally varying 
c     regression coefficients (i.e., tf2). In contrast, tt6 data sets use
c     ttt data calculated with latitudinally fixed regression coefficients
c     (i.e., tf1).
c
c     15-NOV-22.
c     Modified program to allow selection of "paired" data sets:
c      1. STAR SSU3/2/1; STAR TLS/TTT; RSS TLT.
c      2. STAR SSU3/2/1; STAR TLS/TTT; UAH TLT.
c
c     22-FEB-23.
c     Modified program to allow selection of satellite datasets. The
c     updating is through to the end of 2022.
c
c     18-APR-23.
c     Modified program to allow output of sampling distribution of control
c     run noise trends. These are the trends in the projection of the control
c     run annual-mean zonal-mean sm6 or tt6 data onto the fingerprint (i.e.,
c     the trends in the N(t) projection time series). Note that the sampling
c     distribution of trends in N(t) is only output for trend lengths of 37
c     years, which corresponds to the 1986 to 2022 analysis period for the SSU
c     data.
c
c     27-SEP-24.
c     Modified program to allow selection of updated satellite datasets. The
c     updating is through to the end of 2023. Also modified the output of the
c     sampling distribution of noise trends (see note above for 18-APR-23). 
c     The output is now for trends over 38 years.
c
c     28-SEP-24.
c     Made changes to Section 11.6 of program. The fractional detection time
c     was being computed incorrectly in cases where there was a sign change
c     in the S/N ratio between trend lengths with non-significant and with
c     significant values. Note that the frasctional detection time was not used
c     in the 2023 Santer et al. PNAS paper.
c
c     01-OCT-24.
c     Made changes to Section 11.6 of program. The fractional detection time
c     was being computed incorrectly in cases where there was a transition
c     from significant to non-significant S/N ratios at the trend length L
c     increased.
c
c     20-OCT-24.
c     Made changes to subroutines RD3D_YHT_NOMASK and RDMASK. The changes 
c     involve the addition of new domains available for user selection (e.g.,
c     40N-55N).
c
c     23-OCT-24.
c     Made changes to subroutines RD3D_YHT_NOMASK and RDMASK. The changes 
c     involve the addition of new domains available for user selection (e.g.,
c     0-90N and 90S-0).
c
c     24-OCT-24.
c     Fixed output of header lines for selected latitude boundaries.
c
c     21-JAN-25.
c     Changed dimension n2 from 100 to 200 years to accommodate input model
c     data for 1860 to 2050 (191 years).
c
c     22-JAN-25.
c     Updated observed "sixspace" temperature data.
c
c     10-FEB-25.
c     For plotting cartoons, changed the output of distributions of trends
c     in similarity between the fingerprint and control run data. The output
c     trends are now for timescales of 39 years rather than for 38 years.
c
c     REFERENCES:
c
c     R.W. Preisendorfer.
c     Principal Component Analysis in Meteorology and Oceanography.
c     Developments in Atmospheric Science 17, (ed C.D. Mobley). 
c     Elsevier, Amsterdam, 1988.
c
c     Santer, B.D., et al. 2009: Incorporating model quality information in
c     climate change detection and attribution studies. Proceedings of the 
c     National Academy of Sciences, 106, 14778-14783.
c
c     Santer, B.D., et al. 2013: Human and natural influences on the 
c     changing thermal structure of the atmosphere. Proceedings of the 
c     National Academy of Sciences, 110, 17235–17240.
c
c     EXTERNALS:
c
c     subroutine GETDATA - Reads input data
c     subroutine MMULT1  - Matrix multiplication
c     subroutine TRED2   - Householder reduction to tridiagonal matrix
c     subroutine TQLI    - Computes eigenvectors and eigenvalues
c     subroutine EIGSRT  - Sorts eigenvalues and eigenvectors
c     subroutine EIGVAR  - Explained variance of eigenvalues
c     subroutine EXVAR   - Calculates explained spatial variance
c
c
c
c-------------------------------------------------------------------------------
c-------------------------------------------------------------------------------
c
c
c
c
c
c
      program vert1
      implicit real (o-z, a-h)
c
c
c
c
c
c
c-------------------------------------------------------------------------------
c     ** parameters **
c
      parameter(ny     =    36)
      parameter(nh     =     6)
      parameter(n1     =  5000)             
      parameter(n2     =   200)             
      parameter(n3     =   200)             
      parameter(n4     =  5000)             
      parameter(npc    =   216)
      parameter(npcout =    10)
      parameter(nyh    = ny * nh)
c
c     ny        Retrieved latitude dimension
c     nh        Retrieved pressure dimension (no. of atmospheric layers)
c     nyh       Combined spatial dimension (latitude * no. of layers)
c     n1        Number of years for SET1
c     n2        Number of years for SET2
c     n3        Number of years for SET3
c     n4        Number of years for SET4
c     npc       Maximum number of PCs
c     npcout    Maximum number of PCs output
c
      include '/usr/local/uvcdat/Ben_drs/Externals/include/drsdef.h'
c
c
c
c-------------------------------------------------------------------------------
c     ** User options **
c
c     character*24 cexp1     1. Selected expt. (SET1)
c     character*24 cexp2     2. Selected expt. (SET2)
c     character*24 cexp3     3. Selected expt. (SET3)
c     character*24 cexp4     4. Selected expt. (SET4)
c     character*24 cmod1     5. Selected model (SET1)
c     character*24 cmod2     6. Selected model (SET2)
c     character*24 cmod3     7. Selected model (SET3)
c     character*24 cmod4     8. Selected model (SET4)
c     integer itrunc         9. No. of noise EOFs used for fingerprint rotation
c     integer iy1, iy2      10. Years retrieved, used in D&A analysis (SET1)    
c     integer iy3, iy4      11. Years retrieved, used in D&A analysis (SET2)
c     integer iy5, iy6      12. Years retrieved, used in D&A analysis (SET3)
c     integer iy7, iy8      13. Years retrieved, used in D&A analysis (SET4)
c     integer iref1, iref2  14. Reference period for computing anomalies (SET1)
c     integer iref3, iref4  15. Reference period for computing anomalies (SET2)
c     integer iref5, iref6  16. Reference period for computing anomalies (SET3)
c     integer iref7, iref8  17. Reference period for computing anomalies (SET4)
c     integer iover         18. 1=no chunk overlap / else user supplied overlap
c     integer iuserlap      19. User supplied chunk overlap
c     integer istart        20. Start date for computing linear trends in PCS
c     integer imean         21. 1 = retain spatial means / 2 = remove means   
c     integer nfinger       22. Number of EOFs used for defining fingerprint
c                               nfinger = 1             Fingerprint is EOF1 of SET3.      
c                               nfinger = > 1 and < 10  Fingerprint is weighted linear
c                               combination of first nfinger EOFs of SET3
c
c     integer idomain       23. Selected spatial domain
c     integer icasewt       24. Controls weighting of residual fingerprint
c                               0 = No residual fingerprint included
c                               1 = Inclusion of non-optimized residual  
c                               2 = Inclusion of residual with average optimization
c     integer ifiltin       25. Controls temporal filtering of anomaly data
c                               0 = no filter
c                               1 = 13-term Gaussian filter 
c                               2 = 13-term moving average filter 
c                               3 = 5-term binomial filter 
c     integer icritsn       26. Controls S/N threshold for detection time   
c                               0 = 1.641 (5%)              
c                               1 = 2.326 (1%)              
c                               2 = 2.000
c                               3 = 3.000
c                               5 = 5.000
c     integer ivweight      27. Controls vertical weighting
c                               0 = 1/1/1/1/1/1 (s3/s2/s1/m4/m2/mlt; ALL)              
c                               1 = 0/0/0/0/1/1 (s3/s2/s1/m4/m2/mlt; TROP only)              
c                               2 = 1/1/1/0/0/0 (s3/s2/s1/m4/m2/mlt; SSU only)              
c                               3 = 0/0/0/1/1/1 (s3/s2/s1/m4/m2/mlt; MSU only)              
c                               4 = Mass weights
c                               5 = 1/0/0/1/0/1 (s3/s2/s1/m4/m2/mlt; MIN_olap)              
c     real critwt           28. Threshold fraction of weights required
c     logical debug         29. .true. = diagnostics / .false. = no diagnostics
c     character*16 varnam1  30. Name of input variable
c
      integer itrunc, iover, iuserlap, istart, imean, nfinger, idomain
      integer iy1, iy2, iy3, iy4, iy5, iy6, iy7, iy8
      integer iref1, iref2, iref3, iref4, iref5, iref6, iref7, iref8
      integer icasewt, ifiltin, icritsn, ivweight
      real critwt
      logical debug
      character*24 cexp1, cexp2, cexp3, cexp4
      character*24 cmod1, cmod2, cmod3, cmod4
      character varnam1*16
c
c-------------------------------------------------------------------------------
c     ** Input data **
c
c     dimension bas1(ny, nh, n1)    Input data and anomalies, SET1         
c     dimension pro1(ny, nh, n2)    Input data and anomalies, SET2
c     dimension gss1(ny, nh, n3)    Input data and anomalies, SET3         
c     dimension dtt1(ny, nh, n4)    Input data and anomalies, SET4         
c     dimension awts1(ny, nh, n1)   Retrieved area weights, SET1
c     dimension awts2(ny, nh, n2)   Retrieved area weights, SET2
c     dimension awts3(ny, nh, n3)   Retrieved area weights, SET3
c     dimension awts4(ny, nh, n4)   Retrieved area weights, SET4
c     dimension fmask1(ny, nh)      Retrieved SET1 land area fraction (%)
c     dimension fmask2(ny, nh)      Retrieved SET2 land area fraction (%)
c     dimension fmask3(ny, nh)      Retrieved SET3 land area fraction (%)
c     dimension fmask4(ny, nh)      Retrieved SET4 land area fraction (%)
c     dimension cmask(ny, nh)       Common area weight mask of SET1/2/3/4
c     dimension gmask(ny, nh)       Geometrical mean land/sea mask of SET1/2/3/4
c     dimension rtime1(n1)          Time axis, SET1
c     dimension rtime2(n2)          Time axis, SET2
c     dimension rtime3(n3)          Time axis, SET3
c     dimension rtime4(n4)          Time axis, SET4
c     dimension rlat(ny)            Latitude coordinates
c     dimension rlev1(nh)           Pressure coordinates, SET1
c     dimension rlev2(nh)           Pressure coordinates, SET2
c     dimension rlev3(nh)           Pressure coordinates, SET3
c     dimension rlev4(nh)           Pressure coordinates, SET4
c     dimension rmass(nh)           Mass weights
c     dimension rmassn(nh)          Normalized mass weights
c
      dimension bas1(ny, nh, n1), pro1(ny, nh, n2)
      dimension gss1(ny, nh, n3), dtt1(ny, nh, n4)
      dimension awts1(ny, nh, n1), awts2(ny, nh, n2)
      dimension awts3(ny, nh, n3), awts4(ny, nh, n4)
      dimension fmask1(ny, nh), fmask2(ny, nh)
      dimension fmask3(ny, nh), fmask4(ny, nh)
      dimension cmask(ny, nh), gmask(ny, nh)
      dimension rtime1(n1), rtime2(n2), rtime3(n3), rtime4(n4)
      dimension rlat(ny) 
      dimension rlev1(nh), rlev2(nh), rlev3(nh), rlev4(nh)
      dimension rmass(nh), rmassn(nh)
c
c-------------------------------------------------------------------------------
c     ** EOF analysis **
c
c     integer ispace                1 = sample space / 2 = state space
c     integer irank1                Rank of SET1 scatter matrix
c     integer irank2                Rank of SET3 scatter matrix
c     dimension z1(n1, nyh)         Anomaly matrix (SET1) for EOF analysis
c     dimension p1(n2, nyh)         Anomaly matrix (SET2) for projection
c     dimension y1(n3, nyh)         Anomaly matrix (SET3) for EOF analysis
c     dimension d1(n4, nyh)         Anomaly matrix (SET4) for projection
c     dimension zt(nyh, n1)         Transposed anomaly matrix (SET1)
c     dimension yt(nyh, n3)         Transposed anomaly matrix (SET3)
c     dimension s1(n1, n1)          SET1 scatter matrix (state space)
c     dimension s2(n3, n3)          SET3 scatter matrix (state space)
c     dimension t1(n1, n1)          SET1 scatter matrix (sample space)
c     dimension t2(n3, n3)          SET3 scatter matrix (sample space)
c     dimension e1(n1, n1)          SET1 eigenvector matrix (state space)
c     dimension e2(n3, n3)          SET3 eigenvector matrix (state space)
c     dimension f1(n1, n1)          SET1 eigenvector matrix (sample space)
c     dimension f2(n3, n3)          SET3 eigenvector matrix (sample space)
c     dimension phat(nyh, n1)       SET1 EOFs (weighted)
c     dimension ptru(nyh, n1)       SET1 EOFs (unweighted)
c     dimension ptru1(ny, nh, n1)   SET1 EOFs (unweighted, reconstructed)
c     dimension ghat(nyh, n3)       SET3 EOFs (weighted)
c     dimension gtru(nyh, n3)       SET3 EOFs (unweighted)
c     dimension gtru1(ny, nh, n3)   SET3 EOFs (unweighted, reconstructed)
c     dimension val1(n1)            SET1 eigenvalues
c     dimension var1(n1)            SET1 explained variance
c     dimension cum1(n1)            SET1 cumulative variance
c     dimension val2(n3)            SET3 eigenvalues
c     dimension var2(n3)            SET3 explained variance
c     dimension cum2(n3)            SET3 cumulative variance
c     dimension alpha(n1)           Coefficients for projection of fingerprint
c     dimension greof(nyh)          Rotated fingerprint EOF1 (weighted)
c     dimension grtru(nyh)          Rotated fingerprint EOF1 (unweighted)
c     dimension grtru1(ny, nh)      Rotated fingerprint EOF1 (unweighted)
c     dimension fraw(nyh)           Raw fingerprint                   
c     dimension fpro(nyh)           Projection of raw fingerprint onto 
c                                     first itrunc SET1 EOFs
c     dimension fprotru(nyh)        Unweighted form of fpro 
c     dimension fprotru1(ny, nh)    Unweighted form of fpro 
c     dimension fres(nyh)           Residual fingerprint
c     dimension frestru(nyh)        Residual fingerprint (unweighted)
c     dimension frestru1(ny, nh)    Residual fingerprint (unweighted)
c     dimension frot(nyh)           Rotated projected fingerprint
c                                     (intermediate product)
c
      integer ispace, irank1, irank2
      dimension z1(n1, nyh), p1(n2, nyh), y1(n3, nyh), d1(n4, nyh)
      dimension zt(nyh, n1), yt(nyh, n3)
      dimension s1(n1, n1), s2(n3, n3), t1(n1, n1), t2(n3, n3)
      dimension e1(n1, n1), e2(n3, n3), f1(n1, n1), f2(n3, n3)
      dimension phat(nyh, n1), ghat(nyh, n3)
      dimension ptru(nyh, n1), gtru(nyh, n3)
      dimension ptru1(ny, nh, n1), gtru1(ny, nh, n3)
      dimension val1(n1), var1(n1), cum1(n1)
      dimension val2(n3), var2(n3), cum2(n3)
      dimension alpha(n1), greof(nyh), grtru(nyh), grtru1(ny, nh)
      dimension fraw(nyh), fpro(nyh), fres(nyh), frot(nyh)
      dimension frestru(nyh), frestru1(ny, nh)
      dimension fprotru(nyh), fprotru1(ny, nh)
c
c-------------------------------------------------------------------------------
c     ** PCS **
c
c     dimension bpc1(n1, npc)     Projection SET1 on EOFs of s(x,t)
c     dimension ppc1(n2, npc)     Projection SET2 on EOFs of s(x,t)
c     dimension gpc1(n3, npc)     Projection SET3 on EOFs of s(x,t)
c     dimension dpc1(n4, npc)     Projection SET4 on EOFs of s(x,t)
c     dimension bpc2(n1, npc)     Projection SET1 on control run EOFs, c(x,t)   
c     dimension ppc2(n2, npc)     Projection SET2 on control run EOFs, c(x,t)    
c     dimension gpc2(n3, npc)     Projection SET3 on control run EOFs, c(x,t)    
c     dimension dpc2(n4, npc)     Projection SET4 on control run EOFs, c(x,t)    
c     dimension bpc3(n1)          Projection SET1 on rotated fingerprint, f*m(x
c     dimension ppc3(n2)          Projection SET2 on rotated fingerprint, f*m(x)    
c     dimension gpc3(n3)          Projection SET3 on rotated fingerprint, f*m(x)
c     dimension dpc3(n4)          Projection SET4 on rotated fingerprint, f*m(x)
c     dimension bpc4(n1)          Projection SET1 on raw fingerprint, f(x)
c     dimension ppc4(n2)          Projection SET2 on raw fingerprint, f(x)
c     dimension gpc4(n3)          Projection SET3 on raw fingerprint, f(x)
c     dimension dpc4(n4)          Projection SET4 on raw fingerprint, f(x)
c     dimension var1con(npc)      Sum of squares for individual PCs, bpc2
c     dimension var2con(npc)      Sum of squares for individual PCs, ppc2
c     dimension var3con(npc)      Sum of squares for individual PCs, gpc2
c     dimension var4con(npc)      Sum of squares for individual PCs, dpc2
c     dimension var1consum(npc)   Cumulative sum of squares, bpc2
c     dimension var2consum(npc)   Cumulative sum of squares, ppc2
c     dimension var3consum(npc)   Cumulative sum of squares, gpc2
c     dimension var4consum(npc)   Cumulative sum of squares, dpc2
c
      dimension bpc1(n1, npc), ppc1(n2, npc)
      dimension gpc1(n3, npc), dpc1(n4, npc)
      dimension bpc2(n1, npc), ppc2(n2, npc)
      dimension gpc2(n3, npc), dpc2(n4, npc)
      dimension bpc3(n1), ppc3(n2), gpc3(n3), dpc3(n4)
      dimension bpc4(n1), ppc4(n2), gpc4(n3), dpc4(n4)
      dimension var1con(npc), var2con(npc), var3con(npc), var4con(npc)
      dimension var1consum(npc), var2consum(npc)
      dimension var3consum(npc), var4consum(npc)
      dimension diagnostic1(npc), diagnostic2(npc), diagnostic3(npc)
      dimension diagnostic4(npc), diagnostic5(npc), diagnostic6(npc)
c
c-------------------------------------------------------------------------------
c     ** S/N Statistics **
c
c     integer ntr(4)                     Number of trend lengths processed 
c     integer isize1(n1)                 Length of trends (SET1)
c     integer isize2(n2)                 Length of trends (SET2)
c     integer isize3(n3)                 Length of trends (SET3)
c     integer isize4(n4)                 Length of trends (SET4)
c     integer inumb1(n1)                 No. of trends (SET1)
c     integer inumb2(n4)                 No. of trends (SET4)
c     dimension corr1(npcout, npcout)    Pattern correlations (centered)
c     dimension corr2(npcout, npcout)    Pattern correlations (uncentered)
c     dimension ptrend1(n2)              Trends (proj. SET2 on raw fingerprint)
c     dimension gtrend1(n3)              Trends (proj. SET3 on raw fingerprint)
c     dimension ptrend2(n2)              Trends (proj. SET2 on f*m)
c     dimension gtrend2(n3)              Trends (proj. SET3 on f*m)
c     dimension bserr1(n1)               Std. error (proj. SET1 on f)
c     dimension dserr1(n4)               Std. error (proj. SET4 on f)
c     dimension bserr2(n1)               Std. error (proj. SET1 on f*m)
c     dimension dserr2(n4)               Std. error (proj. SET4 on f*m)
c     dimension psn1(n2)                 S/N: S=SET2 PC1, N=SET1
c     dimension psn2(n2)                 S/N: S=SET2 PC1, N=SET1 (f*m)
c     dimension psn3(n2)                 S/N: S=SET2 PC1, N=SET4
c     dimension psn4(n2)                 S/N: S=SET2 PC1, N=SET4 (f*m)
c     dimension gsn1(n3)                 S/N: S=SET3 PC1, N=SET1
c     dimension gsn2(n3)                 S/N: S=SET3 PC1, N=SET1 (f*m)
c     dimension gsn3(n3)                 S/N: S=SET3 PC1, N=SET4
c     dimension gsn4(n3)                 S/N: S=SET3 PC1, N=SET4 (f*m)
c     dimension gmean(nh, n1)            Spatial means
c
      integer ntr(4)
      integer isize1(n1), isize2(n2), isize3(n3), isize4(n4)
      integer inumb1(n1), inumb2(n4)
      dimension corr1(npcout, npcout), corr2(npcout, npcout)
      dimension ptrend1(n2), ptrend2(n2), gtrend1(n3), gtrend2(n3)
      dimension bserr1(n1), dserr1(n4)
      dimension bserr2(n1), dserr2(n4)
      dimension psn1(n2), psn2(n2), psn3(n2), psn4(n2)
      dimension gsn1(n3), gsn2(n3), gsn3(n3), gsn4(n3), gmean(nh, n1)
c
c-------------------------------------------------------------------------------
c     ** Workspace **
c
c     dimension xx(n1)          Workspace
c     dimension zz(n1)          Workspace
c     dimension yy(n1)          Workspace
c     dimension tt(n1)          Workspace
c     dimension work2a(n1)      Workspace (TRED2, TQLI)
c     dimension work2b(n3)      Workspace (TRED2, TQLI)
c     dimension adist(n1)       Workspace (ERR1)
c     dimension bdist(n1)       Workspace (ERR1)
c     dimension cdist(n1)       Workspace (ERR1)
c     dimension stderr(n1)      Workspace (ERR1)
c     dimension eofout1(nyh)    Workspace
c     dimension eofout2(nyh)    Workspace
c     dimension wmask(nyh)      1-D area weights
c     dimension store1(ny, nh)  Workspace       
c
      dimension xx(n1), zz(n1), yy(n1), tt(n1)
      dimension work2a(n1), work2b(n3)
      dimension adist(n1), bdist(n1), cdist(n1), stderr(n1)
      dimension eofout1(nyh), eofout2(nyh), wmask(nyh), store1(ny, nh)
c
c-------------------------------------------------------------------------------
c     ** DRS variables **
c
c     integer aslun             Assigns logical unit for .dat and .nc files
c     integer cluvdb            Clears buffer
c     integer setname           Sets experiment name, variable type, units
c     integer setvdim           Sets irregular dimensions
c     integer putvdim           Writes irregular dimensions
c     integer putdat            Writes data
c     integer cllun             Closes logical unit for .dat and .nc files
c
      integer aslun, cluvdb, setname, setvdim
      integer setdim, putvdim, putdat, cllun
c
c-------------------------------------------------------------------------------
c     ** Character variables **
c
c     character*120 infil1              Input file (SET1 data)
c     character*120 infil2              Input file (SET2 data)
c     character*120 infil3              Input file (SET3 data)
c     character*120 infil4              Input file (SET4 data)
c     character*120 pmask1              Input file (SET1 mask)
c     character*120 pmask2              Input file (SET2 mask)
c     character*120 pmask3              Input file (SET3 mask)
c     character*120 pmask4              Input file (SET4 mask)
c     character  *8 yweights1           Area weights (SET1)
c     character  *8 yweights2           Area weights (SET2)
c     character  *8 yweights3           Area weights (SET3)
c     character  *8 yweights4           Area weights (SET4)
c     character*120 drsdic1             Output DRS data file (EOFS)
c     character*120 drsdat1             Output DRS dictionary file (EOFS)
c     character*120 outfil1             Output file (PCs)
c     character *16 varnam2             Variable name for EOFs (on output)
c     character *40 units1              Units of variable
c     character *80 title1              Title   
c     character*120 source1(npcout)     Text              
c     character*120 source2             Text              
c     character  *9 str0                Text
c     character *60 str1                Text
c     character *40 str2                Text
c     character *24 str3                Text
c     character*133 str4                Text
c     character *10 str5                Text
c     character*142 str6                Text
c     character *20 string2             Text
c     character *20 string3             Text
c     character  *1 cc1                 String
c     character  *2 cc2                 String
c     character  *1 cc3                 String
c     character *10 cc4                 String
c     character  *6 sta                 String
c     character  *5 stb                 String
c     character *24 stc                 String
c     character *32 std                 String
c     character *70 st1, st2, st3, st4  String
c     character *70 st5, st6, st7, st8  String
c     character  *4 cs1, cs2, cs3, cs4  String
c
      character infil1*120, infil2*120, infil3*120, infil4*120
      character pmask1*120, pmask2*120, pmask3*120, pmask4*120
      character yweights1*8, yweights2*8, yweights3*8, yweights4*8
      character drsdic1*120, drsdat1*120, outfil1*120
      character varnam2*16, units1*40, title1*80 
      character source1(npcout)*120, source2*120
      character str0*9, str1*60, str2*40, str3*24, str4*133, str5*10
      character str6*142, string2*20, string3*20
      character cc1*1, cc2*2, cc3*1, cc4*10
      character sta*6, stb*5, stc*24, std*32
      character st1*70, st2*70, st3*70, st4*70
      character st5*70, st6*70, st7*70, st8*70
      character cs1*4, cs2*4, cs3*4, cs4*4
c
c-------------------------------------------------------------------------------
c     ** data statements **
c
      data debug / .true. / ! .true. = debugging on.
c
      data bogus / 1.e20 /  ! Missing data code.
c
      data iover / 1 /      ! 1=No chunk overlap.
cBDS  data iover / 0 /      ! 1=User supplied chunk overlap.
c
      data iuserlap / 1 /   ! User-supplied chunk overlap (not used if iover=1)
c
      data ifact / 3 /      ! Annual-mean input data.
c
      data critwt / 0.25 /
c
c
c
c
c
c
c-------------------------------------------------------------------------------
c     ** 1.0 **
c
c     PREPROCESSING: BATCH INPUTS, CHECKS, DATA INPUT.
c
c     Read batch inputs on unit 1.
c-------------------------------------------------------------------------------
c
      open(unit=1, file='GUESS1', form='formatted', status='old')
c
      read(1, 2) cmod1
      read(1, 2) cexp1
      read(1, 2) cmod2
      read(1, 2) cexp2
      read(1, 2) cmod3
      read(1, 2) cexp3
      read(1, 2) cmod4
      read(1, 2) cexp4
      read(1, *) iy1, iy2
      read(1, *) iy3, iy4
      read(1, *) iy5, iy6
      read(1, *) iy7, iy8
      read(1, 4) varnam1
      read(1, *) iref1, iref2
      read(1, *) iref3, iref4
      read(1, *) iref5, iref6
      read(1, *) iref7, iref8
      read(1, *) itrunc
      read(1, *) istart
      read(1, *) imean 
      read(1, *) nfinger
      read(1, *) idomain
      read(1, *) icasewt
      read(1, *) ifiltin
      read(1, *) icritsn
      read(1, *) ivweight
    2 format(a24)
    3 format(a8)
    4 format(a16)
      close(unit=1)                      
c
c-------------------------------------------------------------------------------
c     ** 1.1 **
c
c     Generate full path to input data.
c-------------------------------------------------------------------------------
c
      call create_path1(cmod1, cexp1, varnam1, infil1, pmask1)
      call create_path1(cmod2, cexp2, varnam1, infil2, pmask2)
      call create_path1(cmod3, cexp3, varnam1, infil3, pmask3)
      call create_path1(cmod4, cexp4, varnam1, infil4, pmask4)
c
c-------------------------------------------------------------------------------
c     ** 1.2 **
c
c     Echo some inputs.
c-------------------------------------------------------------------------------
c
      print *, ' ******************************************************'
      print *, ' ************ OPTIMAL_MSU_MODEL_QUALITY1.F ************'
      print *, ' ******************************************************'
      print *, ' '
      print *, ' '
      print *, ' EOF AND PC ANALYSIS'
      print *, ' '
      print *, ' '
      print *, ' YOU HAVE SELECTED THE FOLLOWING OPTIONS: '
      print *, ' '
      print *, ' '
      print *, ' 1. STANDARD OUTPUT VARIABLE:  ', varnam1
      print *, ' 2. NOISE EOFS:                ', infil1
      print *, ' 3. PROJECTIONS:               ', infil2
      print *, ' 4. GUESS EOFS:                ', infil3
      print *, ' 5. STATISTICS:                ', infil4
      print *, ' '
      print *, ' ******************************************************'
      print *, ' ******************************************************'
c
c-------------------------------------------------------------------------------
c     ** 1.3 **
c
c     Calculate lengths of various character strings. 
c-------------------------------------------------------------------------------
c
      ia1 = index(varnam1, ' ') - 1
c
      ib1 = index(cmod1, ' ') - 1
      ib2 = index(cmod2, ' ') - 1
      ib3 = index(cmod3, ' ') - 1
      ib4 = index(cmod4, ' ') - 1
c
      ic1 = index(cexp1, ' ') - 1
      ic2 = index(cexp2, ' ') - 1
      ic3 = index(cexp3, ' ') - 1
      ic4 = index(cexp4, ' ') - 1
c
c-------------------------------------------------------------------------------
c     ** 1.4 **
c
c     Echo paths to masks containing information on percentage land area
c     fraction (sftlf) for SET1, SET2, SET3, and SET4 data.
c-------------------------------------------------------------------------------
c
      print *, ' ******************************************************'
      print *, ' YOU HAVE SELECTED THE FOLLOWING LAND/SEA MASKS: '
      print *, ' '
      print *, ' pmask1 = ', pmask1
      print *, ' pmask2 = ', pmask2
      print *, ' pmask3 = ', pmask3
      print *, ' pmask4 = ', pmask4
      print *, ' ******************************************************'
c
c-------------------------------------------------------------------------------
c     ** 1.5 **
c
c     Set various character strings needed for CMIP5 MMA files. 
c-------------------------------------------------------------------------------
c
      if (cmod1(1:ib1) .eq. 'MMA_CMIP5_28models') cs1 = 'c5s-'
      if (cmod2(1:ib2) .eq. 'MMA_CMIP5_28models') cs2 = 'c5s-'
      if (cmod3(1:ib3) .eq. 'MMA_CMIP5_28models') cs3 = 'c5s-'
      if (cmod4(1:ib4) .eq. 'MMA_CMIP5_28models') cs4 = 'c5s-'
c
      if (cmod1(1:ib1) .eq. 'MMA_CMIP5_36models') cs1 = 'c5n-'
      if (cmod2(1:ib2) .eq. 'MMA_CMIP5_36models') cs2 = 'c5n-'
      if (cmod3(1:ib3) .eq. 'MMA_CMIP5_36models') cs3 = 'c5n-'
      if (cmod4(1:ib4) .eq. 'MMA_CMIP5_36models') cs4 = 'c5n-'
c
c-------------------------------------------------------------------------------
c     ** 1.6 **
c
c     Set various character strings needed for CMIP6 MMA files. 
c-------------------------------------------------------------------------------
c
      if (cmod1(1:ib1) .eq. 'MMA_CMIP6_21models') cs1 = 'c6s-'
      if (cmod2(1:ib2) .eq. 'MMA_CMIP6_21models') cs2 = 'c6s-'
      if (cmod3(1:ib3) .eq. 'MMA_CMIP6_21models') cs3 = 'c6s-'
      if (cmod4(1:ib4) .eq. 'MMA_CMIP6_21models') cs4 = 'c6s-'
c
      if (cmod1(1:ib1) .eq. 'MMA_CMIP6_30models') cs1 = 'c6n-'
      if (cmod2(1:ib2) .eq. 'MMA_CMIP6_30models') cs2 = 'c6n-'
      if (cmod3(1:ib3) .eq. 'MMA_CMIP6_30models') cs3 = 'c6n-'
      if (cmod4(1:ib4) .eq. 'MMA_CMIP6_30models') cs4 = 'c6n-'
c
      if (cmod1(1:ib1) .eq. 'MMA_CMIP6_09models') then
        if (cexp1(1:ic1) .eq. 'piControl')  cs1 = 'C6n-'
        if (cexp1(1:ic1) .eq. 'hst_ssp585') cs1 = 'C6s-'
        if (cexp2(1:ic2) .eq. 'piControl')  cs2 = 'C6n-'
        if (cexp2(1:ic2) .eq. 'hst_ssp585') cs2 = 'C6s-'
        if (cexp3(1:ic3) .eq. 'piControl')  cs3 = 'C6n-'
        if (cexp3(1:ic3) .eq. 'hst_ssp585') cs3 = 'C6s-'
        if (cexp4(1:ic4) .eq. 'piControl')  cs4 = 'C6n-'
        if (cexp4(1:ic4) .eq. 'hst_ssp585') cs4 = 'C6s-'
      end if
c
c-------------------------------------------------------------------------------
c     ** 1.6 **
c
c     Set various character strings needed for observed files. 
c-------------------------------------------------------------------------------
c
      if (cmod1(1:ib1) .eq. 'SSSRRR_jan25_05x05') cs1 = 'ssr-'
      if (cmod2(1:ib2) .eq. 'SSSRRR_jan25_05x05') cs2 = 'ssr-'
      if (cmod3(1:ib3) .eq. 'SSSRRR_jan25_05x05') cs3 = 'ssr-'
      if (cmod4(1:ib4) .eq. 'SSSRRR_jan25_05x05') cs4 = 'ssr-'
c
      if (cmod1(1:ib1) .eq. 'SSSUUU_jan25_05x05') cs1 = 'ssu-'
      if (cmod2(1:ib2) .eq. 'SSSUUU_jan25_05x05') cs2 = 'ssu-'
      if (cmod3(1:ib3) .eq. 'SSSUUU_jan25_05x05') cs3 = 'ssu-'
      if (cmod4(1:ib4) .eq. 'SSSUUU_jan25_05x05') cs4 = 'ssu-'
c
      if (cmod1(1:ib1) .eq. 'SSSSSS_jan25_05x05') cs1 = 'sss-'
      if (cmod2(1:ib2) .eq. 'SSSSSS_jan25_05x05') cs2 = 'sss-'
      if (cmod3(1:ib3) .eq. 'SSSSSS_jan25_05x05') cs3 = 'sss-'
      if (cmod4(1:ib4) .eq. 'SSSSSS_jan25_05x05') cs4 = 'sss-'
c
c-------------------------------------------------------------------------------
c     ** 1.7 **
c
c     Set various character strings needed for individual CMIP5 models.
c-------------------------------------------------------------------------------
c
      if (cmod1(1:ib1) .eq. 'access1_0_05x05') cs1 = 'A10-'
      if (cmod2(1:ib2) .eq. 'access1_0_05x05') cs2 = 'A10-'
      if (cmod3(1:ib3) .eq. 'access1_0_05x05') cs3 = 'A10-'
      if (cmod4(1:ib4) .eq. 'access1_0_05x05') cs4 = 'A10-'
c
      if (cmod1(1:ib1) .eq. 'access1_3_05x05') cs1 = 'A13-'
      if (cmod2(1:ib2) .eq. 'access1_3_05x05') cs2 = 'A13-'
      if (cmod3(1:ib3) .eq. 'access1_3_05x05') cs3 = 'A13-'
      if (cmod4(1:ib4) .eq. 'access1_3_05x05') cs4 = 'A13-'
c
      if (cmod1(1:ib1) .eq. 'bcc_csm1_1_05x05') cs1 = 'BCC-'
      if (cmod2(1:ib2) .eq. 'bcc_csm1_1_05x05') cs2 = 'BCC-'
      if (cmod3(1:ib3) .eq. 'bcc_csm1_1_05x05') cs3 = 'BCC-'
      if (cmod4(1:ib4) .eq. 'bcc_csm1_1_05x05') cs4 = 'BCC-'
c
      if (cmod1(1:ib1) .eq. 'bcc_csm1_1m_05x05') cs1 = 'BCm-'
      if (cmod2(1:ib2) .eq. 'bcc_csm1_1m_05x05') cs2 = 'BCm-'
      if (cmod3(1:ib3) .eq. 'bcc_csm1_1m_05x05') cs3 = 'BCm-'
      if (cmod4(1:ib4) .eq. 'bcc_csm1_1m_05x05') cs4 = 'BCm-'
c
      if (cmod1(1:ib1) .eq. 'can_esm2_05x05') cs1 = 'CAN-'
      if (cmod2(1:ib2) .eq. 'can_esm2_05x05') cs2 = 'CAN-'
      if (cmod3(1:ib3) .eq. 'can_esm2_05x05') cs3 = 'CAN-'
      if (cmod4(1:ib4) .eq. 'can_esm2_05x05') cs4 = 'CAN-'
c
      if (cmod1(1:ib1) .eq. 'ccsm4_05x05') cs1 = 'CCM-'
      if (cmod2(1:ib2) .eq. 'ccsm4_05x05') cs2 = 'CCM-'
      if (cmod3(1:ib3) .eq. 'ccsm4_05x05') cs3 = 'CCM-'
      if (cmod4(1:ib4) .eq. 'ccsm4_05x05') cs4 = 'CCM-'
c
      if (cmod1(1:ib1) .eq. 'cesm1_bgc_05x05') cs1 = 'CEB-'
      if (cmod2(1:ib2) .eq. 'cesm1_bgc_05x05') cs2 = 'CEB-'
      if (cmod3(1:ib3) .eq. 'cesm1_bgc_05x05') cs3 = 'CEB-'
      if (cmod4(1:ib4) .eq. 'cesm1_bgc_05x05') cs4 = 'CEB-'
c
      if (cmod1(1:ib1) .eq. 'cesm1_cam5_05x05') cs1 = 'CM5-'
      if (cmod2(1:ib2) .eq. 'cesm1_cam5_05x05') cs2 = 'CM5-'
      if (cmod3(1:ib3) .eq. 'cesm1_cam5_05x05') cs3 = 'CM5-'
      if (cmod4(1:ib4) .eq. 'cesm1_cam5_05x05') cs4 = 'CM5-'
c
      if (cmod1(1:ib1) .eq. 'csiro_mk3_6_0_05x05') cs1 = 'CSI-'
      if (cmod2(1:ib2) .eq. 'csiro_mk3_6_0_05x05') cs2 = 'CSI-'
      if (cmod3(1:ib3) .eq. 'csiro_mk3_6_0_05x05') cs3 = 'CSI-'
      if (cmod4(1:ib4) .eq. 'csiro_mk3_6_0_05x05') cs4 = 'CSI-'
c
      if (cmod1(1:ib1) .eq. 'ec_earth_05x05') cs1 = 'ECE-'
      if (cmod2(1:ib2) .eq. 'ec_earth_05x05') cs2 = 'ECE-'
      if (cmod3(1:ib3) .eq. 'ec_earth_05x05') cs3 = 'ECE-'
      if (cmod4(1:ib4) .eq. 'ec_earth_05x05') cs4 = 'ECE-'
c
      if (cmod1(1:ib1) .eq. 'gfdl_cm3_05x05') cs1 = 'GF3-'
      if (cmod2(1:ib2) .eq. 'gfdl_cm3_05x05') cs2 = 'GF3-'
      if (cmod3(1:ib3) .eq. 'gfdl_cm3_05x05') cs3 = 'GF3-'
      if (cmod4(1:ib4) .eq. 'gfdl_cm3_05x05') cs4 = 'GF3-'
c
      if (cmod1(1:ib1) .eq. 'gfdl_esm2g_05x05') cs1 = 'GFG-'
      if (cmod2(1:ib2) .eq. 'gfdl_esm2g_05x05') cs2 = 'GFG-'
      if (cmod3(1:ib3) .eq. 'gfdl_esm2g_05x05') cs3 = 'GFG-'
      if (cmod4(1:ib4) .eq. 'gfdl_esm2g_05x05') cs4 = 'GFG-'
c
      if (cmod1(1:ib1) .eq. 'gfdl_esm2m_05x05') cs1 = 'GFM-'
      if (cmod2(1:ib2) .eq. 'gfdl_esm2m_05x05') cs2 = 'GFM-'
      if (cmod3(1:ib3) .eq. 'gfdl_esm2m_05x05') cs3 = 'GFM-'
      if (cmod4(1:ib4) .eq. 'gfdl_esm2m_05x05') cs4 = 'GFM-'
c
      if (cmod1(1:ib1) .eq. 'giss_e2_h_p1_05x05') cs1 = 'GH1-'
      if (cmod2(1:ib2) .eq. 'giss_e2_h_p1_05x05') cs2 = 'GH1-'
      if (cmod3(1:ib3) .eq. 'giss_e2_h_p1_05x05') cs3 = 'GH1-'
      if (cmod4(1:ib4) .eq. 'giss_e2_h_p1_05x05') cs4 = 'GH1-'
c
      if (cmod1(1:ib1) .eq. 'giss_e2_h_p3_05x05') cs1 = 'GH3-'
      if (cmod2(1:ib2) .eq. 'giss_e2_h_p3_05x05') cs2 = 'GH3-'
      if (cmod3(1:ib3) .eq. 'giss_e2_h_p3_05x05') cs3 = 'GH3-'
      if (cmod4(1:ib4) .eq. 'giss_e2_h_p3_05x05') cs4 = 'GH3-'
c
      if (cmod1(1:ib1) .eq. 'giss_e2_r_p1_05x05') cs1 = 'GR1-'
      if (cmod2(1:ib2) .eq. 'giss_e2_r_p1_05x05') cs2 = 'GR1-'
      if (cmod3(1:ib3) .eq. 'giss_e2_r_p1_05x05') cs3 = 'GR1-'
      if (cmod4(1:ib4) .eq. 'giss_e2_r_p1_05x05') cs4 = 'GR1-'
c
      if (cmod1(1:ib1) .eq. 'giss_e2_r_p2_05x05') cs1 = 'GR2-'
      if (cmod2(1:ib2) .eq. 'giss_e2_r_p2_05x05') cs2 = 'GR2-'
      if (cmod3(1:ib3) .eq. 'giss_e2_r_p2_05x05') cs3 = 'GR2-'
      if (cmod4(1:ib4) .eq. 'giss_e2_r_p2_05x05') cs4 = 'GR2-'
c
      if (cmod1(1:ib1) .eq. 'giss_e2_r_p3_05x05') cs1 = 'GR3-'
      if (cmod2(1:ib2) .eq. 'giss_e2_r_p3_05x05') cs2 = 'GR3-'
      if (cmod3(1:ib3) .eq. 'giss_e2_r_p3_05x05') cs3 = 'GR3-'
      if (cmod4(1:ib4) .eq. 'giss_e2_r_p3_05x05') cs4 = 'GR3-'
c
      if (cmod1(1:ib1) .eq. 'hadgem2_cc_05x05') cs1 = 'HDC-'
      if (cmod2(1:ib2) .eq. 'hadgem2_cc_05x05') cs2 = 'HDC-'
      if (cmod3(1:ib3) .eq. 'hadgem2_cc_05x05') cs3 = 'HDC-'
      if (cmod4(1:ib4) .eq. 'hadgem2_cc_05x05') cs4 = 'HDC-'
c
      if (cmod1(1:ib1) .eq. 'hadgem2_es_05x05') cs1 = 'HDE-'
      if (cmod2(1:ib2) .eq. 'hadgem2_es_05x05') cs2 = 'HDE-'
      if (cmod3(1:ib3) .eq. 'hadgem2_es_05x05') cs3 = 'HDE-'
      if (cmod4(1:ib4) .eq. 'hadgem2_es_05x05') cs4 = 'HDE-'
c
      if (cmod1(1:ib1) .eq. 'miroc5_05x05') cs1 = 'MI5-'
      if (cmod2(1:ib2) .eq. 'miroc5_05x05') cs2 = 'MI5-'
      if (cmod3(1:ib3) .eq. 'miroc5_05x05') cs3 = 'MI5-'
      if (cmod4(1:ib4) .eq. 'miroc5_05x05') cs4 = 'MI5-'
c
      if (cmod1(1:ib1) .eq. 'miroc_esm_05x05') cs1 = 'MIE-'
      if (cmod2(1:ib2) .eq. 'miroc_esm_05x05') cs2 = 'MIE-'
      if (cmod3(1:ib3) .eq. 'miroc_esm_05x05') cs3 = 'MIE-'
      if (cmod4(1:ib4) .eq. 'miroc_esm_05x05') cs4 = 'MIE-'
c
      if (cmod1(1:ib1) .eq. 'miroc-esm-chem_05x05') cs1 = 'MIC-'
      if (cmod2(1:ib2) .eq. 'miroc-esm-chem_05x05') cs2 = 'MIC-'
      if (cmod3(1:ib3) .eq. 'miroc-esm-chem_05x05') cs3 = 'MIC-'
      if (cmod4(1:ib4) .eq. 'miroc-esm-chem_05x05') cs4 = 'MIC-'
c
      if (cmod1(1:ib1) .eq. 'mpi_esm_lr_05x05') cs1 = 'MPL-'
      if (cmod2(1:ib2) .eq. 'mpi_esm_lr_05x05') cs2 = 'MPL-'
      if (cmod3(1:ib3) .eq. 'mpi_esm_lr_05x05') cs3 = 'MPL-'
      if (cmod4(1:ib4) .eq. 'mpi_esm_lr_05x05') cs4 = 'MPL-'
c
      if (cmod1(1:ib1) .eq. 'mpi_esm_mr_05x05') cs1 = 'MPM-'
      if (cmod2(1:ib2) .eq. 'mpi_esm_mr_05x05') cs2 = 'MPM-'
      if (cmod3(1:ib3) .eq. 'mpi_esm_mr_05x05') cs3 = 'MPM-'
      if (cmod4(1:ib4) .eq. 'mpi_esm_mr_05x05') cs4 = 'MPM-'
c
      if (cmod1(1:ib1) .eq. 'mri_cgcm3_05x05') cs1 = 'MRI-'
      if (cmod2(1:ib2) .eq. 'mri_cgcm3_05x05') cs2 = 'MRI-'
      if (cmod3(1:ib3) .eq. 'mri_cgcm3_05x05') cs3 = 'MRI-'
      if (cmod4(1:ib4) .eq. 'mri_cgcm3_05x05') cs4 = 'MRI-'
c
      if (cmod1(1:ib1) .eq. 'noresm1-m_05x05') cs1 = 'NOM-'
      if (cmod2(1:ib2) .eq. 'noresm1-m_05x05') cs2 = 'NOM-'
      if (cmod3(1:ib3) .eq. 'noresm1-m_05x05') cs3 = 'NOM-'
      if (cmod4(1:ib4) .eq. 'noresm1-m_05x05') cs4 = 'NOM-'
c
      if (cmod1(1:ib1) .eq. 'noresm1_me_05x05') cs1 = 'NOE-'
      if (cmod2(1:ib2) .eq. 'noresm1_me_05x05') cs2 = 'NOE-'
      if (cmod3(1:ib3) .eq. 'noresm1_me_05x05') cs3 = 'NOE-'
      if (cmod4(1:ib4) .eq. 'noresm1_me_05x05') cs4 = 'NOE-'
c
c-------------------------------------------------------------------------------
c     ** 1.8 **
c
c     Set various character strings needed for individual CMIP6 models.
c-------------------------------------------------------------------------------
c
      if (cmod1(1:ib1) .eq. 'access_cm2_05x05') cs1 = 'AC2-'
      if (cmod2(1:ib2) .eq. 'access_cm2_05x05') cs2 = 'AC2-'
      if (cmod3(1:ib3) .eq. 'access_cm2_05x05') cs3 = 'AC2-'
      if (cmod4(1:ib4) .eq. 'access_cm2_05x05') cs4 = 'AC2-'
c
      if (cmod1(1:ib1) .eq. 'access_esm1_5_05x05') cs1 = 'AES-'
      if (cmod2(1:ib2) .eq. 'access_esm1_5_05x05') cs2 = 'AES-'
      if (cmod3(1:ib3) .eq. 'access_esm1_5_05x05') cs3 = 'AES-'
      if (cmod4(1:ib4) .eq. 'access_esm1_5_05x05') cs4 = 'AES-'
c
      if (cmod1(1:ib1) .eq. 'cesm2_05x05') cs1 = 'CES-'
      if (cmod2(1:ib2) .eq. 'cesm2_05x05') cs2 = 'CES-'
      if (cmod3(1:ib3) .eq. 'cesm2_05x05') cs3 = 'CES-'
      if (cmod4(1:ib4) .eq. 'cesm2_05x05') cs4 = 'CES-'
c
      if (cmod1(1:ib1) .eq. 'cesm2_waccm_05x05') cs1 = 'CEW-'
      if (cmod2(1:ib2) .eq. 'cesm2_waccm_05x05') cs2 = 'CEW-'
      if (cmod3(1:ib3) .eq. 'cesm2_waccm_05x05') cs3 = 'CEW-'
      if (cmod4(1:ib4) .eq. 'cesm2_waccm_05x05') cs4 = 'CEW-'
c
      if (cmod1(1:ib1) .eq. 'ciesm_05x05') cs1 = 'CIE-'
      if (cmod2(1:ib2) .eq. 'ciesm_05x05') cs2 = 'CIE-'
      if (cmod3(1:ib3) .eq. 'ciesm_05x05') cs3 = 'CIE-'
      if (cmod4(1:ib4) .eq. 'ciesm_05x05') cs4 = 'CIE-'
c
      if (cmod1(1:ib1) .eq. 'cnrm_cm6_1_05x05') cs1 = 'CNR-'
      if (cmod2(1:ib2) .eq. 'cnrm_cm6_1_05x05') cs2 = 'CNR-'
      if (cmod3(1:ib3) .eq. 'cnrm_cm6_1_05x05') cs3 = 'CNR-'
      if (cmod4(1:ib4) .eq. 'cnrm_cm6_1_05x05') cs4 = 'CNR-'
c
      if (cmod1(1:ib1) .eq. 'ec_earth3_05x05') cs1 = 'ECE-'
      if (cmod2(1:ib2) .eq. 'ec_earth3_05x05') cs2 = 'ECE-'
      if (cmod3(1:ib3) .eq. 'ec_earth3_05x05') cs3 = 'ECE-'
      if (cmod4(1:ib4) .eq. 'ec_earth3_05x05') cs4 = 'ECE-'
c
      if (cmod1(1:ib1) .eq. 'ec_earth3_veg_05x05') cs1 = 'ECV-'
      if (cmod2(1:ib2) .eq. 'ec_earth3_veg_05x05') cs2 = 'ECV-'
      if (cmod3(1:ib3) .eq. 'ec_earth3_veg_05x05') cs3 = 'ECV-'
      if (cmod4(1:ib4) .eq. 'ec_earth3_veg_05x05') cs4 = 'ECv-'
c
      if (cmod1(1:ib1) .eq. 'gfdl_cm4_05x05') cs1 = 'GFC-'
      if (cmod2(1:ib2) .eq. 'gfdl_cm4_05x05') cs2 = 'GFC-'
      if (cmod3(1:ib3) .eq. 'gfdl_cm4_05x05') cs3 = 'GFC-'
      if (cmod4(1:ib4) .eq. 'gfdl_cm4_05x05') cs4 = 'GFC-'
c
      if (cmod1(1:ib1) .eq. 'gfdl_esm4_05x05') cs1 = 'GFE-'
      if (cmod2(1:ib2) .eq. 'gfdl_esm4_05x05') cs2 = 'GFE-'
      if (cmod3(1:ib3) .eq. 'gfdl_esm4_05x05') cs3 = 'GFE-'
      if (cmod4(1:ib4) .eq. 'gfdl_esm4_05x05') cs4 = 'GFE-'
c
      if (cmod1(1:ib1) .eq. 'hadgem3_gc31_ll_05x05') cs1 = 'HDL-'
      if (cmod2(1:ib2) .eq. 'hadgem3_gc31_ll_05x05') cs2 = 'HDL-'
      if (cmod3(1:ib3) .eq. 'hadgem3_gc31_ll_05x05') cs3 = 'HDL-'
      if (cmod4(1:ib4) .eq. 'hadgem3_gc31_ll_05x05') cs4 = 'HDL-'
c
      if (cmod1(1:ib1) .eq. 'hadgem3_gc31_mm_05x05') cs1 = 'HDM-'
      if (cmod2(1:ib2) .eq. 'hadgem3_gc31_mm_05x05') cs2 = 'HDM-'
      if (cmod3(1:ib3) .eq. 'hadgem3_gc31_mm_05x05') cs3 = 'HDM-'
      if (cmod4(1:ib4) .eq. 'hadgem3_gc31_mm_05x05') cs4 = 'HDM-'
c
      if (cmod1(1:ib1) .eq. 'ipsl_cm6a_lr_05x05') cs1 = 'IPS-'
      if (cmod2(1:ib2) .eq. 'ipsl_cm6a_lr_05x05') cs2 = 'IPS-'
      if (cmod3(1:ib3) .eq. 'ipsl_cm6a_lr_05x05') cs3 = 'IPS-'
      if (cmod4(1:ib4) .eq. 'ipsl_cm6a_lr_05x05') cs4 = 'IPS-'
c
      if (cmod1(1:ib1) .eq. 'miroc6_05x05') cs1 = 'MI6-'
      if (cmod2(1:ib2) .eq. 'miroc6_05x05') cs2 = 'MI6-'
      if (cmod3(1:ib3) .eq. 'miroc6_05x05') cs3 = 'MI6-'
      if (cmod4(1:ib4) .eq. 'miroc6_05x05') cs4 = 'MI6-'
c
      if (cmod1(1:ib1) .eq. 'miroc_es2l_05x05') cs1 = 'MIE-'
      if (cmod2(1:ib2) .eq. 'miroc_es2l_05x05') cs2 = 'MIE-'
      if (cmod3(1:ib3) .eq. 'miroc_es2l_05x05') cs3 = 'MIE-'
      if (cmod4(1:ib4) .eq. 'miroc_es2l_05x05') cs4 = 'MIE-'
c
      if (cmod1(1:ib1) .eq. 'mpi_esm1_2_hr_05x05') cs1 = 'MPH-'
      if (cmod2(1:ib2) .eq. 'mpi_esm1_2_hr_05x05') cs2 = 'MPH-'
      if (cmod3(1:ib3) .eq. 'mpi_esm1_2_hr_05x05') cs3 = 'MPH-'
      if (cmod4(1:ib4) .eq. 'mpi_esm1_2_hr_05x05') cs4 = 'MPH-'
c
      if (cmod1(1:ib1) .eq. 'mpi_esm1_2_lr_05x05') cs1 = 'MPL-'
      if (cmod2(1:ib2) .eq. 'mpi_esm1_2_lr_05x05') cs2 = 'MPL-'
      if (cmod3(1:ib3) .eq. 'mpi_esm1_2_lr_05x05') cs3 = 'MPL-'
      if (cmod4(1:ib4) .eq. 'mpi_esm1_2_lr_05x05') cs4 = 'MPL-'
c
      if (cmod1(1:ib1) .eq. 'mri_esm2_0_05x05') cs1 = 'MRI-'
      if (cmod2(1:ib2) .eq. 'mri_esm2_0_05x05') cs2 = 'MRI-'
      if (cmod3(1:ib3) .eq. 'mri_esm2_0_05x05') cs3 = 'MRI-'
      if (cmod4(1:ib4) .eq. 'mri_esm2_0_05x05') cs4 = 'MRI-'
c
      if (cmod1(1:ib1) .eq. 'nesm3_05x05') cs1 = 'NES-'
      if (cmod2(1:ib2) .eq. 'nesm3_05x05') cs2 = 'NES-'
      if (cmod3(1:ib3) .eq. 'nesm3_05x05') cs3 = 'NES-'
      if (cmod4(1:ib4) .eq. 'nesm3_05x05') cs4 = 'NES-'
c
      if (cmod1(1:ib1) .eq. 'noresm2_mm_05x05') cs1 = 'NOR-'
      if (cmod2(1:ib2) .eq. 'noresm2_mm_05x05') cs2 = 'NOR-'
      if (cmod3(1:ib3) .eq. 'noresm2_mm_05x05') cs3 = 'NOR-'
      if (cmod4(1:ib4) .eq. 'noresm2_mm_05x05') cs4 = 'NOR-'
c
      if (cmod1(1:ib1) .eq. 'ukesm1_0_ll_05x05') cs1 = 'UKE-'
      if (cmod2(1:ib2) .eq. 'ukesm1_0_ll_05x05') cs2 = 'UKE-'
      if (cmod3(1:ib3) .eq. 'ukesm1_0_ll_05x05') cs3 = 'UKE-'
      if (cmod4(1:ib4) .eq. 'ukesm1_0_ll_05x05') cs4 = 'UKE-'
c
c-------------------------------------------------------------------------------
c     ** 1.9 **
c
c     Using the selected character strings, generate the substrings required for
c     the full output file names.
c-------------------------------------------------------------------------------
c
      str0 = 'evs_NOPT_'
      str1 = cs1//cexp1(1:ic1)//'_OBS_'//cs2//cexp2(1:ic2)
      str2 = '_FING_'//cs3//cexp3(1:ic3)//'_NSTAT_'
      str3 = cs4//cexp4(1:ic4)
c
      print *, 'str0 = ', str0
      print *, 'str1 = ', str1
      print *, 'str2 = ', str2
      print *, 'str3 = ', str3
c
      id1 = index(str1, ' ') - 1
      id2 = index(str2, ' ') - 1
      id3 = index(str3, ' ') - 1
c
      str4 = str0//str1(1:id1)//str2(1:id2)//str3(1:id3)
c
      id4 = index(str4, ' ') - 1
c
      print *, 'str4 = ', str4
c
      str5 = '_'//varnam1(1:ia1)
      id5  = index(str5, ' ') - 1
c
      print *, 'str5 = ', str5
c
      str6 = str4(1:id4)//str5(1:id5)
      id6  = index(str6, ' ') - 1
c
      print *, 'str6 = ', str6
c
      if (itrunc .lt. 10) then
        write(cc1, '(i1)') itrunc
      else if (itrunc .ge. 10) then
        write(cc2, '(i2)') itrunc
      end if
c
      if (nfinger .lt. 10) then
        write(cc3, '(i1)') nfinger
      else
        print *, 'Error in RESIDUAL_FINGER1!'
        print *, 'Number of EOFs for fingerprint cannot exceed 9!'
        print *, 'Selected NFINGER = ', nfinger
        stop
      end if
c
      write(cc4, '(a4, i1, a4, i1)') '_sig', icritsn, '_vwt', ivweight
c
      print *, 'cc4 = ', cc4
c
c-------------------------------------------------------------------------------
c     ** 1.10 **
c
c     Now create full names for output files. There are two files:
c
c     NAME         TYPE
c
c     drsdic1, drsdat1     Empirical Orthogonal Functions (EOFs)
c     outfil1              Principal Components (PCs)
c-------------------------------------------------------------------------------
c
      if (itrunc .lt. 10) then
        drsdic1 = str6(1:id6)//'_itr0'//cc1//'_nf0'//cc3//cc4//'.dic'
      else if (itrunc .ge. 10) then
        drsdic1 = str6(1:id6)//'_itr'//cc2//'_nf0'//cc3//cc4//'.dic'
      end if
c
      id7            = index(drsdic1, ' ') - 5
      drsdat1(1:id7) = drsdic1(1:id7)
      drsdat1        = drsdat1(1:id7)//'.dat'
c 
      if (varnam1(1:3) .eq. 'tt6') then
        outfil1(1:3)   = 'pcs'
      else if (varnam1(1:3) .eq. 'sm6') then
        outfil1(1:3)   = 'pc1'
      end if
c
      outfil1(4:id7) = drsdic1(4:id7)
      outfil1        = outfil1(1:id7)//'.d'
c
      print *, ' ******************************************************'
      print *, ' YOU WILL GENERATE THE FOLLOWING OUTPUT FILES: '
      print *, ' '
      print *, ' drsdic1 = ', drsdic1
      print *, ' drsdat1 = ', drsdat1
      print *, ' outfil1 = ', outfil1
      print *, ' ******************************************************'
c
c
c
c
c
c
c-------------------------------------------------------------------------------
c     ** 2.0 **
c
c     READ DATA.
c
c     Set the weights appropriate for the retrieved data. After regridding, all
c     input datasets should be in the form of zonal averages on a regular 5 
c     degree latitude grid, for three different atmospheric layers (tls, tmt,
c     and tlt). All retrieve area weights (for individual latitude bands) are 
c     cosine weights.
c-------------------------------------------------------------------------------
c
      yweights1 = 'cosine'
      yweights2 = 'cosine'
      yweights3 = 'cosine'
      yweights4 = 'cosine'
c
      print *, ' ******************************************************'
      print *, ' START SETTING WEIGHTS!'
      print *, ' ******************************************************'
      print *, ' '
      print *, ' yweights1 = ', yweights1
      print *, ' yweights2 = ', yweights2
      print *, ' yweights3 = ', yweights3
      print *, ' yweights4 = ', yweights4
      print *, ' '
      print *, ' ******************************************************'
      print *, ' FINISH SETTING WEIGHTS!'
      print *, ' ******************************************************'
      print *, ' '
c
c-------------------------------------------------------------------------------
c     ** 2.1 **
c
c     Retrieve data.
c
c     iopen  = 0   Initialization of EzGet outside subroutine RD3D_YHT_NOMASK.
c     iopen  = 1   Initialization of EzGet within subroutine RD3D_YHT_NOMASK.
c
c     NOTE: Since retrieved data are on a common uniform grid, the latitude  
c     and pressure coordinates (rlat, rlev) and their actual dimensions (iyout, 
c     ihout) should be identical for SET1, SET2, SET3, and SET4! Test that the
c     pressure coordinates are identical.
c-------------------------------------------------------------------------------
c
      print *, ' ******************************************************'
      print *, ' START READING INPUT DATA!'
      print *, ' ******************************************************'
      print *, ' '
c
      iopen  = 1
c
      rt1 = iy1 * 1.0
      rt2 = iy2 * 1.0
      rt3 = iy3 * 1.0
      rt4 = iy4 * 1.0
      rt5 = iy5 * 1.0
      rt6 = iy6 * 1.0
      rt7 = iy7 * 1.0
      rt8 = iy8 * 1.0
c
      call rd3d_yht_nomask(iopen, idomain, ny, nh, n1, rt1, rt2, 
     &infil1, rlat, rlev1, rtime1, awts1, bas1, varnam1, yweights1, 
     &iyout, ihout, nt1, yy1, yy2)       
c
      print *, 'SET1! iyout, ihout, nt1 = ', iyout, ihout, nt1
c
      call rd3d_yht_nomask(iopen, idomain, ny, nh, n2, rt3, rt4,
     &infil2, rlat, rlev2, rtime2, awts2, pro1, varnam1, yweights2, 
     &iyout, ihout, nt2, yy1, yy2)       
c
      print *, 'SET2! iyout, ihout, nt2 = ', iyout, ihout, nt2
c
      call rd3d_yht_nomask(iopen, idomain, ny, nh, n3, rt5, rt6,
     &infil3, rlat, rlev3, rtime3, awts3, gss1, varnam1, yweights3, 
     &iyout, ihout, nt3, yy1, yy2)       
c
      print *, 'SET3! iyout, ihout, nt3 = ', iyout, ihout, nt3
c
      call rd3d_yht_nomask(iopen, idomain, ny, nh, n4, rt7, rt8,
     &infil4, rlat, rlev4, rtime4, awts4, dtt1, varnam1, yweights4, 
     &iyout, ihout, nt4, yy1, yy2)       
c
      print *, 'SET4! iyout, ihout, nt4 = ', iyout, ihout, nt4
      print *, 'Latitude boundaries yy1, yy2 = ', yy1, yy2
c
      print *, ' ******************************************************'
      print *, ' FINISH READING INPUT DATA!'
      print *, ' ******************************************************'
      print *, ' '
c
      do 212 k = 1, ihout ! Loop over number of layers.
        do 211 j2 = 1, 4    ! Loop over number of input data sets.
          do 210 j1 = 1, 4    ! Loop over number of input data sets.
            if (j1. eq. 1) term1 = rlev1(k) 
            if (j1. eq. 2) term1 = rlev2(k) 
            if (j1 .eq. 3) term1 = rlev3(k) 
            if (j1 .eq. 4) term1 = rlev4(k) 
            if (j2. eq. 1) term2 = rlev1(k) 
            if (j2. eq. 2) term2 = rlev2(k) 
            if (j2 .eq. 3) term2 = rlev3(k) 
            if (j2 .eq. 4) term2 = rlev4(k) 
            if (term1 .ne. term2) then
              print *, 'Pressure coordinates differ!'
              print *, 'term1, term2 = ', term1, term2
cBDS          stop
            end if
  210     continue
  211   continue
  212 continue

c-------------------------------------------------------------------------------
c     ** 2.2 **
c
c     Check for missing data. Since all datasets are in the form of zonal
c     averages on a common (regular) 5 degree latitude grid (for three different
c     atmospheric layers), the area weights are identical for SET1, SET2, SET3,
c     and SET4. The area weights are stored for SET1 only.
c-------------------------------------------------------------------------------
c
      print *, ' ******************************************************'
      print *, ' START CHECKING FOR MISSING DATA!'
      print *, ' ******************************************************'
      print *, ' '
c
      do 223 nn = 1, 4
        icnt1 = 0
        sum1  = 0.0
        sum2  = 0.0
c
        if (nn .eq. 1) kyear = nt1
        if (nn .eq. 2) kyear = nt2
        if (nn .eq. 3) kyear = nt3
        if (nn .eq. 4) kyear = nt4
c
        do 222 k = 1, kyear
          do 221 j = 1, ihout
            do 220 i = 1, iyout
c
              if (nn .eq. 1) term =  bas1(i, j, k)
              if (nn .eq. 2) term =  pro1(i, j, k)
              if (nn .eq. 3) term =  gss1(i, j, k)
              if (nn .eq. 4) term =  dtt1(i, j, k)
              if (nn .eq. 1) aa   = awts1(i, j, k)
              if (nn .eq. 2) aa   = awts2(i, j, k)
              if (nn .eq. 3) aa   = awts3(i, j, k)
              if (nn .eq. 4) aa   = awts4(i, j, k)
c
              if (k .eq. 1) sum2 = sum2 + aa
c
              if (term .ne. bogus) then
                sum1 = sum1 + term * aa
              else
                icnt1 = icnt1 + 1
              end if
  220       continue
  221     continue
  222   continue
c
        if (icnt1 .gt. 0) then
          if (nn .eq. 1) print *, 'SET1 has ', icnt1, ' missing values'
          if (nn .eq. 2) print *, 'SET2 has ', icnt1, ' missing values'
          if (nn .eq. 3) print *, 'SET3 has ', icnt1, ' missing values'
          if (nn .eq. 4) print *, 'SET4 has ', icnt1, ' missing values'
        else
          if (nn .eq. 1) print *, 'SET1 has no missing values'
          if (nn .eq. 2) print *, 'SET2 has no missing values'
          if (nn .eq. 3) print *, 'SET3 has no missing values'
          if (nn .eq. 4) print *, 'SET4 has no missing values'
        end if
c
        if (nn .eq. 1) print *, 'Grand mean of SET1: ', sum1 / kyear  
        if (nn .eq. 2) print *, 'Grand mean of SET2: ', sum1 / kyear  
        if (nn .eq. 3) print *, 'Grand mean of SET3: ', sum1 / kyear   
        if (nn .eq. 4) print *, 'Grand mean of SET4: ', sum1 / kyear   
c
        if (nn .eq. 1) print *, 'Sum of weights for SET1: ', sum2
        if (nn .eq. 2) print *, 'Sum of weights for SET2: ', sum2 
        if (nn .eq. 3) print *, 'Sum of weights for SET3: ', sum2
        if (nn .eq. 4) print *, 'Sum of weights for SET4: ', sum2
  223 continue
c
      print *, ' ******************************************************'
      print *, ' FINISH CHECKING FOR MISSING DATA!'
      print *, ' ******************************************************'
      print *, ' '
c
c-------------------------------------------------------------------------------
c     ** 2.3 **
c
c     Echo values of latitude and pressure coordinates.
c-------------------------------------------------------------------------------
c
      print *, ' '
      print *, ' ******************************************************'
      print *, ' ECHOING COORDINATES FOR 1st DIMENSION!'
      print *,   (rlat(i), i = 1, iyout)
      print *, ' ******************************************************'
      print *, ' '
      print *, ' '
      print *, ' ******************************************************'
      print *, ' ECHOING COORDINATES FOR 2nd DIMENSION!'
      print *,   (rlev1(j), j = 1, ihout)
      print *, ' ******************************************************'
      print *, ' '
c
c-------------------------------------------------------------------------------
c     ** 2.4 **
c
c     Now read zonally-averaged land area fraction for SET1, SET2, SET3, and 
c     SET4 data.
c-------------------------------------------------------------------------------
c
      print *, ' ******************************************************'
      print *, ' START READING LAND AREA FRACTION DATA!'
      print *, ' ******************************************************'
c
      call rdmask(ny, nh, iyout, ihout, idomain, pmask1, store1, 
     &fmask1, yweights1)

      call rdmask(ny, nh, iyout, ihout, idomain, pmask2, store1, 
     &fmask2, yweights2)

      call rdmask(ny, nh, iyout, ihout, idomain, pmask3, store1, 
     &fmask3, yweights3)

      call rdmask(ny, nh, iyout, ihout, idomain, pmask4, store1, 
     &fmask4, yweights4)
c
      print *, ' ******************************************************'
      print *, ' FINISH READING LAND AREA FRACTION DATA!'
      print *, ' ******************************************************'
c
c-------------------------------------------------------------------------------
c     ** 2.5 **
c
c     For each atmospheric layer, calculate the global-mean land percentage in 
c     each of the four land area fraction datasets.
c-------------------------------------------------------------------------------
c
      do 252 nn = 1, 4
c
        do 251 j = 1, ihout
          sum1 = 0.0
          do 250 i = 1, iyout
c
            if (nn .eq. 1) term = fmask1(i, j)
            if (nn .eq. 2) term = fmask2(i, j)
            if (nn .eq. 3) term = fmask3(i, j)
            if (nn .eq. 4) term = fmask4(i, j)
            sum1 = sum1 + (term * store1(i, j))
c
  250     continue
          print *, 'Land %, SET', nn, ' LEV', j, ' = ', sum1 
  251   continue
c
  252 continue
c
c-------------------------------------------------------------------------------
c     ** 2.6 **
c
c     Calculate the geometrical mean land area fraction (LAF). This is done for
c     each individual atmospheric layer. Note that:
c
c     "ALL LAND"  = 100.0%
c     "ALL OCEAN" =   0.0%
c
c     Note that in calculating the geometrical mean, a grid-point with an
c     ocean fraction of less than 1% in any of the four input datasets will 
c     be assigned the value zero. This automatically sets that grid-point's 
c     value to zero in the geometrical mean mask.
c-------------------------------------------------------------------------------
c
      print *, ' ******************************************************'
      print *, ' START CALCULATING GEOMETRICAL MEAN LAF!'
      print *, ' ******************************************************'
c
      epsilon = 1.0E-2
      do 261 j = 1, ihout
        print *, 'Atmospheric layer = ', rlev1(j), ' hPa'
        isum1 = 0
        isum2 = 0
        isum3 = 0
        do 260 i = 1, iyout
          term1 = 1.0 - (fmask1(i, j) / 100.0)
          term2 = 1.0 - (fmask2(i, j) / 100.0)
          term3 = 1.0 - (fmask3(i, j) / 100.0)
          term4 = 1.0 - (fmask4(i, j) / 100.0)
c
          if (term1 .lt. epsilon) term1 = 0.0
          if (term2 .lt. epsilon) term2 = 0.0
          if (term3 .lt. epsilon) term3 = 0.0
          if (term4 .lt. epsilon) term4 = 0.0
c
          term5 = (term1 * term2 * term3 * term4) 
          term6 = term5 ** (1.0 / 4.0)
          print *, 'terms = ', term1, term2, term3, term4, term6
c
          if (term6 .lt. epsilon) then  ! Case of "pure" land grid-point.
            isum1       = isum1 + 1
          else if (term6 .eq. 1.0) then ! Case of "pure" ocean grid-point.
            isum2       = isum2 + 1
          else                          ! Fractional land.
            isum3       = isum3 + 1
          end if
c
          gmask(i, j) = term6 ! Assign the geometrical mean LAF to gmask.
c
  260   continue
        print *, 'Pure land points, geometrical mean mask    = ', isum1
        print *, 'Pure ocean points, geometrical mean mask   = ', isum2
        print *, 'Fractl. land points, geometrical mean mask = ', isum3
c
  261 continue
c
      print *, ' ******************************************************'
      print *, ' FINISH CALCULATING GEOMETRICAL MEAN LAF!'
      print *, ' ******************************************************'
c
c-------------------------------------------------------------------------------
c     ** 2.7 **
c
c     Check each dataset's coverage for compatability with other datasets. At 
c     each latitude band where SET1 has missing data, assign a missing data 
c     value to SET2, SET3, and SET4. Do this for all time points (this assumes 
c     that coverage does not change over time). Then repeat for SET2, etc.
c-------------------------------------------------------------------------------
c
      print *, ' ******************************************************'
      print *, ' START ENSURING DATA COVERAGE COMPATABILITY!'
      print *, ' ******************************************************'
      print *, ' '
c
      do 272 nn = 1, 4
        do 271 j = 1, ihout
          do 270 i = 1, iyout
c
            if (nn .eq. 1) term1 = bas1(i, j, 1)
            if (nn .eq. 2) term1 = pro1(i, j, 1)
            if (nn .eq. 3) term1 = gss1(i, j, 1)
            if (nn .eq. 4) term1 = dtt1(i, j, 1)
c
            if (nn .eq. 1) then ! SET1 coverage sets coverage for 2/3/4.
              if (term1 .eq. bogus) then 
                do 2700 k = 1, nt2
                  pro1(i, j, k) = bogus
 2700           continue
                do 2701 k = 1, nt3
                  gss1(i, j, k) = bogus
 2701           continue
                do 2702 k = 1, nt4
                  dtt1(i, j, k) = bogus
 2702           continue
              end if
c
            else if (nn .eq. 2) then ! SET2 coverage sets coverage for 1/3/4.
              if (term1 .eq. bogus) then 
                do 2703 k = 1, nt1
                  bas1(i, j, k) = bogus
 2703           continue
                do 2704 k = 1, nt3
                  gss1(i, j, k) = bogus
 2704           continue
                do 2705 k = 1, nt4
                  dtt1(i, j, k) = bogus
 2705           continue
              end if
c
            else if (nn .eq. 3) then ! SET3 coverage sets coverage for 1/2/4.
              if (term1 .eq. bogus) then 
                do 2706 k = 1, nt1
                  bas1(i, j, k) = bogus
 2706           continue
                do 2707 k = 1, nt2
                  pro1(i, j, k) = bogus
 2707           continue
                do 2708 k = 1, nt4
                  dtt1(i, j, k) = bogus
 2708           continue
              end if
c
            else if (nn .eq. 4) then ! SET4 coverage sets coverage for 1/2/3.
              if (term1 .eq. bogus) then
                do 2709 k = 1, nt1
                  bas1(i, j, k) = bogus
 2709           continue
                do 2710 k = 1, nt2
                  pro1(i, j, k) = bogus
 2710           continue
                do 2711 k = 1, nt3
                  gss1(i, j, k) = bogus
 2711           continue
              end if
            end if
c
  270     continue
  271   continue
  272 continue
c
      print *, ' ******************************************************'
      print *, ' FINISH ENSURING DATA COVERAGE COMPATABILITY!'
      print *, ' ******************************************************'
      print *, ' '
c
c-------------------------------------------------------------------------------
c     ** 2.8 **
c
c     First normalize the mass weights calculated by Cheng-Zhi on 10/31/2022.
c     Normalization is by the sum of the mass weights for the six individual
c     atmospheric layers. The weights are the convolution of each weighting
c     function with the standard atmosphere profile of atmospheric mass. The
c     normalized mass weights are stored in rmassn. 
c
c     For ttt and tlt, there are separate weighting functions for land and 
c     ocean; the overall ttt and tlt mass weights are the weighted averages
c     of the individual mass weights for the land and ocean weighting functions, 
c     with weighting by the global land and ocean fractions. 
c
c     Normalized mass weights are only applied in the ivweight=4 case, For all 
c     other currently implemented cases (ivweight=0/1/2/3/5), only area weights
c     are applied. For ivweight=4, the weight applied to each latitude band of
c     each layer is the product of the area weight and mass weight for that
c     latitude band and layer. 
c
c     After ensuring data coverage compatibility in all four input datasets,
c     any one dataset can be used to set the common area weight mask (cmask). 
c     Assign cmask(i, j) values to bogus at corresponding bas1(i, j, 1) values 
c     with missing data. For bas1(i, j, 1) values with valid data, assign the 
c     area weight of the latitude band to cmask(i, j).
c-------------------------------------------------------------------------------
c
      tttwt = (4.775149E-01 * 0.29) + (4.805052E-01 * 0.71)
      tltwt = (6.132900E-01 * 0.29) + (6.156208E-01 * 0.71)
c
      print *, 'tttwt = ', tttwt
      print *, 'tltwt = ', tltwt
c
      rmass(1) = 4.877587E-03  ! Mass weight for SSU3.
      rmass(2) = 1.075615E-02  ! Mass weight for SSU2.
      rmass(3) = 2.556198E-02  ! Mass weight for SSU1.
      rmass(4) = 8.081517E-02  ! Mass weight for TLS. 
      rmass(5) = tttwt         ! Mass weight for TTT (land+ocean). 
      rmass(6) = tltwt         ! Mass weight for TLT (land+ocean). 
c
      sumvweight1 = 0.0
      do 280 j = 1, ihout ! Loop over atmospheric layers.
       sumvweight1 = sumvweight1 + rmass(j) ! Calculate sum of mass weights.
  280 continue            
c 
      sumvweight2 = 0.0
      do 281 j = 1, ihout ! Loop over atmospheric layers.
        rmassn(j)   = rmass(j) / sumvweight1  ! Normalize by sum of mass weights.
        sumvweight2 = sumvweight2 + rmassn(j) ! Calculate sum of mass weights.
        print *, 'rmassn * 100, j = ', (rmassn(j) *100.0), j
  281 continue            
c 
      print *, 'sumvweight1 = ', sumvweight1
      print *, 'sumvweight2 = ', sumvweight2
c
      do 283 j = 1, ihout ! Loop over atmospheric layers.
        if (ivweight .eq. 0) then      ! Case of no weighting.
          if (j .eq. 1) vwt = 1.0  ! Weight for ssu3.
          if (j .eq. 2) vwt = 1.0  ! Weight for ssu2.
          if (j .eq. 3) vwt = 1.0  ! Weight for ssu1.
          if (j .eq. 4) vwt = 1.0  ! Weight for tls.
          if (j .eq. 5) vwt = 1.0  ! Weight for tf1.
          if (j .eq. 6) vwt = 1.0  ! Weight for tlt.
        else if (ivweight .eq. 1) then ! Case of tropospheric layers only.
          if (j .eq. 1) vwt = 0.0  ! Weight for ssu3.
          if (j .eq. 2) vwt = 0.0  ! Weight for ssu2.
          if (j .eq. 3) vwt = 0.0  ! Weight for ssu1.
          if (j .eq. 4) vwt = 0.0  ! Weight for tls.
          if (j .eq. 5) vwt = 1.0  ! Weight for tf1.
          if (j .eq. 6) vwt = 1.0  ! Weight for tlt.
        else if (ivweight .eq. 2) then ! Case of SSU only.
          if (j .eq. 1) vwt = 1.0  ! Weight for ssu3.
          if (j .eq. 2) vwt = 1.0  ! Weight for ssu2.
          if (j .eq. 3) vwt = 1.0  ! Weight for ssu1.
          if (j .eq. 4) vwt = 0.0  ! Weight for tls.
          if (j .eq. 5) vwt = 0.0  ! Weight for tf1.
          if (j .eq. 6) vwt = 0.0  ! Weight for tlt.
        else if (ivweight .eq. 3) then ! Case of MSU only.                    
          if (j .eq. 1) vwt = 0.0  ! Weight for ssu3.
          if (j .eq. 2) vwt = 0.0  ! Weight for ssu2.
          if (j .eq. 3) vwt = 0.0  ! Weight for ssu1.
          if (j .eq. 4) vwt = 1.0  ! Weight for tls.
          if (j .eq. 5) vwt = 1.0  ! Weight for tf1.
          if (j .eq. 6) vwt = 1.0  ! Weight for tlt.
        else if (ivweight .eq. 4) then ! Case of mass weighting.              
          vwt = rmassn(j)                   
        else if (ivweight .eq. 5) then ! Case of MIN_olap.                    
          if (j .eq. 1) vwt = 1.0  ! Weight for ssu3.
          if (j .eq. 2) vwt = 0.0  ! Weight for ssu2.
          if (j .eq. 3) vwt = 0.0  ! Weight for ssu1.
          if (j .eq. 4) vwt = 1.0  ! Weight for tls.
          if (j .eq. 5) vwt = 0.0  ! Weight for tf1.
          if (j .eq. 6) vwt = 1.0  ! Weight for tlt.
        end if
c
        do 282 i = 1, iyout ! Loop over latitude bands. 
          if (bas1(i, j, 1) .eq. bogus) then
            cmask(i, j) = bogus        
          else
            if (ivweight .lt. 4) then ! Use area weights only.
              if (vwt .eq. 1.0) cmask(i, j) = awts1(i, j, 1)
              if (vwt .eq. 0.0) cmask(i, j) = bogus
            else if (ivweight .eq. 4) then ! Product of area and mass weights.
              cmask(i, j) = awts1(i, j, 1) * vwt
            else if (ivweight .eq. 5) then ! MIN_olap case.                   
              if (vwt .eq. 1.0) cmask(i, j) = awts1(i, j, 1)
              if (vwt .eq. 0.0) cmask(i, j) = bogus
            end if
          end if
c
        print *, 'LEV, LAT, CMASK= ', j, i, cmask(i, j)
c
  282   continue ! End loop over latitude bands.
  283 continue   ! End loop over atmospheric layers. 
c
c
c
c
c
c
c-------------------------------------------------------------------------------
c     ** 3.0 **
c
c     COMPUTE CLIMATOLOGICAL MEANS AND ANOMALIES.
c
c     First set start and end years for reference periods. 
c-------------------------------------------------------------------------------
c
      print *, 'rtime1(1) = ', rtime1(1)
      print *, 'rtime2(1) = ', rtime2(1)
      print *, 'rtime3(1) = ', rtime3(1)
      print *, 'rtime4(1) = ', rtime4(1)
c
      print *, 'rtime1(1) = ', int(rtime1(1))
      print *, 'rtime2(1) = ', int(rtime2(1))
      print *, 'rtime3(1) = ', int(rtime3(1))
      print *, 'rtime4(1) = ', int(rtime4(1))
c
      jyr1 = iref1 - int(rtime1(1)) + 1
      jyr2 = iref2 - int(rtime1(1)) + 1
      jyr3 = iref3 - int(rtime2(1)) + 1
      jyr4 = iref4 - int(rtime2(1)) + 1
      jyr5 = iref5 - int(rtime3(1)) + 1
      jyr6 = iref6 - int(rtime3(1)) + 1
      jyr7 = iref7 - int(rtime4(1)) + 1
      jyr8 = iref8 - int(rtime4(1)) + 1
c
c-------------------------------------------------------------------------------
c     ** 3.1 **
c
c     Compute climatological annual means and anomalies for input data.
c-------------------------------------------------------------------------------
c
      print *, ' ******************************************************'
      print *, ' START COMPUTING ANOMALIES!'
      print *, ' ******************************************************'
      print *, ' '
c
      mcrit = 0
c
      do 314 nn = 1, 4
        if (nn .eq. 1) then
          kyear = nt1
          jj1   = jyr1
          jj2   = jyr2
          if ((iref1 .eq. 0) .or. (iref2 .eq. 0)) goto 314
        else if (nn .eq. 2) then
          kyear = nt2
          jj1   = jyr3
          jj2   = jyr4
          if ((iref3 .eq. 0) .or. (iref4 .eq. 0)) goto 314
        else if (nn .eq. 3) then
          kyear = nt3
          jj1   = jyr5
          jj2   = jyr6
          if ((iref5 .eq. 0) .or. (iref6 .eq. 0)) goto 314
        else if (nn .eq. 4) then
          kyear = nt4
          jj1   = jyr7
          jj2   = jyr8
          if ((iref7 .eq. 0) .or. (iref8 .eq. 0)) goto 315
        end if
c
        if (kyear .gt. n1) then
          print *, 'Parameter n1 is smaller than n2, n3, or n4!'
          print *, 'Increase value of n1 in parameter statement!'
          print *, 'n1, n2, n3, n4 = ', n1, n2, n3, n4
          stop
        end if
c
        print *, 'jj1, jj2, kyear = ', jj1, jj2, kyear
c
        do 313 j = 1, ihout
          do 312 i = 1, iyout
            do 310 k = 1, kyear
              if (nn .eq. 1) xx(k) = bas1(i, j, k)
              if (nn .eq. 2) xx(k) = pro1(i, j, k)
              if (nn .eq. 3) xx(k) = gss1(i, j, k)
              if (nn .eq. 4) xx(k) = dtt1(i, j, k)
  310       continue
c
            call climanom1(n1, kyear, jj1, jj2, mcrit, bogus,
     &      xx, xm)
c
            do 311 k = 1, kyear
              if (nn .eq. 1) bas1(i, j, k) = xx(k)
              if (nn .eq. 2) pro1(i, j, k) = xx(k)
              if (nn .eq. 3) gss1(i, j, k) = xx(k)
              if (nn .eq. 4) dtt1(i, j, k) = xx(k)
              if ((k .eq. kyear) .and. (nn .eq. 1)) then
                print *, 'bas1(i, j, k) = ', bas1(i, j, k), i, j
              else if ((k .eq. kyear) .and. (nn .eq. 2)) then
                print *, 'pro1(i, j, k) = ', pro1(i, j, k), i, j
              else if ((k .eq. kyear) .and. (nn .eq. 3)) then
                print *, 'gss1(i, j, k) = ', gss1(i, j, k), i, j
              else if ((k .eq. kyear) .and. (nn .eq. 4)) then
                print *, 'dtt1(i, j, k) = ', dtt1(i, j, k), i, j
              end if
  311       continue
  312     continue
  313   continue
  314 continue
  315 continue
c
      print *, ' ******************************************************'
      print *, ' FINISH COMPUTING ANOMALIES!'
      print *, ' ******************************************************'
      print *, ' '
c
c-------------------------------------------------------------------------------
c     ** 3.2 **
c
c     Check the sum of the anomalies for each input data set.
c-------------------------------------------------------------------------------
c
      print *, ' ******************************************************'
      print *, ' START CHECKING SUM OF ANOMALIES!'
      print *, ' ******************************************************'
      print *, ' '
c
      do 323 nn = 1, 4
        icnt1 = 0
        sum1  = 0.0
        sum2  = 0.0
c
        if (nn .eq. 1) kyear = nt1
        if (nn .eq. 2) kyear = nt2
        if (nn .eq. 3) kyear = nt3
        if (nn .eq. 4) kyear = nt4
c
        do 322 k = 1, kyear
          do 321 j = 1, ihout
            do 320 i = 1, iyout
c
              if (nn .eq. 1) term =  bas1(i, j, k)
              if (nn .eq. 2) term =  pro1(i, j, k)
              if (nn .eq. 3) term =  gss1(i, j, k)
              if (nn .eq. 4) term =  dtt1(i, j, k)
              if (nn .eq. 1) aa   = awts1(i, j, k)
              if (nn .eq. 2) aa   = awts2(i, j, k)
              if (nn .eq. 3) aa   = awts3(i, j, k)
              if (nn .eq. 4) aa   = awts4(i, j, k)
c
              if (k .eq. 1) sum2 = sum2 + aa
c
              if (term .ne. bogus) then
                sum1 = sum1 + term * aa
              else
                icnt1 = icnt1 + 1
              end if
  320       continue
  321     continue
  322   continue
c
        if (nn .eq. 1) print *, 'CHECK! Grand mean SET1: ', sum1 / kyear
        if (nn .eq. 2) print *, 'CHECK! Grand mean SET2: ', sum1 / kyear
        if (nn .eq. 3) print *, 'CHECK! Grand mean SET3: ', sum1 / kyear
        if (nn .eq. 4) print *, 'CHECK! Grand mean SET4: ', sum1 / kyear
c
  323 continue
c
c
      print *, ' ******************************************************'
      print *, ' END CHECKING SUM OF ANOMALIES!'
      print *, ' ******************************************************'
      print *, ' '
c
c-------------------------------------------------------------------------------
c     ** 3.3 **
c
c     For each atmospheric layer processed, compute spatial averages, and then 
c     subtract these from each latitude band. Do this only for non-missing data.
c     This is the imean=2 case. In the imean=3 case, calculate the overall 
c     global mean for the four stratospheric layers (ssu3, ssu2, ssu1, and tls)
c     and subtract this overall global mean from each latitude band of each of
c     the four stratospheric layers. Also calculate the overall global mean of
c     the tropospheric layers (ttt and tlt) and subtract this overall global 
c     mean from each latitude band of each of the two tropospheric layers.
c
c     Note: If mass weighting has been selected (the ivweight=4 case), do not
c     subtract any spatial averages. In this case, cmask consists of the product
c     of area weights and mass weights (see Section 2.8).
c-------------------------------------------------------------------------------
c
      if ((imean .ne. 1) .and. (ivweight .eq. 4)) then ! Case of mean removal.
        print *, 'User cannot select mass weighting and mean removal!'
        print *, 'Current selections: imean    = ', imean
        print *, 'Current selections: ivweight = ', ivweight
        stop
      end if
c
      if (imean .eq. 1) goto 336 ! Case of no spatial mean subtraction.
c
      print *, ' ******************************************************'
      print *, ' START SUBTRACTION OF SPATIAL MEANS!'
      print *, ' ******************************************************'
      print *, ' '
c
      do 335 nn = 1, 4 ! Loop over number of input data sets.
c
        if (nn .eq. 1) kyear = nt1
        if (nn .eq. 2) kyear = nt2
        if (nn .eq. 3) kyear = nt3
        if (nn .eq. 4) kyear = nt4
c
        sum4 = 0.0
        sum5 = 0.0
        sum6 = 0.0
        do 334 k = 1, kyear ! Loop over years.
          do 331 j = 1, ihout ! Loop over vertical layers.
            sum1 = 0.0
            sum2 = 0.0
            sums = 0.0
            sumt = 0.0
            do 330 i = 1, iyout ! Loop over latitude bands.
              if (nn .eq. 1) term1 = bas1(i, j, k)
              if (nn .eq. 2) term1 = pro1(i, j, k)
              if (nn .eq. 3) term1 = gss1(i, j, k)
              if (nn .eq. 4) term1 = dtt1(i, j, k)
              term2 = cmask(i, j)
c
              if (term1 .ne. bogus) then
                sum1  = sum1 + (term1 * term2)
                sum2  = sum2 + term2
              end if
  330       continue
            gmean(j, k) = sum1 / sum2 ! Spatial average for current layer/year.
            print *, 'SET', nn, ' gmean(j, k) = ', gmean(j, k), k
  331     continue
c
          sums = gmean(1, k) + gmean(2, k) + gmean(3, k) + gmean(4, k)
          sumu = gmean(1, k) + gmean(2, k) + gmean(3, k)
          sumt = gmean(5, k) + gmean(6, k)
c
          gms = sums / 4.0 ! Overall stratospheric global mean.
          gmu = sumu / 3.0 ! Overall SSU global mean.
          gmt = sumt / 2.0 ! Overall tropospheric global mean.
c
          print *, 'SET', nn, ' gms = ', gms, k
          print *, 'SET', nn, ' gmt = ', gmt, k
c
          do 333 j = 1, ihout ! Loop over vertical layers.
            sum3 = 0.0
            do 332 i = 1, iyout ! Loop over latitude bands.
              if (nn .eq. 1) term1 = bas1(i, j, k)
              if (nn .eq. 2) term1 = pro1(i, j, k)
              if (nn .eq. 3) term1 = gss1(i, j, k)
              if (nn .eq. 4) term1 = dtt1(i, j, k)
c
              if (imean .eq. 2) then ! Case of global means for each layer.
                if (term1 .ne. bogus) then
                  term3 = term1 - gmean(j, k)
                else
                  term3 = bogus
                end if
c
              else if (imean .eq. 3) then ! Case of STRAT and TROP global means.
                if (term1 .ne. bogus) then
                  if (ivweight .eq. 0) then ! SSU+MSU case.
                    if (j .ge. 1) term3 = term1 - gms ! ssu3.
                    if (j .ge. 2) term3 = term1 - gms ! ssu2.
                    if (j .ge. 3) term3 = term1 - gms ! ssu1.
                    if (j .ge. 4) term3 = term1 - gms ! tls.  
                    if (j .ge. 5) term3 = term1 - gmt ! ttt.  
                    if (j .ge. 6) term3 = term1 - gmt ! tlt.  
                  else if (ivweight .eq. 1) then ! TROP only case.
                    if (j .ge. 5) term3 = term1 - gmt ! ttt.  
                    if (j .ge. 6) term3 = term1 - gmt ! tlt.  
                  else if (ivweight .eq. 2) then ! SSU only case.
                    if (j .ge. 1) term3 = term1 - gmu ! ssu3.
                    if (j .ge. 2) term3 = term1 - gmu ! ssu2.
                    if (j .ge. 3) term3 = term1 - gmu ! ssu1.
                  else if (ivweight .eq. 3) then ! MSU case.
                    if (j .ge. 4) term3 = term1 - gms ! ttt.  
                    if (j .ge. 5) term3 = term1 - gmt ! ttt.  
                    if (j .ge. 6) term3 = term1 - gmt ! tlt.  
                  else if (ivweight .eq. 5) then ! MIN_olap case.
                    if (j .ge. 1) term3 = term1 - gms ! ssu3  
                    if (j .ge. 4) term3 = term1 - gms ! tls.  
                    if (j .ge. 6) term3 = term1 - gmt ! tlt.  
                  end if
                else
                  term3 = bogus
                end if
              end if
c
              if (nn .eq. 1) bas1(i, j, k) = term3                     
              if (nn .eq. 2) pro1(i, j, k) = term3                     
              if (nn .eq. 3) gss1(i, j, k) = term3                     
              if (nn .eq. 4) dtt1(i, j, k) = term3                     
              if (term1 .ne. bogus) then
                sum3 = sum3 + (term3 * cmask(i, j))
              end if
  332       continue
            print *, 'SET=', nn, ' YR=', k, ' LYR=', j, 'sum3=', sum3
            if (j. eq. 1) sum4 = sum4 + sum3
            if (j. eq. 2) sum5 = sum5 + sum3
            if (j. eq. 3) sum6 = sum6 + sum3
  333     continue
  334   continue
        print *, 'AFTER SPATIAL MEAN SUBTRACTION!'
        print *, 'SUM, SPACE-TIME ANOMALIES IN SET', nn, ' LYR1=', sum4
        print *, 'SUM, SPACE-TIME ANOMALIES IN SET', nn, ' LYR2=', sum5
        print *, 'SUM, SPACE-TIME ANOMALIES IN SET', nn, ' LYR3=', sum6
  335 continue
  336 continue
c
      print *, ' ******************************************************'
      print *, ' FINISH SUBTRACTION OF SPATIAL MEANS!'
      print *, ' ******************************************************'
      print *, ' '
c
c-------------------------------------------------------------------------------
c     ** 3.4 **
c
c     If selected, perform spatial filtering of anomaly data.
c
c     IFILTIN = 0          No filtering
c     IFILTIN = 1          13-term Gaussian filter
c     IFILTIN = 2          13-term moving average filter
c     IFILTIN = 3           5-term binomial filter        
c-------------------------------------------------------------------------------
c
      print *, ' ******************************************************'
      print *, ' BEGIN TEMPORAL FILTERING FOR IFILTIN .GE. 1!'
      print *, ' ******************************************************'
      print *, ' '
c
      do 345 nn = 1, 4
c
        if (nn .eq. 1) kyear = nt1
        if (nn .eq. 2) kyear = nt2
        if (nn .eq. 3) kyear = nt3
        if (nn .eq. 4) kyear = nt4
c
        do 344 j = 1, ihout
          do 343 i = 1, iyout
            do 340 k = 1, kyear
              if (nn .eq. 1) xx(k) = bas1(i, j, k)
              if (nn .eq. 2) xx(k) = pro1(i, j, k)
              if (nn .eq. 3) xx(k) = gss1(i, j, k)
              if (nn .eq. 4) xx(k) = dtt1(i, j, k)
  340       continue
c
            if (ifiltin .eq. 0) then
              do 341 k = 1, kyear
                zz(k) = xx(k)
  341         continue
            else if (ifiltin .eq. 1) then
               call gauss13_pad(n1, kyear, bogus, xx, zz, critwt)
            else if (ifiltin .eq. 2) then
               call linear13_pad(n1, kyear, bogus, xx, zz, critwt)
            else if (ifiltin .eq. 3) then
               call binom5_pad(n1, kyear, bogus, xx, zz, critwt)
            end if
c
            do 342 k = 1, kyear
              if (nn .eq. 1) bas1(i, j, k) = zz(k)
              if (nn .eq. 2) pro1(i, j, k) = zz(k)
              if (nn .eq. 3) gss1(i, j, k) = zz(k)
              if (nn .eq. 4) dtt1(i, j, k) = zz(k)
              if ((k .eq. kyear) .and. (nn .eq. 1)) then
                print *, 'bas1(i, j, k) = ', bas1(i, j, k), i, j
              else if ((k .eq. kyear) .and. (nn .eq. 2)) then
                print *, 'pro1(i, j, k) = ', pro1(i, j, k), i, j
              else if ((k .eq. kyear) .and. (nn .eq. 3)) then
                print *, 'gss1(i, j, k) = ', gss1(i, j, k), i, j
              else if ((k .eq. kyear) .and. (nn .eq. 4)) then
                print *, 'dtt1(i, j, k) = ', dtt1(i, j, k), i, j
              end if
  342       continue
  343     continue
  344   continue
  345 continue
c
      print *, ' ******************************************************'
      print *, ' FINISH TEMPORAL FILTERING FOR IFILTIN .GE. 1!'
      print *, ' ******************************************************'
      print *, ' '
c
c
c
c
c
c
c-------------------------------------------------------------------------------
c     ** 4.0 **
c
c     MATRIX TRANSPOSITION.
c
c     Reduce the spatial dimension of the anomaly data (from a 2-D field of 
c     zonally-averaged atmospheric temperature at individual pressure levels to
c     a 1-D spatial field; the latitude and pressure dimensions are concatenated
c     into a single dimension). The reduced spatial dimension of the field is 
c     nyh1. 
c
c     This section also weights the anomaly data. Weighting is based on the
c     user-selected value of ivweight. For ivweight=0/1/2/3/5, the weight used
c     is the square root of the area weight of the current latitude band. For 
c     ivweight=4 the weight is the square root of the product of the area weight
c     of the current latitude band and the mass weight of the atmospheric layer.
c
c     Note that values of the weights were calculated in Section 2.8 (before
c     taking the square root of each weight) and stored in the 2-D cmask array.
c     After taking the square root of each weight, the result is stored in the
c     1-D wmask array. In Sections 6.8, 6.9, and 7.9, wmask is subsequently
c     used for "unweighting" EOFs for visual display purposes. The wmask array 
c     is time-invariant, and is therefore set for a single year only.
c-------------------------------------------------------------------------------
c
      print *, ' ******************************************************'
      print *, ' START REDUCING SPATIAL DIMENSION!'            
      print *, ' ******************************************************'
      print *, ' '
c
      do 403 nn = 1, 4
c
        if (nn .eq. 1) kyear = nt1
        if (nn .eq. 2) kyear = nt2
        if (nn .eq. 3) kyear = nt3
        if (nn .eq. 4) kyear = nt4
c
        do 402 k = 1, kyear
          ijcnt = 0
          do 401 j = 1, ihout
            do 400 i = 1, iyout
              if (cmask(i, j) .ne. bogus) then
                ijcnt     = ijcnt + 1
                cmbweight = cmask(i, j)
c
                if (nn .eq. 1) then
                  term1        = bas1(i, j, k)
                  z1(k, ijcnt) = term1 * sqrt(cmbweight)
                  if (k .eq. 1) wmask(ijcnt) = cmbweight 
                else if (nn .eq. 2) then
                  term1        = pro1(i, j, k)
                  p1(k, ijcnt) = term1 * sqrt(cmbweight)
                else if (nn .eq. 3) then
                  term1        = gss1(i, j, k)
                  y1(k, ijcnt) = term1 * sqrt(cmbweight)
                else if (nn .eq. 4) then
                  term1        = dtt1(i, j, k)
                  d1(k, ijcnt) = term1 * sqrt(cmbweight)
                end if
c
                if ((k .eq. kyear) .and. (nn .eq. 1)) then
                  print *, 'z1(kyear, ijcnt) = ', z1(kyear, ijcnt), 
     &            ijcnt
                end if
              end if
c
  400       continue
  401     continue
  402   continue
  403 continue
c
      nyh1 = ijcnt
c
      print *, ' ******************************************************'
      print *, ' END REDUCING SPATIAL DIMENSION!'            
      print *, ' REDUCED DIMENSION NYH1 = ', nyh1
      print *, ' ******************************************************'
      print *, ' '
c
c-------------------------------------------------------------------------------
c     ** 4.1 **
c
c     Transpose SET1 and SET3 anomaly matrices.
c-------------------------------------------------------------------------------
c
      print *, ' ******************************************************'
      print *, ' START TRANSPOSING ANOMALY MATRICES!'            
      print *, ' ******************************************************'
      print *, ' '
c
      do 412 nn = 1, 2
        if (nn .eq. 1) kyear = nt1
        if (nn .eq. 2) kyear = nt3
        do 411 ix = 1, nyh1
          do 410 k = 1, kyear
            if (nn .eq. 1) zt(ix, k) = z1(k, ix)
            if (nn .eq. 2) yt(ix, k) = y1(k, ix)
  410     continue
  411   continue
  412 continue
c
      print *, ' ******************************************************'
      print *, ' END TRANSPOSING ANOMALY MATRICES!'            
      print *, ' ******************************************************'
      print *, ' '
c
c
c
c
c
c
c-------------------------------------------------------------------------------
c     ** 5.0 **
c
c     SCATTER MATRIX COMPUTATION.
c
c     Compute scatter matrix s1 for SET1 anomaly data (z1 and zt). First 
c     determine the setting for PCA analysis: state space or sample space. 
c     Definitions of sample and state space are from Preisendorfer, page 64.
c
c     In the terminology used here:
c     it = 1, nt1   (time index)
c     ix = 1, nyh1  (space index)
c
c          Noise data         Scatter matrix               Setting
c
c     z1(it,ix)  *  zt(ix,it)  =   t1(it,it)    Sample space  (nyh1 > nt1 - 1)
c     zt(ix,it)  *  z1(it,ix)  =   s1(ix,ix)    State space   (nyh1 < nt1 - 1)
c
c     In the terminology used by Preisendorfer:
c
c          Noise data         Scatter matrix               Setting
c
c      z(n,p)  *  zt(p,n)  =    t(n,n)           Sample space  (p > n - 1)
c     zt(p,n)  *   z(n,p)  =    s(p,p)           State space   (p < n - 1)
c
c     In dealing with the full spatial field the setting is generally sample
c     space. Note that the maximum (physical) dimension of the scatter matrix
c     (for both sample and state space) is nt1.
c-------------------------------------------------------------------------------
c
      print *, ' ******************************************************'
      print *, ' START COMPUTING SCATTER MATRICES!'              
      print *, ' ******************************************************'
      print *, ' '
c
      if (nyh1 .gt. nt1 - 1) then
        ispace1 = 1
        print *, ' Sample space. nyh1 > nt1 - 1 '
        print *, ' Input matrix is z1(it,ix).'
        call mmult1(n1, nyh, nyh, n1, n1, n1,
     &              nt1, nyh1, nt1, z1, zt, t1)
c
      else if (nyh1 .le. nt1 - 1) then
        ispace1 = 2
        print *, ' State space. nt1 - 1 > nyh1 '
        print *, ' Input matrix is zt(it,ix).'
        call mmult1(nyh, n1, n1, nyh, n1, n1,
     &              nyh1, nt1, nyh1, zt, z1, s1)
      end if
c
      print *, ' SET1 physical dimensions: n1, nyh   = ', n1, nyh
      print *, ' SET1 logical dimensions:  nt1, nyh1 = ', nt1, nyh1
c
c-------------------------------------------------------------------------------
c     ** 5.1 **
c
c     Now compute scatter matrix s2 for SET3 anomaly data (y1 and yt). 
c-------------------------------------------------------------------------------
c
      if (nyh1 .gt. nt3 - 1) then
        ispace2 = 1
        print *, ' Sample space. nyh1 > nt3 - 1 '
        print *, ' Input matrix is y1(it,ix).'
        call mmult1(n3, nyh, nyh, n3, n3, n3,
     &              nt3, nyh1, nt3, y1, yt, t2)
c
      else if (nyh1 .le. nt3 - 1) then
        ispace2 = 2
        print *, ' State space. nt3 - 1 > nyh1 '
        print *, ' Input matrix is yt(it,ix).'
        call mmult1(nyh, n3, n3, nyh, n3, n3,
     &              nyh1, nt3, nyh1, yt, y1, s2)
      end if

      print *, ' SET3 physical dimensions: n3, nyh   = ', n3, nyh
      print *, ' SET3 logical dimensions:  nt3, nyh1 = ', nt3, nyh1

      print *, ' ******************************************************'
      print *, ' FINISH COMPUTING SCATTER MATRICES!'              
      print *, ' ******************************************************'
      print *, ' '
c
c     do 511 it1 = 1, nt1
c       do 510 it2 = 1, nt1
c         t1(it1, it2) = t1(it1, it2) / nt1
c 510   continue
c 511 continue
c
c     do 513 it1 = 1, nt3
c       do 512 it2 = 1, nt3
c         t2(it1, it2) = t2(it1, it2) / nt3
c 512   continue
c 513 continue
c
c
c
c
c
c-------------------------------------------------------------------------------
c     ** 6.0 **
c
c     COMPUTATION OF SET1 EOFS.
c
c     Compute eigenvectors and eigenvalues. The solution is by Householder
c     reduction of a real, symmetric matrix to tridiagonal form (using TRED2)
c     followed by use of tridiagonal ql implicit algorithm (TQLI) to determine
c     eigenvectors and eigenvalues. Routines TRED2 and TQLI are from Press et
c     al., 1987 (see REFERENCES).
c
c     On output from TQLI, eigenvalues are stored in f (sample space) or 
c     e (state space) and eigenvalues in val1. Terminology follows Preisen-
c     dorfer, 1988. Note that both TRED2 and TQLI reference the array t for 
c     the sample space case and s for state space!
c
c     Use EIGSRT (also from Press et al., 1987) to sort the eigenvalues in 
c     ascending order and rearrange the eigenvectors accordingly. Use eigvar
c     to compute the explained variance (var1) and cumulative explained 
c     variance (cum1) for each eigenvalue.
c
c     First operate on SET1 data.
c-------------------------------------------------------------------------------
c
      print *, ' ******************************************************'
      print *, ' BEGIN COMPUTING EOFS!'              
      print *, ' ******************************************************'
      print *, ' '
c
      if (ispace1 .eq. 1) then
        call tred2(nt1, n1, t1, val1, work2a)
        call tqli(nt1, n1, val1, work2a, t1)
c
        do 601 it1 = 1, nt1   
          do 600 it2 = 1, nt1   
            f1(it1, it2) = t1(it1, it2)
  600     continue
  601   continue
c
        call eigsrt(val1, f1, nt1, n1)
        call eigvar(nt1, n1, val1, var1, cum1)
c
      else if (ispace1 .eq. 2) then
        call tred2(nyh1, n1, s1, val1, work2a)
        call tqli(nyh1, n1, val1, work2a, s1)
c
        do 604 ix1 = 1, nyh1
          do 603 ix2 = 1, nyh1
            e1(ix1, ix2) = s1(ix1, ix2)
  603     continue
  604   continue
c
        call eigsrt(val1, e1, nyh1, n1)
        call eigvar(nyh1, n1, val1, var1, cum1)
c
      end if                                        
c
c-------------------------------------------------------------------------------
c     ** 6.1 **
c
c     Compute EOFs for SET3 data.
c-------------------------------------------------------------------------------
c
      if (ispace2 .eq. 1) then
        call tred2(nt3, n3, t2, val2, work2b)
        call tqli(nt3, n3, val2, work2b, t2)
c
        do 611 it1 = 1, nt3  
          do 610 it2 = 1, nt3  
            f2(it1, it2) = t2(it1, it2)
  610     continue
  611   continue
c
        call eigsrt(val2, f2, nt3, n3)
        call eigvar(nt3, n3, val2, var2, cum2)
c
      else if (ispace2 .eq. 2) then
        call tred2(nyh1, n3, s2, val2, work2b)
        call tqli(nyh1, n3, val2, work2b, s2)
c
        do 614 ix1 = 1, nyh1
          do 613 ix2 = 1, nyh1
            e2(ix1, ix2) = s2(ix1, ix2)
  613     continue
  614   continue
c
        call eigsrt(val2, e2, nyh1, n3)
        call eigvar(nyh1, n3, val2, var2, cum2)
c
      end if                                               
c
      print *, ' ******************************************************'
      print *, ' FINISH COMPUTING EOFS!'              
      print *, ' ******************************************************'
      print *, ' '
c
c-------------------------------------------------------------------------------
c     ** 6.2 **
c
c     It is computationally more efficient to work from the scatter matrix
c     with the smallest dimension. If the setting is the sample space, we
c     have computed the eigenvectors f and eigenvalues m of scatter matrix t
c     (notation follows Preisendorfer, pp. 64-66).
c
c     Eigenvectors       Eigenvalues     Defining equation
c
c     f(it,it)           m(it)           t * f = f * m          Sample space
c     e(ix,ix)           l(ix)           s * e = e * l          State space
c
c     Number of non-zero eigenvalues l(ix) = rank z = min(nt1-1, nyh1)
c     Number of non-zero eigenvalues m(it) = rank z = min(nt1-1, nyh1)
c
c     The rank r of z1 is the number of linearly independent rows or columns of
c     z1. Here are some examples:
c
c       Time          Space        Rank (r)     Non-zero eigenvalues
c
c     nt1 = 100     nyh1 = 2048    rank = 99      non-zero m = 99     Sample
c     nt1 = 100     nyh1 =   32    rank = 32      non-zero l = 32     State
c
c     To obtain the eigenvectors:
c
c     phat(ix,j), ix = 1, nyh1; j = 1, irank1
c     ghat(ix,j), ix = 1, nyh1; j = 1, irank2
c
c     from the sample space eigenvectors f(it,it), transform back to the state
c     space setting using a simple matrix multiplication:
c
c     phat = zt * f1
c     ghat = yt * f2
c
c     For either case (state or sample space), it is useful to normalize the
c     computed eigenvectors. This is done automatically by tred2 and tqli. To
c     compensate for the effect of the back-transformation (sample space only),
c     use:
c
c     phat(ix,it) = phat(ix,it) / sqrt (m1(it))
c     ghat(ix,it) = ghat(ix,it) / sqrt (m2(it))
c
c     First process SET1 data.
c-------------------------------------------------------------------------------
c
      irank1 = min(nt1 - 1, nyh1)
      irank2 = min(nt3 - 1, nyh1)
      print *, ' Rank of SET1 matrix = ', irank1
      print *, ' Rank of SET3 matrix = ', irank2
c
      if (ispace1 .eq. 1) then
        print *, ' SET1. Back transformation to state space!'
        call mmult1(nyh, n1, n1, n1, nyh, n1,
     &              nyh1, nt1, nt1, zt, f1, phat)
c
        do 621 j = 1, irank1
          root1 = sqrt(val1(j))
          do 620 ix = 1, nyh1
            phat(ix, j) = phat(ix, j) / root1
  620     continue
  621   continue
c
      else if (ispace1 .eq. 2) then
c
        do 623 j = 1, irank1
          do 622 ix = 1, nyh1
            phat(ix, j) = e1(ix, j)
  622     continue
  623   continue
      end if
c
c-------------------------------------------------------------------------------
c     ** 6.3 **
c
c     Now process SET3 data.
c-------------------------------------------------------------------------------
c
      if (ispace2 .eq. 1) then
        print *, ' SET3. Back transformation to state space!'
        call mmult1(nyh, n3, n3, n3, nyh, n3,
     &              nyh1, nt3, nt3, yt, f2, ghat)
c
        do 631 j = 1, irank2
          root2 = sqrt(val2(j))
          do 630 ix = 1, nyh1
            ghat(ix, j) = ghat(ix, j) / root2
  630     continue
  631   continue
c
      else if (ispace2 .eq. 2) then
c
        do 633 j = 1, irank2
          do 632 ix = 1, nyh1
            ghat(ix, j) = e2(ix, j)
  632     continue
  633   continue
      end if
c
c-------------------------------------------------------------------------------
c     ** 6.4 **
c
c     The normalized eigenvectors [phat(ix,j)] are made unique by choosing, for
c     real-valued components of phat, the first non-zero component's sign as
c     positive (Preisendorfer, pp. 63-64). If the first component of the jth 
c     eigenvector is exactly zero (i.e., phat(1,j) = 0), move to component 
c     phat(2,j), etc., until a non-zero component is located. If the located 
c     component is greater than zero, process the next eigenvector. If the 
c     located component is less than zero, multiply all subsequent components 
c     of the eigenvector by -1.0 before proceeding to the next eigenvector.
c
c     Process SET1 first.
c-------------------------------------------------------------------------------
c
      do 642 j = 1, irank1
        do 641 ix1 = 1, nyh1
          if (phat(ix1, j) .ne. 0.0) then
            if (phat(ix1, j) .lt. 0.0) then
              do 640 ix2 = 1, nyh1
                 phat(ix2, j) = -1.0 * phat(ix2, j)
  640         continue
            end if
            goto 642
          end if
  641   continue
  642 continue
c
c-------------------------------------------------------------------------------
c     ** 6.5 **
c
c     As for (6.4), but for SET3 EOFs.
c-------------------------------------------------------------------------------
c
      do 652 j = 1, irank2
        do 651 ix1 = 1, nyh1
          if (ghat(ix1, j) .ne. 0.0) then
            if (ghat(ix1, j) .lt. 0.0) then
              do 650 ix2 = 1, nyh1
                 ghat(ix2, j) = -1.0 * ghat(ix2, j)
  650         continue
            end if
            goto 652
          end if
  651   continue
  652 continue
c
c-------------------------------------------------------------------------------
c     ** 6.6 **
c
c     Perform orthogonality check for computed eigenvectors. See Preisendorfer, 
c     pg. 34. Compute the integral over space for each EOF. Process SET1 first.
c-------------------------------------------------------------------------------
c
      print *, ' ******************************************************'
      print *, ' START CHECKING ORTHOGONALITY OF UNWEIGHTED EOFS!'      
      print *, ' ******************************************************'
      print *, ' '
c
      if (debug) then
c
        print *, ' Orthogonality check for unweighted NOISE EOFs!'
c
        do 662 j = 1, irank1
          sum1 = 0.0
          do 661 k = 1, irank1
            sum2 = 0.0
            do 660 ix = 1, nyh1
              if (k .eq. 1) sum1 = sum1 + phat(ix, j)
              sum2 = sum2 + phat(ix, j) * phat(ix, k)
  660       continue
            if (j .eq. k) print *, ' sum2 = ', sum2,
     &      '  j = ', j, '  k = ', k
  661     continue
c         print *, ' Integral over space (phat) = ', sum1,
c    &    '  EOF = ', j, k
  662   continue
      end if
c
c-------------------------------------------------------------------------------
c     ** 6.7 **
c
c     As for (6.6), but for EOFs of s(x,t).
c-------------------------------------------------------------------------------
c
      if (debug) then
c
        print *, ' Orthogonality check for unweighted EOFs of s(x,t)!'
c
        do 672 j = 1, irank2
          sum1 = 0.0
          do 671 k = 1, irank2
            sum2 = 0.0
            do 670 ix = 1, nyh1
              if (k .eq. 1) sum1 = sum1 + ghat(ix, j)
              sum2 = sum2 + ghat(ix, j) * ghat(ix, k)
  670       continue
            if (j .eq. k) print *, ' sum2 = ', sum2,
     &      '  j = ', j, '  k = ', k
  671     continue
c         print *, ' Integral over space (ghat) = ', sum1,
c    &    '  EOF = ', j, k
  672   continue
      end if
c
      print *, ' ******************************************************'
      print *, ' FINISH CHECKING ORTHOGONALITY OF UNWEIGHTED EOFS!'     
      print *, ' ******************************************************'
      print *, ' '
c
c-------------------------------------------------------------------------------
c     ** 6.8 **
c
c     "Unweight" the weighted SET1 EOFs (phat) by dividing by the square root
c     of wmask in order to generate the "true" EOFs (ptru). Note that ptru is  
c     used for display purposes only, not for computation of the PC time series.
c     After "unweighting", reconstruct the full spatial field, restoring 
c     missing data in the correct spatial locations. The reconstructed SET1
c     EOFs (ptru1) are output.
c-------------------------------------------------------------------------------
c              
      print *, ' ******************************************************'
      print *, ' START UNWEIGHTING EOFS!'            
      print *, ' ******************************************************'
      print *, ' '
c
      do 681 k = 1, irank1 
        do 680 ix = 1, nyh1  
          ptru(ix, k) = phat(ix, k) / sqrt(wmask(ix))
  680   continue
  681 continue
c
      do 684 k = 1, irank1 
        ix = 0
        do 683 j = 1, ihout
          do 682 i = 1, iyout
            if (cmask(i, j) .ne. bogus) then
              ix = ix + 1
              ptru1(i, j, k) = ptru(ix, k)
            else
              ptru1(i, j, k) = bogus
            end if
  682     continue
  683   continue
  684 continue
c
c-------------------------------------------------------------------------------
c     ** 6.9 **
c
c     As for (6.8), but for SET3 EOFs (ghat).
c-------------------------------------------------------------------------------
c
      do 691 k = 1, irank2
        do 690 ix = 1, nyh1 
          gtru(ix, k) = ghat(ix, k) / sqrt(wmask(ix))
  690   continue
  691 continue
c
      do 694 k = 1, irank2
        ix = 0
        do 693 j = 1, ihout
          do 692 i = 1, iyout
            if (cmask(i, j) .ne. bogus) then
              ix = ix + 1
              gtru1(i, j, k) = gtru(ix, k)
            else
              gtru1(i, j, k) = bogus
            end if
  692     continue
  693   continue
  694 continue
c
      print *, ' ******************************************************'
      print *, ' FINISH UNWEIGHTING EOFS!'            
      print *, ' ******************************************************'
      print *, ' '
c
c
c
c
c
c
c-------------------------------------------------------------------------------
c     ** 7.0 **
c
c     DEFINITION OF RAW FINGERPRINT AND COMPUTATION OF ROTATED FINGERPRINT.
c
c     First define the raw fingerprint (fraw). The fingerprint is the expected
c     signal direction. The user can define the raw fingerprint in the following 
c     ways:
c
c     NFINGER = 1          Fingerprint is EOF#1 of s(x,t).    
c
c     NFINGER > 1 < 10     Fingerprint is weighted linear combination of first
c                          NFINGER EOFs of s(x,t). Weights are the explained
c                          variances of the s(x,t) EOFs. 
c
c     NFINGER > or = 10    Fingerprint is undefined. Program aborts (the
c                          assumption here is that the user will never require
c                          more than 9 s(x,t) EOFs for defining the 
c                          fingerprint).
c-------------------------------------------------------------------------------
c
      print *, ' ******************************************************'
      print *, ' BEGIN COMPUTING RAW FINGERPRINT!'
      print *, ' ******************************************************'
      print *, ' '
c
      if ((nfinger .gt. 1) .and. (nfinger .lt. 10)) then
c
        bnorm = 0.0
        sum2  = 0.0
        do 701 ix = 1, nyh1
          sum1 = 0.0
          do 700 mm = 1, nfinger
            sum1 = sum1 + (ghat(ix, mm) * var2(mm))
  700     continue
          fraw(ix) = sum1 
          bnorm    = bnorm + (fraw(ix) ** 2)
  701   continue
c
c-------------------------------------------------------------------------------
c     ** 7.1 **
c  
c     Renormalize raw fingerprint after forming weighted linear combination.
c-------------------------------------------------------------------------------
c
        do 710 ix = 1, nyh1
          fraw(ix) = fraw(ix) / sqrt(bnorm)
          sum2     = sum2 + fraw(ix) ** 2
  710   continue
        print *, 'NORM = ', sum2
c
c-------------------------------------------------------------------------------
c     ** 7.2 **
c  
c     NFINGER = 1 and NFINGER = 0 cases.
c-------------------------------------------------------------------------------
c
      else if (nfinger .eq. 1) then
c
        do 720 ix = 1, nyh1
          fraw(ix) = ghat(ix, 1)
  720   continue
c
      else if (nfinger .eq. 0) then
c
        do 721 ix = 1, nyh1
          fraw(ix) = ghat(ix, 2)
  721   continue
c
      else  
c
        print *, 'Error in RESIDUAL_FINGER1!'
        print *, 'Number of EOFs for fingerprint cannot exceed 9!'
        print *, 'Selected NFINGER = ', nfinger
        stop
      end if
c
      print *, ' ******************************************************'
      print *, ' END COMPUTING RAW FINGERPRINT!'
      print *, ' ******************************************************'
      print *, ' '
c
c-------------------------------------------------------------------------------
c     ** 7.3 **
c
c     First project the raw fingerprint onto the control EOFs. This yields the 
c     projection coefficients (alpha). Also store the projection vector (fpro).
c     The residual fingerprint is then computed and stored in fres.
c-------------------------------------------------------------------------------
c
      print *, ' ******************************************************'
      print *, ' BEGIN COMPUTING FINGERPRINT PROJECTION!'
      print *, ' ******************************************************'
      print *, ' '
c
      term2 = 0.0
c
      do 731 m = 1, itrunc
        term1 = 0.0
c
        do 730 ix = 1, nyh1
          term1  = term1 + fraw(ix) * phat(ix, m)
  730   continue
c
        alpha(m) = term1
        term2    = term2 + (alpha(m) ** 2)
        print *, 'ALPHA = ', alpha(m), m
c
  731 continue
c
      print *, 'ALPHA**2 = ', term2
c
      cnorm = 0.0
      dnorm = 0.0
c
      do 733 ix = 1, nyh1
        fpro(ix) = 0.0
c
        do 732 m = 1, itrunc
          fpro(ix) = fpro(ix) + (alpha(m) * phat(ix, m))
  732   continue
        cnorm = cnorm + fpro(ix) ** 2
  733 continue
c
      do 734 ix = 1, nyh1
        fres(ix) = fraw(ix) - fpro(ix)
        dnorm    = dnorm + fres(ix) ** 2
  734 continue
c
      print *, 'CNORM = ', cnorm
      print *, 'DNORM = ', dnorm
c
c-------------------------------------------------------------------------------
c     ** 7.4 **
c
c     Normalize the residual fingerprint (fres).
c-------------------------------------------------------------------------------
c
      termd = sqrt(dnorm)
      do 740 ix = 1, nyh1
        fres(ix) = fres(ix) / termd
  740 continue
c
      print *, ' ******************************************************'
      print *, ' END COMPUTING FINGERPRINT PROJECTION!'
      print *, ' ******************************************************'
      print *, ' '
c
c-------------------------------------------------------------------------------
c     ** 7.5 **
c
c     Rotate the projection of the fingerprint. Note that the scaling of the
c     projection coefficients (alpha) is by the inverse of the percentage
c     explained variance of the individual control run eigenvalues.
c-------------------------------------------------------------------------------
c
      print *, ' ******************************************************'
      print *, ' BEGIN ROTATING FINGERPRINT!'
      print *, ' ******************************************************'
      print *, ' '
c
      sum1   = 0.0
      anorm1 = 0.0
c
      do 751 ix = 1, nyh1
        frot(ix) = 0.0
c
        do 750 m = 1, itrunc
          frot(ix) = frot(ix) + (alpha(m) * phat(ix, m) / var1(m))
  750   continue
        anorm1 = anorm1 + frot(ix) ** 2
  751 continue
c
c-------------------------------------------------------------------------------
c     ** 7.6 **
c
c     Normalize the rotated fingerprint (frot).
c-------------------------------------------------------------------------------
c
      terma = sqrt(anorm1)
      do 760 ix = 1, nyh1
        frot(ix) = frot(ix) / terma
  760 continue

      print *, 'Squared length of rotated fingerprint component ',
     &'(ANORM1) = ', anorm1
      print *, 'sqrt{ANORM1} = ', sqrt(anorm1)
c
      print *, ' ******************************************************'
      print *, ' FINISH COMPUTING ROTATED FINGERPRINT!'
      print *, ' ******************************************************'
      print *, ' '
c
c-------------------------------------------------------------------------------
c     ** 7.7 **
c
c     Formation of the optimized fingerprint. 
c
c     CASE                 RESULT
c     icasewt = 0          No residual included         
c     icasewt = 1          Inclusion of non-optimized residual        
c     icasewt = 2          Inclusion of residual with average optimization
c
c     Note: For icasewt = 2, alphaw is used to weight the residual fingerprint 
c     by the inverse of the cumulative variance of the control run (SET1) EOFs 
c     not used in the rotation of the projected fingerprint (i.e., EOFs 
c     itrunc + 1 to z, the rank of the SET1 covariance matrix).
c-------------------------------------------------------------------------------
c
      fact1 = (var1(1) / val1(1)) ** 2
      anorm2 = anorm1 * (fact1)
c
      enorm = 0.0
      if (icasewt .eq. 0) then
        alphaw = 0.0
        print *, 'For no residual included case (icasewt = 0) ',
     &  'alphaw = ', alphaw
      else if (icasewt .eq. 1) then
        alphaw = sqrt(anorm1 / cnorm)
        print *, 'For inclusion of non-optimized residual ',
     &  '(icasewt = 1) alphaw = ', alphaw 
      else if (icasewt .eq. 2) then
        alphaw = 1.0 / (1.0 - (cum1(itrunc) / 100.0))
        print *, 'For inclusion of residual with average ',
     &  'optimization (icasewt = 2)  alphaw = ', alphaw
        print *, 'This is the cumulative variance of SET1 EOFs ',
     &  'not used in rotation.'
      end if
c
      beta = (alphaw * termd) / (terma + (alphaw * termd))
c
      do 770 ix = 1, nyh1
        greof(ix) = ((1.0 - beta) * frot(ix)) + (beta * fres(ix))
        enorm = enorm + greof(ix) ** 2
  770 continue
c
c-------------------------------------------------------------------------------
c     ** 7.8 **
c
c     Normalize the optimized fingerprint.
c-------------------------------------------------------------------------------
c
      terme = sqrt(enorm)
      do 780 ix = 1, nyh1
        greof(ix) = greof(ix) / terme
  780 continue
c
      print *, 'ENORM = ', enorm
c
c-------------------------------------------------------------------------------
c     ** 7.9 **
c
c     "Unweight" the optimized fingerprint (greof) the residual fingerprint
c     (fres), and the projected fingerprint (fpro) by dividing by the square
c     root of wmask. This is for display purposes only, not for computation of 
c     the PC time series. After "unweighting", reconstruct the full spatial 
c     field, restoring missing data in the correct spatial locations.
c-------------------------------------------------------------------------------
c
      ix = 0
      do 791 j = 1, ihout
        do 790 i = 1, iyout
          ix          = ix + 1
          grtru(ix)   = greof(ix) / sqrt(wmask(ix))
          frestru(ix) =  fres(ix) / sqrt(wmask(ix))
          fprotru(ix) =  fpro(ix) / sqrt(wmask(ix))
  790   continue
  791 continue
c
      ix = 0
      do 793 j = 1, ihout
        do 792 i = 1, iyout
          if (cmask(i, j) .ne. bogus) then
            ix = ix + 1
            grtru1(i, j)   =   grtru(ix)
            frestru1(i, j) = frestru(ix)
            fprotru1(i, j) = fprotru(ix)
          else
            grtru1(i, j)   = bogus
            frestru1(i, j) = bogus         
            fprotru1(i, j) = bogus         
          end if
  792   continue
  793 continue
c
c
c
c
c
c
c-------------------------------------------------------------------------------
c     ** 8.0 **
c
c     COMPUTATION OF PRINCIPAL COMPONENT TIME SERIES.                     
c
c     First project the SET1, SET2, SET3, and SET4 anomaly data onto the first
c     irank2 EOFs of the s(x,t) data.
c
c     DATA         Anomalies    PCs          EOFs
c
c     SET1         z1           bpc1         ghat
c     SET2         p1           ppc1         ghat
c     SET3         y1           gpc1         ghat
c     SET4         d1           dpc1         ghat
c-------------------------------------------------------------------------------
c
      print *, ' ********************************************* '
      print *, ' BEGIN PROJECTION ONTO EOFS of s(x,t) DATA!'
      print *, ' ********************************************* '
      print *, ' '
c
      do 803 nn = 1, 4
        if (nn .eq. 1) kyear = nt1
        if (nn .eq. 2) kyear = nt2
        if (nn .eq. 3) kyear = nt3
        if (nn .eq. 4) kyear = nt4
c
        do 802 m = 1, irank2
          do 801 it = 1, kyear
            term1 = 0.0
            do 800 ix = 1, nyh1
              if (nn .eq. 1) term1 = term1 + (ghat(ix, m) * z1(it, ix))
              if (nn .eq. 2) term1 = term1 + (ghat(ix, m) * p1(it, ix))
              if (nn .eq. 3) term1 = term1 + (ghat(ix, m) * y1(it, ix))
              if (nn .eq. 4) term1 = term1 + (ghat(ix, m) * d1(it, ix))
  800       continue
            if (nn .eq. 1) bpc1(it, m) = term1
            if (nn .eq. 2) ppc1(it, m) = term1
            if (nn .eq. 3) gpc1(it, m) = term1
            if (nn .eq. 4) dpc1(it, m) = term1
  801     continue
  802   continue
  803 continue
c
      print *, ' ********************************************* '
      print *, ' FINISH PROJECTION ONTO EOFS OF s(x,t) DATA!'
      print *, ' ********************************************* '
      print *, ' '
c
c-------------------------------------------------------------------------------
c     ** 8.1 **
c
c     Now project the SET1, SET2, SET3, and SET4 anomaly data onto the first
c     irank1 control run EOFs.
c
c     DATA         Anomalies    PCs          EOFs
c
c     SET1         z1           bpc2         phat
c     SET2         p1           ppc2         phat
c     SET3         y1           gpc2         phat
c     SET4         d1           dpc2         phat
c-------------------------------------------------------------------------------
c
      print *, ' ********************************************* '
      print *, ' BEGIN PROJECTION ONTO CONTROL EOFS!'
      print *, ' ********************************************* '
      print *, ' '
c
      do 813 nn = 1, 4
        if (nn .eq. 1) kyear = nt1
        if (nn .eq. 2) kyear = nt2
        if (nn .eq. 3) kyear = nt3
        if (nn .eq. 4) kyear = nt4
c
        do 812 m = 1, irank1
          do 811 it = 1, kyear
            term1 = 0.0
            do 810 ix = 1, nyh1
              if (nn .eq. 1) term1 = term1 + (phat(ix, m) * z1(it, ix))
              if (nn .eq. 2) term1 = term1 + (phat(ix, m) * p1(it, ix))
              if (nn .eq. 3) term1 = term1 + (phat(ix, m) * y1(it, ix))
              if (nn .eq. 4) term1 = term1 + (phat(ix, m) * d1(it, ix))
  810       continue
            if (nn .eq. 1) bpc2(it, m) = term1
            if (nn .eq. 2) ppc2(it, m) = term1
            if (nn .eq. 3) gpc2(it, m) = term1
            if (nn .eq. 4) dpc2(it, m) = term1
  811     continue
  812   continue
  813 continue
c
      print *, ' ********************************************* '
      print *, ' FINISH PROJECTION ONTO CONTROL EOFS!'
      print *, ' ********************************************* '
      print *, ' '
c
c-------------------------------------------------------------------------------
c     ** 8.2 **
c
c     Determine how well the SET1/2/3/4 anomaly data can be represented in the 
c     space of the control EOFs. First compute the total space-time variance
c     of the (area-weighted) anomaly data.
c
c     Variable      Description
c     var1raw       Sum of squares of dataset used in rotation
c     var2raw       Sum of squares of dataset used for projection onto guess
c     var3raw       Sum of squares of dataset used for defining guess
c     var4raw       Sum of squares of dataset used for computing statistics
c-------------------------------------------------------------------------------
c
      do 822 nn = 1, 4
        if (nn .eq. 1) kyear = nt1
        if (nn .eq. 2) kyear = nt2
        if (nn .eq. 3) kyear = nt3
        if (nn .eq. 4) kyear = nt4
c
        term1   = 0.0
c
        do 821 ix = 1, nyh1
          do 820 it = 1, kyear
            if (nn .eq. 1) term1 = term1 + (z1(it, ix) ** 2)
            if (nn .eq. 2) term1 = term1 + (p1(it, ix) ** 2)
            if (nn .eq. 3) term1 = term1 + (y1(it, ix) ** 2)
            if (nn .eq. 4) term1 = term1 + (d1(it, ix) ** 2)
  820     continue
  821   continue
c
        if (nn .eq. 1) var1raw = term1 / kyear
        if (nn .eq. 2) var2raw = term1 / kyear
        if (nn .eq. 3) var3raw = term1 / kyear
        if (nn .eq. 4) var4raw = term1 / kyear
c
  822 continue
c
c-------------------------------------------------------------------------------
c     ** 8.3 **
c
c     Now compute the space-time variance of the SET1/2/3/4 PC time series. Do
c     this both individually and cumulatively.
c
c     Variable      Description
c     var1con       Sum of squares for projection of SET1 onto control EOFs
c     var2con       Sum of squares for projection of SET2 onto control EOFs
c     var3con       Sum of squares for projection of SET3 onto control EOFs
c     var4con       Sum of squares for projection of SET4 onto control EOFs
c
c     var1consum    As for var1b, but cumulative version                   
c     var2consum    As for var2b, but cumulative version                   
c     var3consum    As for var3b, but cumulative version                   
c     var4consum    As for var4b, but cumulative version                   
c-------------------------------------------------------------------------------
c
      do 832 nn = 1, 4
        if (nn .eq. 1) kyear = nt1
        if (nn .eq. 2) kyear = nt2
        if (nn .eq. 3) kyear = nt3
        if (nn .eq. 4) kyear = nt4
c
        term2 = 0.0
        do 831 m = 1, irank1
          term1 = 0.0
          do 830 it = 1, kyear
            if (nn .eq. 1) term1 = term1 + (bpc2(it, m) ** 2)
            if (nn .eq. 2) term1 = term1 + (ppc2(it, m) ** 2)
            if (nn .eq. 3) term1 = term1 + (gpc2(it, m) ** 2)
            if (nn .eq. 4) term1 = term1 + (dpc2(it, m) ** 2)
  830     continue
c
          if (nn .eq. 1) var1con(m)    = term1 / kyear 
          if (nn .eq. 2) var2con(m)    = term1 / kyear        
          if (nn .eq. 3) var3con(m)    = term1 / kyear         
          if (nn .eq. 4) var4con(m)    = term1 / kyear         
c
          term2 = term2 + (term1 / kyear)
c
          if (nn .eq. 1) var1consum(m) = term2         
          if (nn .eq. 2) var2consum(m) = term2        
          if (nn .eq. 3) var3consum(m) = term2         
          if (nn .eq. 4) var4consum(m) = term2         
c
  831   continue
  832 continue
c                                                                                     
c-------------------------------------------------------------------------------
c     ** 8.4 **
c
c     Now project the SET1, SET2, SET3, and SET4 anomaly data onto the rotated
c     fingerprint.
c
c     DATA         Anomalies    PCs          EOFs
c
c     SET1         z1           bpc3         greof
c     SET2         p1           ppc3         greof
c     SET3         y1           gpc3         greof
c     SET4         d1           dpc3         greof
c-------------------------------------------------------------------------------
c
      print *, ' ********************************************* '
      print *, ' BEGIN PROJECTION ONTO ROTATED FINGERPRINT!'
      print *, ' ********************************************* '
      print *, ' '
c
      do 842 nn = 1, 4
        if (nn .eq. 1) kyear = nt1
        if (nn .eq. 2) kyear = nt2
        if (nn .eq. 3) kyear = nt3
        if (nn .eq. 4) kyear = nt4
c
        do 841 it = 1, kyear
          term1 = 0.0
          do 840 ix = 1, nyh1
            if (nn .eq. 1) term1 = term1 + (greof(ix) * z1(it, ix))
            if (nn .eq. 2) term1 = term1 + (greof(ix) * p1(it, ix))
            if (nn .eq. 3) term1 = term1 + (greof(ix) * y1(it, ix))
            if (nn .eq. 4) term1 = term1 + (greof(ix) * d1(it, ix))
  840     continue
          if (nn .eq. 1) bpc3(it) = term1
          if (nn .eq. 2) ppc3(it) = term1
          if (nn .eq. 3) gpc3(it) = term1
          if (nn .eq. 4) dpc3(it) = term1
  841   continue
  842 continue
c
      print *, ' ********************************************* '
      print *, ' FINISH PROJECTION ONTO ROTATED FINGERPRINT!'
      print *, ' ********************************************* '
      print *, ' '
c
c-------------------------------------------------------------------------------
c     ** 8.5 **
c
c     Now project the SET1, SET2, SET3, and SET4 anomaly data onto the raw     
c     fingerprint.
c
c     DATA         Anomalies    PCs          EOFs
c
c     SET1         z1           bpc4         fraw 
c     SET2         p1           ppc4         fraw 
c     SET3         y1           gpc4         fraw 
c     SET4         d1           dpc4         fraw 
c-------------------------------------------------------------------------------
c
      print *, ' ********************************************* '
      print *, ' BEGIN PROJECTION ONTO RAW FINGERPRINT!'
      print *, ' ********************************************* '
      print *, ' '
c
      do 852 nn = 1, 4
        if (nn .eq. 1) kyear = nt1
        if (nn .eq. 2) kyear = nt2
        if (nn .eq. 3) kyear = nt3
        if (nn .eq. 4) kyear = nt4
c
        do 851 it = 1, kyear
          term1 = 0.0
          do 850 ix = 1, nyh1
            if (nn .eq. 1) term1 = term1 + (fraw(ix) * z1(it, ix))
            if (nn .eq. 2) term1 = term1 + (fraw(ix) * p1(it, ix))
            if (nn .eq. 3) term1 = term1 + (fraw(ix) * y1(it, ix))
            if (nn .eq. 4) term1 = term1 + (fraw(ix) * d1(it, ix))
  850     continue
          if (nn .eq. 1) bpc4(it) = term1
          if (nn .eq. 2) ppc4(it) = term1
          if (nn .eq. 3) gpc4(it) = term1
          if (nn .eq. 4) dpc4(it) = term1
  851   continue
  852 continue
c
      print *, ' ********************************************* '
      print *, ' FINISH PROJECTION ONTO RAW FINGERPRINT!'
      print *, ' ********************************************* '
      print *, ' '
c
c-------------------------------------------------------------------------------
c     ** 8.6 **
c
c     Compute the total variance of the projections onto the rotated and        
c     raw fingerprints.
c
c     Variable      Description
c     var1rot       Sum of squares for projection of SET1 on rotated fingerprint
c     var2rot       Sum of squares for projection of SET2 on rotated fingerprint
c     var3rot       Sum of squares for projection of SET3 on rotated fingerprint
c     var4rot       Sum of squares for projection of SET4 on rotated fingerprint
c
c     var1norot     Sum of squares for projection of SET1 on raw fingerprint
c     var2norot     Sum of squares for projection of SET2 on raw fingerprint
c     var3norot     Sum of squares for projection of SET3 on raw fingerprint
c     var4norot     Sum of squares for projection of SET4 on raw fingerprint
c-------------------------------------------------------------------------------
c
      do 862 ll = 1, 2
        do 861 nn = 1, 4
          if (nn .eq. 1) kyear = nt1
          if (nn .eq. 2) kyear = nt2
          if (nn .eq. 3) kyear = nt3
          if (nn .eq. 4) kyear = nt4
c
          term1 = 0.0
          do 860 it = 1, kyear
            if (ll .eq. 1) then
              if (nn .eq. 1) term1 = term1 + (bpc3(it) ** 2)
              if (nn .eq. 2) term1 = term1 + (ppc3(it) ** 2)
              if (nn .eq. 3) term1 = term1 + (gpc3(it) ** 2)
              if (nn .eq. 4) term1 = term1 + (dpc3(it) ** 2)
            else if (ll .eq. 2) then
              if (nn .eq. 1) term1 = term1 + (bpc4(it) ** 2)
              if (nn .eq. 2) term1 = term1 + (ppc4(it) ** 2)
              if (nn .eq. 3) term1 = term1 + (gpc4(it) ** 2)
              if (nn .eq. 4) term1 = term1 + (dpc4(it) ** 2)
            end if
  860     continue
c
          if (ll .eq. 1) then
            if (nn .eq. 1) var1rot   = term1 / kyear
            if (nn .eq. 2) var2rot   = term1 / kyear
            if (nn .eq. 3) var3rot   = term1 / kyear
            if (nn .eq. 4) var4rot   = term1 / kyear
          else if (ll .eq. 2) then
            if (nn .eq. 1) var1norot = term1 / kyear
            if (nn .eq. 2) var2norot = term1 / kyear
            if (nn .eq. 3) var3norot = term1 / kyear
            if (nn .eq. 4) var4norot = term1 / kyear
          end if
  861   continue
  862 continue
c
c-------------------------------------------------------------------------------
c     ** 8.7 **
c
c     Diagnostics for PCs. If the SET1 anomalies were time-centered, then their 
c     projections onto SET1 EOFs should:
c
c     1.  Have the same dimensions as the original SET1 anomaly data.
c     2.  Sum to zero over the time index it. 
c
c     As a check, compute the sum over the time index for each PC time series.
c     Time series have the same dimensions as the entries of z1 (Preisendorfer, 
c     1988, pg. 32).
c-------------------------------------------------------------------------------
c
      if (debug) then
        print *, ' ********************************************* '
        print *, ' START CHECKING THAT SET1 PCS SUM TO ZERO!'
        print *, ' ********************************************* '
        print *, ' '
c
        do 872 m2 = 1, npcout
          do 871 m1 = 1, npcout
            sum1 = 0.0
            sum2 = 0.0
            do 870 it = 1, nt1
              if (m2. eq. 1) sum1 = sum1 + bpc2(it, m1)
              sum2 = sum2 + (bpc2(it, m1) * bpc2(it, m2))
  870       continue
            avv1 = sum1 / float(nt1)
c
            if (m2 .eq. 1) then
              print *, ' Mean of bpc2(it, m1) = ', avv1, ' EOF # ', m1
            end if
            print *, ' sum2 = ', sum2, ' m1, m2 ', m1, m2
  871     continue
  872   continue
c
        print *, ' ********************************************* '
        print *, ' FINISH CHECKING THAT SET1 PCS SUM TO ZERO!'
        print *, ' ********************************************* '
        print *, ' '
c
      end if
c
c-------------------------------------------------------------------------------
c     ** 8.8 **
c
c     As for (8.7) but for SET3 PCs.
c-------------------------------------------------------------------------------
c
      if (debug) then
        print *, ' ********************************************* '
        print *, ' START CHECKING THAT SET3 PCS SUM TO ZERO!'
        print *, ' ********************************************* '
        print *, ' '
c
        do 882 m2 = 1, npcout
          do 881 m1 = 1, npcout
            sum1 = 0.0
            sum2 = 0.0
            do 880 it = 1, nt3
              if (m2. eq. 1) sum1 = sum1 + gpc1(it, m1)
              sum2 = sum2 + (gpc1(it, m1) * gpc1(it, m2))
  880       continue
            avv1 = sum1 / float(nt3)
c
            if (m2 .eq. 1) then
              print *, ' Mean of gpc1(it, m1) = ', avv1, ' EOF # ', m1
            end if
            print *, ' sum2 = ', sum2, ' m1, m2 ', m1, m2
  881     continue
  882   continue
c
        print *, ' ********************************************* '
        print *, ' FINISH CHECKING THAT SET3 PCS SUM TO ZERO!'
        print *, ' ********************************************* '
        print *, ' '
c
      end if
c
c
c
c
c
c
c-------------------------------------------------------------------------------
c     ** 9.0 **
c
c     EXPLAINED VARIANCE DIAGNOSTICS.               
c
c     Variable        Description
c     DIAGNOSTIC1     This gives the fraction of the total space-time variance 
c                     of the signal dataset (SET2) that can be explained by 
c                     projection of the SET2 anomaly data onto the EOFs of SET1
c                     (the control run used for optimization).
c
c     DIAGNOSTIC2     This gives the fraction of the total space-time variance 
c                     of the "guess" dataset (SET3) that can be explained 
c                     by projection of the SET3 anomaly data onto the EOFs of 
c                     SET1 (the control run used for optimization).
c
c     DIAGNOSTIC3     This gives the fraction of the total space-time variance 
c                     of the dataset used for computing natural variability
c                     statistics (SET4) that can be explained by projection of 
c                     the SET4 anomaly data onto the EOFs of SET1 (the control 
c                     run used for optimization).
c
c     DIAGNOSTIC4     This gives the fraction of the total space-time variance 
c                     of the dataset used for optimization (SET1) that can be 
c                     explained by projection of the SET1 anomaly data onto the 
c                     EOFs of SET1. This is a simple check; the diagnostic 
c                     should give the fractional variance explained by each SET1
c                     EOF, so that DIAGNOSTIC3 = var1.
c
c     DIAGNOSTIC5     This is a ratio of DIAGNOSTIC1 and DIAGNOSTIC3.
c
c     DIAGNOSTIC6     This is a ratio of DIAGNOSTIC2 and DIAGNOSTIC3.
c-------------------------------------------------------------------------------
c
      ss1 = 0.0
      ss2 = 0.0
      ss3 = 0.0
      ss4 = 0.0
      ss5 = 0.0
      ss6 = 0.0
      do 900 m = 1, itrunc
        diagnostic1(m) = (var2con(m) / var2raw) * 100.0
        diagnostic2(m) = (var3con(m) / var3raw) * 100.0
        diagnostic3(m) = (var4con(m) / var4raw) * 100.0
        diagnostic4(m) = (var1con(m) / var1raw) * 100.0
        diagnostic5(m) = diagnostic1(m) / diagnostic3(m)
        diagnostic6(m) = diagnostic2(m) / diagnostic3(m)
        ss1 = ss1 + diagnostic1(m)
        ss2 = ss2 + diagnostic2(m)
        ss3 = ss3 + diagnostic3(m)
        ss4 = ss4 + diagnostic4(m)
        ss5 = ss5 + diagnostic5(m)
        ss6 = ss6 + diagnostic6(m)
  900 continue
c
c-------------------------------------------------------------------------------
c     ** 9.1 **
c
c     Variable        Description
c     DIAGNOSTIC7     This gives the fraction of the total space-time variance 
c                     of the signal dataset (SET2) that can be explained by 
c                     projection of the SET2 anomaly data onto the unrotated
c                     EOF1 of SET3 (the first guess EOF).
c
c     DIAGNOSTIC8     This gives the fraction of the total space-time variance 
c                     of the "guess" dataset (SET3) that can be explained 
c                     by projection of the SET3 anomaly data onto the unrotated
c                     EOF1 of SET3 (the first guess EOF).
c
c     DIAGNOSTIC9     This gives the fraction of the total space-time variance 
c                     of the dataset used for computing natural variability
c                     statistics (SET4) that can be explained by projection of 
c                     the SET4 anomaly data onto the unrotated EOF1 of SET3    
c                     (the first guess EOF).         
c
c     DIAGNOSTIC10    This gives the fraction of the total space-time variance 
c                     of the dataset used for optimization (SET1) that can be 
c                     explained by projection of the SET1 anomaly data onto the 
c                     unrotated EOF1 of SET3 (the first guess EOF).         
c
c     DIAGNOSTIC11    Ratio of DIAGNOSTIC7 / DIAGNOSTIC9
c
c     DIAGNOSTIC12    Ratio of DIAGNOSTIC8 / DIAGNOSTIC9
c-------------------------------------------------------------------------------
c
      diagnostic7  = (var2norot / var2raw) * 100.0
      diagnostic8  = (var3norot / var3raw) * 100.0
      diagnostic9  = (var4norot / var4raw) * 100.0
      diagnostic10 = (var1norot / var1raw) * 100.0
      diagnostic11 = diagnostic7 / diagnostic9
      diagnostic12 = diagnostic8 / diagnostic9
c
c-------------------------------------------------------------------------------
c     ** 9.4 **
c
c     Variable        Description
c     DIAGNOSTIC13    This gives the fraction of the total space-time variance
c                     of the signal dataset (SET2) that can be explained by
c                     projection of the SET2 anomaly data onto the rotated
c                     guess.
c
c     DIAGNOSTIC14    This gives the fraction of the total space-time variance 
c                     of the "guess" dataset (SET3) that can be explained 
c                     by projection of the SET3 anomaly data onto the rotated
c                     guess.
c
c     DIAGNOSTIC15    This gives the fraction of the total space-time variance
c                     of the dataset used for computing natural variability
c                     statistics (SET4) that can be explained by projection of
c                     the SET4 anomaly data onto the rotated guess.      
c
c     DIAGNOSTIC16    This gives the fraction of the total space-time variance 
c                     of the dataset used for optimization (SET1) that can be 
c                     explained by projection of the SET1 anomaly data onto the 
c                     rotated guess.         
c
c     DIAGNOSTIC17    Ratio of DIAGNOSTIC13 / DIAGNOSTIC15
c
c     DIAGNOSTIC18    Ratio of DIAGNOSTIC14 / DIAGNOSTIC15
c-------------------------------------------------------------------------------
c
      diagnostic13 = (var2rot / var2raw) * 100.0
      diagnostic14 = (var3rot / var3raw) * 100.0
      diagnostic15 = (var4rot / var4raw) * 100.0
      diagnostic16 = (var1rot / var1raw) * 100.0
      diagnostic17 = diagnostic13 / diagnostic15
      diagnostic18 = diagnostic14 / diagnostic15
c
c
c
c
c
c
c-------------------------------------------------------------------------------
c     ** 10.0 **
c
c     COMPUTE SIGNAL-TO-NOISE STATISTICS.
c
c     First perform spatial pattern correlations: control EOFs versus guess 
c     EOF.
c
c     rr1         Pattern correlation (spatial means subtracted)
c     rr2         Pattern correlation (spatial means included)  
c-------------------------------------------------------------------------------
c
      do 1002 m2 = 1, npcout
        do 1001 m1 = 1, npcout
          do 1000 ix = 1, nyh1   
            eofout1(ix) = ptru(ix, m1)
            eofout2(ix) = gtru(ix, m2)
 1000     continue
c
          call wtrt1(eofout1, eofout2, wmask, nyh, nyh1, rr1)
          call wtct1(eofout1, eofout2, wmask, nyh, nyh1, rr2)
c
          corr1(m1, m2) = rr1
          corr2(m1, m2) = rr2
          if ((m1 .eq. 5) .and. (m2 .eq. 1)) then
            print *, 'NEOF5 vs GEOF1. rr1 = ', rr1
          end if
 1001   continue
 1002 continue
c
c-------------------------------------------------------------------------------
c     ** 10.1 **
c
c     Determine the number of control and guess EOFs required to explain .ge.
c     95% of the total space-time variance of the data.
c-------------------------------------------------------------------------------
c
      do 1010 j = 1, irank1 
        if (cum1(j) .ge. 95.0) then
          icrit1 = j
          goto 1011
        end if
 1010 continue
 1011 continue
c
      do 1012 j = 1, irank2 
        if (cum2(j) .ge. 95.0) then
          icrit2 = j
          goto 1013
        end if
 1012 continue
 1013 continue
c
c
c
c
c
c
c-------------------------------------------------------------------------------
c     ** 11.0  **
c
c     COMPUTE SAMPLING DISTRIBUTION OF LINEAR TRENDS.
c
c     Compute sampling distribution of linear trends from projections of 
c     control run data onto unrotated fingerprint [fraw(ix)] and onto 
c     rotated guess EOF [greof(ix)]. Then compute standard error of the
c     sampling distribution.
c
c     DATA         EOFs         PCs          Standard error of linear trend
c
c     SET1         fraw         bpc4         bserr1
c     SET1         greof        bpc3         bserr2
c
c     Note: If the initial analysis period is 10 years, then set lenchunk = 9. 
c     If the initial analysis period is five years, set lenchunk = 4.
c-------------------------------------------------------------------------------
c
      print *, ' ********************************************* '
      print *, ' START COMPUTING TREND STATISTICS!'         
      print *, ' ********************************************* '
      print *, ' '
c
      icount   =   0
      lenchunk =   4
cBDS  lenchunk =   9
      ilast    =  nt1 / 2
c
      do 1100 it = 1, nt1     
        xx(it) = bpc4(it)
        yy(it) = bpc3(it)
 1100 continue
c
 2000 continue
c
      icount   = icount + 1
      lenchunk = lenchunk + 1
c
      call err1(n1, nt1, lenchunk, iover, iuserlap, ifact,
     &bogus, ntrends, xx, work2a, work2b, adist, bdist, stderr1)
c
      call err1(n1, nt1, lenchunk, iover, iuserlap, ifact,
     &bogus, ntrends, xx, work2a, work2b, adist, cdist, stderr1)
c
c BDS
c
      if (lenchunk .eq. 39) then ! Only for 39-year trends.
        open(unit=2, file='test.d', form='formatted', status='new')
        write(2, 1103) ntrends, lenchunk
        do 1101 ii = 1, ntrends
          write(2, 1102) ii, lenchunk * 1.0, cdist(ii)
 1101   continue
        close(unit=2)
      end if
c
      call err1(n1, nt1, lenchunk, iover, iuserlap, ifact,
     &bogus, ntrends, yy, work2a, work2b, adist, bdist, stderr2)
c
      bserr1(icount) = stderr1
      bserr2(icount) = stderr2
      isize1(icount) = lenchunk
      inumb1(icount) = ntrends 
c
      if (lenchunk .lt. ilast) then
        go to 2000
      else
        go to 2001
      end if
c
 2001 continue
c
 1102 format(i7, 2x, 2f10.4)
 1103 format(2i7)
c
      ntr(1) = icount
c
c-------------------------------------------------------------------------------
c     ** 11.1  **
c
c     Compute sampling distribution of linear trends from projections of 
c     "statistics" data onto unrotated fingerprint [fraw(ix)] and onto 
c     rotated guess EOF [greof(ix)]. Then compute standard error of the
c     sampling distribution.
c
c     DATA         EOFs         PCs          Trend distribution
c
c     SET4         fraw         dpc4         dserr1
c     SET4         greof        dpc3         dserr2
c
c     Note: If the initial analysis period is 10 years, then set lenchunk = 9. 
c     If the initial analysis period is five years, set lenchunk = 4.
c-------------------------------------------------------------------------------
c
      icount   =   0
      lenchunk =   4
cBDS  lenchunk =   9
c
      ilast    =  nt4 / 2
c
      do 1110 it = 1, nt4     
        xx(it) = dpc4(it)
        yy(it) = dpc3(it)
 1110 continue
c
 3000 continue
c
      icount   = icount + 1
      lenchunk = lenchunk + 1
c
      call err1(n1, nt4, lenchunk, iover, iuserlap, ifact,
     &bogus, ntrends, xx, work2a, work2b, adist, bdist, stderr1)
c
      call err1(n1, nt4, lenchunk, iover, iuserlap, ifact,
     &bogus, ntrends, yy, work2a, work2b, adist, bdist, stderr2)
c
      dserr1(icount) = stderr1
      dserr2(icount) = stderr2
      isize4(icount) = lenchunk
      inumb2(icount) = ntrends 
c                                               
      if (lenchunk .lt. ilast) then
        go to 3000
      else
        go to 3001
      end if
c
 3001 continue
c
      ntr(4) = icount
c
c-------------------------------------------------------------------------------
c     ** 11.2  **
c
c     Compute linear trends from projections of SET2 data onto unrotated 
c     fingerprint [fraw(ix)] and onto rotated guess EOF [greof(ix)].
c
c     DATA         EOFs         PCs          Trends
c
c     SET2         fraw         ppc4         ptrend1
c     SET2         greof        ppc3         ptrend2
c
c     Note: If the initial analysis period is 10 years, then set lenchunk = 9. 
c     If the initial analysis period is five years, set lenchunk = 4.
c-------------------------------------------------------------------------------
c
      icount   =    0
      lenchunk =    4   
cBDS  lenchunk =    9   
c
      idetect1 = istart ! Start date for calculation of signal trends (1986)
      idetect2 = idetect1 - int(rtime2(1)) + 1 ! rtime2(1) = start year of OBS
      ilast    = iy4 - idetect1 + 1  ! iy4 = last year of OBS
      print *, 'idetect1 = ', idetect1
      print *, 'idetect2 = ', idetect2
      print *, 'ilast    = ', ilast
c
      icc = 0
      do 1120 it = idetect2, nt2     
        icc     = icc + 1
        tt(icc) = it * 1.0
        xx(icc) = ppc4(it)
        yy(icc) = ppc3(it)
 1120 continue
c
 4000 continue
c
      icount   = icount + 1
      lenchunk = lenchunk + 1
c
      call reg2(n1, lenchunk, tt, xx, aa1, bb1)
      call reg2(n1, lenchunk, tt, yy, aa2, bb2)
c
      ptrend1(icount) = bb1 * 10.0
      ptrend2(icount) = bb2 * 10.0
      isize2(icount)  = lenchunk
c
      if (lenchunk .lt. ilast) then
        go to 4000
      else
        go to 4001
      end if
c
 4001 continue
c
      ntr(2) = icount                 
c
c-------------------------------------------------------------------------------
c     ** 11.3  **
c
c     Compute linear trends from projections of SET3 data onto unrotated
c     fingerprint [fraw(ix)] and onto rotated guess EOF [greof(ix)].
c
c     DATA         EOFs         PCs          Trends
c
c     SET3         fraw         gpc4         gtrend1
c     SET3         greof        gpc3         gtrend2
c
c     NOTE: The multi-model average change in synthetic MSU temperature (which
c     was calculated from the regridded hist_rcp85 experiments) typically starts 
c     in 1860 (not in 1850, because the hadgem2_cc hist_rcp85 simulation starts 
c     in December 1859). This means that rtime3(1), the first retrieved year
c     of the SET3 file containing the multi-model average change in synthetic
c     MSU temperature, is typically 1860.0, so idetect2 = 120 (1979 - 1860 + 1). 
c     If the last year of the retrieved SET3 data (iy6) is 2011, then the total
c     number of signal trends processed (ilast) =  33 (2011 - 1979 + 1).
c-------------------------------------------------------------------------------
c
      icount   =    0
      lenchunk =    4
cBDS  lenchunk =    9
c
      idetect1 = istart ! Start date for calculation of signal trends (1979)
      idetect2 = idetect1 - int(rtime3(1)) + 1 ! rtime3(1) = start year of GSS
      ilast    = iy6 - idetect1 + 1 ! iy6 = last year of GSS
      print *, 'idetect2 = ', idetect2
      print *, 'ilast    = ', ilast
c
      icc = 0
      do 1130 it = idetect2, nt3      
        icc = icc + 1
        tt(icc) = it * 1.0
        xx(icc) = gpc4(it)
        yy(icc) = gpc3(it)
 1130 continue
c
 5000 continue
c
      icount   = icount + 1
      lenchunk = lenchunk + 1
c
      call reg2(n1, lenchunk, tt, xx, aa1, bb1)
      call reg2(n1, lenchunk, tt, yy, aa2, bb2)
c
      gtrend1(icount) = bb1 * 10.0
      gtrend2(icount) = bb2 * 10.0
      isize3(icount)  = lenchunk
c
      if (lenchunk .lt. ilast) then
        go to 5000
      else
        go to 5001
      end if
c
 5001 continue
c
      ntr(3) = icount
c                                           
c-------------------------------------------------------------------------------
c     ** 11.4  **
c
c     Compute S/N ratios for SET2 data. The number of chunk lengths processed
c     is ntr(2).
c-------------------------------------------------------------------------------
c                                           
      isn2 = ntr(2)
c
      do 1140 kk = 1, isn2
        psn1(kk) = ptrend1(kk) / bserr1(kk)
        psn2(kk) = ptrend2(kk) / bserr2(kk)
        psn3(kk) = ptrend1(kk) / dserr1(kk)
        psn4(kk) = ptrend2(kk) / dserr2(kk)
        print *, 'SECTION 11.4! kk, psn2(kk) = ', kk, psn2(kk)
 1140 continue
c
c-------------------------------------------------------------------------------
c     ** 11.5  **
c
c     Compute S/N ratios for SET3 data. The number of chunk lengths processed
c     is ntr(3).
c-------------------------------------------------------------------------------
c
      isn3 = ntr(3)
c
      do 1150 kk = 1, isn3
        gsn1(kk) = gtrend1(kk) / bserr1(kk)
        gsn2(kk) = gtrend2(kk) / bserr2(kk)
        gsn3(kk) = gtrend1(kk) / dserr1(kk)
        gsn4(kk) = gtrend2(kk) / dserr2(kk)
 1150 continue
c                                                  
      print *, ' ********************************************* '
      print *, ' FINISH COMPUTING TREND STATISTICS!'         
      print *, ' ********************************************* '
      print *, ' '
c
c-------------------------------------------------------------------------------
c     ** 11.6  **
c                                                  
c     Compute detection times for SET2 data.
c
c     NOTE: All significance tests performed in Section 11.6 are one-tailed.
c     We assume that the SET3 data used for calculation of the fingerprint
c     provide "directional information" regarding the sign of the expected
c     signal trends. For example, in cases where spatial means are included in
c     the detection analysis, the sign of the anthropogenically-forced
c     temperature-change signal is typically negative for T4 (stratospheric
c     cooling) and positive for T2 (tropospheric warming). The sign of the
c     expected T4 and T2 signals is not obvious a priori in the case where
c     spatial means are removed. To apply one-tailed tests in a consistent way
c     to different variables (T4 and T2) and different cases (mean removed and
c     mean included), we use the final chunk length of the signal trends stored
c     in gtrend1 and gtrend2 to obtain information on the sign of the expected
c     signal.
c-------------------------------------------------------------------------------
c                                           
      print *, ' ********************************************* '
      print *, ' START COMPUTING DETECTION TIMES!'         
      print *, ' ********************************************* '
      print *, ' '
c
      if (icritsn .eq. 0) then
        critsn = 1.641 ! 5% significance threshold, one-tailed test.          
      else if (icritsn .eq. 1) then
        critsn = 2.326 ! 1% significance threshold, one-tailed test.          
      else if (icritsn .eq. 2) then
        critsn = 2.0   ! 2-sigma significance threshold.          
      else if (icritsn .eq. 3) then
        critsn = 3.0   ! 3-sigma significance threshold.          
      else if (icritsn .eq. 5) then
        critsn = 5.0   ! 5-sigma significance threshold.          
      else
        print *, 'Invalid value of icitsn!'
        print *, 'You selected icitsn = ', icritsn
        stop
      end if
c
      print *, 'HERE!'           
      print *, 'isn2 = ', isn2
      print *, 'isn3 = ', isn3
      print *, 'gsn1(isn3) = ', gsn1(isn3)
      print *, 'gsn2(isn3) = ', gsn2(isn3)
      print *, 'gsn3(isn3) = ', gsn3(isn3)
      print *, 'gsn4(isn3) = ', gsn4(isn3)
c
      do 1163 mm = 1, 4
        itimedet = 0
        do 1160 kk = 1, isn2 ! Loop over number of observed trends processed.
          if (mm .eq. 1) term1 = psn1(kk)
          if (mm .eq. 1) term2 = gsn1(isn3) ! isn3 is the number of guess trends.
c
          if (mm .eq. 2) term1 = psn2(kk)
          if (mm .eq. 2) term2 = gsn2(isn3)
c
          if (mm .eq. 3) term1 = psn3(kk)
          if (mm .eq. 3) term2 = gsn3(isn3)
c
          if (mm .eq. 4) term1 = psn4(kk)
          if (mm .eq. 4) term2 = gsn4(isn3)
c
          if ((term1 .ge. 0.0) .and. (term2 .ge. 0.0)) then
            isign = 1
          else if ((term1 .lt. 0.0) .and. (term2 .lt. 0.0)) then
            isign = 2
          else
            isign = 0
          end if
c
          if (mm .eq. 1) then
            print *, 'isign, kk, psn1(kk), gsn1(isn3) = ', isign, kk,
     &      psn1(kk), gsn1(isn3)
          else if (mm .eq. 2) then
c           print *, 'isign, kk, psn2(kk) = ', isign, kk, psn2(kk)
          else if (mm .eq. 3) then
c           print *, 'isign, kk, psn3(kk) = ', isign, kk, psn3(kk)
          else if (mm .eq. 4) then
c           print *, 'isign, kk, psn4(kk) = ', isign, kk, psn4(kk)
          end if
c
          if (isign .eq. 1) then
            if ((term1 .gt. critsn) .and. (itimedet .eq. 0)) then
cBDS          itimedet = kk + idetect1 + 8 ! idetect1 is typically 1986.
              itimedet = kk + idetect1 + 3 ! idetect1 is typically 1986.
              kindex1  = kk
            end if
c
            if (term1 .le. critsn) then
              itimedet = 0
            end if
c
          else if (isign .eq. 2) then
            term3 = abs(term1)
            if ((term3 .gt. critsn) .and. (itimedet .eq. 0)) then
cBDS          itimedet = kk + idetect1 + 8
              itimedet = kk + idetect1 + 3
              kindex1  = kk
            end if
c
            if (term3 .le. critsn) then
              itimedet = 0
            end if
c
          else if (isign .eq. 0) then
            itimedet = 0
          end if
 1160   continue
c
        linterp  = 1000
        rlinterp = 1000.0
c
        if (kindex1 .eq. 1) then ! Case of detection time for first trend.
          eeinc = itimedet * 1.0
c
        else if ((kindex1 .gt. 1) .and. (itimedet .ne. 0)) then ! Case of non-zero detection time.
          if (mm .eq. 1) then
            print *, 'kindex1 = ', kindex1
            dtime1 = psn1(kindex1-1)
            dtime2 = psn1(kindex1)
          else if (mm .eq. 2) then
            dtime1 = psn2(kindex1-1)
            dtime2 = psn2(kindex1)
          else if (mm .eq. 3) then
            dtime1 = psn3(kindex1-1)
            dtime2 = psn3(kindex1)
          else if (mm .eq. 4) then
            dtime1 = psn4(kindex1-1)
            dtime2 = psn4(kindex1)
          end if
c
          sninc1 = (dtime2 - dtime1) / rlinterp
          tninc1 = 1.0 / rlinterp
c
          do 1161 ii = 1, linterp - 1
            ddinc = abs(dtime1 + (ii * sninc1))
            if (ddinc .gt. critsn) then
              eeinc = (itimedet - 1.0) + (ii * tninc1)
              goto 1162
            else
              eeinc = itimedet * 1.0
            end if
 1161     continue
c
 1162     continue
c
cBDS    else  ! Case of no detection time.
        else if (kindex1 .eq. isn2) then ! Case of no detection time.
          eeinc = 0.0
        end if ! End processing non-zero detection times.
c
        print *, 'dtime1 = ', dtime1
        print *, 'dtime2 = ', dtime2
        print *, 'sninc1 = ', sninc1
        print *, 'tninc1 = ', tninc1
        print *, 'ddinc  = ', ddinc
        print *, 'eeinc  = ', eeinc
        print *, 'critsn = ', critsn
c
        if (mm .eq. 1) then
          ipsn1 = itimedet
          rpsn1 = eeinc
        else if (mm .eq. 2) then
          ipsn2 = itimedet
          rpsn2 = eeinc
        else if (mm .eq. 3) then
          ipsn3 = itimedet
          rpsn3 = eeinc
        else if (mm .eq. 4) then
          ipsn4 = itimedet
          rpsn4 = eeinc
        end if
c
 1163 continue
c                                           
c-------------------------------------------------------------------------------
c     ** 11.7  **
c
c     Compute detection times for SET3 data.
c-------------------------------------------------------------------------------
c
      do 1171 mm = 1, 4
        itimedet = 0
        do 1170 kk = 1, isn3
          if (mm .eq. 1) term1 = abs(gsn1(kk))
          if (mm .eq. 2) term1 = abs(gsn2(kk))
          if (mm .eq. 3) term1 = abs(gsn3(kk))
          if (mm .eq. 4) term1 = abs(gsn4(kk))
          if ((term1 .gt. critsn) .and. (itimedet .eq. 0)) then
cBDS        itimedet = kk + idetect1 + 8
            itimedet = kk + idetect1 + 3
          end if
c
          if (term1 .le. critsn) then
            itimedet = 0
          end if
 1170   continue
c
        if (mm .eq. 1) igsn1 = itimedet
        if (mm .eq. 2) igsn2 = itimedet
        if (mm .eq. 3) igsn3 = itimedet
        if (mm .eq. 4) igsn4 = itimedet
 1171 continue
c
      print *, ' ********************************************* '
      print *, ' FINISH COMPUTING DETECTION TIMES!'         
      print *, ' ********************************************* '
      print *, ' '
c
c
c
c                                           
c
c
c-------------------------------------------------------------------------------
c     ** 15.0 **
c
c     DATA OUTPUT: EOFS AND PCS.
c
c     Write the file containing the first npcout control EOFs.
c-------------------------------------------------------------------------------
c
      write(title1, '(1a)') 'Empirical Orthogonal Functions'
      write(units1, '(1a)') 'Dimensionless'
c
      print *, 'title1 = ', title1
      print *, 'units1 = ', units1
      print *, 'drsdic1 = ', drsdic1
      print *, 'drsdat1 = ', drsdat1
      print *, 'CHECK POINT 1A'
      ierr = aslun(7, drsdic1, 8, drsdat1, IDRS_CREATE)
      ierr = setname(' ', 'latitude', ' ', 'Degrees', ' ' )
      ierr = putvdim(7, iyout, rlat, idum1, idum2)
      ierr = setname(' ', 'pressure', ' ', 'hPa', ' ' )
      ierr = putvdim(7, ihout, rlev1, idum3, idum4)
      print *, 'CHECK POINT 1B'
c
c-------------------------------------------------------------------------------
c     ** 15.1 **
c
c     Loop over the number of EOFs output. Control and guess EOFs are output
c     to the same file.
c-------------------------------------------------------------------------------
c
      do 1513 nn = 1, 2
c
        if (nn .eq. 1) then
          string2 = cmod1
          string3 = cexp1
          ky1     = iy1
          ky2     = iy2
          kt      = nt1
          vv      = var1(nn)
        else
          string2 = cmod3
          string3 = cexp3
          ky1     = iy5
          ky2     = iy6
          kt      = nt3
          vv      = var2(m)
        end if
c
        iss2    = index(string2, ' ') - 1
        iss3    = index(string3, ' ') - 1
c
        print *, 'string2 = ', string2
        print *, 'string3 = ', string3
        print *, 'CHECK POINT 2'
c
        do 1512 k = 1, npcout 
c
          ix = 0
          if (nn .eq. 1) vv = var1(k)
          if (nn .eq. 2) vv = var2(k)
c
          do 1511 j = 1, ihout
            do 1510 i = 1, iyout
              ix = ix + 1
              if (nn .eq. 1) eofout1(ix) = ptru1(i, j, k)
              if (nn .eq. 2) eofout1(ix) = gtru1(i, j, k)
 1510       continue
 1511     continue
          print *, 'CHECK POINT 3'
c
c-------------------------------------------------------------------------------
c     ** 15.2 **
c
c     Encode various text strings for output to the DRS file.
c-------------------------------------------------------------------------------
c
          write(source1(k), '(3a, i2, 3a, f5.2, 1a, i4, 1a)')
     &    'MODEL: ', string2(1:iss2), '  EOF: ', k, '  EXPT: ',
     &    string3(1:iss3), '  EV: ', vv, '%  NT: ', kt, ' YEARS'
c
          print *, 'source1(k) = ', source1(k)
c
          if (nn .eq. 1) then
            if (k .lt. 10) write(varnam2, '(1a, i1)') 'ceof', k
            if (k .ge. 10) write(varnam2, '(1a, i2)') 'ceof', k
          else if (nn .eq. 2) then
            if (k .lt. 10) write(varnam2, '(1a, i1)') 'geof', k
            if (k .ge. 10) write(varnam2, '(1a, i2)') 'geof', k
          end if
c
          print *, 'varnam2 = ', varnam2
c
          rnum = k * 1.0
          ierr = cluvdb()
          ierr = setname(source1(k), varnam2, title1, units1, 'R*4')
          ierr = setvdim(1, ' ', 'latitude', ' ', ' ', rlat(1),
     &    rlat(iyout))
          ierr = setvdim(2, ' ', 'pressure', ' ', ' ', rlev1(1),
     &    rlev1(ihout))
c
          ierr = putdat(7, eofout1)
 1512   continue
 1513 continue
      print *, 'CHECK POINT 4'
c
c-------------------------------------------------------------------------------
c     ** 15.3 **
c
c     Write the file containing the rotated guess. This is appended to the file
c     containing the control and unrotated guess EOFs.
c-------------------------------------------------------------------------------
c
      write(source2, '(5a, i3, 1a, f6.4)')
     &'MODEL: ', string2(1:iss2), '  EXPT: ',
     &string3(1:iss3), ' OPT. F: ITRUNC = ', itrunc,
     &' BETA = ', beta
c
      print *, 'source2 = ', source2
c
      write(varnam2, '(1a)') 'greof1'
c
      ierr = cluvdb()
      ierr = setname(source2, varnam2, title1, units1, 'R*4')
      ierr = setvdim(1, ' ', 'latitude', ' ', ' ', rlat(1),
     &rlat(iyout))
      ierr = setvdim(2, ' ', 'pressure', ' ', ' ', rlev1(1),
     &rlev1(ihout))
c
      ierr = putdat(7, grtru1)
      print *, 'CHECK POINT 5'
c
c-------------------------------------------------------------------------------
c     ** 15.4 **
c
c     Write the file containing the residual fingerprint. This is appended to 
c     the file containing the control and unrotated guess EOFs, optimized 
c     fingerprint, etc.
c-------------------------------------------------------------------------------
c
      write(source2, '(3a, f6.4, 3a, i3)')
     &'MODEL: ', string2(1:iss2), ' TERMD = ', termd, '  EXPT: ',
     &string3(1:iss3), ' RESID. F: TRUNCATION = ', itrunc
c
      print *, 'source2 = ', source2
c
      write(varnam2, '(1a)') 'frestru'
c
      ierr = cluvdb()
      ierr = setname(source2, varnam2, title1, units1, 'R*4')
      ierr = setvdim(1, ' ', 'latitude', ' ', ' ', rlat(1),
     &rlat(iyout))
      ierr = setvdim(2, ' ', 'pressure', ' ', ' ', rlev1(1),
     &rlev1(ihout))
c
      ierr = putdat(7, frestru1)
c
      ierr = cllun(7)
      print *, 'CHECK POINT 6'
c
c
c
c
c
c
c-------------------------------------------------------------------------------
c     ** 16.0 **
c
c     OUTPUT S/N RESULTS AND EIGENVALUES.
c-------------------------------------------------------------------------------
c
      print *, ' Writing output to unit 2!'
      open(unit=2, file=outfil1, form='formatted', status='new')
c
      write(2, 9000)
      write(2, 9001) varnam1(1:8)
      write(2, 9010) iy1, iy2, iref1, iref2
      write(2, 9011) iy3, iy4, iref3, iref4
      write(2, 9012) iy5, iy6, iref5, iref6
      write(2, 9013) iy7, iy8, iref7, iref8
      write(2, 9020) infil1
      write(2, 9021) infil2
      write(2, 9022) infil3
      write(2, 9023) infil4
      write(2, 9030) nyh, nyh1, n1, nt1
      write(2, 9031) nyh, nyh1, n2, nt2
      write(2, 9032) nyh, nyh1, n3, nt3
      write(2, 9033) nyh, nyh1, n4, nt4
      if (ispace1 .eq. 1) write(2, 9040) irank1
      if (ispace1 .eq. 2) write(2, 9041) irank1
      if (ispace2 .eq. 1) write(2, 9050) irank2
      if (ispace2 .eq. 2) write(2, 9051) irank2
      write(2, 9060) itrunc
      if (imean .eq. 1) write(2, 9070)
      if (imean .eq. 2) write(2, 9071)
      if (imean .eq. 3) write(2, 9072)
      if ((yy1 .lt. 0.0) .and. (yy2 . gt. 0.0)) then
        write(2, 9073) yy1, yy2
      else if ((yy1 .ge. 0.0) .and. (yy2 .ge. 0.0)) then
        write(2, 9074) yy1, yy2
      else if ((yy1 .le. 0.0) .and. (yy2 .le. 0.0)) then
        write(2, 9075) yy1, yy2
      end if
      write(2, 9080) yweights1
      write(2, 9081) yweights2
      write(2, 9082) yweights3
      write(2, 9083) yweights4
c
      if (nfinger .eq. 1) write(2, 9084)
      if (nfinger .eq. 0) write(2, 9085)
      if (nfinger .gt. 1) write(2, 9086) nfinger
c
      if (icasewt .eq. 0) write(2, 9600)
      if (icasewt .eq. 1) write(2, 9601)
      if (icasewt .eq. 2) write(2, 9602)
c
      if (ifiltin .eq. 0) write(2,9605)
      if (ifiltin .eq. 1) write(2,9606)
      if (ifiltin .eq. 2) write(2,9607)
      if (ifiltin .eq. 3) write(2,9608)
c
      write(2, 9700)
      write(2, 9701)
      write(2, 9702)
      write(2, 9703)
      write(2, 9704)
      write(2, 9705)
c
      write(2, 9710) cmod1(1:ib1), cexp1(1:ic1)
      write(2, 9711) cmod4(1:ib4), cexp4(1:ic4)
      write(2, 9712) cmod3(1:ib3), cexp3(1:ic3)
      write(2, 9713) cmod2(1:ib2), cexp2(1:ic2)
c 
      write(2, 9750) cnorm
      write(2, 9751) dnorm
      write(2, 9752) anorm1
      write(2, 9753) anorm2
      write(2, 9754) enorm
      write(2, 9755) sqrt(anorm1 / cnorm)
      write(2, 9756) 1.0 / (1.0 - (cum1(itrunc) / 100.0))
      write(2, 9757) beta
c
      write(2, 9720)
      write(2, 9721)
      write(2, 9722)
      write(2, 9723)
      write(2, 9724)
      write(2, 9725)
c
      write(2, 9610)
      do 1600 m = 1, itrunc
        write(2, 9615) m, diagnostic1(m), diagnostic2(m),
     &  diagnostic3(m), diagnostic4(m), diagnostic5(m),
     &  diagnostic6(m)
 1600 continue
      write(2, 9616) ss1, ss2, ss3, ss4, ss5, ss6
      write(2, 9999)
c
      write(2, 9730)
      write(2, 9731)
      write(2, 9732)
      write(2, 9733)
      write(2, 9734)
      write(2, 9735)
c
      write(2, 9611)
      write(2, 9616) diagnostic7, diagnostic8,     
     &diagnostic9, diagnostic10, diagnostic11, diagnostic12
      write(2, 9999)
c
      write(2, 9740)
      write(2, 9741)
      write(2, 9742)
      write(2, 9743)
      write(2, 9744)
      write(2, 9745)
c
      write(2, 9612)
      write(2, 9616) diagnostic13, diagnostic14,     
     &diagnostic15, diagnostic16, diagnostic17, diagnostic18
      write(2, 9999)
c
c-------------------------------------------------------------------------------
c     ** 16.1 **
c
c     First write critical number of control and guess EOFs required to explain
c     .ge. 95% of the total space-time variance. Then write pattern correlations 
c     (centered) between control and guess PCs.
c-------------------------------------------------------------------------------
c
      write(2, 9090) icrit1
      write(2, 9091) icrit2
      write(2, 9092)
      write(2, 9094)
      do 1610 m1 = 1, npcout
        write(2, 9095) m1, (corr1(m1, m2), m2 = 1, npcout)
 1610 continue
c
c-------------------------------------------------------------------------------
c     ** 16.2 **
c
c     Write pattern correlations (uncentered) between control and guess PCs.
c-------------------------------------------------------------------------------
c
      write(2, 9093)
      write(2, 9094)
      do 1620 m1 = 1, npcout
        write(2, 9095) m1, (corr2(m1, m2), m2 = 1, npcout)
 1620 continue
c
c-------------------------------------------------------------------------------
c     ** 16.3 **
c
c     Write S/N results for SET2 data.
c-------------------------------------------------------------------------------
c
      sta = 'SIGNAL'                                   
      stb = 'NOISE'                                   
      stc = ' projected onto EOF1 of '         
      std = ' projected onto rotated EOF1 of ' 
c
      st1 = sta//'1: '//cs2//cexp2(1:ic2)//stc//cs3//cexp3(1:ic3)
      st2 = sta//'2: '//cs2//cexp2(1:ic2)//std//cs3//cexp3(1:ic3)
      st3 = sta//'3: '//cs3//cexp3(1:ic3)//stc//cs3//cexp3(1:ic3)
      st4 = sta//'4: '//cs3//cexp3(1:ic3)//std//cs3//cexp3(1:ic3)
c
      st5 = stb//'1:  '//cs1//cexp1(1:ic1)//stc//cs3//cexp3(1:ic3)
      st6 = stb//'2:  '//cs1//cexp1(1:ic1)//std//cs3//cexp3(1:ic3)
      st7 = stb//'3:  '//cs4//cexp4(1:ic4)//stc//cs3//cexp3(1:ic3)
      st8 = stb//'4:  '//cs4//cexp4(1:ic4)//std//cs3//cexp3(1:ic3)
c
      iloop = ntr(2)
      ky1   = idetect1
      ky2   =    0
c
      write(2, 9100) critsn
      write(2, 9101) st1  
      write(2, 9101) st2  
      write(2, 9101) st5  
      write(2, 9101) st6  
      write(2, 9101) st7  
      write(2, 9101) st8  
      write(2, 9110)
      write(2, 9111)
      write(2, 9112)
      write(2, 9113)
      write(2, 9114)
      write(2, 9115)
c
      write(2, 9130)
      do 1630 kk = 1, iloop
        itt = isize2(kk)
        ky2 = ky1 + itt - 1
        write(2, 9132) ky1, ky2, itt, psn1(kk), psn2(kk), 
     &  psn3(kk), psn4(kk), ptrend1(kk), ptrend2(kk), 
     $  bserr1(kk), bserr2(kk), dserr1(kk), dserr2(kk),
     &  inumb1(kk), inumb2(kk)
 1630 continue
      write(2, 9133) ipsn1, ipsn2, ipsn3, ipsn4
      write(2, 9134) rpsn1, rpsn2, rpsn3, rpsn4
c
c-------------------------------------------------------------------------------
c     ** 16.4 **
c
c     Write S/N results for SET3 data.
c-------------------------------------------------------------------------------
c
      iloop = ntr(3)
      ky1   = idetect1
      ky2   =    0
c
      write(2, 9100) critsn
      write(2, 9101) st3  
      write(2, 9101) st4  
      write(2, 9101) st5  
      write(2, 9101) st6  
      write(2, 9101) st7  
      write(2, 9101) st8  
      write(2, 9120)
      write(2, 9121)
      write(2, 9122)
      write(2, 9123)
      write(2, 9114)
      write(2, 9115)
c
      write(2, 9131)
      do 1640 kk = 1, iloop
        itt = isize3(kk)
        ky2 = ky1 + itt - 1
        write(2, 9132) ky1, ky2, itt, gsn1(kk), gsn2(kk),
     &  gsn3(kk), gsn4(kk), gtrend1(kk), gtrend2(kk),
     $  bserr1(kk), bserr2(kk), dserr1(kk), dserr2(kk),
     &  inumb1(kk), inumb2(kk)
 1640 continue
      write(2, 9133) igsn1, igsn2, igsn3, igsn4
c                                                                                  
c-------------------------------------------------------------------------------
c     ** 16.5 **
c
c     For EOFs analysis of SET1 and SET3 anomaly data:
c
c     Write eigenvalues, explained variance, and cumulative explained variance.
c-------------------------------------------------------------------------------
c
      write(2, 9200) cmod1, cexp1
      write(2, 9201)
      do 1650 ix = 1, irank1
        write(2, 9202) ix, val1(ix), var1(ix), cum1(ix)
 1650 continue
      write(2, 9999)
c
      write(2, 9300) cmod3, cexp3
      write(2, 9301)
      do 1651 ix = 1, irank2
        write(2, 9302) ix, val2(ix), var2(ix), cum2(ix)
 1651 continue
      write(2, 9999)
c
c
c
c
c
c
c-------------------------------------------------------------------------------
c     ** 17.0 **
c
c     OUTPUT PRINCIPAL COMPONENTS.
c
c     First write loadings of anomaly data on control run EOFs.    
c-------------------------------------------------------------------------------
c
      do 1701 nn = 1, 4
        if (nn .eq. 1) then
          write(2, 9400) cmod1, cexp1
          kyear = nt1
        else if (nn .eq. 2) then
          write(2, 9400) cmod2, cexp2
          kyear = nt2
        else if (nn .eq. 3) then
          write(2, 9400) cmod3, cexp3
          kyear = nt3
        else if (nn .eq. 4) then
          write(2, 9400) cmod4, cexp4
          kyear = nt4
        end if
c
        write(2, 9401) cmod1, cexp1
        write(2, 9402)
c
        do 1700 it = 1, kyear
          if (nn .eq. 1) then
            if (kyear .eq. 1) print *, 'itt = ', itt
            itt = int(rtime1(it))
            write(2, 9403) itt, (bpc2(it, m), m = 1, npcout)
          else if (nn .eq. 2) then
            itt = int(rtime2(it))
            write(2, 9403) itt, (ppc2(it, m), m = 1, npcout)
          else if (nn .eq. 3) then
            itt = int(rtime3(it))
            write(2, 9403) itt, (gpc2(it, m), m = 1, npcout)
          else if (nn .eq. 4) then
            itt = int(rtime4(it))
            write(2, 9403) itt, (dpc2(it, m), m = 1, npcout)
          end if
 1700   continue
        write(2, 9999)
 1701 continue
c
c-------------------------------------------------------------------------------
c     ** 17.1 **
c
c     Write loadings of anomaly data on guess EOFs.
c-------------------------------------------------------------------------------
c
      do 1711 nn = 1, 4
        if (nn .eq. 1) then
          write(2, 9500) cmod1, cexp1
          kyear = nt1
        else if (nn .eq. 2) then
          write(2, 9500) cmod2, cexp2
          kyear = nt2
        else if (nn .eq. 3) then
          write(2, 9500) cmod3, cexp3
          kyear = nt3
        else if (nn .eq. 4) then
          write(2, 9500) cmod4, cexp4
          kyear = nt4
        end if
c
        write(2, 9501) cmod3, cexp3
        write(2, 9502)
c
        do 1710 it = 1, kyear
          if (nn .eq. 1) then
            itt = int(rtime1(it))
            write(2, 9503) itt, (bpc1(it, m), m = 1, npcout), bpc3(it)
          else if (nn .eq. 2) then
            itt = int(rtime2(it))
            write(2, 9503) itt, (ppc1(it, m), m = 1, npcout), ppc3(it)
          else if (nn .eq. 3) then
            itt = int(rtime3(it))
            write(2, 9503) itt, (gpc1(it, m), m = 1, npcout), gpc3(it)
          else if (nn .eq. 4) then
            itt = int(rtime4(it))
            write(2, 9503) itt, (dpc1(it, m), m = 1, npcout), dpc3(it)
          end if
 1710   continue
        write(2, 9999)
 1711 continue                                          
c
c
c
c
c
c
c-------------------------------------------------------------------------------
c     ** 18.0 **
c
c     FORMAT STATEMENTS.
c-------------------------------------------------------------------------------
c
 9000 format(t2, '***** Results from VERTICAL_MSU3_ZM1.F ',
     &'*****' / )
c
 9001 format(t2, ' 1. Eigenvector and PC analysis of IPCC AR5 ',
     &'historical/rcp8.5 and piControl data. Variable processed: ', a8)
c
 9010 format(t2, ' 2. SET1: NOISE FOR OPTIMIZATION. Years extracted: ',
     &i4, '-', i4, '. Reference period: ', i4, '-', i4)
 9011 format(t2, ' 3. SET2: OBSERVATIONAL DATA.     Years extracted: ',
     &i4, '-', i4, '. Reference period: ', i4, '-', i4)
 9012 format(t2, ' 4. SET3: DATA FOR FINGERPRINT.   Years extracted: ',
     &i4, '-', i4, '. Reference period: ', i4, '-', i4)
 9013 format(t2, ' 5. SET4: NOISE FOR SIG. TESTING. Years extracted: ',
     &i4, '-', i4, '. Reference period: ', i4, '-', i4)
c
 9020 format(t2, ' 6. SET1: ', a125)
 9021 format(t2, ' 7. SET2: ', a125)
 9022 format(t2, ' 8. SET3: ', a125)
 9023 format(t2, ' 9. SET4: ', a125)
c
 9030 format(t2, '10. SET1 maximum/actual dimensions:  ',
     &'nyh = ', i4, ' nyh1 = ', i4, ' n1 = ', i4, ' nt1 = ', i4)
 9031 format(t2, '11. SET2 maximum/actual dimensions:  ',
     &'nyh = ', i4, ' nyh1 = ', i4, ' n2 = ', i4, ' nt2 = ', i4)
 9032 format(t2, '12. SET3 maximum/actual dimensions:  ',
     &'nyh = ', i4, ' nyh1 = ', i4, ' n3 = ', i4, ' nt3 = ', i4)
 9033 format(t2, '13. SET4 maximum/actual dimensions:  ',
     &'nyh = ', i4, ' nyh1 = ', i4, ' n4 = ', i4, ' nt4 = ', i4)
c
 9040 format(t2, '14. SET1 EOFs: sample space: nyh1 > nt - 1. ',
     &'Rank: ', i4)
 9041 format(t2, '14. SET1 EOFs: state space:  nyh1 < nt - 1. ', 
     &'Rank: ', i4)
c
 9050 format(t2, '15. SET3 EOFs: sample space: nyh1 > nt - 1. ', 
     &'Rank: ', i4)
 9051 format(t2, '15. SET3 EOFs: state space:  nyh1 < nt - 1. ', 
     &'Rank: ', i4)
c
 9060 format(t2, '16. Rotated Guess. Truncation dimension: ', i3)
c
 9070 format(t2, '17. Input data available on a 5 x 5 latitude/',
     &'pressure grid. No regridding. Spatial means included')
 9071 format(t2, '17. Input data available on a 5 x 5 latitude/',
     &'pressure grid. No regridding. Spatial means subtracted')
 9072 format(t2, '17. Input data available on a 5 x 5 latitude/',
     &'pressure grid. No regridding. STRAT and TROP spatial ',
     &'means subtracted')
 9073 format(t2, '18. Latitude boundaries: ', f6.2, 'S -', f6.2, 'N')
 9074 format(t2, '18. Latitude boundaries: ', f6.2, 'N -', f6.2, 'N')
 9075 format(t2, '18. Latitude boundaries: ', f6.2, 'S -', f6.2, 'S')
c
 9080 format(t2, '19. SET1 area weights on input: ', a8)       
 9081 format(t2, '20. SET2 area weights on input: ', a8)       
 9082 format(t2, '21. SET3 area weights on input: ', a8)       
 9083 format(t2, '22. SET4 area weights on input: ', a8)       
c
 9084 format(t2, '23. Raw fingerprint is EOF 1 of SET3 data')
 9085 format(t2, '23. Raw fingerprint is EOF 2 of SET3 data')
 9086 format(t2, '23. Raw fingerprint is weighted linear ',
     &'combination of first ', i1, ' SET3 EOFs')
c
 9600 format(t2, '24. No residual fingerprint included')     
 9601 format(t2, '24. Inclusion of non-optimized residual')  
 9602 format(t2, '24. Inclusion of residual with average ',
     &'optimization')
c
 9605 format(t2, '25. No temporal filtering of anomaly data'/)         
 9606 format(t2, '25. Temporal smoothing of anomaly data ',
     &'with 13-term Gaussian filter'/)         
 9607 format(t2, '25. Temporal smoothing of anomaly data ',
     &'with 13-term moving average filter'/)         
 9608 format(t2, '25. Temporal smoothing of anomaly data ',
     &'with 5-term binomial filter'/)         
c
 9700 format(t2, 'OPTIMIZATION DIAGNOSTICS')
 9701 format(t2, 'Notation')                     
 9702 format(t2, ' c(x,t): Control run data used for optimizing ',
     &'fingerprint')
 9703 format(t2, 'c1(x,t): Control run data used for estimating ',
     &'natural variability statistics')
 9704 format(t2, ' s(x,t): Data set used for estimating ',
     &'raw fingerprint, f(x), and optimized fingerprint, f*m(x)')
 9705 format(t2, ' o(x,t): Data set in which f(x) and f*m(x) ',
     &'are searched for')
c
 9710 format(t2, ' c(x,t): ', a14, 2x, a14)
 9711 format(t2, 'c1(x,t): ', a14, 2x, a14)
 9712 format(t2, ' s(x,t): ', a14, 2x, a14)
 9713 format(t2, ' o(x,t): ', a14, 2x, a14/)
c
 9750 format(t2, 'CNORM:  Length **2 of fpro ', f10.4)
 9751 format(t2, 'DNORM:  Length **2 of fres ', f10.4)
 9752 format(t2, 'ANORM1: Length **2 of frot before ',
     &'normalization (rotation by var1)', f10.4) 
 9753 format(t2, 'ANORM2: Length **2 of frot before ',
     &'normalization (rotation by val1)', f10.4) 
 9754 format(t2, 'ENORM:  Length **2 of greof before ',
     &'normalization', f10.4)                                
 9755 format(t2, 'sqrt{ANORM1/CNORM}:     Weight of residual ',
     &'in no residual optimization case', 5x, f10.4, ' (OLD)')
 9756 format(t2, '1/(1-CUM1(ITRUNC)/100): Weight of residual ',
     &'in average residual optimization case', f10.4, ' (OLD)')
 9757 format(t2, 'BETA: Weighting of the normalized ',
     &'residual fingerprint ', f6.4, ' (NEW)'/)
c
 9720 format(t2, 'D01: Percentage of total o(x,t)  variance ',
     &'explained by projection onto jth EOF of c(x,t)')      
 9721 format(t2, 'D02: Percentage of total s(x,t)  variance ',
     &'explained by projection onto jth EOF of c(x,t)')      
 9722 format(t2, 'D03: Percentage of total c1(x,t) variance ',
     &'explained by projection onto jth EOF of c(x,t)')        
 9723 format(t2, 'D04: Percentage variance explained by ',
     &'jth EOF of c(x,t)')
 9724 format(t2, 'D05: Ratio D01 / D03') 
 9725 format(t2, 'D06: Ratio D02 / D03') 
c
 9730 format(t2, 'D07: Percentage of total o(x,t)  variance ',
     &'explained by projection onto raw fingerprint, f(x)')  
 9731 format(t2, 'D08: Percentage of total s(x,t)  variance ',
     &'explained by projection onto raw fingerprint, f(x)')  
 9732 format(t2, 'D09: Percentage of total c1(x,t) variance ',
     &'explained by projection onto raw fingerprint, f(x)')   
 9733 format(t2, 'D10: Percentage of total c(x,t)  variance ',
     &'explained by projection onto raw fingerprint, f(x)')   
 9734 format(t2, 'D11: Ratio D07 / D09')
 9735 format(t2, 'D12: Ratio D08 / D09')
c
 9740 format(t2, 'D13: Percentage of total o(x,t)  variance ',
     &'explained by projection onto rotated fingerprint, f*m(x)')
 9741 format(t2, 'D14: Percentage of total s(x,t)  variance ',
     &'explained by projection onto rotated fingerprint, f*m(x)')
 9742 format(t2, 'D15: Percentage of total c1(x,t) variance ',
     &'explained by projection onto rotated fingerprint, f*m(x)')
 9743 format(t2, 'D16: Percentage of total c(x,t)  variance ',
     &'explained by projection onto rotated fingerprint, f*m(x)')
 9744 format(t2, 'D17: Ratio D13 / D15')
 9745 format(t2, 'D18: Ratio D14 / D15')

 9610 format(t5, 'J', t16, 'DIAGNOSTIC 1', 3x, 'DIAGNOSTIC 2', 3x,
     &'DIAGNOSTIC 3', 3x, 'DIAGNOSTIC 4', 3x, 'DIAGNOSTIC 5', 3x,
     &'DIAGNOSTIC 6' )
 9611 format(t5, ' ', t16, 'DIAGNOSTIC 7', 3x, 'DIAGNOSTIC 8', 3x,
     &'DIAGNOSTIC 9', 3x, 'DIAGNOSTIC10', 3x, 'DIAGNOSTIC11', 3x,
     &'DIAGNOSTIC12' )
 9612 format(t5, ' ', t16, 'DIAGNOSTIC13', 3x, 'DIAGNOSTIC14', 3x,
     &'DIAGNOSTIC15', 3x, 'DIAGNOSTIC16', 3x, 'DIAGNOSTIC17', 3x,
     &'DIAGNOSTIC18' )
c
 9615 format(t2, i4, 7x, 6f15.6)
 9616 format(t2, 11x, 6f15.6)
c
 9090 format(t2, '    Number of noise EOFs required to explain ',
     &'GE 95% of total variance of SET1 data: ', i3)
 9091 format(t2, '    Number of guess EOFs required to explain ',
     &'GE 95% of total variance of SET3 data: ', i3)
 9092 format(t2, '    Pattern correlations (centered) for ',
     &'noise (NEOF) versus guess (GEOF) EOFs' )
 9093 format(t2, '    Pattern correlations (uncentered) for ',
     &'noise (NEOF) versus guess (GEOF) EOFs' )
 9094 format(t12, 'GEOF1', t22, 'GEOF2', t32, 
     &'GEOF3', t42, 'GEOF4', t52, 'GEOF5', t62, 'GEOF6', t72, 
     &'GEOF7', t82, 'GEOF8', t92, 'GEOF9', t102, 'GEOF10' )
 9095 format(t2, 'NEOF', i2, 10f10.4) 
c
 9100 format(t2, '    Signal-to-noise analysis', 2x, 
     &'CRITSN = ',f7.4)
 9101 format(t2, a70)                           
c
 9110 format(t2, '    PSN1 = SIGNAL1 / NOISE1')
 9111 format(t2, '    PSN2 = SIGNAL2 / NOISE2')
 9112 format(t2, '    PSN3 = SIGNAL1 / NOISE3')
 9113 format(t2, '    PSN4 = SIGNAL2 / NOISE4')
 9114 format(t2, '    NUM1 = Number of chunks in ',
     &'NOISE1/2')
 9115 format(t2, '    NUM2 = Number of chunks in ',
     &'NOISE3/4')
c
 9120 format(t2, '    GSN1 = SIGNAL3 / NOISE1')
 9121 format(t2, '    GSN2 = SIGNAL4 / NOISE2')
 9122 format(t2, '    GSN3 = SIGNAL3 / NOISE3')
 9123 format(t2, '    GSN4 = SIGNAL4 / NOISE4')
c
 9130 format(t2, 'START', 2x, 'END', 4x, 'LENGTH',
     &t28, 'PSN1', t38, 'PSN2', t48, 'PSN3', t58, 'PSN4', 
     &t68, 'SIG1', t78, 'SIG2', t89, 'NSE1', t99, 'NSE2',
     &t109, 'NSE3', t119, 'NSE4', t129, 'NUM1', t139, 'NUM2')
 9131 format(t2, 'START', 2x, 'END', 4x, 'LENGTH',
     &t28, 'GSN1', t38, 'GSN2', t48, 'GSN3', t58, 'GSN4',
     &t68, 'SIG3', t78, 'SIG4', t89, 'NSE1', t99, 'NSE2',
     &t109, 'NSE3', t119, 'NSE4', t129, 'NUM1', t139, 'NUM2')
 9132 format(t2, i4, 3x, i4, 3x, i4, 3x, 10f10.4, 2i10)
 9133 format(t2, 'DETECTION TIME', 7x, 4i10)
 9134 format(t2, 'DETECTION TIME', 7x, 4f10.3)
c
 9200 format(t2, '    Eigenvalues and variance for ', 2a20 )
 9201 format(t2, 'Number', t10, 'Eigenvalue', t25, '% Variance',
     &t40, 'Cumulative variance' )
 9202 format(t2, i4, 3f15.5)
c
c
 9300 format(t2, '    Eigenvalues and variance for ', 2a20 )
 9301 format(t2, 'Number', t10, 'Eigenvalue', t25, '% Variance',
     &t40, 'Cumulative variance' )
 9302 format(t2, i4, 3f15.5)
c
 9400 format(t2, '    PC time series. Anomaly data: ', 2a20 )
 9401 format(t2, '    Projected onto EOFs of: ', 2a20 )
 9402 format(t2, 'Time', t11, 'PC1', t21, 'PC2', t31, 'PC3', t41,
     &'PC4', t51, 'PC5', t61, 'PC6', t71, 'PC7', t81, 'PC8', t91,
     &'PC9', t100, 'PC10' )
 9403 format(t2, i4, 1x, 10f10.4) 
c
 9500 format(t2, '    PC time series. Anomaly data: ', 2a20 )
 9501 format(t2, '    Projected onto EOFs of: ', 2a20 )
 9502 format(t2, 'Time', t11, 'PC1', t21, 'PC2', t31, 'PC3', t41,
     &'PC4', t51, 'PC5', t61, 'PC6', t71, 'PC7', t81, 'PC8', t91,
     &'PC9', t100, 'PC10', t110, 'ROT' )
 9503 format(t2, i4, 1x, 11f10.4) 
c
 9999 format(/)
c
c
c
      stop
      end
c-------------------------------------------------------------------------------
c-------------------------------------------------------------------------------
c     SUBROUTINE CREATE_PATH1
c
c     Date:    07-NOV-2022.
c     Author:  Ben Santer, LLNL, Livermore.
c     Version: # 1.0 / Single precision.
c
c
c
c     This subroutine reads creates the full path to an input dataset on
c     /pcmdi/bsanter1.
c
c
c
      subroutine create_path1(cmod, cexp, varnam, fpath1, fpath2)
c
c
c
      implicit real (o-z, a-h)
c
c     ** INPUT **
c     character cmod*24          Model identifier
c     character cexp*24          Experiment identifier
c     character varnam*16        Variable identifier
c
c     ** OUTPUT **
c     character fpath1*120       Full path to input tt6 or sm6 data
c     character fpath2*120       Full path to input sftlf6 data
c
c     ** INTERNAL **
c     character pp1a*29          First part of path to CMIP5 data
c     character pp1b*29          First part of path to CMIP6 data
c     character pp2*23           First part of path to OBS       
c     character p0*90            String
c     character p2*80            String
c     character p3*40            String
c     character p4*30            String
c
c
c
      character cmod*24, cexp*24, varnam*16
      character pp1a*29, pp1b*29, pp2*23
      character p0*90, p2*80, p3*40, p4*30, fpath1*120, fpath2*120
c
      data pp1a / '/export/bsanter1/cmip5/mma05/' /
c
      data pp1b / '/export/bsanter1/cmip6/mma05/' /
c
      data pp2 / '/export/bsanter1/obs05/'       /
c
c
c
c-------------------------------------------------------------------------------
c     ** 1.0 ** 
c
c     Create path.
c-------------------------------------------------------------------------------
c
      ib1 = index(cmod, ' ') - 1
c
      ib2      = index(cexp, ' ')
      ib2      = ib2 - 1      
      ib3      = index(varnam, ' ')
      ib3      = ib3 - 1      
      print *, 'In CREATE_PATH1! pp1a = ', pp1a
      print *, 'In CREATE_PATH1! pp1b = ', pp1b
      print *, 'In CREATE_PATH1! pp2  = ', pp2 
c
      if (ib1 .lt. 0) then
        print *, 'In CREATE_PATH1!'
        print *, 'No blanks in model name (cmod)!'
      end if
      if (ib2 .lt. 0) then
        print *, 'In CREATE_PATH1!'
        print *, 'No blanks in expt name (cexp)!'
      end if
      if (ib3 .lt. 0) then
        print *, 'In CREATE_PATH1!'
        print *, 'No blanks in variable name (varnam)!'
      end if
c
c-------------------------------------------------------------------------------
c     ** 1.1 ** 
c
c     Set path to files containing multi-model averages from CMIP5 extended 
c     historical simulations, and from CMIP5 concatenated piControl runs.
c-------------------------------------------------------------------------------
c
      if (cmod(1:ib1) .eq. 'MMA_CMIP5_28models') then
        p0 = pp1a//cmod(1:ib1)//'/Masks/sftlf6_'//cmod(1:ib1)
        p2 = pp1a//cmod(1:ib1)//'/'//cexp(1:ib2)//'/'//'Yz/'
        if (cexp(1:ib2) .eq. 'hist_rcp85') p4= 'a0_yz_EM_r0000_0000.nc'
c
      else if (cmod(1:ib1) .eq. 'MMA_CMIP5_36models') then
        p0 = pp1a//cmod(1:ib1)//'/Masks/sftlf6_'//cmod(1:ib1)
        p2 = pp1a//cmod(1:ib1)//'/'//cexp(1:ib2)//'/'//'Yz/'
        if (cexp(1:ib2) .eq. 'piControl')  p4= 'a0_yz_CC_r0000_0000.nc'
c
c-------------------------------------------------------------------------------
c     ** 1.2 ** 
c
c     Set path to files containing multi-model averages from CMIP6 extended 
c     historical simulations, and from CMIP6 concatenated piControl runs.
c-------------------------------------------------------------------------------
c
      else if (cmod(1:ib1) .eq. 'MMA_CMIP6_21models') then
        p0 = pp1b//cmod(1:ib1)//'/Masks/sftlf6_'//cmod(1:ib1)
        p2 = pp1b//cmod(1:ib1)//'/'//cexp(1:ib2)//'/'//'Yz/'
        if (cexp(1:ib2) .eq. 'hst_ssp585') p4= 'a0_yz_EM_r0000_0000.nc'
c
      else if (cmod(1:ib1) .eq. 'MMA_CMIP6_30models') then
        p0 = pp1b//cmod(1:ib1)//'/Masks/sftlf6_'//cmod(1:ib1)
        p2 = pp1b//cmod(1:ib1)//'/'//cexp(1:ib2)//'/'//'Yz/'
        if (cexp(1:ib2) .eq. 'piControl')  p4= 'a0_yz_CC_r0000_0000.nc'
c
      else if (cmod(1:ib1) .eq. 'MMA_CMIP6_09models') then
        p0 = pp1b//cmod(1:ib1)//'/Masks/sftlf6_'//cmod(1:ib1)
        p2 = pp1b//cmod(1:ib1)//'/'//cexp(1:ib2)//'/'//'Yz/'
        if (cexp(1:ib2) .eq. 'piControl')  p4= 'am_zm_CC_r0000_0000.nc'
        if (cexp(1:ib2) .eq. 'hst_ssp585') p4= 'a0_yz_EM_r0000_0000.nc'
c
c-------------------------------------------------------------------------------
c     ** 1.3 ** 
c
c     Set path to files containing observations.
c-------------------------------------------------------------------------------
c
      else if (cmod(1:ib1) .eq. 'SSSRRR_jan25_05x05') then
        p0 = pp2//cmod(1:ib1)//'/Masks/sftlf6_'//cmod(1:ib1)
        p2 = pp2//cmod(1:ib1)//'/'//cexp(1:ib2)//'/'//'Yz/'
        if (cexp(1:ib2) .eq. '20c3m_run1') p4= 'a0_yz_AM_r0000_0000.nc'
c
      else if (cmod(1:ib1) .eq. 'SSSUUU_jan25_05x05') then
        p0 = pp2//cmod(1:ib1)//'/Masks/sftlf6_'//cmod(1:ib1)
        p2 = pp2//cmod(1:ib1)//'/'//cexp(1:ib2)//'/'//'Yz/'
        if (cexp(1:ib2) .eq. '20c3m_run1') p4= 'a0_yz_AM_r0000_0000.nc'
c
      else if (cmod(1:ib1) .eq. 'SSSSSS_jan25_05x05') then
        p0 = pp2//cmod(1:ib1)//'/Masks/sftlf6_'//cmod(1:ib1)
        p2 = pp2//cmod(1:ib1)//'/'//cexp(1:ib2)//'/'//'Yz/'
        if (cexp(1:ib2) .eq. '20c3m_run1') p4= 'a0_yz_AM_r0000_0000.nc'
c
c-------------------------------------------------------------------------------
c     ** 1.4 ** 
c
c     Set path to files containing MSU data from individual CMIP5 models.
c-------------------------------------------------------------------------------
c
      else if (cmod(1:ib1) .eq. 'access1_0_05x05') then
        p0 = pp1a//cmod(1:ib1)//'/Masks/sftlf6_'//cmod(1:ib1)
        p2 = pp1a//cmod(1:ib1)//'/'//cexp(1:10)//'/'//'Yz/'
        if (cexp(1:13) .eq. 'hist_rcp85.r1') p4= 'r1i1p1_yz_AM.nc'
c
      else if (cmod(1:ib1) .eq. 'access1_3_05x05') then
        p0 = pp1a//cmod(1:ib1)//'/Masks/sftlf6_'//cmod(1:ib1)
        p2 = pp1a//cmod(1:ib1)//'/'//cexp(1:10)//'/'//'Yz/'
        if (cexp(1:13) .eq. 'hist_rcp85.r1') p4= 'r1i1p1_yz_AM.nc'
c
      else if (cmod(1:ib1) .eq. 'bcc_csm1_1_05x05') then
        p0 = pp1a//cmod(1:ib1)//'/Masks/sftlf6_'//cmod(1:ib1)
        p2 = pp1a//cmod(1:ib1)//'/'//cexp(1:10)//'/'//'Yz/'
        if (cexp(1:13) .eq. 'hist_rcp85.r1') p4= 'r1i1p1_yz_AM.nc'
c
      else if (cmod(1:ib1) .eq. 'bcc_csm1_1m_05x05') then
        p0 = pp1a//cmod(1:ib1)//'/Masks/sftlf6_'//cmod(1:ib1)
        p2 = pp1a//cmod(1:ib1)//'/'//cexp(1:10)//'/'//'Yz/'
        if (cexp(1:13) .eq. 'hist_rcp85.r1') p4= 'r1i1p1_yz_AM.nc'
c
      else if (cmod(1:ib1) .eq. 'can_esm2_05x05') then
        p0 = pp1a//cmod(1:ib1)//'/Masks/sftlf6_'//cmod(1:ib1)
        p2 = pp1a//cmod(1:ib1)//'/'//cexp(1:10)//'/'//'Yz/'
        if (cexp(1:13) .eq. 'hist_LE_r1.r1')  p4= 'r1i1p1_yz_AM.nc'
        if (cexp(1:13) .eq. 'hist_LE_r1.r2')  p4= 'r2i1p1_yz_AM.nc'
        if (cexp(1:13) .eq. 'hist_LE_r1.r3')  p4= 'r3i1p1_yz_AM.nc'
        if (cexp(1:13) .eq. 'hist_LE_r1.r4')  p4= 'r4i1p1_yz_AM.nc'
        if (cexp(1:13) .eq. 'hist_LE_r1.r5')  p4= 'r5i1p1_yz_AM.nc'
        if (cexp(1:13) .eq. 'hist_LE_r1.r6')  p4= 'r6i1p1_yz_AM.nc'
        if (cexp(1:13) .eq. 'hist_LE_r1.r7')  p4= 'r7i1p1_yz_AM.nc'
        if (cexp(1:13) .eq. 'hist_LE_r1.r8')  p4= 'r8i1p1_yz_AM.nc'
        if (cexp(1:13) .eq. 'hist_LE_r1.r9')  p4= 'r9i1p1_yz_AM.nc'
        if (cexp(1:13) .eq. 'hist_LE_r1.r10') p4= 'r10i1p1_yz_AM.nc'
        if (cexp(1:13) .eq. 'hist_LE_r2.r1')  p4= 'r1i1p1_yz_AM.nc'
        if (cexp(1:13) .eq. 'hist_LE_r2.r2')  p4= 'r2i1p1_yz_AM.nc'
        if (cexp(1:13) .eq. 'hist_LE_r2.r3')  p4= 'r3i1p1_yz_AM.nc'
        if (cexp(1:13) .eq. 'hist_LE_r2.r4')  p4= 'r4i1p1_yz_AM.nc'
        if (cexp(1:13) .eq. 'hist_LE_r2.r5')  p4= 'r5i1p1_yz_AM.nc'
        if (cexp(1:13) .eq. 'hist_LE_r2.r6')  p4= 'r6i1p1_yz_AM.nc'
        if (cexp(1:13) .eq. 'hist_LE_r2.r7')  p4= 'r7i1p1_yz_AM.nc'
        if (cexp(1:13) .eq. 'hist_LE_r2.r8')  p4= 'r8i1p1_yz_AM.nc'
        if (cexp(1:13) .eq. 'hist_LE_r2.r9')  p4= 'r9i1p1_yz_AM.nc'
        if (cexp(1:13) .eq. 'hist_LE_r2.r10') p4= 'r10i1p1_yz_AM.nc'
        if (cexp(1:13) .eq. 'hist_LE_r3.r1')  p4= 'r1i1p1_yz_AM.nc'
        if (cexp(1:13) .eq. 'hist_LE_r3.r2')  p4= 'r2i1p1_yz_AM.nc'
        if (cexp(1:13) .eq. 'hist_LE_r3.r3')  p4= 'r3i1p1_yz_AM.nc'
        if (cexp(1:13) .eq. 'hist_LE_r3.r4')  p4= 'r4i1p1_yz_AM.nc'
        if (cexp(1:13) .eq. 'hist_LE_r3.r5')  p4= 'r5i1p1_yz_AM.nc'
        if (cexp(1:13) .eq. 'hist_LE_r3.r6')  p4= 'r6i1p1_yz_AM.nc'
        if (cexp(1:13) .eq. 'hist_LE_r3.r7')  p4= 'r7i1p1_yz_AM.nc'
        if (cexp(1:13) .eq. 'hist_LE_r3.r8')  p4= 'r8i1p1_yz_AM.nc'
        if (cexp(1:13) .eq. 'hist_LE_r3.r9')  p4= 'r9i1p1_yz_AM.nc'
        if (cexp(1:13) .eq. 'hist_LE_r3.r10') p4= 'r10i1p1_yz_AM.nc'
        if (cexp(1:13) .eq. 'hist_LE_r4.r1')  p4= 'r1i1p1_yz_AM.nc'
        if (cexp(1:13) .eq. 'hist_LE_r4.r2')  p4= 'r2i1p1_yz_AM.nc'
        if (cexp(1:13) .eq. 'hist_LE_r4.r3')  p4= 'r3i1p1_yz_AM.nc'
        if (cexp(1:13) .eq. 'hist_LE_r4.r4')  p4= 'r4i1p1_yz_AM.nc'
        if (cexp(1:13) .eq. 'hist_LE_r4.r5')  p4= 'r5i1p1_yz_AM.nc'
        if (cexp(1:13) .eq. 'hist_LE_r4.r6')  p4= 'r6i1p1_yz_AM.nc'
        if (cexp(1:13) .eq. 'hist_LE_r4.r7')  p4= 'r7i1p1_yz_AM.nc'
        if (cexp(1:13) .eq. 'hist_LE_r4.r8')  p4= 'r8i1p1_yz_AM.nc'
        if (cexp(1:13) .eq. 'hist_LE_r4.r9')  p4= 'r9i1p1_yz_AM.nc'
        if (cexp(1:13) .eq. 'hist_LE_r4.r10') p4= 'r10i1p1_yz_AM.nc'
        if (cexp(1:13) .eq. 'hist_LE_r5.r1')  p4= 'r1i1p1_yz_AM.nc'
        if (cexp(1:13) .eq. 'hist_LE_r5.r2')  p4= 'r2i1p1_yz_AM.nc'
        if (cexp(1:13) .eq. 'hist_LE_r5.r3')  p4= 'r3i1p1_yz_AM.nc'
        if (cexp(1:13) .eq. 'hist_LE_r5.r4')  p4= 'r4i1p1_yz_AM.nc'
        if (cexp(1:13) .eq. 'hist_LE_r5.r5')  p4= 'r5i1p1_yz_AM.nc'
        if (cexp(1:13) .eq. 'hist_LE_r5.r6')  p4= 'r6i1p1_yz_AM.nc'
        if (cexp(1:13) .eq. 'hist_LE_r5.r7')  p4= 'r7i1p1_yz_AM.nc'
        if (cexp(1:13) .eq. 'hist_LE_r5.r8')  p4= 'r8i1p1_yz_AM.nc'
        if (cexp(1:13) .eq. 'hist_LE_r5.r9')  p4= 'r9i1p1_yz_AM.nc'
        if (cexp(1:13) .eq. 'hist_LE_r5.r10') p4= 'r10i1p1_yz_AM.nc'
c
      else if (cmod(1:ib1) .eq. 'ccsm4_05x05') then
        p0 = pp1a//cmod(1:ib1)//'/Masks/sftlf6_'//cmod(1:ib1)
        p2 = pp1a//cmod(1:ib1)//'/'//cexp(1:10)//'/'//'Yz/'
        if (cexp(1:13) .eq. 'hist_rcp85.r1') p4= 'r1i1p1_yz_AM.nc'
        if (cexp(1:13) .eq. 'hist_rcp85.r2') p4= 'r2i1p1_yz_AM.nc'
        if (cexp(1:13) .eq. 'hist_rcp85.r3') p4= 'r3i1p1_yz_AM.nc'
c
      else if (cmod(1:ib1) .eq. 'cesm1_bgc_05x05') then
        p0 = pp1a//cmod(1:ib1)//'/Masks/sftlf6_'//cmod(1:ib1)
        p2 = pp1a//cmod(1:ib1)//'/'//cexp(1:10)//'/'//'Yz/'
        if (cexp(1:13) .eq. 'hist_rcp85.r1') p4= 'r1i1p1_yz_AM.nc'
c
      else if (cmod(1:ib1) .eq. 'cesm1_cam5_05x05') then
        p0 = pp1a//cmod(1:ib1)//'/Masks/sftlf6_'//cmod(1:ib1)
        p2 = pp1a//cmod(1:ib1)//'/'//cexp(1:10)//'/'//'Yz/'
        if (cexp(1:14) .eq. 'hist_LE_r1.r01') p4= 'r01i1p1_yz_AM.nc'
        if (cexp(1:14) .eq. 'hist_LE_r1.r02') p4= 'r02i1p1_yz_AM.nc'
        if (cexp(1:14) .eq. 'hist_LE_r1.r03') p4= 'r03i1p1_yz_AM.nc'
        if (cexp(1:14) .eq. 'hist_LE_r1.r04') p4= 'r04i1p1_yz_AM.nc'
        if (cexp(1:14) .eq. 'hist_LE_r1.r05') p4= 'r05i1p1_yz_AM.nc'
        if (cexp(1:14) .eq. 'hist_LE_r1.r06') p4= 'r06i1p1_yz_AM.nc'
        if (cexp(1:14) .eq. 'hist_LE_r1.r07') p4= 'r07i1p1_yz_AM.nc'
        if (cexp(1:14) .eq. 'hist_LE_r1.r08') p4= 'r08i1p1_yz_AM.nc'
        if (cexp(1:14) .eq. 'hist_LE_r1.r09') p4= 'r09i1p1_yz_AM.nc'
        if (cexp(1:14) .eq. 'hist_LE_r1.r10') p4= 'r10i1p1_yz_AM.nc'
        if (cexp(1:14) .eq. 'hist_LE_r1.r11') p4= 'r11i1p1_yz_AM.nc'
        if (cexp(1:14) .eq. 'hist_LE_r1.r12') p4= 'r12i1p1_yz_AM.nc'
        if (cexp(1:14) .eq. 'hist_LE_r1.r13') p4= 'r13i1p1_yz_AM.nc'
        if (cexp(1:14) .eq. 'hist_LE_r1.r14') p4= 'r14i1p1_yz_AM.nc'
        if (cexp(1:14) .eq. 'hist_LE_r1.r15') p4= 'r15i1p1_yz_AM.nc'
        if (cexp(1:14) .eq. 'hist_LE_r1.r16') p4= 'r16i1p1_yz_AM.nc'
        if (cexp(1:14) .eq. 'hist_LE_r1.r17') p4= 'r17i1p1_yz_AM.nc'
        if (cexp(1:14) .eq. 'hist_LE_r1.r18') p4= 'r18i1p1_yz_AM.nc'
        if (cexp(1:14) .eq. 'hist_LE_r1.r19') p4= 'r19i1p1_yz_AM.nc'
        if (cexp(1:14) .eq. 'hist_LE_r1.r20') p4= 'r20i1p1_yz_AM.nc'
        if (cexp(1:14) .eq. 'hist_LE_r1.r21') p4= 'r21i1p1_yz_AM.nc'
        if (cexp(1:14) .eq. 'hist_LE_r1.r22') p4= 'r22i1p1_yz_AM.nc'
        if (cexp(1:14) .eq. 'hist_LE_r1.r23') p4= 'r23i1p1_yz_AM.nc'
        if (cexp(1:14) .eq. 'hist_LE_r1.r24') p4= 'r24i1p1_yz_AM.nc'
        if (cexp(1:14) .eq. 'hist_LE_r1.r25') p4= 'r25i1p1_yz_AM.nc'
        if (cexp(1:14) .eq. 'hist_LE_r1.r26') p4= 'r26i1p1_yz_AM.nc'
        if (cexp(1:14) .eq. 'hist_LE_r1.r27') p4= 'r27i1p1_yz_AM.nc'
        if (cexp(1:14) .eq. 'hist_LE_r1.r28') p4= 'r28i1p1_yz_AM.nc'
        if (cexp(1:14) .eq. 'hist_LE_r1.r29') p4= 'r29i1p1_yz_AM.nc'
        if (cexp(1:14) .eq. 'hist_LE_r1.r30') p4= 'r30i1p1_yz_AM.nc'
        if (cexp(1:14) .eq. 'hist_LE_r1.r31') p4= 'r31i1p1_yz_AM.nc'
        if (cexp(1:14) .eq. 'hist_LE_r1.r32') p4= 'r32i1p1_yz_AM.nc'
        if (cexp(1:14) .eq. 'hist_LE_r1.r33') p4= 'r33i1p1_yz_AM.nc'
        if (cexp(1:14) .eq. 'hist_LE_r1.r34') p4= 'r34i1p1_yz_AM.nc'
        if (cexp(1:14) .eq. 'hist_LE_r1.r35') p4= 'r35i1p1_yz_AM.nc'
        if (cexp(1:14) .eq. 'hist_LE_r1.r36') p4= 'r36i1p1_yz_AM.nc'
        if (cexp(1:14) .eq. 'hist_LE_r1.r37') p4= 'r37i1p1_yz_AM.nc'
        if (cexp(1:14) .eq. 'hist_LE_r1.r38') p4= 'r38i1p1_yz_AM.nc'
        if (cexp(1:14) .eq. 'hist_LE_r1.r39') p4= 'r39i1p1_yz_AM.nc'
        if (cexp(1:14) .eq. 'hist_LE_r1.r40') p4= 'r40i1p1_yz_AM.nc'
c
      else if (cmod(1:ib1) .eq. 'csiro_mk3_6_0_05x05') then
        p0 = pp1a//cmod(1:ib1)//'/Masks/sftlf6_'//cmod(1:ib1)
        p2 = pp1a//cmod(1:ib1)//'/'//cexp(1:10)//'/'//'Yz/'
        if (cexp(1:14) .eq. 'hist_rcp85.r10') p4= 'r10i1p1_yz_AM.nc'
c
      else if (cmod(1:ib1) .eq. 'ec_earth_05x05') then
        p0 = pp1a//cmod(1:ib1)//'/Masks/sftlf6_'//cmod(1:ib1)
        p2 = pp1a//cmod(1:ib1)//'/'//cexp(1:10)//'/'//'Yz/'
        if (cexp(1:13) .eq. 'hist_rcp85.r8') p4= 'r8i1p1_yz_AM.nc'
c
      else if (cmod(1:ib1) .eq. 'gfdl_cm3_05x05') then
        p0 = pp1a//cmod(1:ib1)//'/Masks/sftlf6_'//cmod(1:ib1)
        p2 = pp1a//cmod(1:ib1)//'/'//cexp(1:10)//'/'//'Yz/'
        if (cexp(1:13) .eq. 'hist_rcp85.r1') p4= 'r1i1p1_yz_AM.nc'
c
      else if (cmod(1:ib1) .eq. 'gfdl_esm2g_05x05') then
        p0 = pp1a//cmod(1:ib1)//'/Masks/sftlf6_'//cmod(1:ib1)
        p2 = pp1a//cmod(1:ib1)//'/'//cexp(1:10)//'/'//'Yz/'
        if (cexp(1:13) .eq. 'hist_rcp85.r1') p4= 'r1i1p1_yz_AM.nc'
c
      else if (cmod(1:ib1) .eq. 'gfdl_esm2m_05x05') then
        p0 = pp1a//cmod(1:ib1)//'/Masks/sftlf6_'//cmod(1:ib1)
        p2 = pp1a//cmod(1:ib1)//'/'//cexp(1:10)//'/'//'Yz/'
        if (cexp(1:13) .eq. 'hist_rcp85.r1') p4= 'r1i1p1_yz_AM.nc'
c
      else if (cmod(1:ib1) .eq. 'giss_e2_h_p1_05x05') then
        p0 = pp1a//cmod(1:ib1)//'/Masks/sftlf6_'//cmod(1:ib1)
        p2 = pp1a//cmod(1:ib1)//'/'//cexp(1:10)//'/'//'Yz/'
        if (cexp(1:13) .eq. 'hist_rcp85.r1') p4= 'r1i1p1_yz_AM.nc'
c
      else if (cmod(1:ib1) .eq. 'giss_e2_h_p3_05x05') then
        p0 = pp1a//cmod(1:ib1)//'/Masks/sftlf6_'//cmod(1:ib1)
        p2 = pp1a//cmod(1:ib1)//'/'//cexp(1:10)//'/'//'Yz/'
        if (cexp(1:13) .eq. 'hist_rcp85.r1') p4= 'r1i1p3_yz_AM.nc'
c
      else if (cmod(1:ib1) .eq. 'giss_e2_r_p1_05x05') then
        p0 = pp1a//cmod(1:ib1)//'/Masks/sftlf6_'//cmod(1:ib1)
        p2 = pp1a//cmod(1:ib1)//'/'//cexp(1:10)//'/'//'Yz/'
        if (cexp(1:13) .eq. 'hist_rcp85.r1') p4= 'r1i1p1_yz_AM.nc'
c
      else if (cmod(1:ib1) .eq. 'giss_e2_r_p2_05x05') then
        p0 = pp1a//cmod(1:ib1)//'/Masks/sftlf6_'//cmod(1:ib1)
        p2 = pp1a//cmod(1:ib1)//'/'//cexp(1:10)//'/'//'Yz/'
        if (cexp(1:13) .eq. 'hist_rcp85.r1') p4= 'r1i1p2_yz_AM.nc'
c
      else if (cmod(1:ib1) .eq. 'giss_e2_r_p3_05x05') then
        p0 = pp1a//cmod(1:ib1)//'/Masks/sftlf6_'//cmod(1:ib1)
        p2 = pp1a//cmod(1:ib1)//'/'//cexp(1:10)//'/'//'Yz/'
        if (cexp(1:13) .eq. 'hist_rcp85.r1') p4= 'r1i1p3_yz_AM.nc'
c
      else if (cmod(1:ib1) .eq. 'hadgem2_cc_05x05') then
        p0 = pp1a//cmod(1:ib1)//'/Masks/sftlf6_'//cmod(1:ib1)
        p2 = pp1a//cmod(1:ib1)//'/'//cexp(1:10)//'/'//'Yz/'
        if (cexp(1:13) .eq. 'hist_rcp85.r1') p4= 'r1i1p1_yz_AM.nc'
        if (cexp(1:13) .eq. 'hist_rcp85.r2') p4= 'r2i1p1_yz_AM.nc'
c
      else if (cmod(1:ib1) .eq. 'hadgem2_es_05x05') then
        p0 = pp1a//cmod(1:ib1)//'/Masks/sftlf6_'//cmod(1:ib1)
        p2 = pp1a//cmod(1:ib1)//'/'//cexp(1:10)//'/'//'Yz/'
        if (cexp(1:13) .eq. 'hist_rcp85.r1') p4= 'r1i1p1_yz_AM.nc'
        if (cexp(1:13) .eq. 'hist_rcp85.r2') p4= 'r2i1p1_yz_AM.nc'
c
      else if (cmod(1:ib1) .eq. 'miroc5_05x05') then
        p0 = pp1a//cmod(1:ib1)//'/Masks/sftlf6_'//cmod(1:ib1)
        p2 = pp1a//cmod(1:ib1)//'/'//cexp(1:10)//'/'//'Yz/'
        if (cexp(1:13) .eq. 'hist_rcp85.r1') p4= 'r1i1p1_yz_AM.nc'
c
      else if (cmod(1:ib1) .eq. 'miroc_esm_05x05') then
        p0 = pp1a//cmod(1:ib1)//'/Masks/sftlf6_'//cmod(1:ib1)
        p2 = pp1a//cmod(1:ib1)//'/'//cexp(1:10)//'/'//'Yz/'
        if (cexp(1:13) .eq. 'hist_rcp85.r1') p4= 'r1i1p1_yz_AM.nc'
c
      else if (cmod(1:ib1) .eq. 'miroc-esm-chem_05x05') then
        p0 = pp1a//cmod(1:ib1)//'/Masks/sftlf6_'//cmod(1:ib1)
        p2 = pp1a//cmod(1:ib1)//'/'//cexp(1:10)//'/'//'Yz/'
        if (cexp(1:13) .eq. 'hist_rcp85.r1') p4= 'r1i1p1_yz_AM.nc'
c
      else if (cmod(1:ib1) .eq. 'mpi_esm_lr_05x05') then
        p0 = pp1a//cmod(1:ib1)//'/Masks/sftlf6_'//cmod(1:ib1)
        p2 = pp1a//cmod(1:ib1)//'/'//cexp(1:10)//'/'//'Yz/'
        if (cexp(1:13) .eq. 'hist_rcp85.r1') p4= 'r1i1p1_yz_AM.nc'
        if (cexp(1:13) .eq. 'hist_rcp85.r2') p4= 'r2i1p1_yz_AM.nc'
        if (cexp(1:13) .eq. 'hist_rcp85.r3') p4= 'r3i1p1_yz_AM.nc'
c
      else if (cmod(1:ib1) .eq. 'mpi_esm_mr_05x05') then
        p0 = pp1a//cmod(1:ib1)//'/Masks/sftlf6_'//cmod(1:ib1)
        p2 = pp1a//cmod(1:ib1)//'/'//cexp(1:10)//'/'//'Yz/'
        if (cexp(1:13) .eq. 'hist_rcp85.r1') p4= 'r1i1p1_yz_AM.nc'
c
      else if (cmod(1:ib1) .eq. 'mri_cgcm3_05x05') then
        p0 = pp1a//cmod(1:ib1)//'/Masks/sftlf6_'//cmod(1:ib1)
        p2 = pp1a//cmod(1:ib1)//'/'//cexp(1:10)//'/'//'Yz/'
        if (cexp(1:13) .eq. 'hist_rcp85.r1') p4= 'r1i1p1_yz_AM.nc'
c
      else if (cmod(1:ib1) .eq. 'noresm1-m_05x05') then
        p0 = pp1a//cmod(1:ib1)//'/Masks/sftlf6_'//cmod(1:ib1)
        p2 = pp1a//cmod(1:ib1)//'/'//cexp(1:10)//'/'//'Yz/'
        if (cexp(1:13) .eq. 'hist_rcp85.r1') p4= 'r1i1p1_yz_AM.nc'
c
      else if (cmod(1:ib1) .eq. 'noresm1_me_05x05') then
        p0 = pp1a//cmod(1:ib1)//'/Masks/sftlf6_'//cmod(1:ib1)
        p2 = pp1a//cmod(1:ib1)//'/'//cexp(1:10)//'/'//'Yz/'
        if (cexp(1:13) .eq. 'hist_rcp85.r1') p4= 'r1i1p1_yz_AM.nc'
c
c-------------------------------------------------------------------------------
c     ** 1.5 ** 
c
c     Set path to files containing MSU data from individual CMIP6 models.
c-------------------------------------------------------------------------------
c
      else if (cmod(1:ib1) .eq. 'access_cm2_05x05') then
        p0 = pp1b//cmod(1:ib1)//'/Masks/sftlf6_'//cmod(1:ib1)
        p2 = pp1b//cmod(1:ib1)//'/'//cexp(1:10)//'/'//'Yz/'
        if (cexp(1:13) .eq. 'hst_ssp585.r1') p4= 'r1i1p1f1_yz_AM.nc'
        if (cexp(1:13) .eq. 'hst_ssp585.r2') p4= 'r2i1p1f1_yz_AM.nc'
        if (cexp(1:13) .eq. 'hst_ssp585.r3') p4= 'r3i1p1f1_yz_AM.nc'
c
      else if (cmod(1:ib1) .eq. 'access_esm1_5_05x05') then
        p0 = pp1b//cmod(1:ib1)//'/Masks/sftlf6_'//cmod(1:ib1)
        p2 = pp1b//cmod(1:ib1)//'/'//cexp(1:10)//'/'//'Yz/'
        if (cexp(1:13) .eq. 'hst_ssp585.r1') p4= 'r1i1p1f1_yz_AM.nc'
        if (cexp(1:13) .eq. 'hst_ssp585.r2') p4= 'r2i1p1f1_yz_AM.nc'
        if (cexp(1:13) .eq. 'hst_ssp585.r3') p4= 'r3i1p1f1_yz_AM.nc'
c
      else if (cmod(1:ib1) .eq. 'cesm2_05x05') then
        p0 = pp1b//cmod(1:ib1)//'/Masks/sftlf6_'//cmod(1:ib1)
        p2 = pp1b//cmod(1:ib1)//'/'//cexp(1:10)//'/'//'Yz/'
        if (cexp(1:14) .eq. 'hst_ssp585.r10') p4= 'r10i1p1f1_yz_AM.nc'
        if (cexp(1:14) .eq. 'hst_ssp585.r11') p4= 'r11i1p1f1_yz_AM.nc'
c
      else if (cmod(1:ib1) .eq. 'cesm2_waccm_05x05') then
        p0 = pp1b//cmod(1:ib1)//'/Masks/sftlf6_'//cmod(1:ib1)
        p2 = pp1b//cmod(1:ib1)//'/'//cexp(1:10)//'/'//'Yz/'
        if (cexp(1:13) .eq. 'hst_ssp585.r1') p4= 'r1i1p1f1_yz_AM.nc'
        if (cexp(1:13) .eq. 'hst_ssp585.r2') p4= 'r2i1p1f1_yz_AM.nc'
        if (cexp(1:13) .eq. 'hst_ssp585.r3') p4= 'r3i1p1f1_yz_AM.nc'
c
      else if (cmod(1:ib1) .eq. 'ciesm_05x05') then
        p0 = pp1b//cmod(1:ib1)//'/Masks/sftlf6_'//cmod(1:ib1)
        p2 = pp1b//cmod(1:ib1)//'/'//cexp(1:10)//'/'//'Yz/'
        if (cexp(1:13) .eq. 'hst_ssp585.r1') p4= 'r1i1p1f1_yz_AM.nc'
c
      else if (cmod(1:ib1) .eq. 'cnrm_cm6_1_05x05') then
        p0 = pp1b//cmod(1:ib1)//'/Masks/sftlf6_'//cmod(1:ib1)
        p2 = pp1b//cmod(1:ib1)//'/'//cexp(1:10)//'/'//'Yz/'
        if (cexp(1:13) .eq. 'hst_ssp585.r1') p4= 'r1i1p1f2_yz_AM.nc'
        if (cexp(1:13) .eq. 'hst_ssp585.r2') p4= 'r2i1p1f2_yz_AM.nc'
        if (cexp(1:13) .eq. 'hst_ssp585.r3') p4= 'r3i1p1f2_yz_AM.nc'
        if (cexp(1:13) .eq. 'hst_ssp585.r4') p4= 'r4i1p1f2_yz_AM.nc'
        if (cexp(1:13) .eq. 'hst_ssp585.r5') p4= 'r5i1p1f2_yz_AM.nc'
        if (cexp(1:13) .eq. 'hst_ssp585.r6') p4= 'r6i1p1f2_yz_AM.nc'
c
      else if (cmod(1:ib1) .eq. 'ec_earth3_05x05') then
        p0 = pp1b//cmod(1:ib1)//'/Masks/sftlf6_'//cmod(1:ib1)
        p2 = pp1b//cmod(1:ib1)//'/'//cexp(1:10)//'/'//'Yz/'
        if (cexp(1:13) .eq. 'hst_ssp585.r1')  p4= 'r1i1p1f1_yz_AM.nc'
        if (cexp(1:13) .eq. 'hst_ssp585.r4')  p4= 'r4i1p1f1_yz_AM.nc'
        if (cexp(1:13) .eq. 'hst_ssp585.r6')  p4= 'r6i1p1f1_yz_AM.nc'
        if (cexp(1:14) .eq. 'hst_ssp585.r11') p4= 'r11i1p1f1_yz_AM.nc'
        if (cexp(1:14) .eq. 'hst_ssp585.r13') p4= 'r13i1p1f1_yz_AM.nc'
        if (cexp(1:14) .eq. 'hst_ssp585.r15') p4= 'r15i1p1f1_yz_AM.nc'
c
      else if (cmod(1:ib1) .eq. 'ec_earth3_veg_05x05') then
        p0 = pp1b//cmod(1:ib1)//'/Masks/sftlf6_'//cmod(1:ib1)
        p2 = pp1b//cmod(1:ib1)//'/'//cexp(1:10)//'/'//'Yz/'
        if (cexp(1:13) .eq. 'hst_ssp585.r1')  p4= 'r1i1p1f1_yz_AM.nc'
        if (cexp(1:13) .eq. 'hst_ssp585.r2')  p4= 'r2i1p1f1_yz_AM.nc'
        if (cexp(1:13) .eq. 'hst_ssp585.r3')  p4= 'r3i1p1f1_yz_AM.nc'
        if (cexp(1:13) .eq. 'hst_ssp585.r4')  p4= 'r4i1p1f1_yz_AM.nc'
        if (cexp(1:13) .eq. 'hst_ssp585.r6')  p4= 'r6i1p1f1_yz_AM.nc'
c
      else if (cmod(1:ib1) .eq. 'gfdl_cm4_05x05') then
        p0 = pp1b//cmod(1:ib1)//'/Masks/sftlf6_'//cmod(1:ib1)
        p2 = pp1b//cmod(1:ib1)//'/'//cexp(1:10)//'/'//'Yz/'
        if (cexp(1:13) .eq. 'hst_ssp585.r1')  p4= 'r1i1p1f1_yz_AM.nc'
c
      else if (cmod(1:ib1) .eq. 'gfdl_esm4_05x05') then
        p0 = pp1b//cmod(1:ib1)//'/Masks/sftlf6_'//cmod(1:ib1)
        p2 = pp1b//cmod(1:ib1)//'/'//cexp(1:10)//'/'//'Yz/'
        if (cexp(1:13) .eq. 'hst_ssp585.r1')  p4= 'r1i1p1f1_yz_AM.nc'
c
      else if (cmod(1:ib1) .eq. 'hadgem3_gc31_ll_05x05') then
        p0 = pp1b//cmod(1:ib1)//'/Masks/sftlf6_'//cmod(1:ib1)
        p2 = pp1b//cmod(1:ib1)//'/'//cexp(1:10)//'/'//'Yz/'
        if (cexp(1:13) .eq. 'hst_ssp585.r1')  p4= 'r1i1p1f3_yz_AM.nc'
        if (cexp(1:13) .eq. 'hst_ssp585.r2')  p4= 'r2i1p1f3_yz_AM.nc'
        if (cexp(1:13) .eq. 'hst_ssp585.r3')  p4= 'r3i1p1f3_yz_AM.nc'
c
      else if (cmod(1:ib1) .eq. 'hadgem3_gc31_mm_05x05') then
        p0 = pp1b//cmod(1:ib1)//'/Masks/sftlf6_'//cmod(1:ib1)
        p2 = pp1b//cmod(1:ib1)//'/'//cexp(1:10)//'/'//'Yz/'
        if (cexp(1:13) .eq. 'hst_ssp585.r1')  p4= 'r1i1p1f3_yz_AM.nc'
        if (cexp(1:13) .eq. 'hst_ssp585.r2')  p4= 'r2i1p1f3_yz_AM.nc'
        if (cexp(1:13) .eq. 'hst_ssp585.r3')  p4= 'r3i1p1f3_yz_AM.nc'
c
      else if (cmod(1:ib1) .eq. 'ipsl_cm6a_lr_05x05') then
        p0 = pp1b//cmod(1:ib1)//'/Masks/sftlf6_'//cmod(1:ib1)
        p2 = pp1b//cmod(1:ib1)//'/'//cexp(1:10)//'/'//'Yz/'
        if (cexp(1:13) .eq. 'hst_ssp585.r1')  p4= 'r1i1p1f1_yz_AM.nc'
        if (cexp(1:13) .eq. 'hst_ssp585.r2')  p4= 'r2i1p1f1_yz_AM.nc'
        if (cexp(1:13) .eq. 'hst_ssp585.r3')  p4= 'r3i1p1f1_yz_AM.nc'
        if (cexp(1:13) .eq. 'hst_ssp585.r4')  p4= 'r4i1p1f1_yz_AM.nc'
        if (cexp(1:13) .eq. 'hst_ssp585.r6')  p4= 'r6i1p1f1_yz_AM.nc'
        if (cexp(1:14) .eq. 'hst_ssp585.r14') p4= 'r14i1p1f1_yz_AM.nc'
c
      else if (cmod(1:ib1) .eq. 'miroc6_05x05') then
        p0 = pp1b//cmod(1:ib1)//'/Masks/sftlf6_'//cmod(1:ib1)
        p2 = pp1b//cmod(1:ib1)//'/'//cexp(1:10)//'/'//'Yz/'
        if (cexp(1:13) .eq. 'hst_ssp585.r1')  p4= 'r1i1p1f1_yz_AM.nc'
        if (cexp(1:13) .eq. 'hst_ssp585.r2')  p4= 'r2i1p1f1_yz_AM.nc'
        if (cexp(1:13) .eq. 'hst_ssp585.r3')  p4= 'r3i1p1f1_yz_AM.nc'
        if (cexp(1:13) .eq. 'hst_ssp585.r4')  p4= 'r4i1p1f1_yz_AM.nc'
        if (cexp(1:13) .eq. 'hst_ssp585.r5')  p4= 'r5i1p1f1_yz_AM.nc'
        if (cexp(1:13) .eq. 'hst_ssp585.r6')  p4= 'r6i1p1f1_yz_AM.nc'
        if (cexp(1:13) .eq. 'hst_ssp585.r7')  p4= 'r7i1p1f1_yz_AM.nc'
        if (cexp(1:13) .eq. 'hst_ssp585.r8')  p4= 'r8i1p1f1_yz_AM.nc'
        if (cexp(1:13) .eq. 'hst_ssp585.r9')  p4= 'r9i1p1f1_yz_AM.nc'
        if (cexp(1:14) .eq. 'hst_ssp585.r10') p4= 'r10i1p1f1_yz_AM.nc'
        if (cexp(1:14) .eq. 'hst_ssp585.r11') p4= 'r11i1p1f1_yz_AM.nc'
        if (cexp(1:14) .eq. 'hst_ssp585.r12') p4= 'r12i1p1f1_yz_AM.nc'
        if (cexp(1:14) .eq. 'hst_ssp585.r13') p4= 'r13i1p1f1_yz_AM.nc'
        if (cexp(1:14) .eq. 'hst_ssp585.r14') p4= 'r14i1p1f1_yz_AM.nc'
        if (cexp(1:14) .eq. 'hst_ssp585.r15') p4= 'r15i1p1f1_yz_AM.nc'
        if (cexp(1:14) .eq. 'hst_ssp585.r16') p4= 'r16i1p1f1_yz_AM.nc'
        if (cexp(1:14) .eq. 'hst_ssp585.r17') p4= 'r17i1p1f1_yz_AM.nc'
        if (cexp(1:14) .eq. 'hst_ssp585.r18') p4= 'r18i1p1f1_yz_AM.nc'
        if (cexp(1:14) .eq. 'hst_ssp585.r19') p4= 'r19i1p1f1_yz_AM.nc'
        if (cexp(1:14) .eq. 'hst_ssp585.r20') p4= 'r20i1p1f1_yz_AM.nc'
        if (cexp(1:14) .eq. 'hst_ssp585.r21') p4= 'r21i1p1f1_yz_AM.nc'
        if (cexp(1:14) .eq. 'hst_ssp585.r22') p4= 'r22i1p1f1_yz_AM.nc'
        if (cexp(1:14) .eq. 'hst_ssp585.r23') p4= 'r23i1p1f1_yz_AM.nc'
        if (cexp(1:14) .eq. 'hst_ssp585.r24') p4= 'r24i1p1f1_yz_AM.nc'
        if (cexp(1:14) .eq. 'hst_ssp585.r25') p4= 'r25i1p1f1_yz_AM.nc'
        if (cexp(1:14) .eq. 'hst_ssp585.r26') p4= 'r26i1p1f1_yz_AM.nc'
        if (cexp(1:14) .eq. 'hst_ssp585.r27') p4= 'r27i1p1f1_yz_AM.nc'
        if (cexp(1:14) .eq. 'hst_ssp585.r28') p4= 'r28i1p1f1_yz_AM.nc'
        if (cexp(1:14) .eq. 'hst_ssp585.r29') p4= 'r29i1p1f1_yz_AM.nc'
        if (cexp(1:14) .eq. 'hst_ssp585.r30') p4= 'r30i1p1f1_yz_AM.nc'
        if (cexp(1:14) .eq. 'hst_ssp585.r31') p4= 'r31i1p1f1_yz_AM.nc'
        if (cexp(1:14) .eq. 'hst_ssp585.r32') p4= 'r32i1p1f1_yz_AM.nc'
        if (cexp(1:14) .eq. 'hst_ssp585.r33') p4= 'r33i1p1f1_yz_AM.nc'
        if (cexp(1:14) .eq. 'hst_ssp585.r34') p4= 'r34i1p1f1_yz_AM.nc'
        if (cexp(1:14) .eq. 'hst_ssp585.r35') p4= 'r35i1p1f1_yz_AM.nc'
        if (cexp(1:14) .eq. 'hst_ssp585.r36') p4= 'r36i1p1f1_yz_AM.nc'
        if (cexp(1:14) .eq. 'hst_ssp585.r37') p4= 'r37i1p1f1_yz_AM.nc'
        if (cexp(1:14) .eq. 'hst_ssp585.r38') p4= 'r38i1p1f1_yz_AM.nc'
        if (cexp(1:14) .eq. 'hst_ssp585.r39') p4= 'r39i1p1f1_yz_AM.nc'
        if (cexp(1:14) .eq. 'hst_ssp585.r40') p4= 'r40i1p1f1_yz_AM.nc'
        if (cexp(1:14) .eq. 'hst_ssp585.r41') p4= 'r41i1p1f1_yz_AM.nc'
        if (cexp(1:14) .eq. 'hst_ssp585.r42') p4= 'r42i1p1f1_yz_AM.nc'
        if (cexp(1:14) .eq. 'hst_ssp585.r43') p4= 'r43i1p1f1_yz_AM.nc'
        if (cexp(1:14) .eq. 'hst_ssp585.r44') p4= 'r44i1p1f1_yz_AM.nc'
        if (cexp(1:14) .eq. 'hst_ssp585.r45') p4= 'r45i1p1f1_yz_AM.nc'
        if (cexp(1:14) .eq. 'hst_ssp585.r46') p4= 'r46i1p1f1_yz_AM.nc'
        if (cexp(1:14) .eq. 'hst_ssp585.r47') p4= 'r47i1p1f1_yz_AM.nc'
        if (cexp(1:14) .eq. 'hst_ssp585.r48') p4= 'r48i1p1f1_yz_AM.nc'
        if (cexp(1:14) .eq. 'hst_ssp585.r49') p4= 'r49i1p1f1_yz_AM.nc'
        if (cexp(1:14) .eq. 'hst_ssp585.r50') p4= 'r50i1p1f1_yz_AM.nc'
c
      else if (cmod(1:ib1) .eq. 'miroc_es2l_05x05') then
        p0 = pp1b//cmod(1:ib1)//'/Masks/sftlf6_'//cmod(1:ib1)
        p2 = pp1b//cmod(1:ib1)//'/'//cexp(1:10)//'/'//'Yz/'
        if (cexp(1:13) .eq. 'hst_ssp585.r1')  p4= 'r1i1p1f2_yz_AM.nc'
c
      else if (cmod(1:ib1) .eq. 'mpi_esm1_2_hr_05x05') then
        p0 = pp1b//cmod(1:ib1)//'/Masks/sftlf6_'//cmod(1:ib1)
        p2 = pp1b//cmod(1:ib1)//'/'//cexp(1:10)//'/'//'Yz/'
        if (cexp(1:13) .eq. 'hst_ssp585.r1')  p4= 'r1i1p1f1_yz_AM.nc'
        if (cexp(1:13) .eq. 'hst_ssp585.r2')  p4= 'r2i1p1f1_yz_AM.nc'
c
      else if (cmod(1:ib1) .eq. 'mpi_esm1_2_lr_05x05') then
        p0 = pp1b//cmod(1:ib1)//'/Masks/sftlf6_'//cmod(1:ib1)
        p2 = pp1b//cmod(1:ib1)//'/'//cexp(1:10)//'/'//'Yz/'
        if (cexp(1:13) .eq. 'hst_ssp585.r1')  p4= 'r1i1p1f1_yz_AM.nc'
        if (cexp(1:13) .eq. 'hst_ssp585.r2')  p4= 'r2i1p1f1_yz_AM.nc'
        if (cexp(1:13) .eq. 'hst_ssp585.r3')  p4= 'r3i1p1f1_yz_AM.nc'
        if (cexp(1:13) .eq. 'hst_ssp585.r4')  p4= 'r4i1p1f1_yz_AM.nc'
        if (cexp(1:13) .eq. 'hst_ssp585.r5')  p4= 'r5i1p1f1_yz_AM.nc'
        if (cexp(1:13) .eq. 'hst_ssp585.r6')  p4= 'r6i1p1f1_yz_AM.nc'
        if (cexp(1:13) .eq. 'hst_ssp585.r7')  p4= 'r7i1p1f1_yz_AM.nc'
        if (cexp(1:13) .eq. 'hst_ssp585.r8')  p4= 'r8i1p1f1_yz_AM.nc'
        if (cexp(1:13) .eq. 'hst_ssp585.r9')  p4= 'r9i1p1f1_yz_AM.nc'
        if (cexp(1:13) .eq. 'hst_ssp585.r10') p4= 'r10i1p1f1_yz_AM.nc'
c
      else if (cmod(1:ib1) .eq. 'mri_esm2_0_05x05') then
        p0 = pp1b//cmod(1:ib1)//'/Masks/sftlf6_'//cmod(1:ib1)
        p2 = pp1b//cmod(1:ib1)//'/'//cexp(1:10)//'/'//'Yz/'
        if (cexp(1:13) .eq. 'hst_ssp585.r1')  p4= 'r1i1p1f1_yz_AM.nc'
c
      else if (cmod(1:ib1) .eq. 'nesm3_05x05') then
        p0 = pp1b//cmod(1:ib1)//'/Masks/sftlf6_'//cmod(1:ib1)
        p2 = pp1b//cmod(1:ib1)//'/'//cexp(1:10)//'/'//'Yz/'
        if (cexp(1:13) .eq. 'hst_ssp585.r1')  p4= 'r1i1p1f1_yz_AM.nc'
        if (cexp(1:13) .eq. 'hst_ssp585.r2')  p4= 'r2i1p1f1_yz_AM.nc'
c
      else if (cmod(1:ib1) .eq. 'noresm2_mm_05x05') then
        p0 = pp1b//cmod(1:ib1)//'/Masks/sftlf6_'//cmod(1:ib1)
        p2 = pp1b//cmod(1:ib1)//'/'//cexp(1:10)//'/'//'Yz/'
        if (cexp(1:13) .eq. 'hst_ssp585.r1')  p4= 'r1i1p1f1_yz_AM.nc'
c
      else if (cmod(1:ib1) .eq. 'ukesm1_0_ll_05x05') then
        p0 = pp1b//cmod(1:ib1)//'/Masks/sftlf6_'//cmod(1:ib1)
        p2 = pp1b//cmod(1:ib1)//'/'//cexp(1:10)//'/'//'Yz/'
        if (cexp(1:13) .eq. 'hst_ssp585.r1')  p4= 'r1i1p1f2_yz_AM.nc'
        if (cexp(1:13) .eq. 'hst_ssp585.r2')  p4= 'r2i1p1f2_yz_AM.nc'
        if (cexp(1:13) .eq. 'hst_ssp585.r3')  p4= 'r3i1p1f2_yz_AM.nc'
        if (cexp(1:13) .eq. 'hst_ssp585.r4')  p4= 'r4i1p1f2_yz_AM.nc'
        if (cexp(1:13) .eq. 'hst_ssp585.r8')  p4= 'r8i1p1f2_yz_AM.nc'
c
c
c
      end if
c
c
c
c-------------------------------------------------------------------------------
c     ** 1.6 ** 
c
c     Based on selections, set overall file name.
c-------------------------------------------------------------------------------
c
      if (cexp(1:10) .eq. 'hist_rcp85') then
        if (cmod(1:3) .eq. 'MMA') then
          p3 = varnam(1:ib3)//'_'//cmod(1:ib1)//'_'//cexp(1:10)//'_'
        else
          p3 = varnam(1:ib3)//'_'//cmod(1:ib1)//'_'//cexp(1:10)//'.'
        end if
      else if (cexp(1:10) .eq. 'hist_LE_r1') then
        p3 = varnam(1:ib3)//'_'//cmod(1:ib1)//'_'//cexp(1:10)//'.'
      else if (cexp(1:10) .eq. 'hist_LE_r2') then
        p3 = varnam(1:ib3)//'_'//cmod(1:ib1)//'_'//cexp(1:10)//'.'
      else if (cexp(1:10) .eq. 'hist_LE_r3') then
        p3 = varnam(1:ib3)//'_'//cmod(1:ib1)//'_'//cexp(1:10)//'.'
      else if (cexp(1:10) .eq. 'hist_LE_r4') then
        p3 = varnam(1:ib3)//'_'//cmod(1:ib1)//'_'//cexp(1:10)//'.'
      else if (cexp(1:10) .eq. 'hist_LE_r5') then
        p3 = varnam(1:ib3)//'_'//cmod(1:ib1)//'_'//cexp(1:10)//'.'
      else if (cexp(1:10) .eq. 'hst_ssp585') then
        if (cmod(1:3) .eq. 'MMA') then
          p3 = varnam(1:ib3)//'_'//cmod(1:ib1)//'_'//cexp(1:10)//'_'
        else
          p3 = varnam(1:ib3)//'_'//cmod(1:ib1)//'_'//cexp(1:10)//'.'
        end if
      else if (cexp(1:10) .eq. 'hst_ssp370') then
        p3 = varnam(1:ib3)//'_'//cmod(1:ib1)//'_'//cexp(1:10)//'.'
      else if (cexp(1:9) .eq. 'piControl') then
        p3 = varnam(1:ib3)//'_'//cmod(1:ib1)//'_'//cexp(1:9)//'_'
      else if (cexp(1:10) .eq. '20c3m_run1') then
        p3 = varnam(1:ib3)//'_'//cmod(1:ib1)//'_'//cexp(1:10)//'_'
      end if
c
      print *, 'In CREATE_PATH1! p2 = ', p2
      print *, 'In CREATE_PATH1! p3 = ', p3
      print *, 'In CREATE_PATH1! p4 = ', p4
c
      ib4 = index(p2, ' ')
      ib4 = ib4 - 1
      ib5 = index(p3, ' ')
      ib5 = ib5 - 1
      ib6 = index(p4, ' ')
      ib6 = ib6 - 1
c
      if (ib4 .lt. 0) then
        print *, 'In CREATE_PATH1!'
        print *, 'No blanks in p2!'
      end if
      if (ib5 .lt. 0) then
        print *, 'In CREATE_PATH1!'
        print *, 'No blanks in p3!'
      end if
      if (ib6 .lt. 0) then
        print *, 'In CREATE_PATH1!'
        print *, 'No blanks in p4!'
      end if
c
      fpath1 = p2(1:ib4)//p3(1:ib5)//p4(1:ib6)
      fpath2 = p0                              
c
      print *, 'In CREATE_PATH1! fpath1 = ', fpath1
      print *, 'In CREATE_PATH1! fpath2 = ', fpath2
c
      return
      end
c-------------------------------------------------------------------------------
c-------------------------------------------------------------------------------
c     SUBROUTINE RD3D_YHT_NOMASK 
c
c     Date:    23-OCT-2024.
c     Author:  Ben Santer, LLNL, Livermore.
c     Version: # 1.0 / Single precision.
c
c
c
c     This subroutine reads an input three-dimensional latitude-pressure-time
c     field. There is no masking of input data. It is assumed that the area
c     weights do not change as a function of time in years.
c
      subroutine rd3d_yht_nomask(iopen, idomain, ny, nh, nt, t1, t2, 
     &infil, rlat, rlev, rtime, amask, field1, varnam, yweights, 
     &iyout, ihout, itout, y1, y2)       
c
c
c
      implicit real (o-z, a-h)
c
c     ** INPUT **
c     integer iopen              1 = initialize EzGet, 0 = no initialization
c     integer idomain            Spatial domain extracted         
c     integer ny                 Maximum spatial dimension (latitude)
c     integer nh                 Maximum spatial dimension (pressure) 
c     integer nt                 Maximum time dimension (years)
c     real t1, t2                Selected time range
c     character*(*) infil        Input file name
c     character*16 varnam        Variable name
c     character* 8 yweights      Weighting type
c
c     ** OUTPUT **
c     integer iyout                   Retrieved latitude dimension
c     integer ihout                   Retrieved pressure dimension
c     integer itout                   Retrieved time dimension
c     real y1, y2                     Retrieved latitude boundaries
c     dimension field1(ny, nh, nt)    Output data
c     dimension amask(ny, nh, nt)     Area weights for data
c     dimension rlat(ny)              Latitude coordinates
c     dimension rlev(nh)              Pressure coordinates
c     dimension rtime(nt)             Time coordinates
c
c
c
      integer iopen, ny, nh, nt, iyout, ihout, itout
      integer i1, i2, i3, i4
      real t1, t2, y1, y2
      dimension field1(ny, nh, nt)
      dimension amask(ny, nh, nt)  
      dimension rlat(ny), rlev(nh), rtime(nt)
      character infil*(*), varnam*16, yweights*8
c
c
c
c-------------------------------------------------------------------------------
c     ** 1.0 **
c
c     Define domain extracted. 
c-------------------------------------------------------------------------------
c
      h1 =    1.0
      h2 =  695.0
c
      print *, 'idomain = ', idomain
c
      if (idomain .eq. 1) then
c     Domain retrieved: 90S-90N
        y1 = -90.0
        y2 =  90.0
c
      else if (idomain .eq. 2) then
c     Domain retrieved: 75S-75N
        y1 = -75.0
        y2 =  75.0
c
      else if (idomain .eq. 3) then
c     Domain retrieved: 60S-90N
        y1 = -60.0
        y2 =  90.0
c
      else if (idomain .eq. 4) then
c     Domain retrieved: 60S-60N
        y1 = -60.0
        y2 =  60.0
c
      else if (idomain .eq. 5) then
c     Domain retrieved: 50S-50N
        y1 = -50.0
        y2 =  50.0
c
      else if (idomain .eq. 6) then
c     Domain retrieved: 80S-80N
        y1 = -80.0
        y2 =  80.0
c
      else if (idomain .eq. 7) then
c     Domain retrieved: 70S-80N
        y1 = -70.0
        y2 =  80.0
c
      else if (idomain .eq. 8) then
c     Domain retrieved: 40N-55N
        y1 =  40.0
        y2 =  55.0
c
      else if (idomain .eq. 9) then
c     Domain retrieved: 0-90N
        y1 =   0.0
        y2 =  90.0
c
      else if (idomain .eq. 10) then
c     Domain retrieved: 90S-0
        y1 = -90.0
        y2 =   0.0
c
      end if
c
      if (iopen .eq. 1) call initget
c
      call defdim(1, 1, 'latitude', yweights, 'range', y1, y2, 0.0)
c
      call defdim(1, 2, 'pressure', 'unit', 'assaved', h1, h2, 0.0)
c
      call defdim(1, 3, 'time', 'unit', 'nearest', t1, t2, 0.0) 
      call defmisc('maximum dimension kept', 'integer', 9000)
c
c-------------------------------------------------------------------------------
c     ** 1.1 **
c
c     Define expected dimensions of data.
c-------------------------------------------------------------------------------
c
      i1 = 0 
      i2 = 0  
      i3 = 0  
      i4 = 0  
c
c-------------------------------------------------------------------------------
c     ** 1.2 **
c
c     Define variable 1.
c-------------------------------------------------------------------------------
c
      call defvar(1, varnam, infil)
c
c-------------------------------------------------------------------------------
c     ** 1.3 **
c
c     Extract data.
c-------------------------------------------------------------------------------
c
      call getdata(1, ny, nh, nt, 0, i1, i2, i3, i4, amask, field1)
      call getcoord(1, 1, rlat)
      call getcoord(1, 2, rlev)
      call getcoord(1, 3, rtime)  
c
      iyout = i1
      ihout = i2
      itout = i3
c
      print *, 'field1(1,1,1) = ', field1(1,1,1)
      print *, 'rlat(1)  = ', rlat(1) 
      print *, 'rlat(i1) = ', rlat(i1) 
c
c-------------------------------------------------------------------------------
c     ** 1.4 **
c
c     Release memory.
c-------------------------------------------------------------------------------
c
      if (iopen .eq. 1) call closeget 
c
      return
      end
c-------------------------------------------------------------------------------
c-------------------------------------------------------------------------------
c     SUBROUTINE RDMASK.
c
c     Date:    23-OCT-2024.
c     Author:  Ben Santer and Karl Taylor, LLNL, Livermore.
c     Version: # 1.0 / Single precision.
c
c
c
c     This subroutine reads an input two-dimensional latitude-pressure land area
c     fraction mask. Data are assumed to available in the form of zonal means
c     (at 5 degree latitudinal resolution) for individual atmospheric layers.
c
c
c
      subroutine rdmask(ny, nh, nyi, nhi, idomain, infil, amask, 
     &fmask, yweights)
c
c
c
      implicit real (o-z, a-h)
c
c     ** INPUT **
c     integer ny                  Spatial dimension (latitude)
c     integer nh                  Spatial dimension (pressure)
c     integer idomain             Spatial domain extracted    
c     character*(*)infil          Input file name
c
c     ** OUTPUT **
c     integer nyi                 Output spatial dimension (latitude)
c     integer nhi                 Output spatial dimension (pressure)
c     dimension amask(ny, nh)     Area weights
c     dimension fmask(ny, nh)     Fractional data coverage mask
c
c
c
      integer ny, nh, nyi, nhi, idomain
      dimension fmask(ny, nh), amask(ny, nh)
      character infil*(*), yweights*8
c
c
c
c-------------------------------------------------------------------------------
c     ** 1.0 **
c
c     Initialize EZGET.
c-------------------------------------------------------------------------------
c
      call initget
c
c-------------------------------------------------------------------------------
c     ** 1.1 **
c
c     Define domain extracted.                           
c-------------------------------------------------------------------------------
c
      h1 =    1.0
      h2 =  695.0
c
      if (idomain .eq. 1) then
c     Domain retrieved: 90S-90N
        y1 = -90.0
        y2 =  90.0
c
      else if (idomain .eq. 2) then
c     Domain retrieved: 75S-75N
        y1 = -75.0
        y2 =  75.0
c
      else if (idomain .eq. 3) then
c     Domain retrieved: 60S-90N
        y1 = -60.0
        y2 =  90.0
c
      else if (idomain .eq. 4) then
c     Domain retrieved: 60S-60N
        y1 = -60.0
        y2 =  60.0
c
      else if (idomain .eq. 5) then
c     Domain retrieved: 50S-50N
        y1 = -50.0
        y2 =  50.0
c
      else if (idomain .eq. 6) then
c     Domain retrieved: 80S-80N
        y1 = -80.0
        y2 =  80.0
c
      else if (idomain .eq. 7) then
c     Domain retrieved: 70S-80N
        y1 = -70.0
        y2 =  80.0
c
      else if (idomain .eq. 8) then
c     Domain retrieved: 40N-55N
        y1 =  40.0
        y2 =  55.0
c
      else if (idomain .eq. 9) then
c     Domain retrieved: 0-90N
        y1 =   0.0
        y2 =  90.0
c
      else if (idomain .eq. 10) then
c     Domain retrieved: 90S-0
        y1 = -90.0
        y2 =   0.0
c
      end if
c
      call defdim(1, 1, 'latitude', yweights, 'range', y1, y2, 0.0)
c
      call defdim(1, 2, 'pressure', 'unit', 'assaved', h1, h2, 0.0)
c
c
c-------------------------------------------------------------------------------
c     ** 1.2 **
c
c     Define expected dimensions of data.
c-------------------------------------------------------------------------------
c
      i1 = 0
      i2 = 0
      i3 = 0
      i4 = 0
      n3 = 0
      n4 = 0
c
c-------------------------------------------------------------------------------
c     ** 1.3 **
c
c     Define variable 1. 
c-------------------------------------------------------------------------------
c
      call defvar(1, 'sftlf6', infil)
c
c-------------------------------------------------------------------------------
c     ** 1.4 **
c
c       Extract data.
c-------------------------------------------------------------------------------
c
      call getdata(1, ny, nh, n3, n4, i1, i2, i3, i4, amask, fmask)
      nyi = i1
      nhi = i2
c
c-------------------------------------------------------------------------------
c     ** 1.5 **
c
c     Release memory.
c-------------------------------------------------------------------------------
c
      call closeget
c
      return
      end
