implicit none integer nlon, nlat, mlat, niofiles character*16 abbrev, abbr parameter (niofiles=1) character*120 pathout, inputsic(niofiles), inputsst(niofiles), & outpath, insst, insic C******************* SPECIFY FOLLOWING PARAMETERS ******************* c nlon = number of longitudes in i/o grid c nlat = number of latitudes in i/o grid c mlat = number of latitudes at a time (i.e. latitudes per chunk) c abbr = a filenaming tag based on the number of grid points c outpath = directory in which to place output files c insst = full file specification of input SST c insic = full file specification of input sea ice c parameter (nlon=72, & nlat=36, & mlat=36, & abbr='72x36', & outpath='/clim_var/hc2100/OBS/marine/HadISST/karl', & insst='/clim_var/hc2100/OBS/marine/HadISST/'// & 'actual/HadISST1.1_5d_sst_18701999.pp', & insic='/clim_var/hc2100/OBS/marine/HadISST/'// & 'ice/HadISST1_ice_18701999_5d.pp') C******************** DO NOT EDIT BELOW THIS LINE ******************** c to create an artificial mid-month sea ice fraction or SST data set c that, when linearly interpolated (in time), produces the observed c monthly means c 21 August 1997 c Karl E. Taylor c PCMDI c taylor13@llnl.gov integer nmon, nzon, nzonp, n1, n2, nmon12, & nlagm, nchunks parameter (nmon=1560, & nzon=1, nzonp=nzon+1, n1=nzonp*21, n2=nzonp*41, & nmon12=12, nlagm=13, nchunks=((nlat-1)/mlat+1)) logical caseindp integer i, j, icntl, maxiter integer lagcalc, iie, ie integer m, mp, n, mm, nn, i1, kk, isrc, isrc1, nfilesi, nfileso integer ismc, jsmc, ismf, jsmf, ik, isea, je, k, ig, j1, jchnk, jn integer mons, m1, m2, mon1spin, iyr1spin, nmonspin real conv, r, wtmin, bbmin, dcbias, fac, tmin, offset, & tmax, histo, varmin, dt, amissin, amissout, amissasc real obsclim(nlon,nlat,nmon12), sst(nlon,mlat,nmon), alons(nlon), & alats(nlat), centmon(nlon,nlat,nmon12), gauss(nlat), ss(nmon) real omax(nlon,nlat), omin(nlon,nlat), cmax(nlon,nlat), & cmin(nlon,nlat), ovarmon(nlon,nlat), cvarmon(nlon,nlat), & sstwts(nlon,nlat), wtfrac(nlon,nlat), obsmean(nmon), & space(nlon,nlat,nmon12), spinup(nlon,nlat,nmon12) real woseas(nzon), wcseas(nzon), wovar(nzonp), wcvar(nzonp), & wocorr(nzonp), wccorr(nzonp), & acvarmon(nzonp), aovarmon(nzonp), aobsclim(nmon12,nzon), & acentmon(nmon12,nzon), bpp(19) real ocorrel(nlagm,nzonp), ccorrel(nlagm,nzonp), correl(nmon12) real elem1(2), elemn(2), cenmon(nmon12,2), vecout(nmon), & vecin(nmon), a(nmon), c(nmon), ac(nmon12), cc(nmon12) integer icount(nzonp,-20:20), jcount(nzonp,-10:10), & kcount(nzonp,-20:20), mcount(nzonp,-10:10), lagm(nlagm) integer monlen(nmon12,2), lpp(45), jyr1out(niofiles), kyr1out(3), & iyr1sst(niofiles), iyr1sic(niofiles), iyr1in(niofiles) character*80 title, grid character*120 fbcinfo, fbc, sftfileo, sftfilei, dum(nchunks), & inputfil(niofiles), outfile, tempbc(nchunks), & tempobs(nchunks), fclimbc, fclimobs, fspinup, fobs, outcopy character*16 vname, varin, model, sftnameo, sftnamei character varout*6, varout1*3 character*40 units, units2 character*120 sourceo, sourceb character*30 suff character*3 suffix, c3 character*4 ayr character*16 src, src1 double precision sum character*16 outftype, calclim, calgreg, gtype, typfil, inftype character*80 center character*256 parmtabl, dfltparm integer idvar, ibasemon, ibaseyr, idvara(6), lbfc, iend integer iout(6) integer lats, mon1rd, iyr1rd, monnrd, iyrnrd, mon1clm, iyr1clm, & monnclm, iyrnclm, mon1, iyr1, monn, iyrn, mon1out, iyr1out, & monnout, iyrnout, mon1in, msklndi, msklndo data lagm/0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12/ data ig/1/ data jyr1out/niofiles*-999/, iyr1sst/niofiles*-999/, & iyr1sic/niofiles*-999/, kyr1out/3*-999/ data lpp/45*0/, bpp/19*0.0/ c data monlen/31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31, c & 31, 29, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31/ c *** HadISST *** might assume months of equal length data monlen/30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, & 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30/ c ***************************************************************** inputsst(1) = insst inputsic(1) = insic abbrev = abbr pathout = outpath print*,'HELLO' c SPECIFY controlling parameters: c icntl = 1 for SST, 2 for sic, and 3 for SST and sic icntl = 3 c SPECIFY output wanted: 1=output 0=no output c iout(1) = boundary conditions for period specified by c mon1out, iyr1out, monnout, iyrnout c iout(2) = boundary conditions for year preceeding period c specified in 1 c iout(3) = boundary conditions based on climatology c iout(4) = observed monthly means (normally =0 if input in c pp format, unless you want to change the units) c iout(5) = climatology of monthly means c iout(6) = diagnostic statistics (normally =0 for HadISST) c iout(4) = 0 for HadISST iout(1) = 1 iout(2) = 1 iout(3) = 1 c iout(4) = 1 iout(4) = 0 iout(5) = 1 iout(6) = 1 c SPECIFY source for boundary condition data (used to control c description stored in output file) c (<16 characters) c usually set to 'AMIP 2' or 'HadISST' c isrc is the length of the significant part of the character c string src = 'HadISST' isrc = 7 c SPECIFY src1, which is an abbreviated lower case version of src c used in generating output filenames (must not contain c blanks) (< 16 characters) src1 = 'hadisst' isrc1 = index(src1, ' ') - 1 c SPECIFY either 'obs' (if not regridding input data and EzGet input c has been specified) or a brief indicator of the target grid. c (< 17 characters) c normal for HadISST is model = 'obs' model = 'obs' c SPECIFY the research center that has requested the data (written c as a comment in ascii and lats output files). For grib and c grads files, the center must appear in the parameter table c lists (see parmtabl and dfltparm). c (< 81 characters) c for HadISST centre = 'MetOffice' center = 'MetOffice' c SPECIFY grid or model i.d. for use in generating file names c (< 17 characters) c usually indicating grid. c examples: ncep, 1x1, 4x5, t42, ccm_t63, .... c HadISST: abbrev = '360x180' c abbrev = '72x36' c c SPECIFY grid type, which controls weights and for lats, controls c output grid type: c 'gaussian' or 'cosine' or standard model acronym c (gtype is used only to control weighting if if both c input and output are in pp format) c c HadISST: gtype = 'cosine' gtype = 'cosine' c SPECIFY complete grid description (for ascii output only) c for ascii files following is output to identify grid type: grid = '1x1 degree uniformly-spaced longitude/latitude grid' c SPECIFY file type for output? ('pp', 'drs', 'coards', 'cdf', c 'grads', 'grib', 'ascii', or 'all' (i.e., pp, drs, c coards, grib, and ascii) c ??? why coards and not cdf? c ??? why grib and not grads? c ??? explain here differences between different lats formats outftype = 'pp' c SPECIFY path to all output from this program c pathout = '/clim_var/hc2100/OBS/marine/HadISST/karl' c SPECIFY file type for input? ('pp' or 'ezget') c inftype = 'ezget' !AMIP c inftype = 'pp' !HadISST inftype = 'pp' if (inftype(1:2) .ne. 'pp') then c land/sea mask for output (not used when reading pp format c 0 = no mask c 1 = 0.0 to 100.0 (% land) c 3 = 0.0, 1.0, 2.0 OR 0, 1, 2 indicating land, ocean, c and sea ice, respectively c 4 = detailed geography mask sftnameo = 'sftlf' msklndo=1 C sftfileo = '/pcmdi/staff/longterm/doutriau/ldseamsk/amipII/' // C & 'pcmdi_sftlf_yonu.nc' c land/sea mask for input (see above for key) msklndi = 3 sftnamei = 'ls' C sftfilei = C & '/pcmdi/staff/longterm/doutriau/ldseamsk/bcamip2/mask.1deg.ctl' c sftnamei = 'sftl' c sftfilei = c & '/pcmdi/staff/longterm/gleckler/gleckler/ktaylor/sftl.nc' c sftfilei = '/pcmdi/share1/obs/amip2/bcs/fnl/amip2.bcs.mask.nc' c sftfilei = '/pcmdi/doutriau/geog/sealand/data/sftl_1x1.nc' else sftnameo = ' ' msklndo = 0 sftfileo = ' ' sftnamei = ' ' msklndi = 0 sftfilei = ' ' endif c SPECIFY first month and year for entire period that will be treated c (period of interest plus buffers at ends) mon1 = 1 iyr1 = 1870 c last month and year for entire period that will be treated (period c of interest plus buffers at ends). Note the entire period c treated should be an integral number of years. monn = 12 iyrn = 1999 c first month and year for period in which observed monthly mean data c will be read (must not preceed mon1, iyr1) mon1rd = 1 iyr1rd = 1871 c last month and year for period in which observed monthly mean data c will be read (must not follow monn, iyrn) monnrd = 12 iyrnrd = 1998 c first month and year that will be included in climatological mean mon1clm = 1 iyr1clm = 1871 c last month and year that will be included in climatological mean monnclm = 12 iyrnclm = 1998 c first month and year written to output file mon1out = 1 iyr1out = 1871 c last month and year written to output file monnout = 12 iyrnout = 1998 c first month found in first input file (these values are c ignored if the 'pp' input format has been selected) c data is assumed to be stored by month and if non-pp c format, data is extracted by time-index (not actual c time-coordinate value). c mon1in = 12 !AMIP c mon1in = 0 !HadISST mon1in = 0 c SPECIFY path to SST concentration data (not used if processing c only sic data) Also if more than 1 SST input file, then c specify year of first month of data in each file. c (Note, last month in each file, except possibly last c file, should be December) c inputsst(1) = c & '/clim_var/hc2100/OBS/marine/HadISST/actual/' // c & 'HadISST1.1_5d_sst_18701999.pp' iyr1sst(1) = 1870 c SPECIFY path to sea ice concentration data (not used if processing c only sst data) Also if more than 1 sic input file, then c specify year of first month of data in each file. c (Note, last month in each file, except possibly last c file, should be December) c inputsic(1) = c & '/clim_var/hc2100/OBS/marine/HadISST/ice/' // c & 'HadISST1_ice_18701999_5d.pp' iyr1sic(1) = 1870 c SPECIFY year of first month of boundary condition and monthly mean c data stored in each output file (i.e., fbc and fobs files) c It is not necessary to define jyr1out(1), because this must be c set equal to iyr1out (defined above). c If you want all data placed in a single file, it is not c necessary to define jyr1out c If you want each year to be written to separate files, then set c (jyr1out(i), i=2,number of years) equal to consecutive years c (e.g., 1975, 1976, 1977) c If you want each decade of data to be written to separate files c then set (jyr1out(i), i=2,number of decades) equal to years c separated by increments of decades (e.g., 1942, 1952, 1962) c note: jyr1out(2) - jyr1out(1) can differ from the other c increments C jyr1out(1) = iyr1out c do 26 i=2,19 !AMIP c jyr1out(i) = jyr1out(1) + i - 1 c 26 continue C commented out because want all output to go to one file C do 26 i=2,13 !GISST/HadISST C jyr1out(i) = 1948 + 10*(i - 1) C 26 continue c ***************************************************************** c diagnostics: write(6,*)'icntl',icntl write(6,*)'iout',iout write(6,*)'src',src write(6,*)'isrc',isrc write(6,*)'src1',src1 write(6,*)'isrc1',isrc1 write(6,*)'model',model write(6,*)'center',center write(6,*)'abbrev',abbrev write(6,*)'gtype',gtype write(6,*)'grid',grid write(6,*)'outftype',outftype write(6,*)'pathout',pathout write(6,*)'inftype',inftype write(6,*)' ' write(6,*)' ' write(6,*)' ' c ***************************************************************** c what fraction of a grid cell area must be unmasked if it is to be c included in area mean calculations? wtmin = 0.01 c parameters controlling nonlinear solver: c maximum number of iterations? maxiter = 100 c smallest diagnal element allowed in Jacobian? bbmin = 0.001 c calculate lag correlations? c 0 = no c 1 = calculate monthly lag correlations lagcalc = 1 c base-time? (used only for lats output) ibasemon = 1 c ibaseyr = 1978 !AMIP c ibaseyr = 1870 !HadISST ibaseyr = 1870 c if (src1(1:isrc1) .eq. 'hadisst') then c ***HadISST*** check source description of HadISST data: sourceo = 'HadISST1.1 data: Rayner et al., Hadley Centre '// & 'Clim. Res. Tech. Note xx' sourceb = 'Based on HadISST1.1 data: Rayner et al., Hadley'// & ' Centre Clim. Res. Tech. Note xx' elseif (src1(1:isrc1) .eq. 'amip') then sourceo = 'Observations on 1x1 deg grid obtained from Fiorino, ' & // '1996: http://kuala.llnl.gov/amip2/amip2_bcs.htm' sourceb = 'Based on monthly mean observations on 1x1 deg grid' & // ' obtained from Fiorino 1996:' & // ' http://kuala.llnl.gov/amip2/amip2_bcs.htm' else print*, "Error: src1 should be either 'hadisst' or 'amip')" stop endif c SPECIFY missing data indicators: c amissin : missing data value in input files c amissout : value assigned to missing data in output c files (except ascii files) c amissasc : value assigned to missing data in ascii c output files c amissin = 1.0e20 !AMIP c amissin = -1.0e30 !HadISST amissin = -1.0e30 c amissout = 1.0e20 !AMIP c amissout = -1.0e30 !HadISST amissout = -1.0e30 amissasc = 0.0 c calendar? ('climatology', 'julian', or 'gregorian') c note: 'gregorian' is only option for 'coards' file type c (parameter ignored if pp, drs, or ascii output selected) c ??? generalize for 30-day months calgreg = 'gregorian' if (typfil(1:6) .eq. 'coards') then calclim = 'gregorian' else calclim = 'climatology' endif 22 if ((icntl .eq. 2) .or. (icntl .eq. 3)) then varout = 'sicbcs' varout1 = 'sic' else varout = 'tosbcs' varout1 = 'tos' endif if (varout(1:3) .eq. 'sic') then c c *************** SEA ICE SPECS ***************** c c convergence criteria conv = 1.e-2 c constant to subtract to reduce truncation errors dcbias = 0.0 c factor to multiply input data by (before any cropping). c Output will be scaled by this factor: c fac = 100.0 !AMIP c fac = 1.0 !GISST/HadISST fac = 1.0 c number added to data (after scaling by fac) offset = 0.0 c output units for sea ice (also variance units) c units = '%' !AMIP c units2 = '%**2' !AMIP c units = ' ' !HadISST c units2 = ' ' !HadISST units = ' ' units2 = ' ' c prescribed limits (in output units) used (if necessary) to crop c data and accounted for in creating boundary condition data c tmin = 0.0 !AMIP c tmax = 100.0 !AMIP c *** GISST/HadISST *** check whether tenths are stored as character c data? c tmin = 0.0 !HadISST c tmax = 1.0 !HadISST tmin = 0.0 tmax = 1.0 c smoothing when values go from mininum to maximum c maximum jump in monthly means allowed is c tmax-tmin-dt (in output units) c dt = 2.0 c dt = (tmax-tmin)/0.5 dt = 0.0001 c varin = 'sica' suffix = 'sic' c suff = 'Sea Ice Concentration (%)' !AMIP c suff = 'Sea Ice Concentration Fraction' !HadISST suff = 'Sea Ice Concentration Fraction' c length of above string iie = 30 c histogram fine increment? c histo = (tmax-tmin)/100.0 !AMIP c histo = (tmax-tmin)/100.0 !HadISST histo = (tmax-tmin)/100.0 c minimum variance for grid cell to be included in area means varmin = 1.0 c field code for pp-format output lbfc = 37 c lag correlations for monthly mean anomalies to be used in c generating artificial data if (model(1:3) .eq. 'obs') then c following are from AMIP II data on original 1x1 degree grid correl(1) = 0.53 correl(2) = 0.23 correl(3) = 0.13 correl(4) = 0.08 correl(5) = 0.05 c following are NOT taken from observations correl(6) = 0.02 correl(7) = 0.00 correl(8) = 0.00 correl(9) = 0.00 correl(10) = 0.00 correl(11) = 0.00 correl(12) = 0.00 c observed values: c correl(6) = 0.05 c correl(7) = 0.05 c correl(8) = 0.06 c correl(9) = 0.06 c correl(10) = 0.08 c correl(11) = 0.09 c correl(12) = 0.09 else c following are from AMIP II observations regridded to 3x3 c degree resolution correl(1) = 0.58 correl(2) = 0.27 correl(3) = 0.15 correl(4) = 0.09 correl(5) = 0.06 c following are NOT taken from observations correl(6) = 0.03 correl(7) = 0.00 correl(8) = 0.00 correl(9) = 0.00 correl(10) = 0.00 correl(11) = 0.00 correl(12) = 0.00 c observed values: c correl(6) = 0.05 c correl(7) = 0.05 c correl(8) = 0.06 c correl(9) = 0.07 c correl(10) = 0.08 c correl(11) = 0.09 c correl(12) = 0.08 endif do 31 i=1,niofiles inputfil(i) = inputsic(i) iyr1in(i) = iyr1sic(i) 31 continue c if output is pp format, but input is ezget, c then generate pp header c if ((outftype(1:2) .eq. 'pp') .and. c & (inftype(1:2) .ne. 'pp')) then c print*, 'at present input file type must be ' c print*, 'in pp format if output is in pp format' c stop c call genpphd c endif else c c ********************* SST specs ************************ c c convergence criteria conv = 2.5e-2 c freezing point in output units c dcbias = 273.16 !AMIP c dcbias = 0.0 !HadISST dcbias = 0.0 c factor to multiply input data by (before any cropping). c Output will be scaled by this factor: c fac = 1.0 !AMIP & GISST/HadISST fac = 1.0 c number added to data (after scaling by fac) c to convert from Celsius to Kelvin, set offset = 273.16 c offset = 0.0 !AMIP c offset = 273.16 !HadISST offset = 273.16 c output units for sea ice (also variance units) c ***GISST/HadISST**** check units c units for sst (also variance units) c units = 'K' !AMIP c units2 = 'K**2' !AMIP c units = 'K' !HadISST c units2 = 'K**2' !HadISST units = 'K' units2 = 'K**2' c prescribed limits (in output units) used (if necessary) to crop c data and accounted for in creating boundary condition data c c tmin = 271.38 !AMIP c tmax = 1000.0 !AMIP c tmin = 271.36 !HadISST (-1.8C in K) c tmax = 1000.0 !HadISST tmin = 271.36 tmax = 1000.0 c smoothing when values go from mininum to maximum c maximum jump in monthly means allowed is c tmax-tmin-dt (in output units) dt = 0.001 varin = 'sst' suffix = 'sst' c suff = 'SST (K)' !AMIP c suff = 'SST (K)' !HadISST suff = 'SST (K)' iie = 7 c histogram fine increment? histo = 0.1 c minimum variance for grid cell to be included in area means varmin = 0.0 c field code for pp-format output lbfc = 16 c lag correlations for monthly mean anomalies to be used in c generating artificial data if (model(1:3) .eq. 'obs') then c following are from AMIP II data on original 1x1 degree grid correl(1) = 0.68 correl(2) = 0.46 correl(3) = 0.33 correl(4) = 0.26 correl(5) = 0.20 correl(6) = 0.17 correl(7) = 0.14 correl(8) = 0.11 c following are NOT taken from observations correl(9) = 0.08 correl(10) = 0.05 correl(11) = 0.02 correl(12) = 0.00 c observed values: c correl(9) = 0.09 c correl(10) = 0.07 c correl(11) = 0.06 c correl(12) = 0.03 else c following are from AMIP II observations regridded to 3x3 c degree resolution correl(1) = 0.69 correl(2) = 0.47 correl(3) = 0.34 correl(4) = 0.27 correl(5) = 0.21 correl(6) = 0.17 correl(7) = 0.14 correl(8) = 0.12 c following are NOT taken from observations correl(9) = 0.09 correl(10) = 0.06 correl(11) = 0.03 correl(12) = 0.00 c observed values: c correl(9) = 0.10 c correl(10) = 0.08 c correl(11) = 0.06 c correl(12) = 0.03 endif do 32 i=1,niofiles inputfil(i) = inputsst(i) iyr1in(i) = iyr1sst(i) 32 continue c if output is pp format, but input is ezget, c then generate pp header c if ((outftype(1:2) .eq. 'pp') .and. c & (inftype(1:2) .ne. 'pp')) then c print*, 'at present input file type must be in' c print*, 'in pp format if output is in pp format' c stop c call genpphd c endif endif tmin = tmin - dcbias tmax = tmax - dcbias do 33 i=1,niofiles if (iyr1in(i) .ne. -999) nfilesi = i if (jyr1out(i) .ne. -999) nfileso = i 33 continue je = index(pathout, ' ') - 1 ie = index(abbrev, ' ') if (ie .eq. 0) then ie = 16 else ie = ie-1 endif outfile = pathout(1:je)// & '/'//src1(1:isrc1)//'bc_'//suffix//'_'//abbrev(1:ie)//'.out' c example: outfile = gisstbc_sst_1x1.out c *** fclimbc = pathout(1:je)//'/climbc_'//suffix//'_'//abbrev(1:ie) fclimbc = pathout(1:je)// & '/'//src1(1:isrc1)//'bc_'//suffix//'_'//abbrev(1:ie) fbcinfo = pathout(1:je)//'/bcinfo_'//suffix//'_'//abbrev(1:ie) fbc = pathout(1:je)// & '/'//src1(1:isrc1)//'bc_'//suffix//'_'//abbrev(1:ie) c *** if (src1(1:4) .eq. 'amip') then c fspinup = pathout(1:je)//'/bc1978_'//suffix//'_'//abbrev(1:ie) c else fspinup = pathout(1:je)//'/spinup_'//suffix//'_'//abbrev(1:ie) c endif fobs = pathout(1:je)// & '/'//src1(1:isrc1)//'obs_'//suffix//'_'//abbrev(1:ie) c *** fclimobs = pathout(1:je)//'/climobs_'//suffix//'_'//abbrev(1:ie) fclimobs = pathout(1:je)// & '/'//src1(1:isrc1)//'obs_'//suffix//'_'//abbrev(1:ie) c the following only used for lats output CC parmtabl = '/pcmdi/ktaylor/pcmdi/util/ketgrib.parms' CC dfltparm = '/usr/local/lats/table/amip2.lats.table' open(9, file=outfile, status='new') print*, ' ' print*, 'Output from monthly boundary condition generator ' print*, 'varin = ', varin print*, 'varout = ', varout print*, 'varout1 = ', varout1 print*, 'src = ', src print*, 'src1 = ', src1 print*, 'model = ', model print*, 'grid = ', grid print*, 'gtype = ', gtype print*, 'inftype = ', inftype print*, 'outftype = ', outftype print*, 'calclim = ', calclim print*, 'center = ', center print*, 'abbrev = ', abbrev print*, 'suffix = ', suffix print*, 'suff = ', suff print*, ' ' print*, 'icntl = ', icntl print*, 'lagcalc = ', lagcalc print*, 'iout = ', iout print*, ' ' print*, 'nlon = ', nlon print*, 'nlat = ', nlat print*, 'mlat = ', mlat print*, 'nchunks = ', nchunks print*, ' ' print*, 'nmon = ', nmon print*, 'mon1rd = ', mon1rd print*, 'iyr1rd = ', iyr1rd print*, 'monnrd = ', monnrd print*, 'iyrnrd = ', iyrnrd print*, 'mon1clm = ', mon1clm print*, 'iyr1clm = ', iyr1clm print*, 'monnclm = ', monnclm print*, 'iyrnclm = ', iyrnclm print*, 'mon1 = ', mon1 print*, 'iyr1 = ', iyr1 print*, 'monn = ', monn print*, 'iyrn = ', iyrn print*, 'mon1out = ', mon1out print*, 'iyr1out = ', iyr1out print*, 'monnout = ', monnout print*, 'iyrnout = ', iyrnout print*, 'mon1in = ', mon1in print*, 'ibasemon = ', ibasemon print*, 'ibaseyr = ', ibaseyr print*, ' ' print*, 'conv = ', conv print*, 'wtmin = ', wtmin print*, 'maxiter = ', maxiter print*, 'bbmin = ', bbmin print*, 'histo = ', histo print*, 'varmin = ', varmin print*, 'amissin = ', amissin print*, 'amissout = ', amissout print*, 'amissasc = ', amissasc print*, ' ' print*, 'correl = ', correl print*, ' ' print*, 'monlen = ', monlen print*, ' ' print*, 'tmin = ', tmin print*, 'tmax = ', tmax print*, 'dt = ', dt print*, 'dcbias = ', dcbias print*, 'fac = ', fac print*, 'offset = ', offset print*, ' ' do 41 i=1,nfilesi print*, 'iyr1in = ', iyr1in(i) print*, 'inputfil = ', inputfil(i) 41 continue print*, ' ' print*, 'msklndi = ', msklndi print*, 'sftnamei = ', sftnamei print*, 'sftfilei = ', sftfilei print*, 'msklndo = ', msklndo print*, 'sftnameo = ', sftnameo print*, 'sftfileo = ', sftfileo print*, ' ' print*, 'fclimbc = ', fclimbc print*, 'fbcinfo = ', fbcinfo print*, 'fspinup = ', fspinup print*, 'fclimobs = ', fclimobs print*, 'outfile = ', outfile print*, ' ' print*, 'jyr1out = ', (jyr1out(i), i=1,nfileso) print*, 'fbc = ', fbc print*, 'fobs = ', fobs write(9,*) ' ' write(9,*) 'Output from monthly boundary condition generator ' write(9,*) 'varin = ', varin write(9,*) 'varout = ', varout write(9,*) 'varout1 = ', varout1 write(9,*) 'src = ', src write(9,*) 'src1 = ', src1 write(9,*) 'model = ', model write(9,*) 'grid = ', grid write(9,*) 'gtype = ', gtype write(9,*) 'inftype = ', inftype write(9,*) 'outftype = ', outftype write(9,*) 'calclim = ', calclim write(9,*) 'center = ', center write(9,*) 'abbrev = ', abbrev write(9,*) 'suffix = ', suffix write(9,*) 'suff = ', suff write(9,*) ' ' write(9,*) 'icntl = ', icntl write(9,*) 'lagcalc = ', lagcalc write(9,*) 'iout = ', iout write(9,*) ' ' write(9,*) 'nlon = ', nlon write(9,*) 'nlat = ', nlat write(9,*) 'mlat = ', mlat write(9,*) 'nchunks = ', nchunks write(9,*) ' ' write(9,*) 'nmon = ', nmon write(9,*) 'mon1rd = ', mon1rd write(9,*) 'iyr1rd = ', iyr1rd write(9,*) 'monnrd = ', monnrd write(9,*) 'iyrnrd = ', iyrnrd write(9,*) 'mon1clm = ', mon1clm write(9,*) 'iyr1clm = ', iyr1clm write(9,*) 'monnclm = ', monnclm write(9,*) 'iyrnclm = ', iyrnclm write(9,*) 'mon1 = ', mon1 write(9,*) 'iyr1 = ', iyr1 write(9,*) 'monn = ', monn write(9,*) 'iyrn = ', iyrn write(9,*) 'mon1out = ', mon1out write(9,*) 'iyr1out = ', iyr1out write(9,*) 'monnout = ', monnout write(9,*) 'iyrnout = ', iyrnout write(9,*) 'mon1in = ', mon1in write(9,*) 'ibasemon = ', ibasemon write(9,*) 'ibaseyr = ', ibaseyr write(9,*) ' ' write(9,*) 'conv = ', conv write(9,*) 'wtmin = ', wtmin write(9,*) 'maxiter = ', maxiter write(9,*) 'bbmin = ', bbmin write(9,*) 'histo = ', histo write(9,*) 'varmin = ', varmin write(9,*) 'amissin = ', amissin write(9,*) 'amissout = ', amissout write(9,*) 'amissasc = ', amissasc write(9,*) ' ' write(9,*) 'correl = ', correl write(9,*) ' ' write(9,*) 'monlen = ', monlen write(9,*) ' ' write(9,*) 'tmin = ', tmin write(9,*) 'tmax = ', tmax write(9,*) 'dt = ', dt write(9,*) 'dcbias = ', dcbias write(9,*) 'fac = ', fac write(9,*) 'offset = ', offset write(9,*) ' ' do 42 i=1,nfilesi write(9,*) 'iyr1in = ', iyr1in(i) write(9,*) 'inputfil = ', inputfil(i) 42 continue write(9,*) ' ' write(9,*) 'msklndi = ', msklndi write(9,*) 'sftnamei = ', sftnamei write(9,*) 'sftfilei = ', sftfilei write(9,*) 'msklndo = ', msklndo write(9,*) 'sftnameo = ', sftnameo write(9,*) 'sftfileo = ', sftfileo write(9,*) ' ' write(9,*) 'fclimbc = ', fclimbc write(9,*) 'fbcinfo = ', fbcinfo write(9,*) 'fspinup = ', fspinup write(9,*) 'fclimobs = ', fclimobs write(9,*) 'outfile = ', outfile write(9,*) ' ' write(9,*) 'jyr1out = ', (jyr1out(i), i=1,nfileso) write(9,*) 'fbc = ', fbc write(9,*) 'fobs = ', fobs write(9,*) ' ' if (mod(nlat,mlat) .ne. 0) then print*, 'error in specification of mlat' print*, 'nlat/mlat must be an integer' stop endif if (inftype .ne. 'pp') then if ((mon1rd+nmon12*iyr1rd) .lt. (mon1in+nmon12*iyr1in(1))) then print*, 'error in time specifications: first month read' print*, 'must not preceed first month in input file.' print*, 'check mon1rd, iyr1rd, mon1in, iyr1in(1)' stop endif endif if ((mon1rd+nmon12*iyr1rd) .gt. (monnrd+nmon12*(iyrnrd-1))) then print*, 'error in time specifications: last month read must' print*, 'follow first month read by at least 1 year' print*, 'check mon1rd, iyr1rd, monnrd, iyrnrd' stop endif if ((mon1+nmon12*iyr1) .gt. (mon1rd+nmon12*iyr1rd)) then print*, 'error in time specifications: first month read must ' print*, 'not precede first month to be processed' print*, 'check mon1, iyr1, mon1rd, iyr1rd' stop endif if ((monn+nmon12*iyrn) .lt. (monnrd+nmon12*iyrnrd)) then print*, 'error in time specifications: last month read must ' print*, 'not follow last month processed' print*, 'check monn, iyrn, monnrd, iyrnrd' stop endif if (mod((monn+nmon12-mon1+1),nmon12) .ne. 0) then print*, 'error in time specifications: only integral number' print*, 'of years allowed' print*, 'check monn, mon1' stop endif if ((monn-mon1+nmon12*(iyrn-iyr1)+1) .ne. nmon) then print*, 'error in time specifications: parameter nmon ' print*, 'must be consistent with first and last months' print*, 'specified' print*, 'check nmon, mon1, iyr1, monn, iyrn' stop endif if ((mon1clm+nmon12*iyr1clm) .lt. (mon1rd+nmon12*iyr1rd)) then print*, 'error in time specifications: first month ' print*, 'contributing to climatology must not precede first' print*, 'month read' print*, 'check mon1clm, iyr1clm, mon1rd, iyr1rd' stop endif if ((monnclm+nmon12*iyrnclm) .gt. (monnrd+nmon12*iyrnrd)) then print*, 'error in time specifications: last month ' print*, 'contributing to climatology must not follow last' print*, 'month read' print*, 'check monnclm, iyrnclm, monnrd, iyrnrd' stop endif do 60 i=1,2 sum = 0.0 do 55 m=1,nmon12 cenmon(m,i) = sum + (monlen(m,i)+1.0)/2.0 sum = sum + monlen(m,i) 55 continue 60 continue c calculate jacobian elements for climatological years m = nmon12 mp = 1 do 63 n=1,nmon12 mm = m m = mp mp = mp + 1 if (mp .eq. (nmon12+1)) mp = 1 ac(n) = 2.0*monlen(m,1)/float(monlen(m,1)+monlen(mm,1)) cc(n) = 2.0*monlen(m,1)/float(monlen(m,1)+monlen(mp,1)) 63 continue c calculate jacobian elements for all years i = iyr1 j = 1 if (((mod(i,4) .eq. 0) .and. (mod(i,100) .ne. 0)) & .or. (mod(i,400) .eq. 0)) j = 2 m = mod((mon1+nmon12-2), nmon12) + 1 mp = mod((mon1-1), nmon12) + 1 do 65 n=1,nmon mm = m m = mp mp = mp + 1 if (mp .eq. (nmon12+1)) then mp = 1 i = i + 1 j = 1 if (((mod(i,4) .eq. 0) .and. (mod(i,100) .ne. 0)) .or. & (mod(i,400) .eq. 0)) j = 2 endif a(n) = 2.0*monlen(m,j)/float(monlen(m,j)+monlen(mm,j)) c(n) = 2.0*monlen(m,j)/float(monlen(m,j)+monlen(mp,j)) 65 continue do 100 i=1,nzonp aovarmon(i) = 0.0 acvarmon(i) = 0.0 wovar(i) = 0.0 wcvar(i) = 0.0 wocorr(i) = 0.0 wccorr(i) = 0.0 if (i .lt. nzonp) then woseas(i) = 0.0 wcseas(i) = 0.0 do 70 m=1,nmon12 acentmon(m,i) = 0.0 aobsclim(m,i) = 0.0 70 continue endif do 75 m=1,nlagm ccorrel(m,i) = 0.0 ocorrel(m,i) = 0.0 75 continue do 120 m=-10,10 jcount(i,m) = 0 mcount(i,m) = 0 120 continue do 130 m=-20,20 icount(i,m) = 0 kcount(i,m) = 0 130 continue 100 continue c test c do 101 m=1,nlon c alons(m) = -180. + (m-0.5)*360./nlon c 101 continue c do 102 m=1,nlat c alats(m) = 90. - (m-0.5)*180./nlat c 102 continue c loop through latitude chunks ismc = 0 jsmc = 0 ismf = 0 jsmf = 0 isea = 0 jn = 0 do 1000 jchnk=1,nchunks j1 = jn + 1 jn = min0((j1+mlat-1), nlat) c obtain monthly means lats = jn - j1 + 1 m1 = (iyr1rd-iyr1)*nmon12 + mon1rd - mon1 + 1 m2 = (iyrnrd-iyr1)*nmon12 + monnrd - mon1 + 1 mons = m2 - m1 + 1 print*, 'Reading input for chunk ', jchnk c test c goto 2345 if (inftype .eq. 'pp') then call readpp(inputfil, nlat, nlon, mlat, j1, lats, iyr1rd, & mon1rd, iyr1in, mons, amissin, amissout, alons, alats, & sst(1,1,m1), sstwts, wtfrac, lpp, bpp, space) lpp(29) = 0 bpp(5) = 0.0 bpp(18) = amissout bpp(19) = 1.0 else call getobs(model, gtype, msklndi, sftfilei, sftnamei, & msklndo, sftfileo, sftnameo, inputfil(1), varin, & nlon, nlat, mlat, j1, jn, iyr1in(1), mon1in, iyr1rd, & mon1rd, iyrnrd, monnrd, nmon12, amissin, amissout, wtmin, & sst(1,1,m1), sstwts(1,j1), wtfrac(1,j1), alons, alats(j1)) lpp(15) = nlat*nlon c the following indicates a regular lat/long grid boxes c (grid points are box centres) c ??? change for Gaussian grid? lpp(16) = 2 lpp(18) = nlat lpp(19) = nlon lpp(22) = 2 lpp(23) = lbfc bpp(11) = 90.0 c ??? following incorrect for Gaussian grid bpp(14) = alats(1) + 2.*alats(1)/(nlat-1) bpp(15) = -2.*alats(1)/(nlat-1) bpp(16) = alons(1) bpp(17) = alons(2) - alons(1) bpp(18) = amissout bpp(19) = 1.0 endif c test c 2345 continue print*, 'Finished reading input for chunk ', jchnk print*, 'latitudes ', alats(j1), ' through ', alats(jn) & , ' read successfully.' print*, ' ' write(9,*) 'latitudes ', alats(j1), ' through ', alats(jn) & , ' read successfully.' write(9,*) ' ' do 19 j=j1,jn do 18 i=1,nlon if (wtfrac(i,j) .gt. wtmin) then do 17 m=m1,m2 sst(i,j-j1+1,m) = fac*sst(i,j-j1+1,m) + offset 17 continue endif 18 continue 19 continue if (iout(4) .gt. 0) then c generate file name for chunk of data write(c3, '(i3.3)') jchnk tempobs(jchnk) = pathout(1:je)//'/tempobs_'//c3 nn = (iyrnout-iyr1out)*nmon12 + monnout-mon1out + 1 kk = (iyr1out-iyr1)*nmon12 + mon1out-mon1 + 1 call wrtpp(tempobs(jchnk), 1, nlon, mlat, alats(j1), lats, & iyr1out, mon1out, nn, lpp, bpp, sst(1,1,kk), fac, & amissout) tempobs(jchnk) = pathout(1:je)//'/tempobs_'//c3//'.pp' endif do 29 j=j1,jn do 28 i=1,nlon if (wtfrac(i,j) .gt. wtmin) then do 27 m=m1,m2 sst(i,j-j1+1,m) = sst(i,j-j1+1,m) - dcbias 27 continue endif 28 continue 29 continue c compute climatological means write(9,*) write(9,*) 'Calculating climatological values' c test c go to 3456 call calcclim(j1, jn, nlon, nlat, mlat, nmon, & ismc, jsmc, nzon, nzonp, nlagm, nmon12, mon1clm, iyr1clm, & iyr1, mon1, iyr1rd, mon1rd, iyrnrd, monnrd, & monnclm, iyrnclm, lagcalc, maxiter, amissout, histo, dcbias, & tmin, tmax, dt, wtmin, conv, bbmin, varmin, lagm, & isea, icount, jcount, kcount, mcount, wtfrac, sstwts, ovarmon, & aovarmon, obsclim, aobsclim, sst, woseas, wovar, & ocorrel, centmon, acentmon, cmax, cmin, omax, omin, & vecin, vecout, alons, alats, ac, cc, wcseas, wocorr) print*, ' ' print*, ' ' print*, ' ' print*, '*********** finished climatology *************' print*, ' ' print*, ' ' print*, ' ' c solve for full mid-month values write(9,*) write(9,*) 'Calculating actual monthly values' call calcfull(j1, jn, nlon, nlat, mlat, nmon, & iyr1rd, mon1rd, iyr1, mon1, iyrn, monn, iyrnrd, monnrd, & mon1clm, monnclm, iyr1clm, iyrnclm, ismf, jsmf, & nzon, nzonp, nmon12, lagcalc, nlagm, & maxiter, wtmin, amissout, tmin, tmax, & dt, conv, bbmin, varmin, lagm, sstwts, wtfrac, & wcseas, cvarmon, acvarmon, correl, ccorrel, & obsclim, vecin, vecout, centmon, sst, alons, alats, a, c, & ovarmon, wcvar, wccorr) call finish(j1, jn, nlon, nlat, mlat, mon1, nmon, nmon12, wtmin, & dcbias, amissout, wtfrac, centmon, obsclim, sst) c test c 3456 continue if (iout(1) .gt. 0) then c generate file name for chunk of data write(c3, '(i3.3)') jchnk tempbc(jchnk) = pathout(1:je)//'/tempbc_'//c3 nn = (iyrnout-iyr1out)*nmon12 + monnout-mon1out + 1 kk = (iyr1out-iyr1)*nmon12 + mon1out-mon1 + 1 call wrtpp(tempbc(jchnk), 0, nlon, mlat, alats(j1), lats, & iyr1out, mon1out, nn, lpp, bpp, sst(1,1,kk), fac, & amissout) tempbc(jchnk) = pathout(1:je)//'/tempbc_'//c3//'.pp' endif if (iout(2) .gt. 0) then c store spin-up data: starts in first January preceeding c (and nearest to) first output month (or if this month is c unavailable, starts at imon1, iyr1); ends in December of c same year c m1 = (iyr1out-iyr1)*nmon12 - mon1 + 2 c m2 = m1 + nmon12 - 1 c mon1spin = 1 c if (m1 .gt. 0) then c mon1spin = 1 c else c m1 = 1 c mon1spin = mon1 c endif c m2 = max0(m2, 1) c nmonspin = m2-m1+1 c iyr1spin = iyr1 + (mon1-2+m1)/nmon12 mon1spin = 1 if (mon1out .eq. 1) then iyr1spin = iyr1out -1 else iyr1spin = iyr1out endif if (iyr1spin .le. iyr1) then iyr1spin = iyr1 mon1spin = mon1 endif kyr1out(1) = iyr1spin m1 = (iyr1spin - iyr1)*nmon12 + (mon1spin - mon1) + 1 m2 = m1 + (nmon12-mon1spin) nmonspin = m2-m1+1 do 325 m=m1,m2 do 320 j=j1,jn do 310 i=1,nlon spinup(i,j,m-m1+1) = sst(i,j-j1+1,m) 310 continue 320 continue 325 continue endif 1000 continue print*, ' ' print*, ismc, ' climatological monthly means were smoothed' print*, jsmc, ' grid cells were affected' print*, ' ' write(9,*) ' ' write(9,*) ismc, ' climatological monthly means were smoothed' write(9,*) jsmc, ' grid cells were affected' write(9,*) ' ' write(9,*) ' ' write(9,*) 'ocean grid cells = ', isea write(9,*) ' ' c print*, ' ' print*, ismf, ' monthly means were smoothed' print*, jsmf, ' grid cells were affected' print*, ' ' write(9,*) ' ' write(9,*) ismf, ' monthly means were smoothed' write(9,*) jsmf, ' grid cells were affected' write(9,*) ' ' do 191 j=1,nzonp if (wovar(j) .gt. 0.0) aovarmon(j) = aovarmon(j)/wovar(j) 191 continue do 195 j=1,nzon if (woseas(j) .gt. 0.0) then do 194 n=1,nmon12 aobsclim(n,j) = aobsclim(n,j)/woseas(j) 194 continue else aovarmon(j) = amissout do 196 n=1,nmon12 aobsclim(n,j) = amissout 196 continue endif 195 continue if (lagcalc .gt. 0) then do 260 ik=1,nzonp if (wocorr(ik) .gt. 0.0) then do 255 k=1,nlagm ocorrel(k,ik) = ocorrel(k,ik)/wocorr(ik) 255 continue c else c do 261 k=1,nlagm c ocorrel(k,ik) = amissout c 261 continue endif 260 continue do 560 ik=1,nzonp if (wccorr(ik) .gt. 0.0) then do 555 k=1,nlagm ccorrel(k,ik) = ccorrel(k,ik)/wccorr(ik) 555 continue c else c do 556 k=1,nlagm c ccorrel(k,ik) = amissout c 556 continue endif 560 continue endif do 505 ik=1,nzon do 502 n=-20,20 icount(nzonp,n) = icount(nzonp,n)+icount(ik,n) kcount(nzonp,n) = kcount(nzonp,n)+kcount(ik,n) 502 continue do 503 n=-10,10 jcount(nzonp,n) = jcount(nzonp,n)+jcount(ik,n) mcount(nzonp,n) = mcount(nzonp,n)+mcount(ik,n) 503 continue do 504 m=1,nmon12 if (wcseas(ik) .gt. 0.0) then acentmon(m,ik) = acentmon(m,ik)/wcseas(ik) else acentmon(m,ik) = amissout endif 504 continue 505 continue do 906 ik=1,nzonp if (wcvar(ik) .gt. 0.0) then acvarmon(ik) = acvarmon(ik)/wcvar(ik) elseif (ik .lt. nzonp) then if (wcseas(ik) .eq. 0.0) acvarmon(ik) = amissout endif 906 continue write(9,*) &'Number of cells where difference in climatological max is ' write(9,*) 'within +- ', histo, ' of 0.0: ' write(9,*) (icount(ik,0), ik=1,nzonp) write(9,*) ' ' write(9,*) 'number of cells where difference between constructed' write(9,*) 'mid-point climatological monthly maximum and observed' write(9,*) 'climatological monthly mean exceeds different values:' write(9,*) ' ' if (alats(1) .lt. alats(2)) then write(9,*) 'value ', & '-90 & -60 -60 & -30 -30 & 0 0 & 30 30 & 60 60 & 90', & ' Global' else write(9,*) 'value ', & ' 90 & 60 60 & 30 30 & 0 0 & -30 -30 & -60 -60 & -90', & ' Global' endif write(9,*) 'exceeded' do 201 i=-20,20 if (i .ne. 0) then r = i*histo write(9,'(f10.3, 7i10)') r, (icount(ik,i), ik=1,nzonp) endif 201 continue write(9,*) ' ' do 301 i=-10,10 if (i .ne. 0) then r = i*10*histo write(9,'(f10.3, 7i10)') r, (jcount(ik,i), ik=1,nzonp) endif 301 continue write(9,*) ' ' write(9,*) &'Number of cells where difference in climatological min is ' write(9,*) 'within +- ', histo, ' of 0.0: ' write(9,*) (kcount(ik,0), ik=1,nzonp) write(9,*) ' ' write(9,*)'number of cells where -(difference between constructed' write(9,*) 'mid-point climatological monthly minimum and observed' write(9,*) & 'climatological monthly mean minimum) exceeds different values:' write(9,*) ' ' if (alats(1) .lt. alats(2)) then write(9,*) 'value ', & '-90 & -60 -60 & -30 -30 & 0 0 & 30 30 & 60 60 & 90', & ' Global' else write(9,*) 'value ', & ' 90 & 60 60 & 30 30 & 0 0 & -30 -30 & -60 -60 & -90', & ' Global' endif write(9,*) 'exceeded' do 202 i=-20,20 if (i .ne. 0) then r = i*histo write(9,'(f10.3, 7i10)') r, (kcount(ik,i), ik=1,nzonp) endif 202 continue write(9,*) ' ' do 302 i=-10,10 if (i .ne. 0) then r = i*10*histo write(9,'(f10.3, 7i10)') r, (mcount(ik,i), ik=1,nzonp) endif 302 continue write(9,*) ' ' write(9,*) ' ' write(9,*) 'Climatological Seasonal Variation by Zone' write(9,*) ' ' if (alats(2) .lt. alats(1)) then write(9,*) 'month ', & ' 90 & 60 60 & 30 30 & 0 0 & -30 -30 & -60', & ' -60 & -90' else write(9,*) 'month ', & ' -90 & -60 -60 & -30 -30 & 0 0 & 30 30 & 60 ', & ' 60 & 90' endif write(9,*)' avg mid avg mid avg mid avg', & ' mid avg mid avg mid' do 304 n=1,nmon12 write(9,'(i2, 6x, 12f6.2)') & n, (aobsclim(n,ik), acentmon(n,ik), ik=1,nzon) 304 continue write(9,'(/ ''area '', 12f6.3)') & (woseas(ik), wcseas(ik), ik=1,nzon) c write(9,*) ' ' write(9,*) ' ' write(9,*) 'Variance of Monthly Anomalies by Zone' write(9,*) ' ' if (alats(2) .lt. alats(1)) then write(9,*) & ' 90 & 60 60 & 30 30 & 0 0 & -30 -30 & -60', & ' -60 & -90 global' else write(9,*) & ' -90 & -60 -60 & -30 -30 & 0 0 & 30 30 & 60 ', & ' 60 & 90 global' endif write(9,'(''mon avg.'', 7f10.3)') (aovarmon(ik), ik=1,nzonp) write(9,'(''mid mon '', 7f10.3)') (acvarmon(ik), ik=1,nzonp) write(9,*) ' ' write(9,'(''mon area'', 7f10.5)') (wovar(ik), ik=1,nzonp) write(9,'(''mid area'', 7f10.5)') (wcvar(ik), ik=1,nzonp) write(9,*) ' ' if (lagcalc .gt. 0) then write(9,*) ' ' write(9,*) 'Variance and Auto Correlation by Zone' write(9,*) ' ' write(9,*) ' Monthly Mean Anomalies ' if (alats(2) .lt. alats(1)) then write(9,*) & ' 90 & 60 60 & 30 30 & 0 0 & -30 -30 & -60', & ' -60 & -90 global' else write(9,*) & ' -90 & -60 -60 & -30 -30 & 0 0 & 30 30 & 60 ', & ' 60 & 90 global' endif if (alats(2) .lt. alats(1)) then write(9,*) & ' 90 & 60 60 & 30 30 & 0 0 & -30 -30 & -60', & ' -60 & -90 global' else write(9,*) & ' -90 & -60 -60 & -30 -30 & 0 0 & 30 30 & 60 ', & ' 60 & 90 global' endif write(9,'(''variance:'', 7f10.3 )') (ocorrel(1,ik), ik=1,nzonp) write(9,*) 'lag Auto Correlation' do 306 k=2,nlagm write(9,'(i8, 7f10.3)') lagm(k), (ocorrel(k,ik), ik=1,nzonp) 306 continue write(9,*) ' ' write(9,'(''area '', 7f10.5)') (wocorr(ik), ik=1,nzonp) write(9,*) ' ' write(9,*) ' Mid-Month Anomalies ' if (alats(2) .lt. alats(1)) then write(9,*) & ' 90 & 60 60 & 30 30 & 0 0 & -30 -30 & -60', & ' -60 & -90 global' else write(9,*) & ' -90 & -60 -60 & -30 -30 & 0 0 & 30 30 & 60 ', & ' 60 & 90 global' endif write(9,'(''variance:'', 7f10.3 )') (ccorrel(1,ik), ik=1,nzonp) write(9,*) 'lag Auto Correlation' do 308 k=2,nlagm write(9,'(i8, 7f10.3)') lagm(k), (ccorrel(k,ik), ik=1,nzonp) 308 continue write(9,*) ' ' write(9,'(''area '', 7f10.5)') (wccorr(ik), ik=1,nzonp) write(9,*) ' ' endif if ((outftype(1:2) .eq. 'pp') .or. (outftype(1:3) .eq. 'all')) & then if (iout(1) .gt. 0) then call ppconcat (tempbc, fbc, nmon12, nlon, nlat, nchunks, & iyr1out, mon1out, iyrnout, monnout, jyr1out, & space) endif if (iout(2) .gt. 0) then c write spin-up data write(ayr, '(i4)') iyr1spin iend = index(fspinup, ' ') - 1 outcopy = fspinup(1:iend) // '_' //ayr//' ' call wrtpp(outcopy, 0, nlon, nlat, alats(1), nlat, & iyr1spin, mon1spin, nmonspin, lpp, bpp, spinup, fac, & amissout) endif if (iout(3) .gt. 0) then c write climatology b.c. iend = index(fclimbc, ' ') - 1 outcopy = fclimbc(1:iend) //'_clim ' call wrtpp(outcopy, 0, nlon, nlat, alats(1), & nlat, 0, 1, nmon12, lpp, bpp, centmon, fac, amissout) endif if (iout(4) .gt. 0) then call ppconcat (tempobs, fobs, nmon12, nlon, nlat, nchunks, & iyr1out, mon1out, iyrnout, monnout, jyr1out, & space) endif if (iout(5) .gt. 0) then c write climatology of obs iend = index(fclimobs, ' ') - 1 outcopy = fclimobs(1:iend) //'_clim ' call wrtpp(outcopy, 1, nlon, nlat, alats(1), nlat, & 0, 1, nmon12, lpp, bpp, obsclim, fac, amissout) endif endif if ((outftype(1:3) .eq. 'cdf') .or. & (outftype(1:6) .eq. 'coards') .or. & (outftype(1:5) .eq. 'grads') .or. & (outftype(1:4) .eq. 'grib') .or. & (outftype(1:3) .eq. 'all')) then if (outftype(1:3) .eq. 'all') then typfil = 'coards' i1 = 0 else typfil = outftype i1 = 1 endif if (typfil(1:6) .eq. 'coards') then calclim = 'gregorian' else calclim = 'climatology' endif 330 if (iout(1) .gt. 0) then title = src(1:isrc)// & ' Boundary Condition Data: Mid-Month '//suff(1:iie) nn = (iyrnout-iyr1out)*nmon12 + monnout-mon1out + 1 call wrtlats(ig, 1, 1, 1, 1, 1, 1, & 'mid', 'instantaneous', varout, fbc, typfil, & dfltparm, gtype, calgreg, title, sourceb, model, center, & amissout, nlon, nlat, nchunks, alons, alats, nn, & mon1out, iyr1out, ibasemon, ibaseyr, idvar, space, & tempbc, jyr1out) ig = 0 endif if (iout(2) .gt. 0) then title = src(1:isrc)// & ' Spin-up Boundary Condition Data: Mid-Month ' & //suff(1:iie) call wrtlats(ig, 1, 1, 1, 1, 1, 0, & 'mid', 'instantaneous', varout, fspinup, typfil, & dfltparm, gtype, calgreg, title, sourceb, model, center, & amissout, nlon, nlat, nchunks, alons, alats, nmonspin, & mon1spin, iyr1spin, ibasemon, ibaseyr, idvar, spinup, & dum, kyr1out) ig = 0 endif if (iout(3) .gt. 0) then title = src(1:isrc)// & ' Climatological Bdry. Condition Data: Mid-Month ' & //suff(1:iie) call wrtlats(ig, 1, 1, 1, 1, 1, 0, & 'mid', 'instantaneous', varout, fclimbc, typfil, & dfltparm, gtype, calclim, title, sourceb, model, center, & amissout, nlon, nlat, nchunks, alons, alats, nmon12, & 1, 2, 1, 2, idvar, centmon, dum, jyr1out) ig = 0 endif if (iout(4) .gt. 0) then title = src(1:isrc)// ' Observed Monthly Mean '//suff(1:iie) nn = (iyrnout-iyr1out)*nmon12 + monnout-mon1out + 1 call wrtlats(ig, 1, 1, 1, 1, 1, 1, & 'mean', 'average', varout1, fobs, typfil, & dfltparm, gtype, calgreg, title, sourceo, model, center, & amissout, nlon, nlat, nchunks, alons, alats, nn, & mon1out, iyr1out, ibasemon, ibaseyr, idvar, space, & tempobs, jyr1out) ig = 0 endif if (iout(5) .gt. 0) then title = src(1:isrc)// & ' Observed Climatology: Monthly Mean '//suff(1:iie) call wrtlats(ig, 1, 1, 1, 1, 1, 0, & 'mean','average', varout1, fclimobs, typfil, & dfltparm, gtype, calclim, title, sourceo, model, center, & amissout, nlon, nlat, nchunks, alons, alats, nmon12, & 1, 2, 1, 2, idvar, obsclim, dum, jyr1out) ig = 0 endif if (iout(6) .gt. 0) then vname = varout1//'max' title = src(1:isrc)// & ' Maximum (Monthly) Mean Climatological '// suff(1:iie) call wrtlats(ig, 1, 1, 1, 0, 0, 0, & 'none', 'average', vname, fbcinfo, typfil, & parmtabl, gtype, calclim, title, sourceo, model, center, & amissout, nlon, nlat, nchunks, alons, alats, 0, & 0, 0, 0, 0, idvara(1), omax, dum, jyr1out) ig = 0 vname = varout//'max' title = src(1:isrc)// & ' Maximum (Mid-Month) Climatological '// suff(1:iie) call wrtlats(0, 0, 1, 1, 0, 0, 0, & 'none', 'average', vname, fbcinfo, typfil, & 'none', gtype, calclim, title, sourceb, model, center, & amissout, nlon, nlat, nchunks, alons, alats, 0, & 0, 0, 0, 0, idvara(2), cmax, dum, jyr1out) vname = varout1//'min' title = src(1:isrc)// & ' Minimum (Monthly) Mean Climatological '// suff(1:iie) call wrtlats(0, 0, 1, 1, 0, 0, 0, & 'none', 'average', vname, fbcinfo, typfil, & 'none', gtype, calclim, title, sourceo, model, center, & amissout, nlon, nlat, nchunks, alons, alats, 0, & 0, 0, 0, 0, idvara(3), omin, dum, jyr1out) vname = varout//'min' title = src(1:isrc)// & ' Minimum (Mid-Month) Climatological '// suff(1:iie) call wrtlats(0, 0, 1, 1, 0, 0, 0, & 'none', 'average', vname, fbcinfo, typfil, & 'none', gtype, calclim, title, sourceb, model, center, & amissout, nlon, nlat, nchunks, alons, alats, 0, & 0, 0, 0, 0, idvara(4), cmin, dum, jyr1out) vname = varout1//'var' m1 = index(suff, '(') if (m1 .gt. 0) then m1 = m1 - 1 else m1 = iie endif title = src(1:isrc)// & ' Variance of Monthly Mean Anomalies: '// suff(1:m1) call wrtlats(0, 0, 1, 1, 0, 0, 0, & 'none', 'average', vname, fbcinfo, typfil, & 'none', gtype, calclim, title, sourceo, model, center, & amissout, nlon, nlat, nchunks, alons, alats, 0, & 0, 0, 0, 0, idvara(5), ovarmon, dum, jyr1out) vname = varout//'var' title = src(1:isrc)// & ' Variance of Mid-Month Anomalies: '// suff(1:m1) call wrtlats(0, 0, 1, 1, 0, 0, 0, & 'none', 'average', vname, fbcinfo, typfil, & 'none', gtype, calclim, title, sourceb, model, center, & amissout, nlon, nlat, nchunks, alons, alats, 0, & 0, 0, 0, 0, idvara(6), cvarmon, dum, jyr1out) vname = varout1//'max' call wrtlats(0, 0, 0, 1, 1, 0, 0, & 'none', 'average', vname, fbcinfo, typfil, & 'none', gtype, calclim, title, sourceo, model, center, & amissout, nlon, nlat, nchunks, alons, alats, 0, & 0, 0, 0, 0, idvara(1), omax, dum, jyr1out) vname = varout//'max' call wrtlats(0, 0, 0, 1, 1, 0, 0, & 'none', 'average', vname, fbcinfo, typfil, & 'none', gtype, calclim, title, sourceb, model, center, & amissout, nlon, nlat, nchunks, alons, alats, 0, & 0, 0, 0, 0, idvara(2), cmax, dum, jyr1out) vname = varout1//'min' call wrtlats(0, 0, 0, 1, 1, 0, 0, & 'none', 'average', vname, fbcinfo, typfil, & 'none', gtype, calclim, title, sourceo, model, center, & amissout, nlon, nlat, nchunks, alons, alats, 0, & 0, 0, 0, 0, idvara(3), omin, dum, jyr1out) vname = varout//'min' call wrtlats(0, 0, 0, 1, 1, 0, 0, & 'none', 'average', vname, fbcinfo, typfil, & 'none', gtype, calclim, title, sourceb, model, center, & amissout, nlon, nlat, nchunks, alons, alats, 0, & 0, 0, 0, 0, idvara(4), cmin, dum, jyr1out) vname = varout1//'var' call wrtlats(0, 0, 0, 1, 1, 0, 0, & 'none', 'average', vname, fbcinfo, typfil, & 'none', gtype, calclim, title, sourceo, model, center, & amissout, nlon, nlat, nchunks, alons, alats, 0, & 0, 0, 0, 0, idvara(5), ovarmon, dum, jyr1out) vname = varout//'var' call wrtlats(0, 0, 0, 1, 1, 1, 0, & 'none', 'average', vname, fbcinfo, typfil, & 'none', gtype, calclim, title, sourceb, model, center, & amissout, nlon, nlat, nchunks, alons, alats, 0, & 0, 0, 0, 0, idvara(6), cvarmon, dum, jyr1out) endif if (i1 .eq. 0) then typfil = 'grib' i1 = 1 go to 330 endif endif if ((outftype(1:3) .eq. 'drs') .or. (outftype(1:3) .eq. 'all')) & then elem1(1) = alons(1) elemn(1) = alons(nlon) elem1(2) = alats(1) elemn(2) = alats(nlat) if (caseindp(gtype, 'gau') + .or. caseindp(gtype, 'bmr') .or. caseindp(gtype, 'ccc') + .or. caseindp(gtype, 'cnr') .or. caseindp(gtype, 'col') + .or. caseindp(gtype, 'csi') .or. caseindp(gtype, 'der') + .or. caseindp(gtype, 'ecm') .or. caseindp(gtype, 'gfd') + .or. caseindp(gtype, 'mgo') .or. caseindp(gtype, 'mpi') + .or. caseindp(gtype, 'nca') .or. caseindp(gtype, 'nmc') + .or. caseindp(gtype, 'rpn') .or. caseindp(gtype, 'sun') + .or. caseindp(gtype, 'uga') .or. caseindp(gtype, 'sng') + .or. caseindp(gtype, 'gen') + .or. caseindp(gtype, 'ccm') .or. caseindp(gtype, 'ccs') + .or. caseindp(gtype, 'ech') .or. caseindp(gtype, 'nce')) then do 108 i=1,nlat gauss(i) = alats(i) 108 continue elseif (caseindp(gtype, 'cos') + .or. caseindp(gtype, 'dnm') .or. caseindp(gtype, 'gis') + .or. caseindp(gtype, 'csu') .or. caseindp(gtype, 'ucl') + .or. caseindp(gtype, 'gla') .or. caseindp(gtype, 'gsf') + .or. caseindp(gtype, 'iap') .or. caseindp(gtype, 'mri') + .or. caseindp(gtype, 'uiu') + .or. caseindp(gtype, 'jma') .or. caseindp(gtype, 'nrl') + .or. caseindp(gtype, 'ukm') .or. caseindp(gtype, 'yon')) + then gauss(1) = 0.0 gauss(2) = 0.0 else print*, 'Grid type not recognized ' print*, 'Grid type = ', gtype write(9,*) 'Grid type not recognized ' write(9,*) 'Grid type = ', gtype stop endif if (iout(1) .gt. 0) then vname = varout title = src(1:isrc)// & ' Boundary Condition Data: Mid-Month '//suff(1:iie) call wrtdrs(1) endif if (iout(2) .gt. 0) then vname = varout title = src(1:isrc)// & ' Spin-up Boundary Condition Data: Mid-Month ' & //suff(1:iie) m2 = mon1spin+nmonspin-1 call wrtdrs(1) endif if (iout(3) .gt. 0) then vname = varout title = src(1:isrc)// & ' Climatological Bdry. Condition Data: Mid-Month ' & //suff(1:iie) call wrtdrs(1) endif if (iout(4) .gt. 0) then vname = varout1 title = src(1:isrc)// ' Observed Monthly Mean '//suff(1:iie) call wrtdrs(1) endif if (iout(5) .gt. 0) then vname = varout1 title = src(1:isrc)// & ' Observed Climatology: Monthly Mean '//suff(1:iie) call wrtdrs(1) endif endif if ((outftype(1:3) .eq. 'drs') .or. (outftype(1:3) .eq. 'asc') & .or. (outftype(1:3) .eq. 'all')) then c & .or. (outftype(1:3) .eq. 'pp')) then if (iout(6) .gt. 0) then vname = varout1//'max' title = src(1:isrc)// & ' Maximum (Monthly) Mean Climatological '// suff(1:iie) call wrtdrs(1) vname = varout//'max' title = src(1:isrc)// & ' Maximum (Mid-Month) Climatological '// suff(1:iie) call wrtdrs(0) vname = varout1//'min' title = src(1:isrc)// & ' Minimum (Monthly) Mean Climatological '// suff(1:iie) call wrtdrs(0) vname = varout//'min' title = src(1:isrc)// & ' Minimum (Mid-Month) Climatological '// suff(1:iie) call wrtdrs(0) vname = varout1//'var' m1 = index(suff, '(') if (m1 .gt. 0) then m1 = m1 - 1 else m1 = iie endif title = src(1:isrc)// & ' Variance of Monthly Mean Anomalies: '// suff(1:m1) call wrtdrs(0) vname = varout//'var' title = src(1:isrc)// & ' Variance of Mid-Month Anomalies: '// suff(1:m1) call wrtdrs(0) endif endif if ((outftype(1:5) .eq. 'ascii') .or. (outftype(1:3) .eq. 'all')) & then c should call wrtasc last (if 'all') because missing data value c is different from other formats and sub wrtasc alters c the output arrays. elem1(1) = alons(1) elemn(1) = alons(nlon) elem1(2) = alats(1) elemn(2) = alats(nlat) if (iout(1) .gt. 0) then vname = varout title = src(1:isrc)// & ' Boundary Condition Data: Mid-Month '//suff(1:iie) call wrtasc(fbc, nlon, nlat, nmon12, iyr1out, mon1out, & iyrnout, monnout, elem1, elemn, model, center, sourceb, & vname, title, units, grid, amissout, amissasc, 1, & nchunks, tempbc, jyr1out, space) endif if (iout(2) .gt. 0) then vname = varout title = src(1:isrc)// & ' Spin-up Boundary Condition Data: Mid-Month ' & //suff(1:iie) m2 = mon1spin+nmonspin-1 call wrtasc(fspinup, nlon, nlat, nmon12, iyr1spin,mon1spin, & iyr1spin, m2, elem1, elemn, model, center, sourceb, & vname, title, units, grid, amissout, amissasc, 0, & nchunks, dum, kyr1out, spinup) endif if (iout(3) .gt. 0) then vname = varout title = src(1:isrc)// & ' Climatological Bdry. Condition Data: Mid-Month ' & //suff(1:iie) call wrtasc(fclimbc, nlon, nlat, nmon12, 0, 1, & 0, nmon12, elem1, elemn, model, center, sourceb, & vname, title, units, grid, amissout, amissasc, 0, & nchunks, dum, jyr1out, centmon) endif if (iout(4) .gt. 0) then vname = varout1 title = src(1:isrc)// ' Observed Monthly Mean '//suff(1:iie) call wrtasc(fobs, nlon, nlat, nmon12, iyr1out, mon1out, & iyrnout, monnout, elem1, elemn, model, center, sourceo, & vname, title, units, grid, amissout, amissasc, 1, & nchunks, tempobs, jyr1out, space) endif if (iout(5) .gt. 0) then vname = varout1 title = src(1:isrc)// & ' Observed Climatology: Monthly Mean '//suff(1:iie) call wrtasc(fclimobs, nlon, nlat, nmon12, 0, 1, & 0, nmon12, elem1, elemn, model, center, sourceo, & vname, title, units, grid, amissout, amissasc, 0, & nchunks, dum, jyr1out, obsclim) endif endif close(9) if (icntl .eq. 3) then icntl = 1 go to 22 endif c syscomm = 'rm '//pathout(1:je)//'/temp*.pp' c call system(syscomm) print*, ' ' print*, ' ' print*, 'Program mkgisst finished execution.' print*, ' ' print*, ' ' stop end subroutine finish(j1, jn, nlon, nlat, mlat, mon1, nmon, nmon12, & wtmin, dcbias, amissout, wtfrac, centmon, obsclim, sst) implicit none integer nlon, nlat, mon1, nmon, nmon12, mlat real wtmin, dcbias, amissout real wtfrac(nlon,nlat), centmon(nlon,nlat,nmon12), & obsclim(nlon,nlat,nmon12), sst(nlon,mlat,nmon) integer i, j, m, n, j1, jn c calculate full mid-month time-series do 990 j=j1,jn do 980 i=1,nlon if (wtfrac(i,j) .le. wtmin) then do 920 n=1,nmon sst(i,j-j1+1,n) = amissout 920 continue else do 930 m=1,nmon12 centmon(i,j,m) = centmon(i,j,m) + dcbias obsclim(i,j,m) = obsclim(i,j,m) + dcbias 930 continue do 970 n=1,nmon m = mod((n+mon1-2), nmon12) + 1 sst(i,j-j1+1,n) = sst(i,j-j1+1,n) + centmon(i,j,m) 970 continue endif 980 continue 990 continue return end subroutine calcclim(j1, jn, nlon, nlat, mlat, nmon, & ismc, jsmc, nzon, nzonp, nlagm, nmon12, mon1clm, iyr1clm, & iyr1, mon1, iyr1rd, mon1rd, iyrnrd, monnrd, & monnclm, iyrnclm, lagcalc, maxiter, amiss, histo, dcbias, & tmin, tmax, dt, wtmin, conv, bbmin, varmin, lagm, & isea, icount, jcount, kcount, mcount, wtfrac, sstwts, ovarmon, & aovarmon, obsclim, aobsclim, sst, woseas, wovar, ocorrel, & centmon, acentmon, cmax, cmin, omax, omin, vecin, & vecout, alons, alats, ac, cc, wcseas, wocorr) implicit none integer nlon, nlat, nmon, nzon, j1, jn, mon1clm, iyr1clm, monnclm, & iyrnclm, nzonp, nlagm, nmon12, & lagcalc, maxiter, isea, mlat, & iyr1, mon1, iyr1rd, mon1rd, iyrnrd, monnrd real amiss, histo, dcbias, tmin, tmax, dt, wtmin, conv, & bbmin, varmin integer icount(nzonp,-20:20), jcount(nzonp,-10:10), & kcount(nzonp,-20:20), mcount(nzonp, -10:10), lagm(nlagm) real wtfrac(nlon,nlat), sstwts(nlon,nlat), ovarmon(nlon,nlat), & aovarmon(nzonp), obsclim(nlon,nlat,nmon12), & aobsclim(nmon12,nzon), sst(nlon,mlat,nmon), & woseas(nzon), wovar(nzonp), ocorrel(nlagm,nzonp), & centmon(nlon,nlat,nmon12), & acentmon(nmon12,nzon), cmax(nlon,nlat), cmin(nlon,nlat), & omax(nlon,nlat), omin(nlon,nlat), vecin(nmon), & vecout(nmon), alons(nlon), alats(nlat), ac(nmon12), & cc(nmon12), wcseas(nzon), wocorr(nzonp) integer ismc, jsmc, i, j, ik, n, m, nn, mn, k, mm, momax, & momin, mmax, mmin, ii, jj, kk, mmnn, ji, m1, m2, m3, m4 real var, ocentmax, ocentmin, centmax, centmin, diffmax, diffmin, & rmin, rmax, rrmin, rrmax double precision sum, avg1, avg2, covar mmnn = 0 rrmin = 0.0 rrmax = 0.0 m1 = nmon12*(iyr1rd-iyr1) + mon1rd-mon1 + 1 m2 = nmon12*(iyrnrd-iyr1) + monnrd-mon1 + 1 m3 = nmon12*(iyr1clm-iyr1) + mon1clm-mon1 + 1 m4 = nmon12*(iyrnclm-iyr1) + monnclm-mon1 + 1 do 190 j=j1,jn ji = j-j1+1 ik = (j-1)/(nlat/nzon) + 1 ik = min0(ik, nzon) print*, 'computing climatology for latitude ', alats(j) do 180 i=1,nlon if (wtfrac(i,j) .le. wtmin) then ovarmon(i,j) = amiss do 148 n=1,nmon12 obsclim(i,j,n) = amiss 148 continue else do 150 n=1,nmon12 obsclim(i,j,n) = 0.0 150 continue mn = 0 rmin = 0.0 rmax = 0.0 do 155 n=m1,m2 if (sst(i,ji,n) .gt. tmax) then if (sst(i,ji,n) .gt. tmax+1.0e-2) then print*, 'tmax exceeded ', i, j, n, sst(i,ji,n) rmax = amax1(rmax, (sst(i,ji,n)-tmax)) mn = mn + 1 endif sst(i,ji,n) = tmax elseif (sst(i,ji,n) .lt. tmin) then if (sst(i,ji,n) .lt. tmin-dt) then c print*, 'tmin exceeded ', i, j, n, sst(i,ji,n) if ((tmin-sst(i,ji,n)) .lt. 900.) then rmin = amin1(rmin, (sst(i,ji,n)-tmin)) mn = mn + 1 endif endif sst(i,ji,n) = tmin endif 155 continue if (mn .gt. 0) then mmnn = mmnn + 1 rrmin = amin1(rrmin, rmin) rrmax = amax1(rrmax, rmax) write(9,*) ' ' write(9,*) & ' WARNING -- observed value exceeds limits at' write(9,*) mn, ' time points at latitude ', alats(j) write(9,*) 'and longitude ', alons(i) write(9,*) 'max error = ', rmax, & ' min error = ', rmin endif c compute climatology woseas(ik) = woseas(ik) + sstwts(i,j) do 170 n=1,nmon12 m = mod((mon1clm+n-2), nmon12) + 1 obsclim(i,j,m)=0.0 nn = 0 do 160 k = m3+n-1, m4, nmon12 nn = nn + 1 obsclim(i,j,m) = obsclim(i,j,m)+sst(i,ji,k) 160 continue obsclim(i,j,m) = obsclim(i,j,m)/nn aobsclim(m,ik) = aobsclim(m,ik) + & obsclim(i,j,m)*sstwts(i,j) 170 continue c remove climatology to generate anomalies do 175 n=m1,m2 mm = mod((mon1rd+n-m1-1), nmon12) + 1 sst(i,ji,n) = sst(i,ji,n) - obsclim(i,j,mm) 175 continue sum = 0.0 avg1 = 0.0 do 177 n=m3,m4 sum = sum + sst(i,ji,n)*sst(i,ji,n) avg1 = avg1 + sst(i,ji,n) 177 continue ovarmon(i,j) = & (sum-(avg1*avg1/(m4-m3+1)))/float(m4-m3) if (ovarmon(i,j) .ge. varmin) then aovarmon(ik) = aovarmon(ik) + ovarmon(i,j)*sstwts(i,j) aovarmon(nzonp) = aovarmon(nzonp) + & ovarmon(i,j)*sstwts(i,j) wovar(ik) = wovar(ik) + sstwts(i,j) wovar(nzonp) = wovar(nzonp) + sstwts(i,j) endif endif 180 continue 190 continue if (mmnn .gt. 0) then write(9,*) ' ' write(9,*) ' WARNING -- observed value exceeds limits at' write(9,*) mmnn, ' grid cells' write(9,*) 'max error = ', rrmax, ' min error = ', rrmin endif if (lagcalc .gt. 0) then do 250 j=j1,jn print*, 'computing lag correlations for latitude ', alats(j) ik = (j-1)/(nlat/nzon) + 1 ik = min0(ik, nzon) ji = j-j1+1 do 240 i=1,nlon if ((wtfrac(i,j) .gt. wtmin) .and. (ovarmon(i,j) .ge. varmin)) & then wocorr(ik) = wocorr(ik) + sstwts(i,j) wocorr(nzonp) = wocorr(nzonp) + sstwts(i,j) do 230 k=1,nlagm covar = 0.0 avg1 = 0.0 avg2 = 0.0 nn = m4 - lagm(k) do 220 n=m3,nn mm = n+lagm(k) covar = covar + sst(i,ji,n)*sst(i,ji,mm) avg1 = avg1 + sst(i,ji,n) avg2 = avg2 + sst(i,ji,mm) 220 continue if (k .eq. 1) then c calculate variance var = (covar - (avg1*avg2/(nn-m3+1)))/(nn-m3) ocorrel(k,ik) = ocorrel(k,ik) + sstwts(i,j)*var ocorrel(k,nzonp) = ocorrel(k,nzonp) + sstwts(i,j)*var elseif (var .gt. 0.) then c calculate correlations ocorrel(k,ik) = ocorrel(k,ik) + sstwts(i,j)* & ((covar - (avg1*avg2/(nn-m3+1)))/(nn-m3))/var ocorrel(k,nzonp) = ocorrel(k,nzonp) + sstwts(i,j)* & ((covar - (avg1*avg2/(nn-m3+1)))/(nn-m3))/var endif 230 continue endif 240 continue 250 continue endif c solve for climatological mid-month values do 500 j=j1,jn print*, & 'Computing climatological mid-month values for latitude ', & alats(j) ik = (j-1)/(nlat/nzon) + 1 ik = min0(ik, nzon) do 400 i=1,nlon if (wtfrac(i,j) .le. wtmin)then do 308 m=1,nmon12 centmon(i,j,m) = amiss 308 continue cmax(i,j) = amiss cmin(i,j) = amiss omax(i,j) = amiss omin(i,j) = amiss else do 310 m=1,nmon12 vecin(m) = obsclim(i,j,m) vecout(m) = obsclim(i,j,m) 310 continue call solvmid(alons(i), alats(j), nmon12, conv, dt, tmin, & tmax, bbmin, maxiter, ac, cc, vecin, vecout, ismc, jsmc) do 320 m=1,nmon12 centmon(i,j,m) = vecout(m) 320 continue ocentmax = -1.e20 ocentmin = 1.e20 centmax = -1.e20 centmin = 1.e20 do 390 m=1,nmon12 c find max and min values and months of max and min sst if (obsclim(i,j,m) .gt. ocentmax) then ocentmax = obsclim(i,j,m) momax = m endif if (obsclim(i,j,m) .lt. ocentmin) then ocentmin = obsclim(i,j,m) momin = m endif if (centmon(i,j,m) .gt. centmax) then centmax = centmon(i,j,m) mmax = m endif if (centmon(i,j,m) .lt. centmin) then centmin = centmon(i,j,m) mmin = m endif 390 continue cmax(i,j) = centmax + dcbias cmin(i,j) = centmin + dcbias omax(i,j) = ocentmax + dcbias omin(i,j) = ocentmin + dcbias wcseas(ik) = wcseas(ik) + sstwts(i,j) do 395 m=1,nmon12 acentmon(m,ik) = acentmon(m,ik) + & centmon(i,j,m)*sstwts(i,j) 395 continue c count ocean grid cells isea = isea + 1 diffmax = centmax - ocentmax diffmin = ocentmin - centmin if ((diffmax .gt. 700.0) .or. (diffmin .gt. 700.0)) then print*, 'exceeds 700 at lat = ', alats(j), & ' alon = ', alons(i) print*, 'kk = ', kk, ' obs = ', (obsclim(i,j,m), m=1,12) print*, 'kk = ', kk, ' mid = ', (centmon(i,j,m), m=1,12) endif ii = diffmax/histo if (ii .gt. 20) ii=20 if (ii .lt. -20) ii=-20 jj = 0.1*diffmax/histo if (jj .gt. 10) jj=10 if (jj .lt. -10) jj=-10 if (ii .gt. 0) then do 84 m=1,ii icount(ik,m) = icount(ik,m) + 1 84 continue elseif (ii .lt. 0) then do 85 m=-1,ii,-1 icount(ik,m) = icount(ik,m) + 1 85 continue else icount(ik,0) = icount(ik,0) + 1 endif if (jj .gt. 0) then do 86 m=1,jj jcount(ik,m) = jcount(ik,m) + 1 86 continue elseif (jj .lt. 0) then do 87 m=-1,jj,-1 jcount(ik,m) = jcount(ik,m) + 1 87 continue else jcount(ik,0) = jcount(ik,0) + 1 endif kk = diffmin/histo if (kk .gt. 20) kk=20 if (kk .lt. -20) kk=-20 mm = 0.1*diffmin/histo if (mm .gt. 10) mm=10 if (mm .lt. -10) mm=-10 if (kk .gt. 0) then do 94 m=1,kk kcount(ik,m) = kcount(ik,m) + 1 94 continue elseif (kk .lt. 0) then do 95 m=-1,kk,-1 kcount(ik,m) = kcount(ik,m) + 1 95 continue else kcount(ik,0) = kcount(ik,0) + 1 endif if (mm .gt. 0) then do 96 m=1,mm mcount(ik,m) = mcount(ik,m) + 1 96 continue elseif (mm .lt. 0) then do 97 m=-1,mm,-1 mcount(ik,m) = mcount(ik,m) + 1 97 continue else mcount(ik,0) = mcount(ik,0) + 1 endif endif 400 continue 500 continue return end subroutine calcfull(j1, jn, nlon, nlat, mlat, nmon, & iyr1rd, mon1rd, iyr1, mon1, iyrn, monn, iyrnrd, monnrd, & mon1clm, monnclm, iyr1clm, iyrnclm, ismf, jsmf, & nzon, nzonp, nmon12, lagcalc, nlagm, & maxiter, wtmin, amiss, tmin, tmax, & dt, conv, bbmin, varmin, lagm, sstwts, wtfrac, & wcseas, cvarmon, acvarmon, correl, ccorrel, & obsclim, vecin, vecout, centmon, sst, alons, alats, a, c, & ovarmon, wcvar, wccorr) implicit none integer nlon, nlat, nmon, iyr1rd, mon1rd, iyr1, mon1, iyrn, monn, & mon1clm, monnclm, iyr1clm, iyrnclm, ismf, jsmf, & iyrnrd, monnrd, nzon, nzonp, nmon12, lagcalc, nlagm, & maxiter, mlat, ji, nprior, nafter, j1, jn real wtmin, amiss, tmin, tmax, dt, conv, bbmin, varmin integer lagm(nlagm) real sstwts(nlon,nlat), wtfrac(nlon,nlat), wcseas(nzon), & cvarmon(nlon,nlat), acvarmon(nzonp), & correl(nmon12), ccorrel(nlagm,nzonp), & obsclim(nlon,nlat,nmon12), vecin(nmon), vecout(nmon), & centmon(nlon,nlat,nmon12), sst(nlon,mlat,nmon), alons(nlon), & alats(nlat), a(nmon), c(nmon), ovarmon(nlon,nlat), & wcvar(nzonp), wccorr(nzonp) integer i, j, k, m, mm, ik, nn, n, m3, m4 real tt, s, var, cc double precision avg1, avg2, covar, sum m3 = nmon12*(iyr1clm-iyr1) + mon1clm-mon1 + 1 m4 = nmon12*(iyrnclm-iyr1) + monnclm-mon1 + 1 do 900 j=j1,jn print*, 'Calculating mid-month values for latitude ', alats(j) ik = (j-1)/(nlat/nzon) + 1 ik = min0(ik, nzon) ji = j-j1+1 do 800 i=1,nlon if (wtfrac(i,j) .le. wtmin) then cvarmon(i,j) = amiss else c fill in some months prior to the 1st observed month and c after the last observed month nprior = nmon12*(iyr1rd-iyr1) + mon1rd - mon1 if (nprior .gt. 0) then do 708 m=1,nprior if ((nprior+1-m) .gt. nmon12) then cc = 0.0 else cc = correl(nprior+1-m) endif mm = mod((mon1-2+m), nmon12) + 1 tt = sst(i,ji,nprior+1)*cc if ((obsclim(i,j,mm)+tt) .gt. tmax) & tt = tmax-obsclim(i,j,mm) if (obsclim(i,j,mm) .ge. tmax) tt = 0.0 if ((obsclim(i,j,mm)+tt) .lt. tmin) & tt = tmin-obsclim(i,j,mm) if (obsclim(i,j,mm) .le. tmin) tt = 0.0 sst(i,ji,m) = tt 708 continue endif nafter = nmon12*(iyrn-iyrnrd) + monn - monnrd if (nafter .gt. 0) then do 709 m=1,nafter if (m .gt. nmon12) then cc = 0.0 else cc = correl(m) endif mm = mod((monnrd-1+m), nmon12) + 1 tt = sst(i,ji,nmon-nafter)*cc if ((obsclim(i,j,mm)+tt) .gt. tmax) & tt = tmax-obsclim(i,j,mm) if (obsclim(i,j,mm) .ge. tmax) tt = 0.0 if ((obsclim(i,j,mm)+tt) .lt. tmin) & tt = tmin-obsclim(i,j,mm) if (obsclim(i,j,mm) .le. tmin) tt = 0.0 sst(i,ji,nmon-nafter+m) = tt 709 continue endif c copy data into vecin, vecout do 710 n=1,nmon mm = mod((mon1+n-2), nmon12) + 1 vecin(n) = sst(i,ji,n) + obsclim(i,j,mm) vecout(n) = vecin(n) 710 continue call solvmid(alons(i), alats(j), nmon, conv, dt, tmin, tmax, & bbmin, maxiter, a, c, vecin, vecout, ismf, jsmf) sum = 0.0 avg1 = 0.0 do 718 n=m3, m4 mm = mod((mon1clm+n-m3-1), nmon12) + 1 s = amax1(tmin, amin1(tmax, vecout(n))) - & amax1(tmin, amin1(tmax, centmon(i,j,mm))) sum = sum + s*s avg1 = avg1 + s 718 continue cvarmon(i,j) = (sum-avg1*avg1/(m4-m3+1))/(m4-m3) do 720 n=1,nmon mm = mod((mon1+n-2), nmon12) + 1 sst(i,ji,n) = vecout(n) - centmon(i,j,mm) 720 continue if (ovarmon(i,j) .ge. varmin) then wcvar(ik) = wcvar(ik) + sstwts(i,j) acvarmon(ik) = acvarmon(ik) + cvarmon(i,j)*sstwts(i,j) wcvar(nzonp) = wcvar(nzonp) + sstwts(i,j) acvarmon(nzonp) = acvarmon(nzonp) + & cvarmon(i,j)*sstwts(i,j) endif endif 800 continue 900 continue if (lagcalc .gt. 0) then do 550 j=j1,jn print*, & 'calculating mid-month lag correlations for latitude ', & alats(j) ji = j-j1+1 ik = (j-1)/(nlat/nzon) + 1 ik = min0(ik, nzon) do 540 i=1,nlon if ((wtfrac(i,j) .gt. wtmin) .and. (ovarmon(i,j) .ge. varmin)) & then wccorr(ik) = wccorr(ik) + sstwts(i,j) wccorr(nzonp) = wccorr(nzonp) + sstwts(i,j) do 530 k=1,nlagm covar = 0.0 avg1 = 0.0 avg2 = 0.0 nn = m4 - lagm(k) do 520 n=m3,nn mm = n+lagm(k) covar = covar + sst(i,ji,n)*sst(i,ji,mm) avg1 = avg1 + sst(i,ji,n) avg2 = avg2 + sst(i,ji,mm) 520 continue if (k .eq. 1) then c calculate variance var = (covar - (avg1*avg2/(nn-m3+1)))/(nn-m3) ccorrel(k,ik) = ccorrel(k,ik) + sstwts(i,j)*var ccorrel(k,nzonp) = ccorrel(k,nzonp) + sstwts(i,j)*var elseif (var .gt. 0.) then c calculate correlations ccorrel(k,ik) = ccorrel(k,ik) + sstwts(i,j)* & ((covar - (avg1*avg2/(nn-m3+1)))/(nn-m3))/var ccorrel(k,nzonp) = ccorrel(k,nzonp) + sstwts(i,j)* & ((covar - (avg1*avg2/(nn-m3+1)))/(nn-m3))/var endif 530 continue endif 540 continue 550 continue endif return end subroutine solvmid(alon, alat, nmon, conv, dt, tmin, tmax, & bbmin, maxiter, a, c, obsmean, ss, icnt, jcnt) implicit none integer nmont, nmon12 parameter (nmont=12*130, nmon12=12) C ^^^ C is this correct???????????????????????????????????????? integer nmon, icnt, jcnt, maxiter real conv, tmin, tmax, dt, bbmin, alon, alat real obsmean(nmon), a(nmon), c(nmon), ss(nmon) integer i, n, imethod, n1, n2, nn, jj, i1, i2, i3, nnn, jend, & kk, j, k, kkk, nm, np integer jbeg(nmont) real relax, residmax, resid, dxm, dxp, s1, s2, addmax, addmin real r(nmont), avg(nmont), aa(nmont), bb(nmont), cc(nmont), & add(nmont) double precision s(nmont), sum imethod = 1 c check following value relax = 1.0 if (nmon .gt. nmont) then print*, 'error-- nmont not declared large enough in ' print*, 'subroutine solvmid' stop endif c check for occurance where obs monthly means are consecutively c at upper and lower limits. If so, smooth data, being careful c to preserve annual mean. do 48 n=1,nmon add(n) = 0.0 48 continue n2 = nmon do 50 n=1,nmon n1 = n2 n2 = n if ((obsmean(n2)-obsmean(n1)) .gt. (tmax-tmin-2.*dt)) then add(n1) = add(n1) + & (0.5*(obsmean(n2)-obsmean(n1)-tmax+tmin) + dt)/c(n1) add(n2) = add(n2) - & (0.5*(obsmean(n2)-obsmean(n1)-tmax+tmin) + dt)/a(n2) elseif ((obsmean(n1)-obsmean(n2)) .gt. (tmax-tmin-2.*dt)) then add(n1) = add(n1) - & (0.5*(obsmean(n1)-obsmean(n2)-tmax+tmin) + dt)/c(n1) add(n2) = add(n2) + & (0.5*(obsmean(n1)-obsmean(n2)-tmax+tmin) + dt)/a(n2) endif 50 continue nn = 0 addmax = 0.0 addmin = 0.0 do 51 n=1,nmon if (add(n) .ne. 0.0) then c print*, c & 'monthly means go from one limit to the other in 1 month' c print*, 'alat = ', alat, ' alon = ', alon, ' n = ', n, c & ' add = ', add(n) addmax = amax1(addmax, add(n)) addmin = amin1(addmin, add(n)) obsmean(n) = obsmean(n) + add(n) nn = nn + 1 endif 51 continue icnt = icnt + nn if (nn .gt. 0) then jcnt = jcnt + 1 c if (nmon .eq. nmon12) then c print*, 'Climatology: ' c write(9,*) 'Climatology: ' c endif if (jcnt .le. 200) then print*, nn, & ' monthly values smoothed at lat= ', alat, ' lon= ', alon print*, 'max added = ', addmax, & ' max subtracted = ', addmin write(9,*) nn, & ' monthly values smoothed at lat= ', alat, ' lon= ', alon write(9,*) 'max added = ', addmax, & ' max subtracted = ', addmin endif if (jcnt .eq. 200) then print*, ' ' if (nmon .eq. nmon12) then print*, & 'No more warnings will be printed concerning smoothing '// & 'of climatological data' else print*, & 'No more warnings will be printed concerning smoothing '// & 'of monthly data' endif print*, ' ' write(9,*) ' ' if (nmon .eq. nmon12) then write(9,*) & 'No more warnings will be printed concerning smoothing '// & 'of climatological data' else write(9,*) & 'No more warnings will be printed concerning smoothing '// & 'of monthly data' endif write(9,*) ' ' endif endif c check if all are le tmin or all are ge tmax if (obsmean(1) .le. (tmin+0.01*dt)) then do 80 i=2,nmon if (obsmean(i) .gt. (tmin+0.01*dt)) go to 99 80 continue do 85 i=1,nmon ss(i) = tmin 85 continue c if (nmon .eq. nmon12) print*, 'Climatology: ' c print*, 'all values were at minimum at this grid cell:' c print*, 'latitude = ', alat, ' longitude = ', alon return elseif (obsmean(1) .ge. (tmax-0.01*dt)) then do 90 i=2,nmon if (obsmean(i) .lt. (tmax-0.01*dt)) go to 99 90 continue do 95 i=1,nmon ss(i) = tmax 95 continue c if (nmon .eq. nmon12) print*, 'Climatology: ' c print*, 'all values were at maximum at this grid cell:' c print*, 'latitude = ', alat, ' longitude = ', alon return endif 99 jj = 0 do 100 i=1,nmon i1 = i i2 = mod(i,nmon) + 1 i3 = mod((i+1), nmon) + 1 if (((obsmean(i1) .le. tmin+0.01*dt) .and. & (obsmean(i2) .le. tmin+0.01*dt) .and. & (obsmean(i3) .gt. tmin+0.01*dt)) .or. & ((obsmean(i1) .ge. tmax-0.01*dt) .and. & (obsmean(i2) .ge. tmax-0.01*dt) .and. & (obsmean(i3) .lt. tmax-0.01*dt))) then jj = jj + 1 jbeg(jj) = i2 endif 100 continue if (jj .eq. 0) then c simple cyclic treatment c latest approximation of means (given mid-month values) nnn = 0 105 nnn = nnn + 1 sum = 0.0 residmax = 0.0 do 110 n = 1, nmon nm = mod((n+nmon-2), nmon) + 1 np = mod(n, nmon) + 1 bb(n) = 0.0 avg(n) = 0.0 if (nnn .lt. imethod) then call approx(tmin, tmax, a(n), c(n), ss(nm), ss(n), & ss(np), aa(n), bb(n), cc(n), avg(n)) else call numer(conv, tmin, tmax, bbmin, a(n), c(n), ss(nm), & ss(n), ss(np), aa(n), bb(n), cc(n), avg(n)) endif r(n) = obsmean(n) - avg(n) sum = sum + r(n)**2 residmax = amax1(residmax, abs(r(n))) 110 continue resid = dsqrt(sum)/nmon if (residmax .gt. conv) then if (nnn .gt. maxiter*0.5) then print*, 'iteration = ', nnn, ' residual = ', resid, & ' maximum residual = ', residmax endif if (nnn .gt. maxiter*0.9) then print*, ' ' print*, 'latitude = ', alat, ' longitude = ', alon do 115 n=1,nmon write(*,'(8(1pe10.2))') obsmean(n), avg(n), r(n), & s(n), ss(n), aa(n), bb(n), cc(n) 115 continue endif if (nnn .gt. maxiter) then print*, 'latitude = ', alat, ' longitude = ', alon print*, 'does not converge' write(9,*) 'latitude = ', alat, ' longitude = ', alon write(9,*) 'does not converge' stop else c solve for new estimate of mid-month values call cyclic(alon, alat, aa, bb, cc, cc(nmon), aa(1), r, s, & nmon) do 120 n=1,nmon ss(n) = ss(n) + relax*s(n) 120 continue c if ss exceeds tmax or tmin, then it should exceed it no c more than absolutely necessary: do 130 n=1,nmon nm = mod((n+nmon-2), nmon) + 1 np = mod(n, nmon) + 1 if (ss(n) .gt. tmax) then if (ss(nm) .le. tmax) then dxm = (ss(n)-tmax)/((ss(n)-ss(nm))*a(n)) else dxm = 0.0 endif if (ss(np) .le. tmax) then dxp = (ss(n)-tmax)/((ss(n)-ss(np))*c(n)) else dxp = 0.0 endif if ((dxm .gt. 0.5) .and. (dxp .gt. 0.5)) then s1 = tmax + (tmax-ss(nm))*a(n)/(2.-a(n)) s2 = tmax + (tmax-ss(np))*c(n)/(2.-c(n)) ss(n) = amin1(s1, s2) elseif ((dxm .eq. 0.0) .and. (dxp .eq. 0.0)) then ss(n) = tmax elseif ((dxp .eq. 0.0) .and. (dxm .gt. 0.5)) then ss(n) = tmax + (tmax-ss(nm))*a(n)/(2.-a(n)) elseif ((dxm .eq. 0.0) .and. (dxp .gt. 0.5)) then ss(n) = tmax + (tmax-ss(np))*c(n)/(2.-c(n)) endif elseif (ss(n) .lt. tmin) then if (ss(nm) .ge. tmin) then dxm = (ss(n)-tmin)/((ss(n)-ss(nm))*a(n)) else dxm = 0.0 endif if (ss(np) .ge. tmin) then dxp = (ss(n)-tmin)/((ss(n)-ss(np))*c(n)) else dxp = 0.0 endif if ((dxm .gt. 0.5) .and. (dxp .gt. 0.5)) then s1 = tmin + (tmin-ss(nm))*a(n)/(2.-a(n)) s2 = tmin + (tmin-ss(np))*c(n)/(2.-c(n)) ss(n) = amin1(s1, s2) elseif ((dxm .eq. 0.0) .and. (dxp .eq. 0.0)) then ss(n) = tmin elseif ((dxp .eq. 0.0) .and. (dxm .gt. 0.5)) then ss(n) = tmin + (tmin-ss(nm))*a(n)/(2.-a(n)) elseif ((dxm .eq. 0.0) .and. (dxp .gt. 0.5)) then ss(n) = tmin + (tmin-ss(np))*c(n)/(2.-c(n)) endif endif 130 continue go to 105 endif endif else c treat independent segments do 300 j=1,jj jend = jbeg(j) 150 jend = jend + 1 i1 = mod((jend-1), nmon) + 1 i2 = mod(jend, nmon) + 1 if (((obsmean(i1) .le. tmin+0.01*dt) .and. & (obsmean(i2) .le. tmin+0.01*dt)) .or. & ((obsmean(i1) .ge. tmax-0.01*dt) .and. & (obsmean(i2) .ge. tmax-0.01*dt))) then c calculate values for interval jbeg(j) to jend c c latest approximation of means (given mid-month values) nnn = 0 205 nnn = nnn + 1 kk = jend - jbeg(j) + 1 n = jbeg(j) avg(1) = obsmean(n) r(1) = 0.0 n = mod((jend-1), nmon) + 1 avg(kk) = obsmean(n) r(kk) = 0.0 sum = 0.0 residmax = 0.0 do 210 k = 2, kk-1 nm = mod((k+jbeg(j)-3), nmon) + 1 n = mod((k+jbeg(j)-2), nmon) + 1 np = mod((k+jbeg(j)-1), nmon) + 1 bb(k) = 0.0 avg(k) = 0.0 if (nnn .lt. imethod) then call approx(tmin, tmax, a(n), c(n), ss(nm), ss(n), & ss(np), aa(k), bb(k), cc(k), avg(k)) else call numer(conv, tmin, tmax, bbmin, a(n), c(n), ss(nm), & ss(n), ss(np), aa(k), bb(k), cc(k), avg(k)) endif r(k) = obsmean(n) - avg(k) sum = sum + r(k)**2 residmax = amax1(residmax, abs(r(k))) 210 continue resid = dsqrt(sum)/(kk-2) if (residmax .gt. conv) then if (nnn .gt. maxiter*0.5) then print*, 'iter = ', nnn, ' kk = ', kk, ' residual = ', resid, & ' maximum residual = ', residmax endif if (nnn .gt. maxiter*0.9) then print*, ' ' print*, 'latitude = ', alat, ' longitude = ', alon do 215 k=1,kk n = mod((k+jbeg(j)-2), nmon) + 1 write(*,'(8(1pe10.2))') obsmean(n), avg(k), r(k), & s(k), ss(n), aa(k), bb(k), cc(k) 215 continue endif if (nnn .gt. maxiter) then print*, 'latitude = ', alat, ' longitude = ', alon print*, 'does not converge' write(9,*) 'latitude = ', alat, ' longitude = ', alon write(9,*) 'does not converge' stop else c solve for new estimate of mid-month values kkk = kk - 2 call tridag(alon, alat, aa(2), bb(2), cc(2), r(2), & s(2), kkk) do 220 k=2,kk-1 n = mod((k+jbeg(j)-2), nmon) + 1 ss(n) = ss(n) + relax*s(k) 220 continue c if ss exceeds tmax or tmin, then it should exceed it no c more than absolutely necessary: n = mod((jbeg(j)-1), nmon) + 1 np = mod(jbeg(j), nmon) + 1 if (obsmean(n) .ge. (tmax-0.01*dt)) then ss(n) = & amax1(tmax, (tmax + (tmax-ss(np))*c(n)/(2.-c(n)))) else ss(n) = & amin1(tmin, (tmin + (tmin-ss(np))*c(n)/(2.-c(n)))) endif nm = mod((jend+nmon-2), nmon) + 1 n = mod((jend-1), nmon) + 1 if (obsmean(n) .ge. (tmax-0.01*dt)) then ss(n) = & amax1(tmax, (tmax + (tmax-ss(nm))*a(n)/(2.-a(n)))) else ss(n) = & amin1(tmin, (tmin + (tmin-ss(nm))*a(n)/(2.-a(n)))) endif do 230 k=2,kk-1 nm = mod((k+jbeg(j)+nmon-3), nmon) + 1 n = mod((k+jbeg(j)-2), nmon) + 1 np = mod((k+jbeg(j)-1), nmon) + 1 if (ss(n) .gt. tmax) then if (ss(nm) .le. tmax) then dxm = (ss(n)-tmax)/((ss(n)-ss(nm))*a(n)) else dxm = 0.0 endif if (ss(np) .le. tmax) then dxp = (ss(n)-tmax)/((ss(n)-ss(np))*c(n)) else dxp = 0.0 endif if ((dxm .gt. 0.5) .and. (dxp .gt. 0.5)) then s1 = tmax + (tmax-ss(nm))*a(n)/(2.-a(n)) s2 = tmax + (tmax-ss(np))*c(n)/(2.-c(n)) ss(n) = amin1(s1, s2) elseif ((dxm .eq. 0.0) .and. (dxp .eq. 0.0)) then ss(n) = tmax elseif ((dxp .eq. 0.0) .and. (dxm .gt. 0.5)) then ss(n) = tmax + (tmax-ss(nm))*a(n)/(2.-a(n)) elseif ((dxm .eq. 0.0) .and. (dxp .gt. 0.5)) then ss(n) = tmax + (tmax-ss(np))*c(n)/(2.-c(n)) endif elseif (ss(n) .lt. tmin) then if (ss(nm) .ge. tmin) then dxm = (ss(n)-tmin)/((ss(n)-ss(nm))*a(n)) else dxm = 0.0 endif if (ss(np) .ge. tmin) then dxp = (ss(n)-tmin)/((ss(n)-ss(np))*c(n)) else dxp = 0.0 endif if ((dxm .gt. 0.5) .and. (dxp .gt. 0.5)) then s1 = tmin + (tmin-ss(nm))*a(n)/(2.-a(n)) s2 = tmin + (tmin-ss(np))*c(n)/(2.-c(n)) ss(n) = amax1(s1, s2) elseif ((dxm .eq. 0.0) .and. (dxp .eq. 0.0)) then ss(n) = tmin elseif ((dxp .eq. 0.0) .and. (dxm .gt. 0.5)) then ss(n) = tmin + (tmin-ss(nm))*a(n)/(2.-a(n)) elseif ((dxm .eq. 0.0) .and. (dxp .gt. 0.5)) then ss(n) = tmin + (tmin-ss(np))*c(n)/(2.-c(n)) endif endif 230 continue go to 205 endif endif go to 300 else go to 150 endif if ((ss(kk) .gt. 700.0) .or. (ss(kk) .lt. -700.0)) then print*, 'exceeds 700 at lat = ', alat, & ' alon = ', alon print*, 'kk = ', kk, ' obs = ', obsmean(kk-1), obsmean(kk), & obsmean(kk+1) print*, 'kk = ', kk, ' ss = ', ss(kk-1), ss(kk), & ss(kk+1) endif 300 continue c fill in values where consecutive means are outside limits do 250 i=1,nmon i1 = mod((i-2+nmon), nmon) + 1 i2 = mod((i-1), nmon) + 1 i3 = mod(i, nmon) + 1 if ((obsmean(i1) .le. (tmin+0.01*dt)) .and. & (obsmean(i2) .le. (tmin+0.01*dt)) .and. & (obsmean(i3) .le. (tmin+0.01*dt))) then ss(i2) = tmin elseif ((obsmean(i1) .ge. (tmax-0.01*dt)) .and. & (obsmean(i2) .ge. (tmax-0.01*dt)) .and. & (obsmean(i3) .ge. (tmax-0.01*dt))) then ss(i2) = tmax endif 250 continue endif return end subroutine numer(conv, tmin, tmax, bbmin, a, c, ssm, ss, ssp, aa, & bb, cc, avg) c ********************************************************************* implicit none real conv, tmin, tmax, bbmin, a, c, ssm, ss, ssp, aa, bb, cc, avg real ssmm, ssmp, sssm, sssp, sspm, sspp, r real amean avg = amean(tmin,tmax,a,c,ssm,ss,ssp) ssmm = ssm - conv ssmp = ssm + conv sssm = ss - conv sssp = ss + conv sspm = ssp - conv sspp = ssp + conv aa = (amean(tmin,tmax,a,c,ssmp,ss,ssp) - & amean(tmin,tmax,a,c,ssmm,ss,ssp)) / (2.*conv) bb = (amean(tmin,tmax,a,c,ssm,sssp,ssp) - & amean(tmin,tmax,a,c,ssm,sssm,ssp)) / (2.*conv) cc = (amean(tmin,tmax,a,c,ssm,ss,sspp) - & amean(tmin,tmax,a,c,ssm,ss,sspm)) / (2.*conv) aa = amin1(aa, bb) cc = amin1(cc, bb) if (bb .lt. bbmin) then bb = bbmin r = 0.2*bbmin aa = amax1(r, aa) cc = amax1(r, cc) endif return end function amean(tmin, tmax, a, c, ssm, ss, ssp) c ********************************************************************* implicit none real amean real tmin, tmax, a, c, ssm, ss, ssp real dx, dy, avg avg = 0.0 if (ss .le. tmin) then if (ssm .le. tmin) then avg = avg + tmin*0.5 elseif (ssm .ge. tmax) then dx = (ss-tmin)/((ss-ssm)*a) dy = (ss-tmax)/((ss-ssm)*a) if (dx .ge. 0.5) then avg = avg + tmin*0.5 elseif (dy .le. 0.5) then avg = avg + tmin*dx + tmax*(.5-dy) + (dy-dx)*.5*(tmin+tmax) else avg = avg + & tmin*dx + (0.5-dx)*0.5*(tmin + ss - 0.5*a*(ss-ssm)) endif else dx = (ss-tmin)/((ss-ssm)*a) if (dx .ge. 0.5) then avg = avg + tmin*0.5 else avg = avg + & tmin*dx + (0.5-dx)*0.5*(tmin + ss - 0.5*a*(ss-ssm)) endif endif if (ssp .le. tmin) then avg = avg + tmin*0.5 elseif (ssp .ge. tmax) then dx = (ss-tmin)/((ss-ssp)*c) dy = (ss-tmax)/((ss-ssp)*c) if (dx .ge. 0.5) then avg = avg + tmin*0.5 elseif (dy .le. 0.5) then avg = avg + tmin*dx + tmax*(.5-dy) + (dy-dx)*.5*(tmin+tmax) else avg = avg + & tmin*dx + (0.5-dx)*0.5*(tmin + ss - 0.5*c*(ss-ssp)) endif else dx = (ss-tmin)/((ss-ssp)*c) if (dx .ge. 0.5) then avg = avg + tmin*0.5 else avg = avg + & tmin*dx + (0.5-dx)*0.5*(tmin + ss - 0.5*c*(ss-ssp)) endif endif elseif (ss .ge. tmax) then if (ssm .ge. tmax) then avg = avg + tmax*0.5 elseif (ssm .le. tmin) then dx = (ss-tmax)/((ss-ssm)*a) dy = (ss-tmin)/((ss-ssm)*a) if (dx .ge. 0.5) then avg = avg + tmax*0.5 elseif (dy .le. 0.5) then avg = avg + tmax*dx + tmin*(.5-dy) + (dy-dx)*.5*(tmin+tmax) else avg = avg + & tmax*dx + (0.5-dx)*0.5*(tmax + ss - 0.5*a*(ss-ssm)) endif else dx = (ss-tmax)/((ss-ssm)*a) if (dx .ge. 0.5) then avg = avg + tmax*0.5 else avg = avg + & tmax*dx + (0.5-dx)*0.5*(tmax + ss - 0.5*a*(ss-ssm)) endif endif if (ssp .ge. tmax) then avg = avg + tmax*0.5 elseif (ssp .le. tmin) then dx = (ss-tmax)/((ss-ssp)*c) dy = (ss-tmin)/((ss-ssp)*c) if (dx .ge. 0.5) then avg = avg + tmax*0.5 elseif (dy .le. 0.5) then avg = avg + tmax*dx + tmin*(.5-dy) + (dy-dx)*.5*(tmin+tmax) else avg = avg + & tmax*dx + (0.5-dx)*0.5*(tmax + ss - 0.5*c*(ss-ssp)) endif else dx = (ss-tmax)/((ss-ssp)*c) if (dx .ge. 0.5) then avg = avg + tmax*0.5 else avg = avg + & tmax*dx + (0.5-dx)*0.5*(tmax + ss - 0.5*c*(ss-ssp)) endif endif else if (ssm .le. tmin) then dx = (ss-tmin)/((ss-ssm)*a) if (dx .ge. 0.5) then avg = avg + 0.5*0.5*(2.*ss - 0.5*(ss-ssm)*a) else avg = avg + tmin*(.5-dx) + dx*0.5*(tmin + ss) endif elseif (ssm .ge. tmax) then dx = (ss-tmax)/((ss-ssm)*a) if (dx .ge. 0.5) then avg = avg + 0.5*0.5*(2.*ss - 0.5*(ss-ssm)*a) else avg = avg + tmax*(.5-dx) + dx*0.5*(tmax + ss) endif else avg = avg + 0.5*0.5*(2.*ss - 0.5*(ss-ssm)*a) endif if (ssp .le. tmin) then dx = (ss-tmin)/((ss-ssp)*c) if (dx .ge. 0.5) then avg = avg + 0.5*0.5*(2.*ss - 0.5*(ss-ssp)*c) else avg = avg + tmin*(.5-dx) + dx*0.5*(tmin + ss) endif elseif (ssp .ge. tmax) then dx = (ss-tmax)/((ss-ssp)*c) if (dx .ge. 0.5) then avg = avg + 0.5*0.5*(2.*ss - 0.5*(ss-ssp)*c) else avg = avg + tmax*(.5-dx) + dx*0.5*(tmax + ss) endif else avg = avg + 0.5*0.5*(2.*ss - 0.5*(ss-ssp)*c) endif endif amean = avg return end subroutine approx(tmin, tmax, a, c, ssm, ss, ssp, aa, bb, cc, avg) c ********************************************************************* implicit none real tmin, tmax, a, c, ssm, ss, ssp, aa, bb, cc, avg real dx, dy if (ss .le. tmin) then if (ssm .le. tmin) then avg = avg + tmin*0.5 aa = a/32. bb = bb + 0.125 - a/32. elseif (ssm .ge. tmax) then dx = (ss-tmin)/((ss-ssm)*a) dy = (ss-tmax)/((ss-ssm)*a) if (dx .ge. 0.5) then avg = avg + tmin*0.5 aa = a/32. bb = bb + 0.125 - a/32. elseif (dy .le. 0.5) then avg = avg + tmin*dx + tmax*(.5-dy) + (dy-dx)*.5*(tmin+tmax) aa = a/16. bb = bb + 0.25 - a/16. else avg = avg + & tmin*dx + (0.5-dx)*0.5*(tmin + ss - 0.5*a*(ss-ssm)) aa = a/16. bb = bb + 0.25 - a/16. endif else dx = (ss-tmin)/((ss-ssm)*a) if (dx .ge. 0.5) then avg = avg + tmin*0.5 aa = a/32. bb = bb + 0.125 - a/32. else avg = avg + & tmin*dx + (0.5-dx)*0.5*(tmin + ss - 0.5*a*(ss-ssm)) aa = a/16. bb = bb + 0.25 - a/16. endif endif if (ssp .le. tmin) then avg = avg + tmin*0.5 cc = c/32. bb = bb + 0.125 - c/32. elseif (ssp .ge. tmax) then dx = (ss-tmin)/((ss-ssp)*c) dy = (ss-tmax)/((ss-ssp)*c) if (dx .ge. 0.5) then avg = avg + tmin*0.5 cc = c/32. bb = bb + 0.125 - c/32. elseif (dy .le. 0.5) then avg = avg + tmin*dx + tmax*(.5-dy) + (dy-dx)*.5*(tmin+tmax) cc = c/16. bb = bb + 0.25 - c/16. else avg = avg + & tmin*dx + (0.5-dx)*0.5*(tmin + ss - 0.5*c*(ss-ssp)) cc = c/16. bb = bb + 0.25 - c/16. endif else dx = (ss-tmin)/((ss-ssp)*c) if (dx .ge. 0.5) then avg = avg + tmin*0.5 cc = c/32. bb = bb + 0.125 - c/32. else avg = avg + & tmin*dx + (0.5-dx)*0.5*(tmin + ss - 0.5*c*(ss-ssp)) cc = c/16. bb = bb + 0.25 - c/16. endif endif elseif (ss .ge. tmax) then if (ssm .ge. tmax) then avg = avg + tmax*0.5 aa = a/32. bb = bb + 0.125 - a/32. elseif (ssm .le. tmin) then dx = (ss-tmax)/((ss-ssm)*a) dy = (ss-tmin)/((ss-ssm)*a) if (dx .ge. 0.5) then avg = avg + tmax*0.5 aa = a/32. bb = bb + 0.125 - a/32. elseif (dy .le. 0.5) then avg = avg + tmax*dx + tmin*(.5-dy) + (dy-dx)*.5*(tmin+tmax) aa = a/16. bb = bb + 0.25 - a/16. else avg = avg + & tmax*dx + (0.5-dx)*0.5*(tmax + ss - 0.5*a*(ss-ssm)) aa = a/16. bb = bb + 0.25 - a/16. endif else dx = (ss-tmax)/((ss-ssm)*a) if (dx .ge. 0.5) then avg = avg + tmax*0.5 aa = a/32. bb = bb + 0.125 - a/32. else avg = avg + & tmax*dx + (0.5-dx)*0.5*(tmax + ss - 0.5*a*(ss-ssm)) aa = a/16. bb = bb + 0.25 - a/16. endif endif if (ssp .ge. tmax) then avg = avg + tmax*0.5 cc = c/32. bb = bb + 0.125 - c/32. elseif (ssp .le. tmin) then dx = (ss-tmax)/((ss-ssp)*c) dy = (ss-tmin)/((ss-ssp)*c) if (dx .ge. 0.5) then cc = c/32. bb = bb + 0.125 - c/32. avg = avg + tmax*0.5 elseif (dy .le. 0.5) then avg = avg + tmax*dx + tmin*(.5-dy) + (dy-dx)*.5*(tmin+tmax) cc = c/16. bb = bb + 0.25 - c/16. else avg = avg + & tmax*dx + (0.5-dx)*0.5*(tmax + ss - 0.5*c*(ss-ssp)) cc = c/16. bb = bb + 0.25 - c/16. endif else dx = (ss-tmax)/((ss-ssp)*c) if (dx .ge. 0.5) then avg = avg + tmax*0.5 cc = c/32. bb = bb + 0.125 - c/32. else avg = avg + & tmax*dx + (0.5-dx)*0.5*(tmax + ss - 0.5*c*(ss-ssp)) cc = c/16. bb = bb + 0.25 - c/16. endif endif else if (ssm .le. tmin) then dx = (ss-tmin)/((ss-ssm)*a) if (dx .ge. 0.5) then avg = avg + 0.5*0.5*(2.*ss - 0.5*(ss-ssm)*a) aa = a/16. bb = bb + 0.25 - a/16. else avg = avg + tmin*(.5-dx) + dx*0.5*(tmin + ss) aa = a/16. bb = bb + 0.25 - a/16. endif elseif (ssm .ge. tmax) then dx = (ss-tmax)/((ss-ssm)*a) if (dx .ge. 0.5) then avg = avg + 0.5*0.5*(2.*ss - 0.5*(ss-ssm)*a) aa = a/16. bb = bb + 0.25 - a/16. else avg = avg + tmax*(.5-dx) + dx*0.5*(tmax + ss) aa = a/16. bb = bb + 0.25 - a/16. endif else avg = avg + 0.5*0.5*(2.*ss - 0.5*(ss-ssm)*a) aa = a/8. bb = bb + 0.5 - a/8. endif if (ssp .le. tmin) then dx = (ss-tmin)/((ss-ssp)*c) if (dx .ge. 0.5) then avg = avg + 0.5*0.5*(2.*ss - 0.5*(ss-ssp)*c) cc = c/16. bb = bb + 0.25 - c/16. else avg = avg + tmin*(.5-dx) + dx*0.5*(tmin + ss) cc = c/16. bb = bb + 0.25 - c/16. endif elseif (ssp .ge. tmax) then dx = (ss-tmax)/((ss-ssp)*c) if (dx .ge. 0.5) then avg = avg + 0.5*0.5*(2.*ss - 0.5*(ss-ssp)*c) cc = c/16. bb = bb + 0.25 - c/16. else avg = avg + tmax*(.5-dx) + dx*0.5*(tmax + ss) cc = c/16. bb = bb + 0.25 - c/16. endif else avg = avg + 0.5*0.5*(2.*ss - 0.5*(ss-ssp)*c) cc = c/8. bb = bb + 0.5 - c/8. endif endif return end SUBROUTINE tridag(alon,alat,a,b,c,r,u,n) implicit none INTEGER n,nmax REAL alon,alat,a(n),b(n),c(n),r(n) double precision u(n) PARAMETER (nmax=12*130) c ^^^ number of years INTEGER j REAL bet, gam(nmax) if (nmax .lt. n) then print*, 'Error nmax not declared large enough' print*, 'in tridag' stop endif if(b(1).eq.0.) then print*, 'longitude = ', alon, ' latitude = ', alat stop 'tridag: rewrite equations' endif bet=b(1) u(1)=r(1)/bet if (n .gt. 1) then do 11 j=2,n gam(j)=c(j-1)/bet bet=b(j)-a(j)*gam(j) if(bet.eq.0.) then print*, 'longitude = ', alon, ' latitude = ', alat stop 'tridag failed' endif u(j)=(r(j)-a(j)*u(j-1))/bet 11 continue do 12 j=n-1,1,-1 u(j)=u(j)-gam(j+1)*u(j+1) 12 continue endif return END SUBROUTINE cyclic(alon,alat,a,b,c,alpha,beta,r,x,n) implicit none INTEGER n,nmax real alon,alat,alpha,beta,a(n),b(n),c(n),r(n) double precision x(n) PARAMETER (nmax=12*130) c ^^^ number of years CU USES tridag INTEGER i REAL fact,gamma,bb(nmax),u(nmax) double precision z(nmax) if(n.le.2) stop 'n too small in cyclic' if(n.gt.nmax) stop 'nmax too small in cyclic' gamma=-b(1) bb(1)=b(1)-gamma bb(n)=b(n)-alpha*beta/gamma do 11 i=2,n-1 bb(i)=b(i) 11 continue call tridag(alon,alat,a,bb,c,r,x,n) u(1)=gamma u(n)=alpha do 12 i=2,n-1 u(i)=0. 12 continue call tridag(alon,alat,a,bb,c,u,z,n) fact=(x(1)+beta*x(n)/gamma)/(1.+z(1)+beta*z(n)/gamma) do 13 i=1,n x(i)=x(i)-fact*z(i) 13 continue return END subroutine wrtlats(ldefgrid, ldefvar, lcreate, ldefmisc, lwrite, & lclose, lconcat, date, stat, varname, outfile, outftype, & parmtabl, gtype, calendar, varcomm, filecomm, model, center, & amiss, nlon, nlat, nchunks, alons, alats, nmon, mon1, iyr1, & ibasemon, ibaseyr, idvar, array, chnkname, jyr1out) implicit none integer ldefgrid, ldefvar, lcreate, ldefmisc, lwrite, lclose, & lconcat, nlon, nlat, nchunks, nmon, mon1, iyr1, ibasemon, & ibaseyr, idvar, jyr1out(*) real alons(*), alats(*), amiss, array(*) character*80 center, varcomm character*120 outfile, chnkname(nchunks), filecomm c character*256 parmtabl character*16 outftype, calendar, gtype, model character*(*) varname, stat, date, parmtabl print*, 'Your version of this code cannot write with the "lats"' print*, 'subroutine. Only pp format or ascii is allowed.' stop return end subroutine wrtdrs(lcreate) implicit none integer lcreate print*, 'Your version of the code is unable to write ' print*, 'DRS-formatted files.' stop return end subroutine wrtasc(outfile,nlon,nlat,nmon12,iyr1,mon1,iyrn, & monn,elem1,elemn,model,center,source,vname,title,units, & grid,amissout,amissasc,lconcat,nchunks,chnkname, & jyr1out, array) implicit none integer nlon, nlat, iyr1, iyrn, mon1, monn, nmon12, & lconcat, nchunks, mm, iendyr, jyr1out(*) real array(*),elem1(2), elemn(2), amissout, amissasc character*16 vname, model character*120 source, fullfile character*80 title, grid, center character*120 outfile, chnkname(nchunks) character*16 ayr character*40 units integer n, nn, nnn, imon, jmon, iyr, iend, n1, n2, ichunk, & m, lpp(45) real bpp(19) print*, 'Preparing to write ascii data ...' print*, title write(9,*) title nnn = 0 if (lconcat .gt. 0) then do 200 ichunk=1,nchunks OPEN (12+ichunk,FILE=chnkname(ichunk),STATUS='OLD', & FORM='UNFORMATTED') 200 continue endif mm = 0 do 1000 iyr = iyr1, iyrn if (iyr .eq. iyr1) then imon = mon1 else imon = 1 endif if (iyr .eq. iyrn) then jmon = monn else jmon = nmon12 endif if ((iyr .eq. iyr1) .or. (iyr .eq. jyr1out(mm+1))) then mm = mm + 1 if (iyr1 .eq. 0) then write(ayr, '(''clim.asc'')') elseif (iyr .eq. iyrn) then write(ayr, '(i4, ''.asc'')') iyr iendyr = iyrn elseif (jyr1out(mm+1) .eq. -999) then write(ayr, '(i4, ''_'', i4,''.asc'')') iyr, iyrn iendyr = iyrn elseif (jyr1out(mm) .eq. (jyr1out(mm+1)-1)) then write(ayr, '(i4,''.asc'')') iyr iendyr = iyr else iendyr = jyr1out(mm+1) - 1 iendyr = min0(iendyr, iyrn) write(ayr, '(i4, ''_'', i4,''.asc'')') iyr, iendyr endif iend = index(outfile, ' ') - 1 fullfile = outfile(1:iend)//'_'//ayr print*, 'Opening ', fullfile open(12, file=fullfile, status='new') write(12,'(a16 / a80 / a80 / a40 / a40 / a16, 2x, a16 / a80)') & vname, title, source(1:80), source(81:120), units, center, & model, grid write(12,'(i5, '' = number of longitudes'')') nlon write(12,'(i5, '' = number of latitudes'')') nlat write(12,'(f10.4, '' = first longitude'')') elem1(1) write(12,'(f10.4, '' = last longitude'')') elemn(1) write(12,'(f10.4, '' = first latitude'')') elem1(2) write(12,'(f10.4, '' = last latitude'')') elemn(2) if (iyr1 .eq. 0) then write(12,'('' 0 - 0 = climatology'')') else write(12,'(i5, '' - '', i5, '' = years written'')') iyr, & iendyr endif write(12,'(i5, & '' = first month (jan, feb, mar, ... = 1, 2, 3, ...)'')') imon write(12,'(i5, & '' = last month (jan, feb, mar, ... = 1, 2, 3, ...)'')') jmon write(12,'(1x)') endif if (lconcat .eq. 0) then nn = nnn + 1 nnn = nnn + nlon*nlat*(jmon-imon+1) else c loop over months; read pp format temporary files c Process each chunk for this month n2 = 0 do 400 m=imon,jmon do 300 ichunk = 1,nchunks c read chunk in pp-format read(12+ichunk)lpp,bpp n1 = n2+1 n2 = n1+lpp(15)-1 read(12+ichunk) (array(n), n=n1,n2) 300 continue 400 continue nn = 1 nnn = n2 endif do 100 n=nn,nnn if (array(n) .eq. amissout) array(n) = amissasc if ((array(n) .gt. 9999.999) .or. & (array(n) .lt. -999.999)) then print*, 'Error: output data contains numbers that' print*, 'are too large to express in f8.3 format' stop endif 100 continue write(12,'(10f8.3)') (array(n), n=nn,nnn) if ((iyr .eq. iyrn) .or. ((iyr+1) .eq. jyr1out(mm+1))) & close(12) 1000 continue if (lconcat .gt. 0) then do 600 ichunk=1,nchunks close(12+ichunk) 600 continue endif c return end SUBROUTINE readpp(filename, nlat, nlon, mlat, lat1, lats, iyr1rd, & mon1rd, iyr1in, mons, amissin, amissout, alons, alats, & array, amask, wtfrac, lpp, bpp, space) c Read in every month to month value from pp format file within c the given chunk c c assumes data stored from north to south c NB Uses unit number 10 c IN: c filename : Filename of pp format file to be read c nlat : total number of latitudes spanning globe c nlon : Size of input data longitude dimension c mlat : Size of input data latitude dimension c lat1 : Index of first latitude to be read c lats : Number of latitudes to be read c iyr1in(1): Year at beginning of period to be read c iyr1in(i): for i>1, 1st year of succeeding input files c iyr1rd : Year at beginning of period to be read c mon1rd : Month at beginning of period to be read c mons : Number of months to be read c amissin : Value that indicates missing data in file read c amissout : value assigned to missing data before returned c space : scratch space needed to read in data c OUT: c alons : Vector containing the values of each longitude c in array c alats : Vector containing the values of each latitude in c array c array : Chunk of data to be returned c amask : proportional to grid-cell area; 0. if missing c wtfrac : 0. or 1. if missing or not. c lpp : integer header c bpp : real header IMPLICIT NONE CHARACTER*120 filename(*) INTEGER nlat,lat1,lats,nlon,mlat,iyr1rd,mon1rd,iyr1in(*),mons REAL array(nlon,mlat,mons),amask(nlon,nlat),amissin, & amissout, wtfrac(nlon,nlat), space(nlon,nlat) REAL alons(nlon),alats(nlat) INTEGER lyear,lmon INTEGER i,j,cnt,ncnt,nn, ij c Part one of pp-header INTEGER lpp(45) c Part two of pp-header REAL bpp(19) REAL radlat1,radlat2,sindiff,pi PI=4.0*atan(1.) c 3.1415926 print*, 'Opening input file: ', filename(1) OPEN (10,FILE=filename(1),STATUS='OLD',FORM='UNFORMATTED') lyear = 0 lmon = 0 ij = 1 c input file counter nn = 1 c number of data points read (excluding missing data) ncnt = 0 c This maintains an index of the current month cnt = 1 c Repeat until mons months have been processed DO WHILE (cnt.LE.mons) if ((lyear .eq. (iyr1in(nn+1)-1)) .and. (lmon .eq. 12)) then ij = cnt - 1 print*, 'found cnt= ', ij,' year= ',lyear,' month=', lmon close(10) nn = nn + 1 print*, 'Opening input file: ', filename(nn) OPEN (10,FILE=filename(nn),STATUS='OLD',FORM='UNFORMATTED') ij = 0 endif c Read in next field from data file READ (10,END=999)lpp,bpp c Extract the month and year for this input data array lyear = lpp(1) lmon = lpp(2) c If the data is for the period that we want IF ((lyear.GT.iyr1rd).OR. & ((lyear.EQ.iyr1rd).AND.(lmon.GE.mon1rd))) THEN if (nlat .ne. lpp(18)) then print*, 'ERROR: For input data in pp format, the number of' print*, ' latitudes stored in input file should equal' print*, ' the number declared in the parameter statement' print*, ' number stored = ', lpp(18) print*, ' nlat = ', nlat stop endif if (nlon .ne. lpp(19)) then print*, 'ERROR: For input data in pp format, the number of' print*, ' longitudes stored in input file should equal' print*, ' the number declared in the parameter statement' print*, ' number stored = ', lpp(19) print*, ' nlon = ', nlon stop endif READ (10,END=999) space if ((ij .eq. 1) .or. (cnt .eq. mons)) & print*, 'found cnt= ', cnt,' year= ',lyear,' month=', lmon ij=0 c if (cnt .eq. 1) then DO j = lat1, lat1+lats-1 c Set vector of latitude values alats(j) = bpp(14) + j*bpp(15) c Find the difference in sines of north/south boundaries c of square (assumes that if data stored from north to c south then bpp(15)<0.) radlat1 = (2.0*pi/360.0) * (bpp(14) + (j-0.5)*bpp(15)) radlat2 = (2.0*pi/360.0) * (bpp(14) + (j+0.5)*bpp(15)) sindiff = SIN(radlat1) - SIN(radlat2) do i = 1, nlon amask(i,j) = abs(sindiff) / (2.*nlon) wtfrac(i,j) = 1.0 enddo enddo c Assign values to alons do i=1,nlon alons(i) = bpp(16) + i*bpp(17) enddo endif do j = lat1, lat1+lats-1 DO i = 1,nlon IF (space(i,j) .NE. amissin) THEN ncnt = ncnt + 1 array(i,j - lat1 + 1,cnt) = space(i,j) ELSE array(i,j - lat1 + 1,cnt) = amissout amask(i,j) = 0.0 wtfrac(i,j) = 0.0 ENDIF ENDDO ENDDO c Increment count of months used cnt = cnt + 1 else read(10) ENDIF ENDDO go to 1000 999 CONTINUE print*, 'Error in subroutine readpp:' print*, 'reached end of file before all data were read' print*, 'Check specification of months and years you ' print*, 'want to read and input file name.' print*, 'Also check specification of nlon and nlat' stop 1000 continue PRINT*, 'number of longitude x latitude x month data read' print*, ' (excluding missing data) = ', ncnt ncnt = ncnt/(cnt-1) PRINT*, 'number of months read = ', cnt-1 CLOSE (10) print*, 'returning to main program' RETURN END subroutine wrtpp(filename, iobs, nlon, mlat, alat1, lats, & iyear1, mon1, nmon, lpp, bpp, array, fac, amissout) c Write the adjusted data to a series of intermediate files, c each containing a chunk of data, which can later be adjusted c to pp-format c Uses logical unit number 11 c IN: c Filename : The name of the file to which the chunk will be c written c iobs : 1 if observed data; 0 otherwise c nlon,mlat : The dimensions of the data array c alat1 : First latitude of chunk c lats : Number of latitudes in each chunk c iyr1 : Year at which chunk begins c mon1 : Month at which chunk begins c nmon : Number of months in chunk c array : Chunk itself c fac : factor to divide output data by (EBH) c amissout : Missing data indicator (EBH) c OUT :No value is returned, but the chunk stored in array c is written to unit number 11, in pp format to c retain info about nlat/nlon coverage c IMPLICIT NONE INTEGER nlon,lats,iyear1,mon1,nmon,mlat,iobs,mon REAL array(nlon,mlat,nmon),alat1 CHARACTER*120 filename, outfile INTEGER ijk REAL bpp(*), bbpp(19), day, fac, amissout INTEGER lpp(*), llpp(45), i,j,iyr,imon, iyr1, imon1, iday1, & ihour1, iyr2, imon2, iday2, ihour2, lastyr, ij, iend integer monlen(12,2) c data monlen/31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31, c & 31, 29, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31/ c *** GISST/HadISST *** might assume months of equal length data monlen/30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, & 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30/ lastyr = (iyear1*12+mon1+nmon-2)/12 c Define pp-header DO i=1,19 bbpp(i)=bpp(i) ENDDO DO i=1,45 llpp(i)=lpp(i) ENDDO c Size of chunk llpp(15) = nlon*lats c Height of chunk llpp(18) = lats c ??? check the following c ??? should change for different calendar if (iobs .gt. 0) then llpp(25) = 128 if (iyear1 .eq. 0) llpp(13) = 32 if (iyear1 .ne. 0) llpp(13) = 22 else llpp(25) = 0 if (iyear1 .eq. 0) llpp(13) = 32 if (iyear1 .ne. 0) llpp(13) = 2 endif c Top of chunk bbpp(14) = alat1 - bbpp(15) iend = index(filename, ' ') - 1 outfile = filename(1:iend)//'.pp' print*, 'opening pp file: ' print*, outfile ijk = index(outfile, 'temp') if (ijk .eq. 0) then OPEN (11,FILE=outfile,STATUS='NEW',FORM='UNFORMATTED') else OPEN (11,FILE=outfile,STATUS='unknown',FORM='UNFORMATTED') endif c Initialise the month and year to the beginning of chunk iyr = iyear1 imon = mon1 - 1 DO mon = 1,nmon imon = imon + 1 IF (imon.EQ.13) THEN iyr = iyr + 1 imon = 1 ENDIF if (iyear1 .eq. 0) then c we have climatological data if (iobs .gt. 0) then c we have monthly mean data ihour1 = 0 iday1 = 1 imon1 = imon iyr1 = iyear1 ihour2 = 0 iday2 = 1 imon2 = imon + 1 iyr2 = lastyr if (imon2 .gt. 12) then imon2 = 1 iyr2 = iyr2 + 1 endif else c we have mid-month data day = monlen(imon,1)/2.0 + 1.01 iday1 = day if ((day - iday1) .gt. 0.1) then ihour1 = 12 else ihour1 = 0 endif iyr1 = iyear1 imon1 = imon ihour2 = ihour1 iday2 = iday1 imon2 = imon iyr2 = lastyr endif else c we don't have climatology if (iobs .gt. 0) then c we have monthly mean data ihour1 = 0 iday1 = 1 imon1 = imon iyr1 = iyr ihour2 = 0 iday2 = 1 imon2 = imon + 1 iyr2 = iyr if (imon2 .eq. 13) then imon2 = 1 iyr2 = iyr2 + 1 endif else c we have mid-month data c ??? generalize for Julian, 360-day, 365-day years ij = 1 if (((mod(iyr,4) .eq. 0) .and. (mod(iyr,100) .ne. 0)) .or. & (mod(iyr,400) .eq. 0)) ij = 2 day = monlen(imon,ij)/2.0 + 1.01 iday1 = day if ((day - iday1) .gt. 0.1) then ihour1 = 12 else ihour1 = 0 endif imon1 = imon iyr1 = iyr ihour2 = ihour1 iday2 = iday1 imon2 = imon1 iyr2 = iyr endif endif llpp(1) = iyr1 llpp(2) = imon1 llpp(3) = iday1 llpp(4) = ihour1 llpp(7) = iyr2 llpp(8) = imon2 llpp(9) = iday2 llpp(10) = ihour2 c Convert back to original units by dividing by fac DO i=1,nlon DO j=1,lats IF (array(i,j,mon).NE.amissout) then array(i,j,mon)=array(i,j,mon)/fac ENDIF ENDDO ENDDO c Write each month of chunk in pp-format WRITE(11) llpp, bbpp WRITE(11) ((array(i,j,mon), i=1,nlon), j=1,lats) ENDDO CLOSE (11) RETURN END SUBROUTINE ppconcat (chnkname, outfile, nmon12, nlon, nlat, & nchunks, iyr1out, mon1out, iyrnout, monnout, jyr1out, & aa) c Concatenate all the intermediary files back into a monthly c pp file c Note, it uses logical units 12 through 12+number of chunks c IN: c chnkname : The names of the files that each of the chunks c are stored in, in N-S order c outfile : Name of output file c nmon12 : Number of months in a year c nlon,nlat : Dimensions of output file c nchunks : The number of chunks c iyr1out : The first year for which output is written c mon1out : The first month for which output is written c iyrnout : The last year for which output is written c monnout : The last month for which output is written c jyr1out(*) : First year of each output file c aa(nlon,nlat): scratch space for reading and writing c OUT: c : No output is returned, but the chunks stored in c the files given by chnkname have been c concatenated together IMPLICIT NONE integer nmon12, nlon, nlat, nchunks, iyr1out, mon1out, & iyrnout, monnout, lpp(45), jyr1out(*) CHARACTER*120 chnkname(nchunks),outfile, fullfile real bpp(19), aa(nlon,nlat) character*16 ayr INTEGER imon, mm, iendyr, iyr, m, jmon, i, j, & i2, j1, j2, ichunk, iend real alat0 DO ichunk = 1,nchunks OPEN (12+ichunk,FILE=chnkname(ichunk),STATUS='OLD', & FORM='UNFORMATTED') ENDDO mm = 0 do 1000 iyr = iyr1out, iyrnout if (iyr .eq. iyr1out) then imon = mon1out else imon = 1 endif if (iyr .eq. iyrnout) then jmon = monnout else jmon = nmon12 endif if ((iyr .eq. iyr1out) .or. (iyr .eq. jyr1out(mm+1))) then mm = mm + 1 if (iyr1out .eq. 0) then write(ayr, '(''clim.pp'')') elseif (iyr .eq. iyrnout) then write(ayr, '(i4, ''.pp'')') iyr elseif (jyr1out(mm+1) .eq. -999) then write(ayr, '(i4, ''_'', i4,''.pp'')') iyr, iyrnout elseif (jyr1out(mm) .eq. (jyr1out(mm+1)-1)) then write(ayr, '(i4,''.pp'')') iyr elseif (mm .eq. 1 .and. jyr1out(mm) .eq. -999) then write(ayr, '(i4, ''_'', i4, ''.pp'')') iyr1out, iyrnout else iendyr = jyr1out(mm+1) - 1 iendyr = min0(iendyr, iyrnout) write(ayr, '(i4, ''_'', i4,''.pp'')') iyr, iendyr endif iend = index(outfile, ' ') - 1 fullfile = outfile(1:iend)//'_'//ayr print*, 'Opening pp file:' print*, fullfile OPEN (12,FILE=fullfile,STATUS='unknown',FORM='UNFORMATTED') endif c Concatenate each month consecutively do 500 m=imon,jmon c Process each chunk for this year j2 = 0 DO 300 ichunk = 1,nchunks c Read chunk in pp-format READ(12+ichunk)lpp,bpp if (ichunk .eq. 1) alat0 = bpp(14) i2 = lpp(19) j1 = j2 + 1 j2 = j1 + lpp(18) - 1 read(12+ichunk) ((aa(i,j), i=1,i2), j=j1,j2) 300 continue c Assign the correct pp-header lpp(18)= j2 lpp(15) = i2*j2 bpp(14)= alat0 c Write the regenerated monthly field to unit 12 WRITE(12)lpp,bpp WRITE(12) ((aa(i,j), i=1,i2), j=1,j2) 500 continue if ((iyr .eq. iyrnout) .or. ((iyr+1) .eq. jyr1out(mm+1))) & close(12) 1000 continue DO ichunk = 1,nchunks close(12+ichunk) ENDDO RETURN END subroutine getobs(model, gtype, msklndi, sftfilei, sftnamei, & msklndo, sftfileo, sftnameo, inputfil, varin, & nlon, nlat, mlat, lat1, latn, iyr1in, mon1in, iyr1rd, & mon1rd, iyrnrd, monnrd, nmon12, amissin, amissout, wtmin, & sst, sstwts, wtfrac, alons, alats) implicit none integer msklndi, msklndo, nlon, nlat, mlat, lat1, latn, iyr1in, & mon1in, iyr1rd, mon1rd, iyrnrd, monnrd, nmon12 real sst(nlon,mlat,*), alons(nlon), alats(*), & sstwts(nlon,*), wtfrac(nlon,*), amissin, amissout, wtmin character*120 sftfilei, sftfileo, inputfil character*16 sftnamei, sftnameo, varin, model, gtype print*, 'Your version of the code cannot read files with' print*, 'the EzGet software (i.e., files other than pp-format).' stop return end c ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c logical function caseindp c ------- c c Description: c ----------- c This function tests for character string equivalence c c Usage: c ------ c c if (caseindp(a, b)) then c c ------ c Date: 9/17/94 c ---- c c ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc logical function caseindp(a, b) c This function tests for character string equivalence c c compare what follows any leading blanks c check for equivalence of strings ignoring case. c only check through the last character of the shorter string c (where "shorter" is computed after removing leading blanks) c The only exception to this is if 1 string is null or completely c blank; then equivalence requires both strings to be null or c completely blank. character*(*) a, b character*1 c, d integer i, j, k, m, n, j1, k1, ii c look for leading blanks and neglect j = len(a) k = len(b) c check whether both strings are length zero if ((j .eq. 0) .or. (k .eq. 0)) then if (j .eq. k) then caseindp = .true. else caseindp = .false. endif return endif j1 = 0 do 10 i=1,j if (a(i:i) .ne. ' ') go to 20 j1 = j1 + 1 10 continue 20 continue k1 = 0 do 30 i=1,k if (b(i:i) .ne. ' ') go to 40 k1 = k1 + 1 30 continue 40 continue c check whether either string is completely blank if ((j1 .eq. j) .or. (k1 .eq. k)) then if ((j1 .eq. j) .and. (k1 .eq. k)) then caseindp = .true. else caseindp = .false. endif return endif c look for trailing blanks and neglect j2 = j do 50 i=1,j if (a(j-i+1:j-i+1) .ne. ' ') go to 60 j2 = j2 - 1 50 continue 60 continue k2 = k do 70 i=1,k if (b(k-i+1:k-i+1) .ne. ' ') go to 80 k2 = k2 - 1 70 continue 80 continue j = j2 - j1 k = k2 - k1 ii = min0(j, k) do 100 i=1,ii m = i + j1 n = i + k1 c = a(m:m) if (lge(c,'A') .and. lle(c,'Z')) c = char(ichar(c)+32) d = b(n:n) if (lge(d,'A') .and. lle(d,'Z')) d = char(ichar(d)+32) if (c .ne. d) then caseindp = .false. return endif 100 continue caseindp = .true. return end