PROGRAM FINALSD c@***h* OCEAN_HEAT_CONTENT/finalsd.f c c NAME c finalsd.f c c PURPOSE c Specifically for OCEAN_HEAT_CONTENT (OHC), the mean, standard c deviation, and number of observations are calculated on a c one-degree latitude/longitude grid for the variable temperature c anomaly (TA), which is the temperature value at a standard depth level c minus the monthly climatological mean at the same standard level. c c DESCRIPTION c finalsd.f is a multipurpose program within the World Ocean c Database (WOD) system. It reads standard level in situ profile c data and calculates mean, standard deviation, and number of c observations for grid boxes of chosen size for each chosen c geographic surface. The program applies quality flags as c designated by user to eliminate data points from the calculation. c The program can run a simple box mean or run iteratively, each time c running parallel five-degree latitude/longitude box mean and c standard deviations. If a data point is outside a chosen number of c standard deviations of the five-degree mean, that data point will c not be used for the next iteration. If a profile (set of c depth/variable measurements taken at the same date/time/position c from the same instrument) has more than a chosen number of c measurements outside the standard deviation criteria, the entire c profile is eliminated from use in the next iteration. c c INPUT c The set of monthly climatological means used is the World Ocean Atlas c 2009. The quality flags used are listed in file 'maskchoice_tanom.d', c which must be copied to 'maskchoice.d'. The variables set for TA c calculation are in file 'devfile_tanom.d', which must be copied to c 'devfile.d'. File subset.sheet_empty must be copied to 'subset.sheet'. c This file lists designated subsets of WOD on which to run c finalsd program, but for TA calculation, entire database is used. c File 'probefile.OHC' must be copied to probefile.d. This file lists c WOD dataset information for running OHC. Some of the data sets c within WOD must be pre-binned as their time frequency is much higher c than other data sets. Those data sets which are pre-binned are c listed in 'tobebinned.d' which is either in the directory in which c the program is run, or in [MAINDIR]/sys.inf. [MAINDIR] is the c location of the main directory of the WOD system. The c expendable bathythermograph (XBT) data set is bias corrected to the c Levitus et al. (2009) correction. The corrections are found in c file [MAINDIR]/sys.inf/antonov_xbtbias. The mechanical c bathythermograph (MBT) data are bias corrected to the c Levitus et al. (2009) correction. These corrections are found in c [MAINDIR]/sys.inf/antonov_mbtbias. Standard levels and corresponding c depths are found in file [MAINDIR]/sys.inf/standard_depths_orig.dat c The directory to which output TA mean, standard deviation, and c number of observations files is found in ./statpath.d or c [MAINDIR]/sys.inf/statpath.d. The output file for TA will be c in the given statpath/1.00/mean, statpath/1.00/sd, statpath/1.00/num c respectively. c cOUTPUT c output will be files of mean, standard deviation, number of c observations of TA for given time period. The files are in c World Ocean Atlas (WOA) gridded binary format. c cAUTHOR c Tim Boyer, National Centers for Environmental Information (NCEI) c cCOPYRIGHT cTHIS SOFTWARE AND ITS DOCUMENTATION ARE CONSIDERED TO BE IN THE cPUBLIC DOMAIN AND THUS ARE AVAILABLE FOR UNRESTRICTED PUBLIC USE. cTHEY ARE FURNISHED "AS IS." THE AUTHORS, THE UNITED STATES GOVERNMENT, cITS INSTRUMENTALITIES, OFFICERS, EMPLOYEES, AND AGENTS MAKE NO cWARRANTY, EXPRESS OR IMPLIED, AS TO THE USEFULNESS OF THE SOFTWARE cAND DOCUMENTATION FOR ANY PURPOSE. THEY ASSUME NO RESPONSIBILITY c(1) FOR THE USE OF THE SOFTWARE AND DOCUMENTATION; OR (2) TO PROVIDE cTECHNICAL SUPPORT TO USERS. c cHISTORY c August 8, 2017 Initial version prepared for Reference Environmental c Data Record (REDR) program c c***** parameter (minmeasures=5) parameter (maxparm=100, maxlevel=6000, maxcalc=200) parameter (maxprobe=100, maxmask=100, bmiss=-1.E10) parameter (jdim=180, idim=360, kdim=102, jyear=1800) parameter (ibinmax=500000,ncruisebinmax=100000) parameter (istanprofile=77,istantotfile=14) parameter (maxcc=200,maxfreq=10) parameter (idim5=72,jdim5=36) parameter (inum=500000000) parameter (kdim2=137, maxpas=5) parameter (ndays=366) parameter (maxsurf=1000, maxpot=100) parameter (xovertop=-1000.,xunderbot=-2000.) character citf(maxlevel)*15,citf2(maxlevel)*6 character mtype(maxlevel)*3 character*2 ncc(maxcc) character cc*2,cfile(maxmask)*80 character*80 fileout,fileoutd character ctype(maxprobe)*11,season(0:40)*9 character anlyfile*80,cyon*1 character csea*2,cyc*4 character*15 actionrec common /countryc/ ncc dimension ip2(0:maxlevel), numlevels(0:1) dimension jp(0:maxprobe),nopas(maxprobe),isur(maxcalc) dimension nprec(maxcalc,maxprobe), itoto(0:maxcalc) dimension irecord(0:2*maxparm), ioldrecord(0:2*maxparm) dimension ifn(maxmask), ipar2(maxmask), ipro2(maxmask) dimension nmask(5),isignif(maxlevel,maxcalc) dimension idt(3) dimension mask(idim5,jdim5,kdim) dimension sav(idim,jdim,kdim) dimension num5(idim5,jdim5,kdim),sntemp(maxlevel,maxcalc) dimension num(idim,jdim,kdim) dimension nprm(maxprobe,maxparm) dimension ndep(2,kdim,maxpas),inume(kdim) dimension nfreq(0:maxfreq+1,maxpas) dimension ittot(maxprobe,maxpas),itty(maxprobe,maxpas) dimension itot(maxprobe,maxpas),ndeptot(2,maxpas) dimension dz(kdim2),nbad(maxprobe,maxpas),nalt(maxprobe,maxpas) dimension npatot(maxparm),depth(maxlevel) dimension icounter(maxlevel,maxcalc) dimension istc(0:5,0:1) dimension istanoutfile(maxprobe,2) dimension npos(ibinmax),snbin(ibinmax,kdim) dimension nbdate(ibinmax) dimension jpotlev(maxpot),densesurf(maxsurf) dimension idensefile(0:maxparm,maxpot),iposdfile(maxpot) dimension idenserecord(0:maxparm),iolddenserecord(0:maxparm) c*********************************************************** c c These arrays are double precision to accomodate the c large numbers generated which need great accuracy c c*********************************************************** double precision tempsq5(idim5,jdim5,kdim) double precision tempsq(idim,jdim,kdim) double precision tempav(idim,jdim,kdim) double precision tempav5(idim5,jdim5,kdim) double precision temp(kdim2,maxcalc) common /thedata/ depth,sntemp common /levelsfull/ nfulldeps,dz c common /levels/ dz common /nopars/ nopas common /reqprobe/ jp common /extra/ ireqpars2,ireqsecond,ireqbio common /surrogate/ isur common /twoparm/ itwo(0:maxparm) common /parmord/ nprec common /eachparm/ numeach(6),numtot2,numtot3,ntax,numsecset common /reqps/ ip2 common /probes/ ctype common /thename/ anlyfile common /parnames/ citf common /parname1/ mtype common /parname4/ citf2 common /fyear/ nfin common /partype/ nprm common /parm2/ ipar2,ifn,ipro2 common /totalparms/ npatot common /maskfile/ cfile common /significant/ isignif common /fracta/ numsuc,nsub common /geop/ igeot,jlev,dsurf common/sdata/ temp common /multy/ mask common /squared/ tempsq common /mean/ tempav common /dist/ num common/anom/ sav common /stanmin/ minstan common/five/ tempav5,tempsq5,num5 common /iMlockACT/ ilockact common /iMlockNUM/ ilocknum common /cMlockTXT/ actionrec c**************************************************************** c c season contains the time period names: c 0 - Annual c 1-12 - the months c 13 -16 - the seasons c 17 -40 - half months c c**************************************************************** data season/'Annual','January','February',' March ', * ' April',' May ',' June ',' July ','August','September', * 'October','November','December', * 'Winter','Spring','Summer','Autumn', * 'Jan1-15','Jan15-31','Feb1-15','Feb15-29','Mar1-15', * 'Mar16-31','Apr1-15','Apr16-30','May1-15','May16-30', * 'Jun1-15','Jun16-30','Jul1-15','Jul16-31','Aug1-15', * 'Aug16-31','Sep1-15','Sep16-30','Oct1-15','Oct16-31', * 'Nov1-15','Nov16-30','Dec1-15','Dec16-31'/ ilockact = 0 ilocknum=0 actionrec = "finalsd" call input(npro,ireqpars,isoor) call setstandepsfromfile(nfulldeps,dz) call idate(idt) jyearmax=idt(3) c Get user input information (from file devfile.d) c write(6,*) 'mtype',mtype(1),mtype(2) call inputf(npass,isea,istc,msk,iprint,namis,itprint,npro, * ik,grids,numo,ndepcheck,nyc,nyc2,nomal,nomaly,nsuct,iylook, * cyon, kdimx, nloop, imn1, imn2, bmiss, ip2, ionegrid,kdim1, * idimgs,idimge,idim5gs,idim5ge,jdimgs,jdimge,jdim5gs,jdim5ge, * iannual,kdima,ibinsurf,ngrid,jquad,idimt,jdimt,idimt5,jdimt5, * i1,j1,i2,j2,i15,j15,i25,j25,mpass1,ireqpars,ik2,isoor) c write(6,*) 'depths',kdimx,kdim1,kdima,iannual c Get density surface information call getdensesurfs(igeot,nsurfs,densesurf,npots, * jpotlev,maxpot,maxsurf,ireflev,dsurf,kdim1,kdimx, * ndsurf,dz,jlev,kdim2,kdim) c Open accounting files (stantots[].d, stanpros[].d) call opensdcountfiles(ctype(jp(1)),npro,ik,nomal, * isea,csea,nyc,nyc2,cyc,istantotfile,istanprofile) c Initialize accounting arrays call initsdcounters(itoto,ndep,nfreq,ittot, * itty,itot,nbad,nalt,icounter,kdim,maxcalc,maxpas, * maxprobe,maxlevel,maxfreq,ndeptot) c c start 1/4 degree quadrant loop. For each time this loop is c run, a different part of the earth is run through the c standard deviation checks and calculations. This is c necessary due to space limitations. Arrays are set c to the one degree size of the earth. Smaller studies c need to run smaller parts of the earth successively. c For 1-degree calculations, jquad=1 c do 90 nquad=1,jquad write(6,*) 'Section ',nquad, ' of ', jquad c Initialize statistics arrays (tempav,tempsq, c num, tempav5, tempsq5, num5) call initialsd(kdim1,kdimx,i1,i2,j1,j2,i15,i25,j15,j25, * bmiss,ndsurf,igeot,msk) c set grid boundaries to run appropriate section of earth call setforquad(nquad,ngrid,lonquadi,latquadi, * idimgs,jdimgs,idimge,jdimge,idim5gs,jdim5gs,idim5ge, * jdim5ge,ionegrid) c write(6,*) 'from lon:',idimgs,idimge c write(6,*) 'from lat:',jdimgs,jdimge c start passes loop do 27 mpass=mpass1,npass c write(6,*) 'mpass',mpass,nloop,npro c nam - number of individual profiles written to file nam=0 c c start anomaly calculating loop. c c do 52 nlo=1,nloop do 52 nlo=1,1 c c start probe loop c do 21 nx=1,npro m=jp(nx) c Set up given probe for use: c --------------------------- c Check that given probe measures desired variable (iyes) c Open cruise files if data it to be binned c Adjust for surface-only data sets c Open desired masks (from maskchoice.d) c Open pentadal anomaly sdout masks if doing yearly anomalies c Reset variable association for each mask for calculated variables c Open data files c write(6,*) 'm',m,nopas(m),nprm(m,1) call setforspecprobe(m,ik,ik2,iyes,ibinyes,icruisefile, * icruisefilep,icruisefile2,ncruisebin,isurf,latindex, * lonindex,itimeindex,isoor,ireqtots,ireqpars,ioldrecord, * ienda,jj,jj1,ireqtotz,nmask,nomal,nyc,nyc2,iend,iuniqn, * cyon,msk) c write(6,*) 'here binyes',iyes,ibinyes if ( iyes .ne. 1 ) goto 21 c Set up given probe for density surfaces c Open density surface files c Reset position arrays call setforspecprobedense(cyon,igeot,npots,jpotlev, * idensefile,maxpot,iposdfile,ireqtotz,maxparm, * idenserecord,ioldenserecord,m,ireqpars,ip2,maxlevel, * densesurf,maxsurf,jjd) c Create automatic output file names, if necessary c Automatic output files are sdout[].s and sdout[].d files c They are created for msk=3 and msk=4 call createoutmasks(msk,mpass,npass,nquad,nomal,ngrid, * csea,minstan,istanoutfile,ik,m,jquad,isdsfile,isddfile, * nmask,cyc,fileout,fileoutd,ifn,maxmask) c Get geographic mask (if quarter-degree, or specific c geographic box requested). Integrate mask with c other requested masks call setsdgeography(nmaskg,ionegrid,jquad, * idimgs,idimge,jdimgs,jdimge,idimt,m,ik,nmask) c Set time loops (imn1,imn2,nyc,nyc2) call settimeloops(imn1,imn2,isea,ibinsurf,numsuc, * nsub,ihalf,isea2,iannual,nomal,nomaly,nlo,nyc,nyc2, * jyear,jyearmax,iylook,imn1x,imn2x) c write(6,*) 'ihalf',ihalf c c Start year loop (used for anomaly calculation, otherwise c nyc=nyc2) c do 48 iyx = nyc, nyc2 c c Start month loop c do 49 imx = imn1, imn2 c Get station sequential numbers for stations within desired c time period call setsddate(m,iyx,imx,nmask,nmaskg,nmasko, * ik,iend,isea2,iylook,nomal,jyear,jyearmax,isurf,nyc, * nprofcount,iyx0f,iyx1f,imx0,imx1) c write(6,*) 'nprofcount',nprofcount,imn1,imn2,nyc,nyc2 if ( nprofcount .eq. -2 ) goto 449 if ( nprofcount .eq. -1 ) goto 49 c read in appropriate analyzed field to calculate c anomalies (if necessary: nomal > 0, nprofcount > 0) c write(6,*) 'nameanomin' call nameanom(anlyfile,nomaly,kdim1,kdimx,imx, * iannual,nomal,nprofcount,kdima,idimgs,idimge, * jdimgs,jdimge,idimt,jdimt) write(6,*) 'nameanomout',sav(360,2,1) c Start the profile loop do 50 jj0=1,inum c Read in necessary data and prepare for calculations c read from database (readdata) c read density surface c check if desired variable is present c special inclusion for NO2+NO3 data c calculate geopotential thickness (if desired) c bin high density data (if desired) c bin surface-only data c mask out levels based on masks from maskchoice.d c set nlvs, kall c VARIABLE DATA WILL BE IN DOUBLE PRECISION TEMP ARRAY c AFTER THIS SUBROUTINE FOR ALL DATA EXCEPT BINNED C DATA AND SURFACE-ONLY DATA c write(6,*) 'jj0',jj0,iannual call preparesddata(jj0,nmask,jj,cc,numlevels, * iyear,month,iday,time,icruise,rlat,rlon,jj1, * ireqtotz,isoor,irecord,ioldrecord,m,ik,ik2,ienda, * ishallow,ideep,xovertop,xunderbot,nsurfs,idimt, * isea,nyc,nyc2,itoto,icounter,grids,latindex,lonindex, * itimeindex,kall,nlvs,iyear0,kdim1,kdimx,ibinyes,kdima, * idensefile,iposdfile,idenserecord,iolddenserecord,jjd, * icruisefile,icruisefile2,icruisefilep,ireqtots,isurf, * nbdate,npos,snbin,iyx0f,iyx1f,imx0,imx1,ibinsurf,iend, * ialready,idimgs,idimge,jdimgs,jdimge,iuniqn,ju,mpass, * nomal,nomaly,ireqpars,iannual,nquad,ipreparestat) c if ( ipreparestat .eq. 2 .or. ipreparestat .eq. 1) c if ( ju .eq. 11781361 ) c write(6,*) 'jjz',jj,iannual,ipreparestat,kall if ( ipreparestat .eq. 2 ) goto 50 if ( ipreparestat .eq. 1 ) goto 4 if ( mpass .eq. 1 ) write(98,*) jj do 500 nk=1,kall c prepare bin data for calculations c surface-only and binned data are placed c in double precision array temp c lat,lon,lat2,lon2 grid indices calculated c bin is checked too see if it is within c date/position criteria c if ( lat .eq. 164 .and. lon .eq. 288 ) then c write(6,*) 'tCASE:',lon,lat,temp(1,1), c * ju,jj,m,iyear,month,iday,tempav(lon,lat,1) c endif call preparesdbin(ialready,ibinyes,nk,rlat,rlon, * lvs1,kdim1,iyear0,iyear,month,iday,time,iannual, * nomal,iax,imn1x,imn2x,ihalf,nbdate,npos,grids,lat, * lon,lat2,lon2,idimt,idimgs,idimge,jdimgs,jdimge, * idim5gs,idim5ge,jdim5gs,jdim5ge,ity,itti,numeach(1), * isurf,lonindex,latindex,itimeindex,inume,snbin, * ik,isea,ngrid,idim,jdim,kdimx,ionegrid,jquad,ibinsurf, * ibinstat,ju) c write(6,*) 'preparesd',ibinstat,lvs1,nlvs,ibinyes c if ( jj .eq. 512008 ) c * write(6,*) 'ibinstat',ibinstat,lat,lon,numlevels(1), c * nlvs,iax,month if ( ibinstat .eq. 1 ) goto 500 c if ( lat .eq. 164 .and. lon .eq. 288 ) then c write(6,*) 'CASE:',lon,lat,temp(1,1), c * ju,jj,m,iyear,month,iday,tempav(lon,lat,1) c endif c for each level: c calculate mean, standard deviation c perform standard deviation check c remove data from statistics if standard deviation c check is failed if ( lvs1 .gt. nlvs ) write(6,*) jj do 42 ijz=lvs1,nlvs c write(6,*) 'ijz',ijz call performsdcheck(ijz,ndsurf,iax,kdima,ik, * ndep,mpass,ndeptot,itti,itot,npro,nx,nomal,msk, * lon,lat,lon2,lat2,ity,isurf,inume,nk,ittot,npass, * nmask,isddfile,numo,igeot,minstan,minmeasures,bmiss, * ju,iuniqn,jquad,jj,iperformstat) 42 continue c update accounting statistics call updatesdstats(iax,ity,ialready,nx,mpass, * npro,numo,isurf,ndepcheck,nalt,nbad,ndep,ndeptot, * itti,nfreq,maxfreq,namis,nam,numlevels,isoor, * igeot,kdim1,kdimx,isea,m,ik,rlon,rlat,iyear,month, * iday,cc,icruise,inume,lon,lat,lon2,lat2,bmiss,season, * jj,ndsurf,npass,nmask,isdsfile,maxpas,kdim,kdim2, * maxprobe,itty,iprint,msk,istanprofile,ju) 500 continue 50 continue 4 continue 449 continue c Reset masks to original maskchoice.d configuration c (remove present geographic/date masks) call resetmasks(isea2,iylook,nomal,nprofcount,nmask, * nmasko) 49 continue 48 continue c Close all files for present dataset call probesdclose(ifn,maxmask,nmask,ireqtotz,isoor, * icruisefile,icruisefile2,icruisefilep,m,iposdfile, * idensefile,jpotlev,npots,ireqpars,maxparm,maxpot, * ip2,maxlevel,densesurf,maxsurf,ibinyes,igeot) 21 continue 52 continue write(6,*) 'standev' c Calculate five-degree statistics if needed c Do this if five-degree statistics are requested,or c if this is not the last pass call calcfivedegstats(mpass,npass,kdim1,kdimx, * igeot,ndsurf,j15,j25,i15,i25,bmiss,istc) write(6,*) 'standev' c Calculate one-degree statistics if needed c Do this if this is the last pass call calconedegstats(mpass,npass,msk,kdim1,kdimx, * igeot,ndsurf,i1,i2,j1,j2) c write out accounting statistics if desired call printot(mpass,npro,itot,citf(ik),jp, * ctype,m,itty,ittot,nfreq,ndep,dz,ndeptot, * nbad,nalt,season(isea),numo,itprint,nquad, * jquad,maxfreq,istantotfile) c write out statistics, if necessary if ( mpass .eq. npass ) then write(6,*) 'final',tempav(360,73,1) call writeoutsdstats(igeot,kdim1,ndsurf,kdimx,istc, * idimgs,jdimgs,idimge,jdimge,idimt,jdimt,idim5gs,jdim5gs, * idim5ge,jdim5ge,idimt5,jdimt5,npass-mpass) endif 27 continue c Sort quarter-degree masks call sortqmasks(nquad,jquad,msk,npro, * istanoutfile,jp,maxprobe,fileout,fileoutd) 90 continue stop end SUBROUTINE ANOMALYMASKS(nyc1,nyc2,nmask,m,ip,nperiod,msk) C ANOMALYMASKS ADDS MASKS OF ANOMALY STANDARD DEVIATON C OUTLIERS FROM THE DESIGNATED PERIOD TO THE SET OF MASKS TO CHECK parameter (maxprobe=100,maxmask=100) character*12 filenames,filenamed,filename dimension nmask(5) common /parm2/ ipar2(maxmask),ifn(maxmask),ipro2(maxmask) filenames='sdout??^^.s'//CHAR(0) filenamed='sdout??^^.d'//CHAR(0) if ( ip .eq. 0 ) return ny1=nyc1-1900 ny2=nyc2-1900 numask=ny2-ny1+1 ncount=0 isame=0 if ( nyc2 - nyc1 +1 .eq. nperiod ) isame=1 if ( msk .eq. 0 ) isame=0 do 30 nn=ny1-nperiod+1,ny2 if ( isame .eq. 1 .and. nn .eq. ny1 ) goto 30 filename=filenames call insertyear(nn,'?',filename,11) call insertyear(nn+nperiod-1,'^',filename,11) call addoldmask(filename,nmask,2,m,ip,0) call uniqtoorder1mask(ifn(nmask(2)),m,1) filename=filenamed call insertyear(nn,'?',filename,11) call insertyear(nn+nperiod-1,'^',filename,11) call addoldmask(filename,nmask,3,m,ip,0) call uniqtoorder1mask(ifn(nmask(3)),m,2) 30 continue return end SUBROUTINE BINPROFS(jj,jj1,irecord,ioldrecord,m,bmiss, * maxlevel,maxcalc,numlevels,itwo,depth, * sntemp,isignif,nprec,nopas,maxparm, * icrcon,ibinmax,kdim,snbin,nbins,npos, * xgrid,idim,icruisefile,icruisefilep,ndate, * icruisefile2,intime,nyc,nyc2,maxprobe, * ibintype,ik,ireqtots,ip2,ireqpars2,iannual, * itoto,icounter,nmask,ifn,ipar2,ipro2, * maxmask) C BINPROFS BINS DATA FROM A CRUISE BY TIME PERIOD SPECIFIED C BY USER c*********************************************************** c c Passed Variables: c O jj - sequential order of station in dataset c O jj1 - present position in dataset header file c O ioldrecord - present position in measured variable file c I m - present dataset code c I bmiss - missing value marker c I maxlevel - maximum number of levels c I maxcalc - maximum number of measured and calculated variables c I numlevels(1) - number of standard levels (0=observed) c I itwo - 2=biological variable, 1=second header, 0=measured var c O depth, sntemp - depth, measured variable data c O isignif - number of significant figuers for corresponding c element in sntemp c I nprec - order for variable for dataset c I nopas - number of measured variables for dataset c I maxparm - maximum number of measured variables c O icrcon - full cruise number (cc+icruise) c I ibinmax - maximun number of data bins c I kdim - maximun number of examined levels c O snbin - output data bins c O nbins - number of output bins c O npos - geographic mean position for each bin c I xgrid - grid size (in degrees) c I idim - number of grids in longitudinal direction c I icruisefile,icruisefile2,icruisefilep - file ID numbers c for cruise index files c O ndate - mean date for bin c I intime - time period for binning c I nyc, nyc2 - years for binning c I maxprobe - maximum number of data sets c I inbintype - 0=whole cruise bin, 3=yearly bin, c 2=monthly bin, 1=daily bin c I ik - desired measured variable code c I ireqtots - number of measured variables to read in c I ip2 - variable codes for variables to read in c I ireqpars2 - number of calculated variables to read in c I iannual - set to one if annual statistics are to be used below c a certain level (kdima) c O itoto, icounter - counters set in depthmask, not used elsewhere c in this program c I nmask - number of each type of mask c I ifn - file ID number for each mask c I ipar2 - measured variable for each mask c I ipro2 - measured variable for each mask c I maxmask - maximum number of masks c c*************************************************************************** parameter (ibinmaxl=100000,kdiml=102,maxincruise=500000) parameter (maxmaskl=100,kdiml2=137,maxcalcl=200) character*2 cc dimension npos(ibinmax),xnumbin(ibinmaxl,kdiml) dimension snbin(ibinmax,kdim),ndate(ibinmax) dimension irecord(0:2*maxparm), ioldrecord(0:2*maxparm) dimension nprec(maxcalc,maxprobe),nopas(maxprobe) dimension numlevels(0:1),itwo(0:maxparm) dimension depth(maxlevel),sntemp(maxlevel,maxcalc) dimension isignif(maxlevel,maxcalc),jjin(maxincruise) dimension ip2(0:maxlevel),itoto(0:maxcalc) dimension ifn(maxmask), ipar2(maxmask), ipro2(maxmask) dimension icounter(maxlevel,maxcalc),nmask(5) dimension maskpos(maxmaskl),jpresent(maxmaskl) dimension igpresent(maxmaskl) common /geop/ igeot,jlev,dsurf common /levelsfull/ nfulldeps,dz(kdiml2) common /surrogate/ isur(maxcalcl) nbins=0 iend=0 ieof=0 ienda=0 c Mark positions of masks for reinitialization do 135 n=1,nmask(2) call filepos(ifn(n),0,1,maskpos(n)) call getjpresent(n,jpresent(n)) 135 continue do 136 n=nmask(2)+1,nmask(3) call filepos(ifn(n),0,2,maskpos(n)) call getjpresent(n,jpresent(n)) call getgpresent(n,igpresent(n)) 136 continue c get stations from cruise call findcruise(icrcon,ncruise,jjin,icruisefile, * icruisefilep,icruisefile2) c read in stations from cruise jjm1=0 iyear1=0 do 50 jj0=1,ncruise if ( icrcon .eq. 99009682 ) write(6,*) 'icr',jj0,jjin(jj0),jj if ( jjin(jj0) .lt. jj ) goto 50 iwrite=0 c if ( jjin(jj0) .ge. 30000 ) iwrite=1 if ( nmask(2) .gt. 0 .and. jjin(jj0) .gt. jj) then if ( jjm1 .lt. jjin(jj0) ) jjm1=jjin(jj0)-1 call seqread(jjm1,nmask(2),nmask(1),ieof,iend,ifn) if (jjm1 .ne. jjin(jj0) ) goto 50 endif call readdata(cc,sntemp,numlevels,jjin(jj0),iyear, * month,iday,time,icruise,rlat,rlon,jj1,depth, * nopas(m),ioldrecord,ireqtots,ip2,1, * maxlevel,irecord,maxcalc,m,nprec,itwo,bmiss, * 1,1,ieof,ireqpars2,isignif) if ( jjin(jj0) .eq. 1085883 ) * write(6,*) 'jjin0',sntemp(1,2),rlat,rlon,icruise,jj0 c modify station based on depth masks iwrite=0 if ( nmask(3) .gt. nmask(2) ) then call depthmask(maxlevel,jjin(jj0),ifn, * nmask(2), nmask(3), * sntemp,ipar2,icounter,itoto,ienda,ireqtots,ip2, * numlevels(1),bmiss,ireqpars2) endif if ( jjin(jj0) .eq. 1085883 ) * write(6,*) 'jjin1',sntemp(16,2) c set gradient if ( igeot .eq. 3 ) call setgradientx(numlevels, * sntemp,maxlevel,maxcalc,bmiss,ik,0) c get initial year as reference if ( iyear1 .eq. 0 ) iyear1=iyear if ( jj0 .eq. 1 ) then c Mark positions of masks for reinitialization do 35 n=1,nmask(2) call filepos(ifn(n),0,1,maskpos(n)) call getjpresent(n,jpresent(n)) 35 continue do 36 n=nmask(2)+1,nmask(3) call filepos(ifn(n),0,2,maskpos(n)) call getjpresent(n,jpresent(n)) call getgpresent(n,igpresent(n)) 36 continue endif c check if the necessary variable is present ipcheck=ip2(1) if ( ip2(1) .gt. ip2(maxlevel) ) ipcheck=isur(ip2(1)) if ( irecord(nopas(m)+nprec(ipcheck,m)) .le. -1 ) goto 50 c get position bin call grid(rlat,rlon,lat,lon,xgrid) ibinpos=idim*(lat-1)+lon if ( jjin(jj0) .eq. 1085883 ) * write(6,*) 'pos',rlat,rlon,lat,lon,idim,ibinpos c get date bin iyearadd=iyear-iyear1 ibindate=0 iseason=0 if ( ibintype .eq. 0 ) then c whole cruise ibindate=0 elseif ( ibintype .eq. 3 ) then c seasonal bin iseason=month/3 if ( mod(month,3) .ne. 0 ) iseason=iseason+1 ibindate=(4*iyearadd)+iseason elseif ( ibintype .eq. 2 ) then c monthly bin ibindate=(12*iyearadd)+month elseif ( ibintype .eq. 1 ) then c daily bin call juliant(iyear1,iyear,month,iday,time,2, * xjulian,jsig) ibindate=xjulian endif c if ( jjin(jj0) .eq. 1085883 ) if ( icrcon .eq. 99009682 ) * write(6,*) 'date',ibindate,iyear,nyc,nyc2, * intime,iseason,iannual,ibintype,iyearadd, * iseason c Make sure station is in correct time period c Year specific if ( nyc .gt. 1 ) then if ( iyear .lt. nyc .or. iyear .gt. nyc2 ) goto 50 endif c Season specific if ( iannual .eq. 0 ) then if ( intime .ge. 13 .and. intime .le. 16 ) then iseason=12+month/3 if ( mod(month,3) .ne. 0 ) iseason=iseason+1 if ( iseason .ne. intime ) goto 50 endif c Monthly specific if ( intime .ge. 1 .and. intime .le. 12 ) then if ( month .ne. intime ) goto 50 endif endif c See if this is a new bin or needs to be added to an old bin inbin=0 do 55 nx=1,nbins if ( npos(nx) .eq. ibinpos .and. * ndate(nx) .eq. ibindate ) then inbin=nx goto 56 endif 55 continue 56 continue if ( jjin(jj0) .eq. 1085883 ) * write(6,*) 'inbin',inbin if ( inbin .eq. 0 ) then nbins=nbins+1 if ( nbins .gt. ibinmax ) then write(6,*) 'NUMBER OF BINS HAS EXCEEDED MAXIMUM:',ibinmax write(6,*) 'PROGRAM WILL BE TERMINATED' stop endif inbin=nbins npos(nbins)=ibinpos ndate(nbins)=ibindate do 57 k=1,kdiml xnumbin(nbins,k)=0. snbin(nbins,k)=0.0 57 continue if ( jjin(jj0) .eq. 1085883 ) * write(6,*) 'nbin',jj,nbins,ibindate endif c Add data to bin do 60 nn=1,numlevels(1) if ( sntemp(nn,ik) .gt. bmiss ) then xnumbin(inbin,nn)=xnumbin(inbin,nn)+1. snbin(inbin,nn)=snbin(inbin,nn)+sntemp(nn,ik) if ( jjin(jj0) .eq. 1085883 ) write(6,*) 'nbin',nn,ik, * xnumbin(inbin,nn),snbin(inbin,nn) endif 60 continue 50 continue c Calculate means for bins do 80 n=1,nbins do 85 k=1,kdim if ( xnumbin(n,k) .gt. 0.0 ) then snbin(n,k)=snbin(n,k)/xnumbin(n,k) else snbin(n,k)=bmiss endif if ( jj .eq. 1085883 ) write(6,*) 'kdim',n,k, * xnumbin(n,k),snbin(n,k) 85 continue 80 continue c Reinitialize masks do 90 n=1,nmask(2) call posfile(maskpos(n)+1,ifn(n)) call putjpresent(n,jpresent(n)) 90 continue do 92 n=nmask(2)+1,nmask(3) call pos2file(maskpos(n)+1,ifn(n)) call putjpresent(n,jpresent(n)) call putgpresent(n,igpresent(n)) 92 continue return end SUBROUTINE BINSURF(depth,sntemp,maxcalc,maxlevel,nlevels, * latindex,lonindex,itimeindex,ip,xgrid, * nparams,idate) C BINSURF BINS SURFACE DATA AS PER USER INSTRUCTIONS FOR C USE IN CALCULATING MEANS c******************************************************* c c Passed Variables: c c O depth,sntemp - depth, measured variables at each level c I maxcalc - maximum number of calculated and measured c variables c I maxlevel - maximun number of depth levels c O nlevels - number of surface-only measurement sets c I latindex,lonindex,itimeindex - measured variable codes for c latitude, longitude, julian day c I ip - chosen measured variable code c I xgrid - grid size, in degrees c I nparams - total number of measured variables c I idate - -1 = no bins (date or geographic) c 0 = no date bins c 1 = bin by day c 2 = bin by month c 3 = bin by season c > 3 = bin only for designated year c c***************************************************** parameter (xdays=366.,bmiss=-1.E10) dimension depth(maxlevel),sntemp(maxlevel,maxcalc) dimension xmon(13) data xmon/0,31,59,90,120,151,181,212,243,273,304,334,365/ c************************************************** c c Bin choices: c c position - grid size, one-degree, quarter-degree, etc. c temporal - daily, monthly,yearly, cruise period c idate - -1 = no bins (date or geographic) c 0 = no date bins c 1 = bin by day c 2 = bin by month c 3 = bin by season c > 3 = bin only for designated year c c************************************************** c write(6,*) 'idate',idate,nl if ( idate .eq. -1 ) return nl=nlevels nlevels=0 iyearx=0 rjulmin=0.0 rjulmax=0.0 xsub=0.0 do 50 n=1,nl if ( n .le. maxlevel ) then d1=depth(n) ic=n np=0 else call overmax(maxlevel,nparams,n,0,ic,np) d1=sntemp(ic,np) endif if ( sntemp(ic,np+ip) .gt. bmiss ) then icheck=0 c Check for lat/lon if ( sntemp(ic,np+latindex) .gt. bmiss * .and. sntemp(ic,np+lonindex) .gt. bmiss ) then call grid(sntemp(ic,np+latindex),sntemp(ic,np+lonindex), * lat,lon,xgrid) c Check for date if ( idate .gt. 0 .and. * sntemp(ic,np+itimeindex) .gt. bmiss ) then icheck=1 rjul=sntemp(ic,np+itimeindex) iyearadd=0 if ( rjul .gt. xdays ) * call reducejulian(iyear,rjul,iyearadd) iyearx=iyear+iyearadd if ( (iyearx/4)*4 .eq. iyearx ) xsub=1. c set restricting julian days for this bin if ( idate .eq. 1 ) then rjulmin=int(rjul) rjulmax=int(rjul)+1. elseif ( idate .eq. 2 .or. idate .eq. 3) then xm=0. do 400 nn=2,13 if ( nn .eq. 3 ) xm=xsub if ( rjul .lt. xmon(nn)+xm ) then rjulmin=xmon(nn-1) if ( nn-1 .gt. 2 ) rjulmin=rjulmin+xsub rjulmax=xmon(nn) if ( nn .ge. 2 ) rjulmax=rjulmax+xsub if ( idate .eq. 3 ) then ijulmin=(nn-1)/3 if ( ijulmin*3 .eq. nn-1) ijulmin=ijulmin-1 ijulmin=(ijulmin*3)+1 rjulmin=xmon(ijulmin) if ( ijulmin .gt. 1 ) rjulmin=rjulmin+xsub rjulmax=xmon(ijulmin+3)+xsub endif goto 401 endif 400 continue 401 continue endif elseif ( idate .eq. 0 ) then icheck=1 endif endif c write(6,*) 'icheck',icheck,nlevels if ( icheck .eq. 1 ) then nlevels=nlevels+1 if ( nlevels .le. maxlevel ) then depth(nlevels)=d1 ic0=nlevels np0=0 else call overmax(maxlevel,nparams,nlevels,0,ic0,np0) sntemp(ic0,np0)=d1 endif sntemp(ic0,np0+ip)=sntemp(ic,np+ip) sntemp(ic0,np0+latindex)=sntemp(ic,np+latindex) sntemp(ic0,np0+lonindex)=sntemp(ic,np+lonindex) if ( sntemp(ic0,np0+lonindex) .gt. -180.001 .and. * sntemp(ic0,np0+lonindex) .lt. -179.999) then sntemp(ic0,np0+lonindex)=180.0 endif sntemp(ic0,np0+itimeindex)=sntemp(ic,np+itimeindex) if ( nlevels .lt. n ) sntemp(ic,np+ip)=bmiss nbin=1 c write(6,*) 'ic0',ic09,np0,ip,sntemp(ic0,np0+ip) do 100 n2=n+1,nl if ( n2 .le. maxlevel ) then d2=depth(n2) ic2=n2 np2=0 else call overmax(maxlevel,nparams,n2,0,ic2,np2) d2=sntemp(ic2,np2) endif if ( sntemp(ic2,np2+ip) .gt. bmiss .and. * d2 .ge. d1 - 5. .and. d2 .le. d1 + 5 ) then rlat2=sntemp(ic2,np2+latindex) rlon2=sntemp(ic2,np2+lonindex) if ( rlat2 .gt. bmiss .and. rlon2 .gt. bmiss ) then call grid(rlat2,rlon2,lat2,lon2,xgrid) if ( lat2 .eq. lat .and. lon2 .eq. lon ) then ibin=0 rjul2=sntemp(ic2,np2+itimeindex) if ( idate .gt. 0 .and. rjul2 .gt. bmiss ) then if ( idate .gt. 3 .and. iyearx .eq. idate ) then ibin=1 elseif ( rjul2 .ge. rjulmin .and. * rjul2 .lt. rjulmax ) then ibin=1 endif elseif ( idate .eq. 0 ) then ibin=1 endif if ( ibin .eq. 1 ) then sntemp(ic0,np0+ip)=sntemp(ic0,np0+ip)+sntemp(ic2,np2+ip) sntemp(ic0,np0+latindex)=sntemp(ic0,np0+latindex)+ * sntemp(ic2,np2+latindex) sntemplon=sntemp(ic2,np2+lonindex) if ( sntemplon .gt. -180.001 .and. * sntemplon .lt. -179.999) then sntemplon=180.0 endif sntemp(ic0,np0+lonindex)=sntemp(ic0,np0+lonindex)+ * sntemplon sntemp(ic0,np0+itimeindex)=sntemp(ic0,np0+itimeindex)+ * sntemp(ic2,np2+itimeindex) sntemp(ic2,np2+ip)=bmiss nbin=nbin+1 endif endif endif endif 100 continue c Calculate mean value for present bin sntemp(ic0,np0+ip)=sntemp(ic0,np0+ip)/nbin c write(6,*) 'sntemp',ic0,np0,ip,sntemp(ic2,np2+ip),nbin sntemp(ic0,np0+latindex)=sntemp(ic0,np0+latindex)/nbin sntemp(ic0,np0+lonindex)=sntemp(ic0,np0+lonindex)/nbin sntemp(ic0,np0+itimeindex)=sntemp(ic0,np0+itimeindex)/nbin endif endif 50 continue return end SUBROUTINE CALCFIVEDEGSTATS(mpass,npass,kdim1,kdimx, * igeot,ndsurf,j15,j25,i15,i25,bmiss,istc) C CALCFIVEDEGSTATS CALCULATES FIVE DEGREE C STATISTICS FROM ONE-DEGREE STATISTICS c**************************************************** c c Passed Variables: c I mpass - present pass through standard deviation check c I npass - total number of passes through standard deviation check c I kdim1 - first level to perform calculations c I kdimx - deepest level to perform calculations c I igeot - 1=calculate geopotential thickness c -1=perform calculations on density surfaces c 0=normal c I ndsurf - index for desired density surface c I i15,j15,i25,j25 - unchanging 5x5 grid boundaries c I bmiss - missing value marker c I istc - information on which statistics to save: c 0. one grid box means c 1. one grid box data distributions c 2. one grid box standard deviations c 3. five grid box means c 4. five grid box data distributions c 5. five grid box standard deviations c each array member is set to file ID number if c statistic are to be saved, zero otherwise. c c*************************************************** parameter (jdim=180, idim=360, kdim=102) parameter (jdim5=36,idim5=72) dimension num5(idim5,jdim5,kdim) dimension num(idim,jdim,kdim) dimension istc(0:5,0:1) double precision tempsq5(idim5,jdim5,kdim) double precision tempsq(idim,jdim,kdim) double precision tempav(idim,jdim,kdim) double precision tempav5(idim5,jdim5,kdim) common /squared/ tempsq common /mean/ tempav common /dist/ num common /five/ tempav5,tempsq5,num5 c******************************************************** c c If this is not the last pass, or five degree statistics c are requested, calculate five degree statistics c c******************************************************** if ( mpass.lt.npass.or. * istc(3,0).gt.0.or.istc(4,0).gt.0.or.istc(5,0).gt.0) then do 222 kz=kdim1,kdimx if ( igeot .ne. -1 ) then k=kz else k=kz-ndsurf+1 endif do 223 j=j15,j25 js=(j-1)*5+1 do 224 i=i15,i25 is=(i-1)*5+1 c************************************************** c c Reset five grid box statistic c c************************************************** tempav5(i,j,k)=bmiss tempsq5(i,j,k)=bmiss num5(i,j,k)=0 c************************************************** c c Add together statistics from each one grid box c in the five grid box area c c************************************************** do 1223 j9=js,js+4 do 1224 i9=is,is+4 if (num(i9,j9,k).gt.0) then num5(i,j,k)=num5(i,j,k)+num(i9,j9,k) if (tempav5(i,j,k).gt.bmiss) then tempav5(i,j,k)=tempav5(i,j,k)+tempav(i9,j9,k) tempsq5(i,j,k)=tempsq5(i,j,k)+tempsq(i9,j9,k) else tempav5(i,j,k)=tempav(i9,j9,k) tempsq5(i,j,k)=tempsq(i9,j9,k) endif c******************************************************* c c Reset one grid box statistics c c******************************************************* if ( mpass .ne. npass ) then tempav(i9,j9,k)=bmiss tempsq(i9,j9,k)=bmiss num(i9,j9,k)=0 endif endif 1224 continue 1223 continue c****************************************************** c c Calculate standard deviation and mean c c****************************************************** if (num5(i,j,k).gt.0) then tempav5(i,j,k)=tempav5(i,j,k)/num5(i,j,k) tempsq5(i,j,k)=tempsq5(i,j,k)/num5(i,j,k) if (tempsq5(i,j,k)-tempav5(i,j,k)**2.le.0.D+00) then if ( mpass.eq.npass ) then tempsq5(i,j,k)=0.D+00 else tempsq5(i,j,k)=50.D+00 endif else tempsq5(i,j,k)=dsqrt(tempsq5(i,j,k)-tempav5(i,j,k)**2) endif endif 224 continue 223 continue 222 continue endif return end SUBROUTINE CALCONEDEGSTATS(mpass,npass,msk,kdim1,kdimx, * igeot,ndsurf,i1,i2,j1,j2) C CALCONEDEGSTATS CALCULATES ONE-DEGREE STATISTICS c**************************************************** c c Passed Variables: c I mpass - present pass through standard deviation check c I npass - total number of passes through standard deviation check c I msk = 0 - simply calculate means c = 1 - calculate means after throwing out c standard deviation outliers c = 2 - calculate means after checking for c standard deviation outliers against c current five-degree statistics c = 3 - same as (2), but standard deviation c outliers masks are automatically created c = 4 - same as (1), but standard deviation c outliers masks are automatically created c I kdim1 - first level to perform calculations c I kdimx - deepest level to perform calculations c I igeot - 1=calculate geopotential thickness c -1=perform calculations on density surfaces c 0=normal c I ndsurf - index for desired density surface c I i1,j1,i2,j2 - unchanging 1x1 grid boundaries c c*************************************************** parameter (jdim=180, idim=360, kdim=102) dimension num(idim,jdim,kdim) double precision tempsq(idim,jdim,kdim) double precision tempav(idim,jdim,kdim) common /squared/ tempsq common /mean/ tempav common /dist/ num c************************************************************* c c If this is the last pass through the standard deviation c check, compute mean and standard deviation for each c one grid box square c c************************************************************* if (mpass.eq.npass.and. (msk.lt.2 .or. msk .eq. 4)) then do 22 kz=kdim1,kdimx if ( igeot .ne. -1 ) then k=kz else k=kz-ndsurf+1 endif do 23 j=j1,j2 do 24 i=i1,i2 c if ( j .eq. 75 .and. i .eq. 357 ) c * write(6,*) 'n',k,kz,tempav(i,j,k),num(i,j,k), c * tempav(i,j,k)/num(i,j,k) if (num(i,j,k).gt.0) then tempav(i,j,k)=tempav(i,j,k)/num(i,j,k) c if ( tempav(i,j,k) .ge. 24.135 .and. c * tempav(i,j,k) .lt. 24.145 ) c * write(6,*)'thisone?',i,j,k,kz,tempav(i,j,k), c * num(i,j,k) tempsq(i,j,k)=tempsq(i,j,k)/num(i,j,k) if (tempsq(i,j,k)-tempav(i,j,k)**2.le.0.D+00) then tempsq(i,j,k)=0.D+00 else tempsq(i,j,k)=dsqrt(tempsq(i,j,k)-tempav(i,j,k)**2) endif endif 24 continue 23 continue 22 continue endif return end SUBROUTINE CREATEOUTMASKS(msk,mpass,npass,nquad,nomal,ngrid, * csea,minstan,istanoutfile,ik,m,jquad,isdsfile,isddfile,nmask, * cyc,fileout,fileoutd,ifn,maxmask) C CREATEOUTMASKS CREATES SDOUT MASKS AUTOMATICALLY c I msk = 0 - simply calculate means c = 1 - calculate means after throwing out c standard deviation outliers c = 2 - calculate means after checking for c standard deviation outliers against c current five-degree statistics c = 3 - same as (2), but standard deviation c outliers masks are automatically created c = 4 - same as (1), but standard deviation c outliers masks are automatically created c I mpass - present pass through data c I npass - number of passes to run through standard deviation c check c I nquad - present section for quarter-degree (1-16) c always 1 for one-degree c I nomal - set to one if anomalies are to be calculated, zero c for full parameter values c I ngrid - the number of grids at the requested grid size c needed to equal one one-degree grid in either c latitude or longitude c I csea - character representation of time period ('00'=annual) c I minstan - minimum number of standard deviations from mean for c outlier check (3 is default) c O istanoutfile - C file ID number for standard deviation outlier c masks for quarter-degree. This is necessary c because these masks must be kept open through c 16 iterations through database c I ik - variable code (1=temperature, 2=salinity, etc.) c I m - probe code c I jquad - number of iterations (1 for 1-degree, 16 for 1/4-degree) c O isdsfile - present sdoutXX.s mask file ID. If user has setup c this output mask, number=ifn(nmask(4)) c O isddfile - present sdoutXX.d mask file ID. If user has set up c this output mask, number=ifn(nmask(4)+1) c I nmask - number of each type of mask c I cyc - character representation of beginning-end year period c 4 characters, first two are last two digits of first c year periods, second two are last two digits of last c year period (for > 1999, use A in place of first zero, c B in place of first one, etc. So 2001=A1, 2011=B1. c O fileout - sdout[].s file name c O fileoutd - sdout[].d file name parameter (maxprobe=100) character*2 csea character*4 cyc character*11 ctype(maxprobe) character*80 title,title2,fileout,fileoutd character*80 filename common /probes/ ctype dimension istanoutfile(maxprobe,2) dimension nmask(5),ifn(maxmask) c initialization nfadd=0 ny1r=0 ny2r=0 call clearstring(fileout,80) call clearstring(fileoutd,80) call clearstring(title,80) call clearstring(title2,80) c msk=3, msk=4 require sdout masks to be opened c automatically if (msk.ge.3) then c create automatic masks on last pass through if ( mpass .eq. npass .and. npass .gt. 1 .and. * nquad .eq. 1) then call clearstring(fileout,80) c start with generic names c sdout00.s, sdout00qt.s for regular run c sdout0000.s for anomaly run if ( nomal .eq. 0 ) then nfadd=0 if ( ngrid .eq. 1 ) then fileout(1:9)='sdout00.s' elseif ( ngrid .eq. 4 ) then fileout(1:12)='sdout00.qt.s' nfadd=3 endif c Insert season to generic name for regular run fileout(6:7)=csea else c Insert years to generic name for anomaly run fileout(1:11)='sdout0000.s' fileout(6:9)=cyc nfadd=2 endif c Add minimum number of standard deviations c if not default value (3) if ( minstan .gt. 9 ) then write(fileout(9+nfadd:10+nfadd),'(i2)') minstan fileout(11+nfadd:14+nfadd)='.s'//CHAR(0) elseif ( minstan .gt. 3 ) then fileout(9+nfadd:9+nfadd)='0' write(fileout(10+nfadd:10+nfadd),'(i1)') minstan fileout(11+nfadd:14+nfadd)='.s'//CHAR(0) else fileout(8+nfadd:11+nfadd)='.s'//CHAR(0) endif c Open actual mask for one-degree c Open temporary mask for q-degree for later sorting istanoutfile(m,1)=-1 if ( jquad .le. 1 ) then call maskname(fileout,ik,filename,m) call fileopen(filename,'w+'//CHAR(0),istanoutfile(m,1)) else call tempfile(istanoutfile(m,1)) endif c create name for .d mask fileoutd=fileout if ( minstan .le. 3 ) then fileoutd(9+nfadd:9+nfadd) = 'd' else fileoutd(12+nfadd:12+nfadd)='d' endif c open .d file istanoutfile(m,2)=-1 if ( jquad .le. 1 ) then call maskname(fileoutd,ik,filename,m) call fileopen(filename,'w+'//CHAR(0),istanoutfile(m,2)) else call tempfile(istanoutfile(m,2)) endif endif endif c Set presently used sd file IDs if ( nmask(4) .gt. nmask(3) ) then isdsfile=ifn(nmask(4)) else isdsfile=istanoutfile(m,1) endif if ( nmask(5) .gt. nmask(4) ) then isddfile=ifn(nmask(5)) else isddfile=istanoutfile(m,2) endif return end SUBROUTINE GETDENSESURFS(igeot,nsurfs,densesurf,npots, * jpotlev,maxpot,maxsurf,ireflev,dsurf,kdim1,kdimx, * ndsurf,dz,jlev,kdim2,kdim) c Passed Variables c c I igeot - calculate geopotential thickness (1) or not (0) c -1 = calculate statistics on given density surface c O nsurfs - number of density surfaces c O densesurf - individual density surfaces c O npots - number of potential surfaces c O jpotlev - level of reference for each potential surface c I maxpot - maximum number of potential surfaces c I maxsurf - maximum number of density surfaces c O ireflev - index of desired reference level c I jlev - desired reference level c I dsurf - desired density surface c O kdim1 - index for desired density surface c O kdimx - also index for desired density surface c kdim1 and kdimx are used instead of normal depth indices c for density surfaces. They may represent a range of c density surfaces c O dz - value(s) of density surface(s) c I kdim2 - number of elements in dz c I kdim - number of levels on which check can be performed dimension dz(kdim2) dimension jpotlev(maxpot),densesurf(maxsurf) if ( igeot .eq. -1 ) then call setdensesurfs(nsurfs,densesurf,npots,jpotlev, * maxpot,maxsurf) write(6,*) 'densesurf',nsurfs,densesurf(1) ireflev=0 do 719 nz=1,npots 719 if ( jpotlev(nz) .eq. jlev ) ireflev=nz if ( ireflev .eq. 0 ) then write(6,*) 'Reference level not found:',jlev stop endif npots=1 jpotlev(1)=jlev kdim1=0 kdimx=0 ndsurf=0 if ( dsurf .gt. -1. ) then do 720 nz=1,nsurfs write(6,*) dsurf,densesurf(nz),nz if ( dsurf .ge. densesurf(nz)-1.E-5 .and. * dsurf .le. densesurf(nz)+1.E-5 ) ndsurf=nz 720 continue kdim1=ndsurf kdimx=kdim1 if ( kdim1 .eq. 0 ) then write(6,*) 'Density surface not found:',dsurf stop endif else open(9,file='./dum.densesurf',status='old') read(9,*) ndsurf close(9) kdim1=ndsurf kdimx=ndsurf+kdim-1 if ( kdimx .gt. nsurfs ) kdimx=nsurfs write(6,*) 'density surfaces:', * densesurf(kdim1),densesurf(kdimx) endif do 744 nz=kdim1,kdimx dz(nz-ndsurf+1)=densesurf(nz) 744 continue endif return end SUBROUTINE INITIALSD(kdim1z,kdimxz,idimgs,idimge,jdimgs,jdimge, * idim5gs,idim5ge,jdim5gs,jdim5ge,bmiss,ndsurf,igeot,msk) C INITIALSD INITILIAZES ARRAYS HOLDING STATISTICS c********************************************************* c c Passed Variables c I kdim1z - first level to perform calculations c I kdimxz - deepest level to perform calculations c I igeot - 1=calculate geopotential thickness c -1=perform calculations on density surfaces c 0=normal c I ndsurf - index for desired density surface c I idimgs - starting longitude to perform checks on c I idimge - ending longitude to perform checks on c I jdimgs - starting latitude to perform checks on c I jdimge - ending latitude to perform checks on c I idim5gs - starting five grid longitude to perform checks on c I idim5gs - ending five grid longitude to perform checks on c I jdim5gs - starting five grid latitude to perform checks on c I jdim5ge - ending five grid latitude to perform checks on c I bmiss - missing value marker c I msk = 0 - simply calculate means c = 1 - calculate means after throwing out c standard deviation outliers c = 2 - calculate means after checking for c standard deviation outliers against c current five-degree statistics c = 3 - same as (2), but standard deviation c outliers masks are automatically created c = 4 - same as (1), but standard deviation c outliers masks are automatically created c c********************************************************** parameter (idim=360,jdim=180,kdim=102,idim5=72,jdim5=36) c double precision tempsq(idim,jdim,kdim) double precision tempav(idim,jdim,kdim) double precision tempav5(idim5,jdim5,kdim) double precision tempsq5(idim5,jdim5,kdim) double precision dmiss c common /squared/ tempsq common /mean/ tempav common /dist/ num(idim,jdim,kdim) common /five/ tempav5,tempsq5, * num5(idim5,jdim5,kdim) dmiss=bmiss*1.D00 if ( igeot .ne. -1 ) then kdim1=kdim1z kdimx=kdimxz else kdim1=kdim1z-ndsurf+1 kdimx=kdimxz-ndsurf+1 endif c c Initialize Counters c do 700 k=kdim1,kdimx c c do 701 j=jdimgs,jdimge do 702 i=idimgs,idimge tempav(i,j,k)=dmiss tempsq(i,j,k)=dmiss num(i,j,k)=0 702 continue 701 continue if ( msk .ne. 2 .and. msk .ne. 3 ) then do 703 j=jdim5gs,jdim5ge do 704 i=idim5gs,idim5ge tempav5(i,j,k)=dmiss tempsq5(i,j,k)=dmiss num5(i,j,k)=0 704 continue 703 continue endif 700 continue c return end SUBROUTINE INITSDCOUNTERS(itoto,ndep,nfreq,ittot, * itty,itot,nbad,nalt,icounter,kdim,maxcalc,maxpas, * maxprobe,maxlevel,maxfreq,ndeptot) C INITSDCOUNTERS INITIALIZES ACCOUNTING C ARRAYS AND VARIABLES FOR FINALSD c**************************************************** c c Passed Variables: c O itty - number of profiles with outliers (per probe/pass) c O nalt - number of altered profiles (per probe/pass) c altered profile is profile with standard deviation c outliers masked out, but still used c O nbad - number of profiles not used due to outliers (per probe/pass) c O ndep - 1=frequency of error (per depth/pass) c 2=total number of measurements (per depth/pass) c O ndeptot - 1=total number of errors (per pass) c 2=total number of measurements (per pass) c O itot - number of profiles with desired variable (per probei/per pass) c O ittot - number of profiles with outliers (per probe/per pass) c O nfreq - number of profiles by the number of outliers c contained (per pass) c O itoto,icounter - counters associated with depthmask. c not used here c I maxpas - maximum number of passes through standard deviation check c I maxfreq - maximum number of outliers per profile listed in statistics c I maxprobe - maximum number of data sets c I maxlevel - maximun number of depth levels for a profile c I kdim - number of depth levels examined in standard deviation check c c*********************************************************************** dimension itoto(0:maxcalc) dimension ndep(2,kdim,maxpas) dimension nfreq(0:maxfreq+1,maxpas) dimension ittot(maxprobe,maxpas),itty(maxprobe,maxpas) dimension itot(maxprobe,maxpas),ndeptot(2,maxpas) dimension nbad(maxprobe,maxpas),nalt(maxprobe,maxpas) dimension icounter(maxlevel,maxcalc) call arrayinit(itoto,0,maxcalc) call arrayinit3(ndep,1,2,1,kdim,1,maxpas) call arrayinit2(nfreq,0,maxfreq+1,1,maxpas) call arrayinit2(ittot,1,maxprobe,1,maxpas) call arrayinit2(itty,1,maxprobe,1,maxpas) call arrayinit2(itot,1,maxprobe,1,maxpas) call arrayinit2(ndeptot,1,2,1,maxpas) call arrayinit2(nbad,1,maxprobe,1,maxpas) call arrayinit2(nalt,1,maxprobe,1,maxpas) call arrayinit2(icounter,1,maxlevel,1,maxcalc) return end SUBROUTINE INITSDVARS( nsuct, isea, nomal, nomaly, * iylook, nsucy, numsuc, nyc, nyc2, istc, grids, msk, numo, * ndepcheck, cyon, cyon2, uninfile, unoutfile, statpath, * anlypath, iprint, namis, itprint, uninparm, kdimx, npass, * ionegrid, latgrid, longrid, kdim1, iannual, kdima, nsub, * igeot,ibinsurf,jlev,dsurf ) C INITSDVARS INITIALIZES USER SET VARIABLES C FOR FINALSD character*1 cyon,cyon2 character*80 unoutfile,uninfile,statpath,anlypath character*80 uninparm dimension latgrid(2),longrid(2) nsuct=0 isea=0 nomal=0 nomaly=0 iylook=0 nsucy=0 istc=0 grids=0.0 msk=0 numo=0 ndepcheck=0 call clearstring(cyon,1) call clearstring(cyon2,1) call clearstring(unoutfile,80) call clearstring(uninfile,80) call clearstring(statpath,80) call clearstring(anlypath,80) iprint=0 namis=0 itprint=0 call clearstring(uninparm,80) kdimx=0 npass=0 ionegrid=0 latgrid(1)=0 longrid(1)=0 latgrid(2)=0 longrid(2)=0 kdim1=0 iannual=0 kdima=0 minstan=0 nsub=0 igeot=0 ibinsurf=0 jlev=0 dsurf=0 return end SUBROUTINE INPUTF(npass,isea,istc,msk,iprint,namis,itprint,npro, * ik,grids,numo,ndepcheck,nyc,nyc2,nomal,nomaly,nsuct,iylook, * cyon, kdimx, nloop, imn1, imn2, bmiss, ip2, ionegrid,kdim1, * idimgs,idimge,idim5gs,idim5ge,jdimgs,jdimge,jdim5gs,jdim5ge, * iannual,kdima,ibinsurf,ngrid,jquad,idimt,jdimt,idimt5,jdimt5, * i1,j1,i2,j2,i15,j15,i25,j25,mpass1,ireqpars,ik2,isoor) C INPUTF READS FILE DEVFILE.D AND ARRANGES USER INPUT INFORMATION C TO BE USED TO RUN STATISTICAL CHECKS AND COMPILATIONS c***************************************************************** c c Passed Variables: c O npass - number of passes through standard deviation check. c One pass is all that is completed for simple c calculations, npass + 1 passes for deviation checks. c O isea - time period of statistical checks: c 00 - annual c 01-12 - months c 13-16 - seasons c 17-40 - half months c O istc - information on which statistics to save: c 0. one grid box means c 1. one grid box data distributions c 2. one grid box standard deviations c 3. five grid box means c 4. five grid box data distributions c 5. five grid box standard deviations c each array member is set to file ID number if c statistic are to be saved, zero otherwise. c O msk = 0 - simply calculate means c = 1 - calculate means after throwing out c standard deviation outliers c = 2 - calculate means after checking for c standard deviation outliers against c current five-degree statistics c = 3 - same as (2), but standard deviation c outliers masks are automatically created c = 4 - same as (1), but standard deviation c outliers masks are automatically created c O iprint - set to one if the user requests individual profile c printouts, zero otherwise c O namis - the number of individual profiles requested c O itprint - set to one if statistics totals are to be printed out c O npro - number of requested probe types c O ik - code for requested measured variable c O grids - grid size of statistical check c O numo - number of standard deviation outliers which invalidate a c profile c O ndepcheck - set to one to discard standard deviation outliers c in useable profiles c O nyc - starting year for statistics checks c O nyc2 - ending year of statistics checks c O nomal - set to one if anomalies are to be calculated, zero c for full variable values c O nomaly - read which type of anomaly: annual - 0, seasonal - 1, c monthly - 2 c O nsuct - set to zero for one time period, one for all time periods c requested in devshell c O iylook -look at only certain years (2), months (1), or all values c O cyon - set to 'n' if a non-standard input file name is used c set to 'y' otherwise c O cyon2 - set to 'n' if a non-standard output file name is used c set to 'y' otherwise c O kdimx - deepest level to perform checks c O nloop - number of times through checks to perform checks on c seasonal anomalies (4 times), monthly anomalies (12 times) c or regular checks (1 time). c O imn1 - first month(season) to check c O imn2 - last month(season) to check c I bmiss - missing value marker c O ip2 - parameter codes c O ionegrid - set to one if only one five grid box is to be checked, c zero otherwise c O kdim1 - first level to perform checks on c O idimgs - starting longitude to perform checks on c O idimge - ending longitude to perform checks on c O jdimgs - starting latitude to perform checks on c O jdimge - ending latitude to perform checks on c O idim5gs - starting five grid longitude to perform checks on c O idim5gs - ending five grid longitude to perform checks on c O jdim5gs - starting five grid latitude to perform checks on c O jdim5ge - ending five grid latitude to perform checks on c O iannual - set to one if annual statistics are to be used below c a certain depth level c O kdima - depth level at which annual statistics are to be used c O ibinsurf - surface data binning:-1=no date or geographic binning c 0=no date binning,gridbox c 1=bin by day,gridbox c 2=bin by month,gridbox c 3=bin by season,gridbox c O ngrid - the number of grids at the requested grid size c needed to equal one one-degree grid in either latitude or longitude c direction c O jquad - number of iterations (1 for 1-degree, 16 for 1/4-degree) c O idimt,jdimt - ngrid*idim, ngrid*jdim c O idimt5,jdimt5 - ngrid*idim5, ngrid*jdimt5 c O i1,j1,i2,j2,i15,j15,i25,j25 - unchanging grid boundaries. c (initially set to idimge,jdimgs, etc. It is necessary to store c these initial values, since idimge, jdimgs, etc. may change c during program run) c O mpass1 - pass start number (if this is a check against five-degree c statistics (msk=2 or 3), mpass1=npass) c O ireqpars - number of variables to read from file. Always one, c except for NO3, in which NO2+NO3 is also read in, c and surface-only data, in which lat,lon,time are also c read in c O ik2 - second measured variable. zero, except when NO3 is first c variable, in which case ik2=23 for NO2+NO3 c O isoor - always 1 for standard levels c c******************************************************************* c*********************************************************** c c PARAMETERS: c c maxlevel - maximum number of depth levels c maxparm - maximum number of measured parameters c maxcalc - maximum number of measured and calculated parameters c maxprobe - maximum number of probes c maxmask - maximum number of masks c jdim - number of latitudinal grid spacings c idim - number of longitudinal grid spacings c kdim - number of standard levels for which calculations are made c idim5 - number of five degree longitudinal grid spacings c jdim5 - number of five degree latitudinal grid spacings c c************************************************************** parameter (maxparm=100, maxlevel=6000, maxcalc=200) parameter (maxprobe=100, maxmask=100) parameter (jdim=180, idim=360, kdim=102) parameter (idim5=72,jdim5=36,kdim2=137) c************************************************************** c c Character Arrays: c c statpath - standard statistics output path c uninparm - name of non-standard input file c meanfile - name of mean statistics file c numfile - name of data distribution file c sgrid - character representation of grid size c csea - character representation of time period c mtype - 1-3 letter parameter character code c anlypath - path for analyzed field c anlyfile - name of analyzed field file c unoutfile - non-standard output file name c uninfile - non-standard input file name c cyon - set to 'n' if non-standard input file name c cyon2 - set to 'n' if non-standard output file name c standfile - standard deviation file name c stand5file - five grid box standard deviation file name c mean5file - five grid box mean file name c num5file - five grid box data distribution file name c filename - constructed extra file name c c*************************************************************** character*4 cdepthset character statpath*80, uninparm*80 character meanfile*80,numfile*80 character sgrid*4,csea*2,mtype*3 character anlypath*80,anlyfile*80,filename*80 character anlyfile2*80 character unoutfile*80,uninfile*80,cyon*1,cyon2*1 character*80 standfile,stand5file,mean5file,num5file character*80 outfile,infile c**************************************************************** c c Arrays: c tempav - mean of data c tempsq - standard deviation of data c tempav5 - mean of five grid point squares c tempsq5 - standard deviation of data in five grid point squares c ip2 - parameter code numbers c istc decides which statistics to save c sntempav - single precision parameter means c c**************************************************************** dimension ip2(0:maxlevel), sntempav(idim,jdim,kdim) dimension istc(0:5,0:1),latgrid(2),longrid(2) double precision tempav(idim,jdim,kdim) double precision tempsq(idim,jdim,kdim) double precision tempav5(idim5,jdim5,kdim) double precision tempsq5(idim5,jdim5,kdim) c*************************************************************** c c Common Blocks: c c extra contains: c ireqpars2 - number of requested calculated parameters c ireqsecond - number of requested second header parameters c ireqbio - number of requested biological parameters c parnames - contains citf, 1-15 letter parameter names c parname4 - contains citf2, 1-6 letter parameter names c used in writing out individual profiles. c parname1 - contains mtype, 1-3 letter parameter code used c in making file names for analyses and other applications. c probes - contains ctype, 1-11 letter probe names c thename - contains anlyfile, the pathname to the analyzed data c c data - contains the double precision parameter data c multy - contains mask, the standard deviation multiplier c squared - contains tempsq, the double precision parameter c values standard deviation c mean - contains tempav, the double precision mean parameter c values at each grid point c dist - contains the data distributios at each grid point c anom - contains sav, the analyzed parameter data c five - contains the five grid point average, standard deviation c and data distributions c totalparm - contains npatot, number of probes which measure c each type of parameter c partype - contains nprm, the parameters which c a probe measures, in probefile.d c nopars - contains nopas, the number of parameters measured c by each probe type c c****************************************************************** common /thename/ anlyfile common /multy/ mask(idim5,jdim5,kdim) common /squared/ tempsq common /mean/ tempav common /dist/ num(idim,jdim,kdim) common /anom/ sav(idim,jdim,kdim) common/parname1/ mtype(maxlevel) common /reqprobe/ jp(0:maxprobe) common /totalparms/ npatot(maxparm) common /nopars/ nopas(maxprobe) common /partype/ nprm(maxprobe,maxparm) common /five/ tempav5, * tempsq5, * num5(idim5,jdim5,kdim) common /stanmin/ minstan common /fracta/ numsuc,nsub common /geop/ igeot,jlev,dsurf common /extra/ ireqpars2,ireqsecond,ireqbio common /parmord/ nprec(maxcalc,maxprobe) common /levelsfull/ nfulldeps,dz(kdim2) c Initializations ik=0 ik2=0 call arrayinit2(istc,0,5,0,1) iwrite3=0 iwrite4=0 iwrite5=0 nycx=0 nycx2=0 c*************************************************************** c c Get user input information from file devfile.d c c usersd: c O nsuct - set to zero for one time period, one for all time periods c requested in devshell c O isea - time period of statistical checks: c 00 - annual c 01-12 - months c 13-16 - seasons c 17-40 - half months c O nomal - set to one if anomalies are to be calculated, zero c for full parameter values c O nomaly - read which type of anomaly: annual - 0, seasonal - 1, c monthly - 2 c O iylook -look at only certain years (2), months (1), or all values c O nsucy - set to one if successive years are desired (using c devshell2) c O numsuc - number of successive years to composite c O nyc - starting year for statistics checks c O nyc2 - ending year of statistics checks c O istc2 - information on which statistics to save. c O grids - grid size of statistical check c O numo - number of standard deviation outliers which invalidate a c profile c O ndepcheck - set to one to discard standard deviation outliers c in useable profiles c O msk = 0 - simply calculate means c = 1 - calculate means after throwing out c standard deviation outliers c = 2 - calculate means after checking for c standard deviation outliers against c current five-degree statistics c = 3 - same as (2), but standard deviation c outliers masks are automatically created c = 4 - same as (1), but standard deviation c outliers masks are automatically created c O cyon - set to 'n' if a non-standard input file name is used c set to 'y' otherwise c O cyon2 - set to 'n' if a non-standard output file name is used c set to 'y' otherwise c O unoutfile - non-standard output file name c O uninfile - non-standard input file name c O statpath - standard statistics output path c O anlypath - path for analyzed field c O iprint - set to one if the user requests individual profile c printouts, zero otherwise c O namis - the number of individual profiles requested c O itprint - set to one if statistics totals are to be printed out c O uninparm - name of non-standard input file c O kdimx - lowest level to perform checks c O npass - number of passes through standard deviation check. c One pass is all that is completed for simple c calculations, npass + 1 passes for deviation checks. c O ionegrid - set to one if only one five grid box is to be checked, c zero otherwise c O latgrid - latitude grid box for one grid check c O longrid - longitude grid box for one longitude check c O kdim1 - first level to perform checks on c O iannual - set to one if annual statistics are to be used below c a certain depth level c O kdima - depth level at which annual statistics are to be used c O nsub - which fraction of year (used only if numsuc < 0 ) c O igeot - calculate geopotential thickness (1) c calculate statistics for mixed layer depth (2) c calculate gradient (3) c -1 = calculate statistics on given density surface c O ibinsurf - surface data binning:-1=no date or geographic binning c 0=no date binning,gridbox c 1=bin by day,gridbox c 2=bin by month,gridbox c 3=bin by season,gridbox c O jlev - reference level for density calculation c O dsurf - desired density surface c c******************************************************************* call usersd( nsuct, isea, nomal, nomaly, * iylook, nsucy, numsuc, nyc, nyc2, istc2, grids, msk, numo, * ndepcheck, cyon, cyon2, uninfile, unoutfile, statpath, * anlypath, iprint, namis, itprint, uninparm, kdimx, npass, * ionegrid,latgrid,longrid, kdim1, iannual, kdima, nsub, * igeot,ibinsurf,jlev,dsurf ) write(6,*) 'xxxx',nyc,nyc2 c Set main values ik=ip2(1) if ( ireqpars2 .gt. 0 ) then ireqpars = 0 call issurf(isurf,jp(1),ilatindex,ilonindex,itimeindex) if ( ireqpars .gt. 0 .and. ip2(1) .eq. ilatindex ) * ik=ip2(ireqpars+1) elseif ( ip2(1) .eq. 8 ) then ireqpars=2 ik2=23 else ireqpars = 1 endif isoor = 1 c ipreyes - if 0 - no pre statistics (second to last run through sd check) c if 1 - record pre statistics (if msk =1,4) or c use pre statistics (if msk=2,3) ipreyes=useprestats(grids,ik,nomal,iylook) c Turn it off for now ipreyes=0 c******************************************************************* c c If this is a special five degree check (msk=3) find c which parameter is to be checked c c******************************************************************* if ( msk .eq. 3 ) then call readinparam(ip2(1)) ik = ip2(1) iyes=0 do 138 nn=1,nopas(jp(1)) 138 if ( ik .eq. nprm(jp(1),nn) ) iyes=1 if ( iyes .ne. 1 ) stop endif c if ( igeot .gt. 0 .and. ik .lt. ip2(maxlevel) ) then if ( igeot .eq. 1 .and. ireqpars2 .lt. 1 ) then write(6,*) 'Geopotential thickness is marked' write(6,*) ik stop endif c***************************************************************** c c decide which statistics to save: c the code entered by the user is a power of two code: c 0 # of Observations 1 c 1 Means 2 c 2 Standard Deviation 4 c 3,4,5 Five degree square statistics (# of observations,means,sd) 8,16,32 c (To save more than one statistic, add numbers together, c for example to save everything, enter 63) c c This code is disassembled into its component parts. If a certain c statistic is requested, a one will be placed in its array member. c c***************************************************************** do 55 n=5,0,-1 if (istc2.ge.2**n) then istc(n,0)=-1 c ipreyes only for 5-degree statistics if ( n .gt. 2 ) istc(n,ipreyes)=-1 istc2=istc2-2**N endif 55 continue cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c c IF SUCCESSIVE TIME PERIODS ARE TO BE PROCESSED, c c READ IN NEXT TIME PERIOD. THIS IS IN CONJUNCTION WITH c C SHELL DEVSHELL c c c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc if (nsuct.gt.0) then open(13,file='./dum3',status='old') read(13,*) isea close(13) endif c If time period is annual and anomalies are not calcualted, c set iannual to zero, as annual statistics are already c used below given depth c If time period is annual and anomalies are not calcualted, c set iannual to zero, as annual statistics are already c used below given depth if ( isea.eq.0 .and. nomal .eq. 0) iannual=0 c******************************************************************* c c If successive years are requested, read in next year from c dum2. This is in conjunction with shell devshell2 c c******************************************************************* if (nsucy.gt.0) then open(12,file='./dum2',status='old') if ( numsuc .lt. 0 ) then read(12,*) nycx,nsub else read(12,*) nycx endif close(12) if ( numsuc .ge. 1 ) then nyc2x = nycx+numsuc-1 else nyc2x = nycx endif else nycx = nyc nyc2x = nyc2 endif write(6,*) 'yyyy',nsucy,numsuc,nyc2x nsucyx=nsucy c************************************************************ c c convert time period to characters for file name c c************************************************************ csea='??' call insertnum('?',isea,csea,2) c************************************************************ c c do the same for grid spacing c c************************************************************ if ( grids .eq. 1.00 ) sgrid = '1.00' if ( grids .eq. .5 ) sgrid = '0.50' if ( grids .eq. .25 ) sgrid = '0.25' if ( grids .le. 0.10 ) sgrid = '0.10' ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c c CREATE FILE NAMES c c c ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c************************************************************* c c if non-standard input file, adjust parameters c only one parameter is allowed in non-standard files c c npro - number of probes c jp - probe code c nopas - number of parameters a probe measures c nprm - code of parameter measured by probe c c************************************************************* if (cyon.eq.'n' .and. nomal .eq. 0) then npro=1 jp(npro)=1 nopas(1)=1 nprm(1,1)=1 c******************************************************************* c c Open non standard files c c fileopen: c c 'xxxxx'//CHAR(0) - file name (without pathname) followed by a c C end of character string marker (\0) c 'x' - file opening instruction: c 'w' - open a new file c 'r' - open an old file to read c 'rb+' - open an old file for reading and writing c ## - file identification number c c******************************************************************* call addCend(uninfile,80,0) call fileopen(uninfile, 'rb+'//CHAR(0), 1) call addCend(uninparm,80,0) call fileopen(uninparm, 'rb+'//CHAR(0), 1) endif c************************************************************** c c if output file is non-standard, construct file names c c************************************************************** if (cyon2.eq.'n') then if ( nomal .eq. 0 ) then call insertnum('?',isea,unoutfile,80) call findlastchar(nqtit,unoutfile,80) else call insertyear(nycx,'?',unoutfile,80) call insertyear(nyc2x,'!',unoutfile,80) endif meanfile=unoutfile(1:nqtit)//'.mean'//CHAR(0) numfile=unoutfile(1:nqtit)//'.dist'//CHAR(0) standfile=unoutfile(1:nqtit)//'.sd'//CHAR(0) meanfile=unoutfile(1:nqtit)//'.mean'//CHAR(0) numfile=unoutfile(1:nqtit)//'.dist'//CHAR(0) standfile=unoutfile(1:nqtit)//'.sd'//CHAR(0) c*************************************************************** c c Read statpath and anlypath from file if requested c c*************************************************************** elseif ( statpath(1:1) .eq. '0' ) then open(14,file='./statpath.d',status='old',iostat=imaster) if ( imaster .ne. 0 ) then call clearstring(filename,80) if ( ik .eq. 1 ) then call extraname('sys.inf'//CHAR(0),'statpath_t.d'//CHAR(0), * filename) open(14,file=filename,status='old',iostat=ios) if ( ios .ne. 0 ) then call extraname('sys.inf'//CHAR(0),'statpath.d'//CHAR(0), * filename) else close(14) endif elseif ( ik .eq. 2 ) then call extraname('sys.inf'//CHAR(0),'statpath_s.d'//CHAR(0), * filename) open(14,file=filename,status='old',iostat=ios) if ( ios .ne. 0 ) then call extraname('sys.inf'//CHAR(0),'statpath.d'//CHAR(0), * filename) else close(14) endif else call extraname('sys.inf'//CHAR(0),'statpath.d'//CHAR(0), * filename) endif open(14,file=filename,status='old') endif read(14,'(a80)') statpath read(14,'(a80)') anlypath close(14) endif call findlastchar(nsps,statpath,80) call findlastchar(nspa,anlypath,80) c************************************************************* c c read in multiplier c c************************************************************* call checkwhichdepthfile(cdepthset) imultifile=-1 infile='multimask_'//cdepthset//'.d'//CHAR(0) call extraname('sys.inf'//CHAR(0),infile,filename) call fileopen(filename,'rb+'//CHAR(0),imultifile) if ( nfulldeps .gt. kdim ) then kread=kdim else kread=nfulldeps endif write(6,*) 'kread',kread call readleveli(1,1,1,1,1,idim5,jdim5,idim5,jdim5, * idim5,jdim5,1,kread,imultifile,mask) call onefilecloseft(imultifile) c********************************************************** c c create loop for calculating anomalies c c********************************************************** if (nomal.eq.1.and.nomaly.eq.2) then nloop=4 elseif (nomal.eq.1.and.nomaly.eq.3) then nloop=12 else nloop=1 endif c********************************************************* c c Figure number of Passes: c 1 pass only if simple calculations are performed (msk=0) c 2 passes if 5X5 check (msk=2,3) c 2-5 passes if statistical checks (msk=1,4) c (input number of passes plus one) c c********************************************************* if (msk.lt.1) then npass=1 elseif (msk.gt.1.and.msk.lt.4) then npass=2 else npass=npass+1 endif c*********************************************************** c c Set dimensional boundaries if only one grid square is c to be evaluated c c*********************************************************** if ( ionegrid .eq. 1 ) then idimgs = (longrid(1)-1) * 5 + 1 idimge = idimgs + 4 idim5gs = longrid(1) idim5ge = longrid(1) jdimgs = (latgrid(1)-1) * 5 + 1 jdimge = jdimgs + 4 jdim5gs = latgrid(1) jdim5ge = latgrid(1) elseif ( ionegrid .eq. 2 ) then idimgs = longrid(1) idimge = longrid(2) idim5gs = (longrid(1)/5) + 1 idim5g3 = (longrid(2)/5) + 1 jdimgs = latgrid(1) jdimge = latgrid(2) jdim5gs = (latgrid(1)/5) + 1 jdim5ge = (latgrid(2)/5) + 1 else idimgs = 1 idimge = idim idim5gs = 1 idim5ge = idim5 jdimgs = 1 jdimge = jdim jdim5gs = 1 jdim5ge =jdim5 endif c set unchanging grid boundaries. (idimgs, jdimgs, etc. may c change in main program) i1=idimgs i2=idimge j1=jdimgs j2=jdimge i15=idim5gs i25=idim5ge j15=jdim5gs j25=jdim5ge c******************************************************************* c c set parameters for small grids c c ngrid is the numbe of grids at the requested grid size c needed to equal one grid in either latitude or longitude c direction c c jquad is the number of grids necessary at the requested c grid size to make a one degree grid c c jquad is also the number of iterations of standard c deviation checks which will be required to cover c the globe c c************************************************************ ngrid=int(1/grids) idimt=ngrid*idim jdimt=ngrid*jdim idimt5= idimt/5 jdimt5= jdimt/5 if (grids.lt.1.0) then jquad=ngrid*ngrid else jquad=1 endif c************************************************************ c c If only one grid is requested, set the number of c iterations of standard deviation checks to one c c************************************************************ if ( ionegrid.gt.0 ) jquad=1 if ( cyon2 .ne. 'n' ) then c********************************************************** c c open standard files c c fileconst: c Conventional file name consists of: c c Directory pathname (statpath) + /grid size (sgrid)/ + c /field type (nt) (no field type for analyzed data)/ + c file name c c File name consists of: c one to three letter parameter code (mtype) + c one character probe codes of all probe types involved (tcar) or c zero if all probes which measure the parameter in question c are involved + c [ 'y' of 'm' if specific years or months only are involved] + c [ two digit beginning year or month (01 -99)] + c [ two digit ending year or month (01-99)] + c two digit time period (00 - 16): 00 is annual, 01 to 12 are c the months, 13-16 are the seasons + c [ an A if anomalies are being investigated] + c [ a Y, S, or M depending if yearly, seasonal, or monthly c analyses are being investigated ] c c [] brackets denote additions used only in the specific c special situation listed. c ik - parameter number c npro - number of probes c nyc, nyc2 - starting and ending years (months) c nt - type of field: c 0 - unnamed (analyzed) c 1 - means c 2 - number of observations c 3 - standard deviations c csea - character representation of time period c statpath - pathname to file c nfile - constructed file name c nomal - denotes anomaly value if set to one c nomaly -type of anomaly: c 0 - not an anomaly c 1 - yearly c 2 - seasonal c 3 - monthly c iylook - look at specific years if set to one c nsucy - successive year marker if greater than zero, c a shell is being run on successive years c sgrid - character representation of grid size c c************************************************************** ifive=0 call fileconstx( ik, npro, nycx, nyc2x, 1,csea , * statpath, meanfile, nomal, nomaly, iylook, nsucyx, sgrid, * ifive) call fileconstx( ik, npro, nycx, nyc2x, 2,csea , * statpath, numfile, nomal, nomaly, iylook, nsucyx, sgrid, * ifive) call fileconstx( ik, npro, nycx, nyc2x, 3,csea , * statpath, standfile, nomal, nomaly, iylook, nsucyx, sgrid, * ifive) endif c*********************************************************** c c Construct a generic analyzed data path c c*********************************************************** if ( nomal .gt. 0 ) then c Enter annual climatology at required levels if ( cyon .ne. 'n' ) then call fileconstx( ik, npro, nyc, nyc2, -1, '00', * anlypath, anlyfile2, 0, 0, 0, 0, sgrid, * ifive) else anlyfile2=uninparm call insertnum('N',00,anlyfile2,80) endif call addCend(anlyfile2,80,0) call fileopen(anlyfile2,'rb+'//CHAR(0),12) call readlevel(1,1,kdima,1,1,idim,jdim,idim,jdim,idim,jdim, * kdima,kdimx,12,sav) call onefileclose(12) c Set up analysis path for monthly means c switch maxprobe back to npro * if probe list desired if ( cyon .ne. 'n' ) then call fileconstx( ik, npro, nyc, nyc2, 0, 'NN', * anlypath, anlyfile, 0, 0, 0, 0, sgrid, * ifive) else anlyfile=uninparm endif endif if ( nomal .gt. 0 .or. iylook .gt. 0 ) then nyc = nycx + 1900 nyc2 = nyc2x + 1900 else nyc = nyc + 1900 nyc2 = nyc2 + 1900 endif c************************************************************** c c If this is a 5X5 check, name five grid box files c c************************************************************** if ( msk .eq. 2 .or. msk .eq. 3 ) then nprox=npatot(ik) else nprox=npro endif c THIS IS TO USE NITRATE CLIMATOLOGY FOR NO2+NO3 DATA ikreal=ik if ( ik .eq. 23 ) ik=8 c call fileconstx( ik, nprox, nycx, nyc2x, 1,csea , * statpath, mean5file, nomal, nomaly, iylook, nsucyx, sgrid, * 1) call fileconstx( ik, nprox, nycx, nyc2x, 2,csea , * statpath, num5file, nomal, nomaly, iylook, nsucyx, sgrid, * 1) call fileconstx( ik, nprox, nycx, nyc2x, 3,csea , * statpath, stand5file, nomal, nomaly, iylook, nsucyx, sgrid, * 1) ik=ikreal c******************************************************* c c open appropriate files c c******************************************************* if (istc(1,0).eq.-1) then call addCend(meanfile,80,0) call fileopen(meanfile,'rb+'//CHAR(0),istc(1,0)) endif if ( istc(0,0).eq.-1) then call addCend(numfile,80,0) call fileopen(numfile,'rb+'//CHAR(0),istc(0,0)) endif if ( istc(2,0).eq.-1) then call addCend(standfile,80,0) call fileopen(standfile,'rb+'//CHAR(0),istc(2,0)) endif if ( istc(3,0).eq.-1 .or. (ipreyes .eq. 0 .and. * (msk .eq. 2 .or. msk .eq. 3)) ) then iwrite3=1 if ( istc(3,0) .ne. -1 ) iwrite3=0 istc(3,0)=-1 call clearstring(outfile,80) outfile=num5file call addCend(outfile,80,0) call fileopen(outfile,'rb+'//CHAR(0),istc(3,0)) endif if ( istc(3,1).eq.-1 .or. (ipreyes .eq. 1 .and. * (msk .eq. 2 .or. msk .eq. 3 )) ) then iwrite3=1 if ( istc(3,1) .ne. -1 ) iwrite3=0 istc(3,1)=-1 call findlastchar(nsp,num5file,80) num5file=num5file(1:nsp)//'.pre'//CHAR(0) call fileopen(num5file,'rb+'//CHAR(0),istc(3,1)) endif if ( istc(4,0).eq.-1 .or. (ipreyes .eq. 0 .and. * (msk .eq. 2 .or. msk .eq. 3)) ) then iwrite4=1 if ( istc(4,0) .ne. -1 ) iwrite4=0 istc(4,0)=-1 call clearstring(outfile,80) outfile=mean5file call addCend(outfile,80,0) call fileopen(outfile,'rb+'//CHAR(0),istc(4,0)) endif if ( istc(4,1).eq.-1 .or. (ipreyes .eq. 1 .and. * (msk .eq. 2 .or. msk .eq. 3 )) ) then iwrite4=1 if ( istc(4,1) .ne. -1 ) iwrite4=0 istc(4,1)=-1 call findlastchar(nsp,mean5file,80) mean5file=mean5file(1:nsp)//'.pre'//CHAR(0) call fileopen(mean5file,'rb+'//CHAR(0),istc(4,1)) endif if ( istc(5,0).eq.-1 .or. (ipreyes .eq. 0 .and. * (msk .eq. 2 .or. msk .eq. 3)) ) then iwrite5=1 if ( istc(5,0) .ne. -1 ) iwrite5=0 istc(5,0)=-1 call clearstring(outfile,80) outfile=stand5file call addCend(outfile,80,0) call fileopen(outfile,'rb+'//CHAR(0),istc(5,0)) endif if ( istc(5,1).eq.-1 .or. (ipreyes .eq. 1 .and. * (msk .eq. 2 .or. msk .eq. 3 )) ) then iwrite5=1 if ( istc(5,1) .ne. -1 ) iwrite5=0 istc(5,1)=-1 call findlastchar(nsp,stand5file,80) stand5file=stand5file(1:nsp)//'.pre'//CHAR(0) call fileopen(stand5file,'rb+'//CHAR(0),istc(5,1)) endif c************************************************************* c c If only one grid is requested, read in data c c************************************************************* if (ionegrid.gt.0) then call readlevel(1,1,kdim1,1,1,idim,jdim,idim,jdim,idim,jdim, * kdim1,kdimx,istc(1,0),sntempav) call readleveli(1,1,kdim1,1,1,idim,jdim,idim,jdim, * idim,jdim,kdim1,kdimx,istc(0,0),num) do Kx = kdim1,kdimx do IN = 1,idim do JN = 1,jdim tempav(IN,JN,Kx) = sntempav(IN,JN,Kx)*1.D+00 enddo enddo enddo endif c*************************************************** c c If five grid data is needed, read in c c*************************************************** if ( msk.eq.2 .or. msk .eq. 3 ) then call readlevel(1,1,kdim1,1,1,idim5,jdim5,idim5,jdim5, * idim,jdim,kdim1,kdimx,istc(4,ipreyes),sntempav) if ( iwrite4 .eq. 0 ) then call onefilecloseft(istc(4,ipreyes)) istc(4,ipreyes)=0 endif call readleveli(1,1,kdim1,1,1,idim5,jdim5,idim5,jdim5, * idim5,jdim5,kdim1,kdimx,istc(3,ipreyes),num5) if ( iwrite3 .eq. 0 ) then call onefilecloseft(istc(3,ipreyes)) istc(3,ipreyes)=0 endif do Kx = kdim1,kdimx do IN = 1,idim5 do JN = 1,jdim5 tempav5(IN,JN,Kx) = sntempav(IN,JN,Kx)*1.D+00 enddo enddo enddo call readlevel(1,1,kdim1,1,1,idim5,jdim5,idim5,jdim5, * idim,jdim,kdim1,kdimx,istc(5,ipreyes),sntempav) if ( iwrite5 .eq. 0 ) then call onefilecloseft(istc(5,ipreyes)) istc(5,ipreyes)=0 endif do Kx = kdim1,kdimx do IN = 1,idim5 do JN = 1,jdim5 tempsq5(IN,JN,Kx) = sntempav(IN,JN,Kx)*1.D+00 enddo enddo enddo endif c******************************************************************* c c Set the start number of passes through check. c If this is a 5X5 check (msk=2) the start pass is the same c number as the end pass. c c******************************************************************* if ( msk.gt.1 .and. msk .lt. 4 ) then mpass1=npass else mpass1=1 endif if ( nomal .ne. 0 .and. cyon .eq. 'n' ) cyon='N' return end SUBROUTINE ISTOBINPROF(ibinyes,ctype) C ISTOBINPROF CHECKS IF THIS DATA WILL BE BINNED parameter (maxbin=100) character*(*) ctype character*11 firsttime character*11 binprobe(maxbin) character*80 filename dimension nspace(maxbin) save firsttime save npro,nspace save binprobe if ( firsttime .ne. 'istobinprof' ) then firsttime='istobinprof' npro=0 open(96,file='./tobebinned.d',status='old', * iostat=ios) if ( ios .ne. 0 ) then call clearstring(filename,80) call extraname('sys.inf'//CHAR(0),'tobebinned.d'//CHAR(0), * filename) open(96,file=filename,status='old',iostat=ios) endif if ( ios .eq. 0 ) then do 50 n=1,maxbin read(96,'(a11)',end=4) binprobe(n) npro=npro+1 call findlastchar(nspace(n),binprobe(n),11) 50 continue 4 continue close(96) endif endif ibinyes=0 call findlastchar(nsp,ctype,11) do 55 n=1,npro 55 if ( nsp .eq. nspace(n) .and. * ctype(1:nsp) .eq. binprobe(n)(1:nsp) ) * ibinyes=1 return end SUBROUTINE NAMEANOM(anlyfile,nomaly,kdim1,kdimxin,nlo, * iannual,nomal,nprofcount,kdima,idimgs,idimge,jdimgs, * jdimge,idimt,jdimt) C NAMEANOM OPENS THE DESIGNATED TIME PERIOD CLIMATOLOGY TO C BE USED FOR CALCULATING ANOMALIES, AND READS IN DATA c********************************************* c c Passed VAriables: c I anlyfile - generic name of analyzed field file c I nomaly - read which type of anomaly: annual - 0, seasonal - 1, c monthly - 2 c I nomal - set to one if anomalies are to be calculated, zero c for full parameter values c I iannual - set to one if annual statistics are to be used below c a certain depth level c I nlo - designated time period c I kdima - depth level at which annual statistics are to be used c I kdimxin - maximum depth level to compute anomalies c I kdim1 - first level to compute anomalies c I nprofcount - set to zero if no data for given time period c O anlyfile2 - full filename generic pathname with time period c O sav - (in common /anom/) climatological data for designated c time period c c**************************************************** parameter (idim=360,jdim=180,kdim=102,icmax=10) c character*8 firsttime character*80 anlyfile,anlyfile2,text,filename dimension ifullnum(kdim),icomma(icmax,2) c common /anom/ sav(idim,jdim,kdim) save firsttime save ilevnum,ifullnum if ( firsttime .ne. 'nameanom' ) then c get full depth listings for subsetting (if requested) firsttime='nameanom' ilevnum=0 open(12,file='./w13_depth_subset',status='old', * iostat=ios) if ( ios .eq. 0 ) then close(12) call clearstring(filename,80) call extraname('sys.inf'//CHAR(0), * 'standepths_levtofull.txt'//CHAR(0), * filename) open(12,file=filename,status='old') read(12,'(a)') text read(12,'(a)') text do 45 nn=1,kdim call clearstring(text,80) read(12,'(a)',end=46) text call findlastchar(nsp,text,80) call commafindx(nsp,icmax,ncomma,icomma,text) read(text(icomma(2,1):icomma(2,1)+icomma(2,2)-1),*) ilevnum read(text(icomma(3,1):icomma(3,1)+icomma(3,2)-1),*) * ifullnum(ilevnum) 45 continue 46 continue close(12) endif endif kdimx=0 if ( nomal .eq. 1.and. nprofcount .gt. 0 ) then if ( iannual .eq. 1 ) then kdimx=kdima-1 else kdimx=kdimxin endif else return endif k1=kdim1 k2=kdimx if ( ilevnum .gt. 0 ) then k1=ifullnum(k1) k2=ifullnum(k2) endif nad=0 anlyfile2=anlyfile do 78 nb=1,80 if (anlyfile(nb:nb).eq.'N'.and.nad.lt.1) then nad=1 if (nomaly.eq.2) then write(anlyfile2(nb:nb+1),'(i2)') nlo+12 elseif (nomaly.eq.3) then if (nlo.lt.10) then anlyfile2(nb:nb)='0' write(anlyfile2(nb+1:nb+1),'(i1)') nlo else write(anlyfile2(nb:nb+1),'(i2)') nlo endif elseif (nomaly.eq.0) then anlyfile2(nb:nb+1)='00' endif endif 78 continue call addCend(anlyfile2,80,0) call fileopen(anlyfile2,'rb+'//CHAR(0),12) write(6,*)'read;ev',idimgs,jdimgs,idimge,jdimge, * idimt,jdimt,idim,jdim call readlevel(1,1,1,idimgs,jdimgs,idimge,jdimge, * idimt,jdimt,idim,jdim,k1,k2,12,sav) c call readlevel(1,1,1,1,1,idim,jdim,idim,jdim,idim,jdim, c * k1,k2,12,sav) call onefileclose(12) if ( ilevnum .gt. 0 ) then k1=kdim1 k2=kdimx do 55 k=k1,k2 c write(6,*) 'k1',k,ifullnum(k) if ( k .ne. ifullnum(k) ) then do 60 j=1,jdim do 65 i=1,idim sav(i,j,k)=sav(i,j,ifullnum(k)) 65 continue 60 continue endif 55 continue endif return end SUBROUTINE OPENSDCOUNTFILES(ctype,npro,ik,nomal, * isea,csea,nyc,nyc2,cyc,istantotfile,istanprofile) C OPENSDCOUNTFILES OPENS STANTOTS.D, ACCOUNTING C STATITSICS FOR FINALSD RUN AND STANPROS.D, C FIRST 20 STATIONS THROWN OUT FOR STANDARD C DEVIATION OUTLIERS c****************************************************** c c Passed Variables: c ctype - dataset name c I npro - number of requested datasets c I ik - variable code for requested variable c I nomal - set to one if anomalies are to be calculated, zero c for full variable values c I isea - time period of statistical checks: c 00 - annual c 01-12 - months c 13-16 - seasons c 17-40 - half months c O csea - character representation of time period c I nyc - starting year for statistics checks c I nyc2 - ending year of statistics checks c O cyc - character representation of beginning-end year period c I istantotfile - FORTRAN file ID number for stantots file c I istanprofile - FORTRAN file ID number for stanpros file c c******************************************************** character*2 csea character*4 cyc character*11 ctype character*80 title,title2,statsdir character*160 filename nsps=0 open(9,file='./statsdirectout.d',status='old', * iostat=ios) call clearstring(statsdir,80) if ( ios .eq. 0 ) then read(9,'(a80)') statsdir close(9) else statsdir="." endif call findlastchar(nsps,statsdir,80) call clearstring(title2,80) if ( npro .eq. 1 ) then call findlastchar(nspace,ctype,11) title2(1:nspace+1)='.'//ctype(1:nspace) nspace=nspace+1 else nspace=0 endif title2(nspace+1:nspace+1)=CHAR(0) c Create program statistics file names call clearstring(title,80) if ( nomal .eq. 0 ) then csea='??' call insertnum('?',isea,csea,2) title = 'stantots.???.'//csea//title2 else cyc='??>>' call insertyear(nyc,'?',cyc,4) call insertyear(nyc2,'>',cyc,4) title = 'stantots.???.'//cyc//title2 endif call insertnum('?',ik,title,80) call clearstring(filename,160) c if ( ios .eq. 0 ) then call findlastchar(nspt,title,80) filename=statsdir(1:nsps)//'/'//title(1:nspt) nsps=nsps+1 c else c filename=title c endif open(istantotfile,file=filename,status='unknown') filename(nsps+1:nsps+8) ='stanpros' open(istanprofile,file=filename,status='unknown') return end SUBROUTINE PERFORMSDCHECK(ijz,ndsurf,iax,kdima,ik, * ndep,mpass,ndeptot,itti,itot,npro,nx,nomal,msk, * lon,lat,lon2,lat2,ity,isurf,inume,nk,ittot,npass, * nmask,isddfile,numo,igeot,minstan,minmeasures,bmiss, * ju,iuniqn,jquad,jj,iperformstat) C PERFORMSDCHECK CALCULATES MEAN AND STANDARD DEVIATION C AND PERFORMS STANDARD DEVIATION OUTLIER CHECK IF C NECESSARY. IF THE CHECK IS FAILED, STATISTICS ARE C UPDATED c******************************************************** c c Passed Variables: c I ijz - present level c I ndsurf - index for desired density surface c I iax - 1=check annual statistics below certain level,0=do not c I kdima - depth level at which annual statistics are to be used c I ik - variable code for requested variable c O ndep - 1=frequency of error (per depth/pass) c 2=total number of measurements (per depth/pass) c I mpass - present pass through standard deviation check c O ndeptot - 1=total number of errors (per pass) c 2=total number of measurements (per pass) c O itti - set to one when first value of variable is found c O itot - number of profiles with desired variable (per probe) c I npro - total number of requested data sets c I nx - present sequential data set number c I nomal - set to one if anomalies are to be calculated, zero c for full parameter values c O msk = 0 - simply calculate means c = 1 - calculate means after throwing out c standard deviation outliers c = 2 - calculate means after checking for c standard deviation outliers against c current five-degree statistics c = 3 - same as (2), but standard deviation c outliers masks are automatically created c = 4 - same as (1), but standard deviation c outliers masks are automatically created c I lon,lat - longitude/latitude grid numbers for 1x1 gridbox c I lon2,lat2 - longitude/latitude grid numbers for 5x5 gridbox c O ity - number of outliers in profile c I isurf - if 1, dataset is surface-only c O inume - levels with standard deviation outliers for profile c I nk - present bin (level for surface-only data) c O ittot - number of profiles with outliers (per probe/per pass) c I npass - number of passes through standard deviation check. c I nmask - number of each type of mask c I isddfile - file ID number for sdout[].d mask c I numo - number of standard deviation outliers which invalidate an c entire profile from use c I igeot - 1=calculate geopotential thickness c -1=perform calculations on density surfaces c 0=normal c I minstan - minimum number of standard deviations for check c I minmeasures - minumum number of measurements to perform c a valid standard deviation check c I bmiss - missing value marker c I ju - unique station number c I iuniqn - 0=no unique stats, <0=temp unique stats c >0 unique stats c I jj - sequential station order c I jquad - number of iterations (1 for 1-degree, 16 for 1/4-degree) c O iperformstat - 0=statistics calculated c 1=no statistics calculated c c*********************************************************** parameter (kdim=102,maxlevel=6000,maxcalc=200) parameter (idim5=72,jdim5=36) parameter (idim=360,jdim=180) parameter (kdim2=137,maxprobe=100,maxpas=5) dimension jp(0:maxprobe),nmask(5) dimension inume(kdim) dimension mask(idim5,jdim5,kdim) dimension sav(idim,jdim,kdim) dimension num5(idim5,jdim5,kdim) dimension num(idim,jdim,kdim) dimension ndep(2,kdim,maxpas) dimension ittot(maxprobe,maxpas) dimension itot(maxprobe,maxpas) dimension ndeptot(2,maxpas) double precision tempsq5(idim5,jdim5,kdim) double precision tempsq(idim,jdim,kdim) double precision tempav(idim,jdim,kdim) double precision tempav5(idim5,jdim5,kdim) double precision temp(kdim2,maxcalc) common/sdata/ temp common /multy/ mask common /squared/ tempsq common /mean/ tempav common /dist/ num common/anom/ sav common/five/ tempav5,tempsq5,num5 common /reqprobe/ jp iperformstat=0 itotwrite=0 if ( igeot .ne. -1 ) then ij=ijz else ij=ijz-ndsurf+1 endif c**************************************************************** c c If this profile is not within the time period and annual c statistics are not yet used at this level, skip to the c next level c c**************************************************************** c write(6,*) 'iax',iax,ijz,kdima if ( iax.eq.1.and.ijz.lt.kdima ) then iperformstat=1 return endif c if ( lon .eq. 1437 .and. lat .eq. 255 .and. ijz .eq. 1) then if ( lon .eq. 357 .and. lat .eq. 75) then c if ( ju .eq. 7806336 ) then write(6,*) 'TESTcase',nx,ju,ijz,ik,temp(ijz,ik) write(6,*) 'tc',iax,kdima endif if (temp(ijz,ik).gt.bmiss) then c**************************************************************** c c If this is the correct time period, add to counters: c c ndep - total number of values at this depth c ndeptot - total number of values c c if this is the first occurrance of a value in this profile: c c itot - number of profiles with desired parameter for c this probe c itot - number of profiles with desired parameter c c*************************************************************** if ( iax.eq.0 ) then ndep(2,ij,mpass)=ndep(2,ij,mpass)+1 ndeptot(2,mpass)=ndeptot(2,mpass)+1 if (itti.eq.0) then itot(nx,mpass)=itot(nx,mpass)+1 itot(npro+1,mpass)=itot(npro+1,mpass)+1 itti=1 endif endif c****************************************************** c c calculate anomaly by subtracting out the analyzed c field. If there is no analyzed field here, there is c no anomaly c c****************************************************** if (nomal.gt.0) then if ( lat .eq. 255 .and. lon .eq. 1437 ) * write(6,*) 'sav1',sav(lon,lat,ij),temp(ijz,ik) if (sav(lon,lat,ij).gt.bmiss) then temp(ijz,ik)=temp(ijz,ik)-sav(lon,lat,ij) c if ( jj .eq. 2044107 ) if ( lon .eq. 1437 .and. lat .eq. 255 ) * write(6,*) ij,sav(lon,lat,ij),temp(ijz,ik) else iperformstat=1 return endif if ( lat .eq. 255 .and. lon .eq. 14e7 ) * write(6,*) 'sav',sav(lon,lat,ij),temp(ijz,ik) endif c************************************************************** c c Perform statistics checks only if this is not the first c pass, simple calculations are not the only check called c for, and there are more than minmeasures data points within the c five grid box area c c************************************************************** call switchfivedegree(lon2,lat2,lon,lat,lon5,lat5) if (mpass.gt.1.and.msk.gt.0.and. * num5(lon5,lat5,ij).gt.minmeasures) then c************************************************************** c c figure out difference between data value and data mean c c************************************************************** tdif=dabs(temp(ijz,ik)-tempav5(lon5,lat5,ij)) else c************************************************************** c c If the criteria were not met, set the difference to a c value outside the checking range. c c************************************************************* tdif=(bmiss*1.0D+00)*(mask(lon5,lat5,ij)+1) endif c************************************************************* c c see if this difference is within standard deviation times c the area multiplier c c************************************************************* if ( mask(lon5,lat5,ij) .lt. minstan) * mask(lon5,lat5,ij)=minstan c if ( lon .eq. 77 .and. lat .eq. 52 .and. ij .eq. 16) then if ( jj .eq. 1085883 ) then write(6,*) 'testcase',tempav(lon,lat,ij) write(6,*) mpass,tdif,lon5,lat5,ij,tempsq5(lon5,lat5,ij), * mask(lon2,lat2,ij) write(6,*) temp(ijz,ik),tempav5(lon5,lat5,ij),num5(lon5,lat5,ij) write(6,*) tdif,tempsq5(lon5,lat5,ij)*mask(lon5,lat5,ij) endif if (mpass.gt.1.and. * tdif-1.D-06.gt.tempsq5(lon5,lat5,ij)*mask(lon5,lat5,ij)) * then if ( jj .eq. 1085883 ) write(6,*) 'flag',ij c************************************************************* c c record outlier in statistics: c c ity - number of outliers in this profile c inume - level of each outlier c ittot - total number of outliers per probe c ittot - total number of outliers c c************************************************************* if (iax.eq.0) then ity=ity+1 if ( isurf .eq. 0 ) then inume(ity)=ij else inume(ity)=nk endif ittot(nx,mpass)=ittot(nx,mpass)+1 ittot(npro+1,mpass)=ittot(npro+1,mpass)+1 c************************************************************** c c If a writing depth mask is present, write outlier information c out to file if this is the last pass through checks. c c writetomask2: c jj- profile number c itotwrite - total number of values written to files c ifn - file identification number for mask to be written to c ij - depth level of outlier c c*************************************************************** if ( mpass .eq. npass .and. * (msk .ge. 3 .or. nmask(5)-nmask(4) .gt. 0) ) then if ( iuniqn .ne. 0 .and. jquad .le. 1) then jout=ju else jout=jj endif if ( isurf .eq. 0 ) then call writetomask2(jout,itotwrite,isddfile,ij) else call writetomask2(jout,itotwrite,isddfile,nk) endif endif endif c********************************************************** c c If the number of standard deviation outliers is equal c to the number which invalidates a profile, remove c all prior levels from calculated values c c********************************************************** if (ity.eq.numo .and. isurf .eq. 0) then do 33 ij2=1,ij-1 if ( igeot .ne. -1 ) then ij2z=ij2 else ij2z=ij2+ndsurf-1 endif c********************************************************** c c If this profile is not from the correct time period c and the level is above the level at which annual c statistics are used, skip this level c c********************************************************** if ( iax.eq.1.and.ij2z.lt.kdima) goto 33 c********************************************************* c c If this level was not recorded due to a standard deviation c outlier, skip this level c c********************************************************* do 34 nx2=1,ity 34 if (inume(nx2).eq.ij2) goto 33 if (temp(ij2z,ik).gt.bmiss) then c******************************************************** c c Remove data from calculated values c c******************************************************** num(lon,lat,ij2)=num(lon,lat,ij2)-1 if (num(lon,lat,ij2).eq.0) then tempav(lon,lat,ij2)=bmiss tempsq(lon,lat,ij2)=bmiss else tempav(lon,lat,ij2)=tempav(lon,lat,ij2)-temp(ij2z,ik) tempsq(lon,lat,ij2)=tempsq(lon,lat,ij2)- * temp(ij2z,ik)*temp(ij2z,ik) endif endif 33 continue endif elseif (ity.lt.numo .or. isurf .eq. 1 ) then c********************************************************** c c Record statistics if the number of outliers in this c profile is less than the number which invalidate a c profile. c c Statistics are compiled by adding one to number c of profiles (num), adding the value at this level to c all other values at this level (tempav) and adding c the square of the value to composited squares of all c values in this grid box. c c c********************************************************** if (mpass.lt.npass.or.iax.eq.0) then num(lon,lat,ij)=num(lon,lat,ij)+1 c if ( lon .eq. 317 .and. lat .eq. 65 .and. c * ij .ge. 1 ) write(6,*) 'numcase',num(lon,lat,ij), c * tempav(lon,lat,ij),temp(ijz,ik),ijz if (tempav(lon,lat,ij).gt.bmiss) then tempav(lon,lat,ij)=tempav(lon,lat,ij)+temp(ijz,ik) tempsq(lon,lat,ij)= * tempsq(lon,lat,ij)+temp(ijz,ik)*temp(ijz,ik) else tempav(lon,lat,ij)=temp(ijz,ik) tempsq(lon,lat,ij)=temp(ijz,ik)*temp(ijz,ik) endif endif endif endif return end SUBROUTINE PREPARESDBIN(ialready,ibinyes,nk,rlat,rlon, * lvs1,kdim1,iyear0,iyear,month,iday,time,iannual, * nomal,iax,imn1x,imn2x,ihalf,nbdate,npos,grids,lat, * lon,lat2,lon2,idimt,idimgs,idimge,jdimgs,jdimge, * idim5gs,idim5ge,jdim5gs,jdim5ge,ity,itti,nsec,isurf, * lonindex,latindex,itimeindex,inume,snbin,ik,isea, * ngrid,idim,jdim,kdimx,ionegrid,jquad,ibinsurf,ibinstat, * ju) C PREPARESDBIN PREPARES INDIVIDUAL BINNED C DATA (SURFACE AND DESIGNATED CRUISE BINS) C FOR CALCULATIONS. IT ALSO MAKES SURE THE C DATA FITS DATE/POSITION CRITERIA c******************************************************** c c Passed Variables: c O ialready - if 1, record profile specific statistics c (initialized here) c I ibinyes - set to one if data from this probe should be binne c I nk - present bin set c I rlat,rlon - position of station, from header c O lvs1 - first level to look at for this station/bin c I kdim1 - first level on which to perform calculations c I iyear0 - base year for calculating julian day c I iax - 1=check annual statistics below certain level,0=do not c I imn1x,imn2x - original monthly values c I ihalf - for biweekly statistics, ihalf is set to 1 for first c half of month, 2 for second half c O nbdate - mean date for bin c O npos - geographic mean position for each bin c O iyear,month,iday,time - date/time for bin c I grids - grid size (1.00, 0.25, 5.00) in degrees c O lon,lat - longitude/latitude grid numbers for 1x1 gridbox c O lon2,lat2 - longitude/latitude grid numbers for 5x5 gridbox c I idimt - number of grid squares in longitude direction c I idimgs - starting longitude to perform checks on c I idimge - ending longitude to perform checks on c I jdimgs - starting latitude to perform checks on c I jdimge - ending latitude to perform checks on c I ity - number of outliers in profile (initialized here) c I itti - set to one when first value of variable is found c (initialized here) c I nsec - measured variable code for second header c I isurf - 1 if surface-only data, 0 otherwise c I lonindex,latindex,itimeindex - measured variable codes for c longitude, latitude, julian day c I inume - levels with standard deviation outliers for profile c (initialized here) c I snbin - measured variable data bins c I kdimx - deepest level to perform calculations c I ionegrid - set to one if only one five grid box is to be checked, c zero otherwise c I jquad - number of iterations (1 for 1-degree, 16 for 1/4-degree) c I idim,jdim - number of one-degree lat/lon spacings c I ngrid - the number of grids at the requested grid size c I isea - time period of statistical checks: c 00 - annual c 01-12 - months c 13-16 - seasons c 17-40 - half months c O ik - code for requested measured variable c O ibinstat - bin status: 0=ok, 1=skip c c******************************************************************** parameter (maxcalc=200, maxlevel=6000) parameter (ibinmax=500000,kdim2=137,kdim=102) dimension depth(maxlevel),sntemp(maxlevel,maxcalc) dimension isignif(maxlevel,maxcalc) dimension nbdate(ibinmax),npos(ibinmax) dimension inume(kdim),snbin(ibinmax,kdim) double precision temp(kdim2,maxcalc) common /thedata/ depth,sntemp common /significant/ isignif common/sdata/ temp ibinstat=0 c write(6,*) 'binin',ibinstat lvs1=kdim1 xidim=idim yjdim=jdim yjdim2=yjdim/2. if ( ibinyes .eq. 1 ) ialready=0 if ( isurf .gt. 0 ) then if ( nk .le. maxlevel ) then np=0 ic=nk dp1=depth(nk) else call overmax(maxlevel,nsec,nk,0,ic,np) dp1=sntemp(ic,np) endif lvs1=1 rlat=sntemp(ic,np+latindex) rlon=sntemp(ic,np+lonindex) xjulian=sntemp(ic,np+itimeindex) jsig=isignif(ic,np+ik) if ( dp1 .ge. 5. ) lvs1=2 temp(lvs1,ik)=sntemp(ic,np+ik) * 1.D+00 call nailujt(iyear0,iyear,month,iday,time,itimesig, * xjulian,jsig) endif c write(6,*) 'nk',nk,sntemp(nk,latindex),latindex,isurf,rlat, c * lvs1,ik,temp(lvs1,ik) c***************************************************************** c c If annual statistics are being used to check partial annual c data, set iax to one, unless the profile is from c the partial annual time period of statistical c compilation c c***************************************************************** imfact=0 imlft=0 if ( ibinyes .eq. 1 .and. isurf .eq. 0 ) then imfact=12 if ( ibinsurf .eq. 3 ) imfact=4 c cant do daily yet c if ( ibinsurf .eq. 1 ) then imlft=mod(nbdate(nk),imfact) if ( imlft .eq. 0 .and. nbdate(nk) .ne. 0 ) imlft=imfact if ( ibinsurf .eq. 3 ) imlft=(imlft-1)*3+1 month=imlft endif iax=0 if ( iannual.eq.1 .and. nomal .eq. 0 ) then iax=1 if ( month.ge.imn1x.and.month.le.imn2x ) iax=0 c if ( ju .eq. 7806336 ) if ( lat .eq. 73 .and. lon .eq. 360 ) * write(6,*) 'LONX',nk,nbdate(nk),imfact,isea,imlft,month, * imn1x,imn2x,iax endif c********************************************************* c c Check for half month data c c********************************************************* if ( ihalf .gt. 0 ) then ihmon=15 if ( month .eq. 2 ) ihmon=14 if ( (ihalf .eq. 1 .and. iday .gt. ihmon) .or. * (ihalf .eq. 2 .and. iday .le. ihmon )) then c write(6,*) 'binstat',ibinstat,ihalf,iday,ihmon,iannual if ( iannual .eq. 1 ) then iax=0 else ibinstat=1 return endif endif endif c***************************************************************** c c Transform Longitude and Latitude to grid numbers (lon and lat) c and to five grid numbers (lon2,lat2) c c For binned data, take the position at the middle of the c grid square c c***************************************************************** if ( ibinyes .eq. 1 .and. isurf .eq. 0 ) then lat=(npos(nk)/idimt)+1 lon=mod(npos(nk),idimt) c write(6,*) 'latl',nk,npos(nk),idimt,lat if ( lon .ne. 0 ) then rlon=(lon*grids)-(grids/2.) if ( rlon .gt. idim/2. ) rlon=rlon-idim rlat=((lat*grids)-(grids/2.))-(yjdim2) else rlon=-(grids/2.) rlat=((lat*grids)-1.-(grids/2.))-(jdim/2.) endif endif izink=0 if ( npos(nk) .eq. 3887140 ) izink=1 c if ( npos(nk) .eq. 3887140 ) c * write(6,*) 'rlat',nk,npos(nk),idimt,grids,ngrid, c * idim,jdim,rlat,rlon,ibinyes,idimgs,idimge, c * jdimgs,jdimge call grid(rlat,rlon,lat,lon,grids) call grid(rlat,rlon,lat2,lon2,grids*5.) if ( lon .eq. 360 .and. lat .eq. 73 ) izink=1 if ( izink .eq. 1) * write(6,*) 'LONX2',nk,nbdate(nk),imfact,isea,imlft,month, * imn1x,imn2x,rlat,rlon,jdimgs,jdimge,idimgs,idimge,lat,lon, * kdim1,kdimx,snbin(nk,3),npos(nk),ju c***************************************************************** c c Check that lat and lon are within prescribed limits c c shift latitude and longitude into 1-360, 1-180 range for c longitude and latitude respectively c c***************************************************************** c write(6,*) 'binstat',ibinstat,lat,lon,jdimgs, c * jdimge,idimgs,idimge,rlat,rlon if ( lat .lt. jdimgs .or. lat .gt. jdimge .or. * lon .lt. idimgs .or. lon .gt. idimge ) then ibinstat=1 return endif if ( ibinyes .eq. 1 .and. isurf .eq. 0 ) then do 425 nb1=kdim1,kdimx temp(nb1,ik)=snbin(nk,nb1) * 1.D+00 425 continue endif if ( ionegrid.lt.1 .and. jquad.gt.1 ) then lat= lat -jdimgs+1 lon= lon -idimgs+1 lat2= lat2 - jdim5gs+1 lon2= lon2 - idim5gs+1 endif c****************************************************************** c c Initialize counters: c ity - number of outliers in profile c itti - set to one when one value of parameter is found c c****************************************************************** ity=0 itti=0 call arrayinit(inume,1,kdim) return end SUBROUTINE PREPARESDDATA(jj0,nmask,jj,cc,numlevels, * iyear,month,iday,time,icruise,rlat,rlon,jj1, * ireqtotz,isoor,irecord,ioldrecord,m,ik,ik2,ienda, * ishallow,ideep,xovertop,xunderbot,nsurfs,idimt, * isea,nyc,nyc2,itoto,icounter,grids,latindex,lonindex, * itimeindex,kall,nlvs,iyear0,kdim1,kdimx,ibinyes,kdima, * idensefile,iposdfile,idenserecord,iolddenserecord,jjd, * icruisefile,icruisefile2,icruisefilep,ireqtots,isurf, * nbdate,npos,snbin,iyx0f,iyx1f,imx0,imx1,ibinsurf,iend, * ialready,idimgs,idimge,jdimgs,jdimge,iuniqn,ju,mpass, * nomal,nomaly,ireqpars,iannual,nquad,ieof) C PREPARESDDATA READS IN DATA AND PREPARES C IT FOR MEAN/STANDARD DEVIATION CALCULATIONS c******************************************************* c c Passed Variables: c I jj0 - number of stations examined for this dataset c I nmask - number of each type of mask c O jj - sequential order of station in dataset c O numlevels - number of levels with measured data in profile c O cc,iyear,month,iday,time,icruise,rlat,rlon - c cruise/date/position information for station c O jj1 - present position in dataset header file c ireqtotz - set to zero for density surfaces, MLD, ireqtots otherwise c I isoor - set to 1 for standard levels (0=observed levels) c O irecord - desired position in measured variable file c set to -1 if no measured variable data for this c station c O ioldrecord - present position in measured variable file c I m - present dataset code c I ik - desired measured variable code c I ik2 - secondary desired measured variable code c O ienda - set to 1 if all depth masks have been read in completely c O ishallow - index for shallowest density surface c O ideep - index for deepest density surface c I xovertop,xunderbot - markers for above ocean surface or c below ocean bottom for density surface c O nsurfs - number of density surfaces between surface and c ocean bottom c I idimt - number of grid squares in longitude direction c I isea - time period of statistical checks: c 00 - annual c 01-12 - months c 13-16 - seasons c 17-40 - half months c I nyc - starting year for statistics checks c I nyc2 - ending year of statistics checks c O itoto, icounter - used for depthmask subroutine, not c used elsewhere in program c I grids - grid size (1.00, 0.25, 5.00) in degrees c I latindex,lonindex,itimeindex - measured variable code c for latitude, longitude, julian day c O kall - number of bins (1 for non-binned profile data, c number of levels for surface-only data c O nlvs - number of levels (1 for surface-only data, c number of levels (numlevels) for profile data c O iyear0 - base year for surface-only data julian day calculation c I kdim1 - first level to perform calculations c I kdimx - deepest level to perform calculations c I ibinyes - set to one if data from this probe should be binned c I kdima - depth level at which annual statistics are to be used c O jjd - present position in density surface position file c I idensefile - file ID numbers for density surface files c I iposdfile - file ID numbers for position density surface files c I idenserecord,iolddenserecord - holds position for each c density surface file opened c I icruisefile,icruisefile2,icruisefilep - file ID number for c cruise index. Used for binned data c I ireqtots - number of measured variables to read c I isurf - 1=surface-only data, 0=no c O nbdate - dates for binned data (julian day) c O npos - positions for binned data c I iyx0f,iyx1f,imx0,imx1 - beginning,ending years/months c for subsetting surface-only data c I ibinsurf - surface data binning:-1=no date or geographic binning c 0=no date binning,gridbox c 1=bin by day,gridbox c 2=bin by month,gridbox c 3=bin by season,gridbox c >3=bin only for specified years c O iend - set to 1 if all profile masks have been read in completely c O ialready - if 1, record profile specific statistics c (initialized here) c I idimgs - starting longitude to perform checks on c I idimge - ending longitude to perform checks on c I jdimgs - starting latitude to perform checks on c I jdimge - ending latitude to perform checks on c I iuniqn - 0=no unique stats, <0=temp unique stats c >0-unique stat c O ju - unique station number c I mpass - pass through mean calculation c I nomaly - read which type of anomaly: annual - 0, seasonal - 1, c monthly - 2 c I nomal - set to one if anomalies are to be calculated, zero c for full parameter values c I nquad - present quadrant (for grid resolution > 1) c O ieof - set to 1 if dataset has been completed c set to 2 to skip present station c c********************************************************************* parameter (maxlevel=6000, maxcalc=200) parameter (maxprobe=100, maxmask=100, bmiss=-1.E10) parameter (kdim2=137,maxpot=100,maxsurf=1000) parameter (maxcc=200,ibinmax=500000) parameter (ncruisebinmax=500000,kdim=102) parameter (maxparm=100,noffset=200) parameter (maskclosenum=300000000) character*2 cc,ncc(maxcc) character*11 ctype(maxprobe) character*13 firsttime common /countryc/ ncc common /probes/ ctype dimension depth(maxlevel),sntemp(maxlevel,maxcalc) dimension ip2(0:maxlevel), numlevels(0:1) dimension nopas(maxprobe),isur(maxcalc) dimension nprec(maxcalc,maxprobe), itoto(0:maxcalc) dimension irecord(0:2*maxparm), ioldrecord(0:2*maxparm) dimension ifn(maxmask), ipar2(maxmask), ipro2(maxmask) dimension nmask(5),isignif(maxlevel,maxcalc) dimension nprm(maxprobe,maxparm),dz(kdim2) dimension icounter(maxlevel,maxcalc) dimension idensefile(0:maxparm,maxpot),iposdfile(maxpot) dimension idenserecord(0:maxparm),iolddenserecord(0:maxparm) dimension npos(ibinmax),snbin(ibinmax,kdim) dimension nbdate(ibinmax),icrbin(ncruisebinmax) dimension itx(maxprobe) double precision temp(kdim2,maxcalc) common /thedata/depth,sntemp common /nopars/ nopas common /extra/ ireqpars2,ireqsecond,ireqbio common /surrogate/ isur common /twoparm/ itwo(0:maxparm) common /parmord/ nprec common /eachparm/ numeach(6),numtot2,numtot3,ntax,numsecset common /reqps/ ip2 common /partype/ nprm common /significant/ isignif common /geop/ igeot,jlev,dsurf common /levelsfull/ nfulldeps,dz common /parm2/ ipar2,ifn,ipro2 common/sdata/ temp common /geopnt/ igpnt,xlon,xlat,xdist save firsttime save icrbin,ncruisebin,itx,ipassnow,imx0now,imx1now, * ixbt_adjust,mnow,ixbtis,ikcheck,ik2check, * levelmaxbias,levelmaxxbt,levelmaxmbt,nquadnow ieof=0 ialready=0 igpnt=0 isecond =1 if ( ireqtotz .eq. 0 ) isecond=0 if ( firsttime .ne. 'preparesddata' ) then firsttime='preparesddata' ncruisebin=0 call arrayinit(itx,1,maxprobe) ipassnow=mpass imx0now=imx0 imx1now=imx1 nquadnow=nquad call finddepthlevelnum(dz,kdim2,nfulldeps,100.,levelmaxbias) call finddepthlevelnum(dz,kdim2,nfulldeps,250.,levelmaxmbt) call finddepthlevelnum(dz,kdim2,nfulldeps,950.,levelmaxxbt) c write(6,*) 'maxlev',levelmaxbias,levelmaxmbt,levelmaxxbt ikcheck=ik if ( ik .gt. ip2(maxlevel) ) ikcheck=isur(ik) ik2check=ik2 if ( ik2 .gt. ip2(maxlevel) ) ik2check=isur(ik) c check for XBT adjustment ixbt_adjust=1 open(9,file='./no_adjust_XBT_data',status='old',iostat=ios) if ( ios .eq. 0 ) then close(9) ixbt_adjust=0 endif mnow=0 ixbtis=0 endif if ( mnow .ne. m ) then mnow=m ixbtis=0 if ( ctype(m)(1:3) .eq. 'XBT' ) ixbtis=1 endif if ( ipassnow .ne. mpass ) then ipassnow=mpass ncruisebin=0 endif if ( nomal .gt. 0 .and. (imx0now .ne. imx0 .or. * imx1now .ne. imx1) ) then imx0now=imx0 imx1now=imx1 ncruisebin=0 endif c reset encountered cruise count if ( nquad .ne. nquadnow ) then nquadnow=nquad ncruisebin=0 endif if ( mod(jj0,10000) .eq. 0 ) write(6,*) jj0 c******************************************************* c c If there are sequential reading masks, c call seqread to find the next proper profile c number to input, otherwise input sequentially c c******************************************************* if (nmask(2).gt.0) then c******************************************************* c c seqread finds the next profile number to input c based on sequential masks.c c c jj - chosen profile number c nmask(2) - number of inclusive and exclusive sequential c masks c nmask(1) - number of inclusive sequential masks c ieof - set to one if all masks files have reached their c end of file. If this is set to one, exit c profile loop. c iend - end of file marker (1) for each separate mask c ifn - file identification numbers for mask files c c********************************************************* iend=0 if ( jj .ge. maskclosenum) jj=0 call seqread(jj,nmask(2),nmask(1),ieof,iend,ifn) if (ieof.eq.1) return else jj=jj0 endif c if ( jj .ge. 1421282 .and. jj .le. 1422200 ) then c write(6,*) 'jjix',jj c endif c********************************************************* c c Call readdata to read in profile. c c cc - NODC country code c temp - parameter data c jj - profile number c iyear,month,iday - date of profile c time - GMT time of profile (starting time) c icruise - NODC cruise number (0 for other formats)d c rlat,rlon - geographic position of profile c jj1 - present position of header file pointer c depth - observed depths c nopas - number of parameters measured by a probe c ioldrecord - present position of data file pointers c ireqtots - number of requested parameters c ip2 - parameter codes c isoor - observed (0) or standard (1) levels c maxlevel - maximum number of levels c irecord - desired file position c maxcalc - maximum number of parameters c m - probe type code c itwo - set to one for second header, two for bio data c bmiss - missing value marker c 0 - do not read second header (1 to read) c 1 - which set of files to read from c ieof - marked with a one if end of file has been reached c ireqpars2 - requested number of calculated parameters c isignif - the number of significant figures in a data value c c******************************************************** call readdata(cc,sntemp,numlevels,jj,iyear, * month,iday,time,icruise,rlat,rlon,jj1,depth, * nopas(m),ioldrecord,ireqtotz,ip2,isoor, * maxlevel,irecord,maxcalc,m,nprec,itwo,bmiss, * isecond,1,ieof,ireqpars2,isignif) c write(6,*) 'readd',jj,ieof if (ieof.eq.1) return c Find unique station number if ( iuniqn .ne. 0 ) then call uniqfind(jj,ju,m,itx(m)) else ju=0 endif if ( icruise .eq. 34062 ) write(6,*) 'jjcc',jj if ( ju .eq. 4683251 ) write(6,*) 'isreadin',numlevels(1),kdimx c******************************************************** c c Get mixed-layer depth information c c******************************************************** if ( igeot .eq. 2 ) then call readmixedlayer(ju,sntemp(1,ik),isignif(1,ik),numlevels(1), * irecord((nopas(m)*isoor)+nprec(ik,m))) c if ( sntemp(1,ik) .gt. bmiss ) write(6,*) 'ju',ju,sntemp(1,ik),ik c write(6,*) 'ju',ju,sntemp(1,ik),ik c********************************************************* c c Calculate gradient c c********************************************************* elseif ( igeot .eq. 3 ) then iwrite=0 if ( ju .eq. 8214623 ) iwrite=1 if ( iwrite .eq. 1 ) write(6,*) 'jj',jj call setgradientx(numlevels,sntemp,maxlevel,maxcalc, * bmiss,ik,iwrite) endif c********************************************************* c c If the desired parameter was not measured in this c profile, skip to the next profile c c********************************************************* c write(6,*) 'irec',m,nopas(m),isoor,ikcheck,nprec(ikcheck,m), c * irecord((nopas(m)*isoor)+nprec(ikcheck,m)),ik2check, c * irecord((nopas(m)*isoor)+nprec(ik2check,m)),ik if ( irecord((nopas(m)*isoor)+nprec(ikcheck,m)) .lt. 0 ) then if ( ik .ne. 8 .and. ik .ne. 25) then ieof=2 return elseif ( ik .ne. 25 ) then if ( irecord((nopas(m)*isoor)+nprec(ik2check,m)) .lt. 0 ) then ieof=2 return endif endif endif c********************************************************* c c Read in density surface variables c c********************************************************* c write(6,*) 'igeot',igeot,ireqpars,ip2(1) if ( igeot .eq. -1 ) then c if ( jj .eq. 106017 ) then call filepos(iposdfile(1),0,1,ipospos) write(6,*) 'position',ipospos c endif call readsurfaces(nopas(m),jj,jjd,iposdfile(1), * ireqpars,ip2,idenserecord,iolddenserecord, * idensefile,1,maxpot,sntemp,maxlevel, * maxparm,maxcalc,m,nprec,ishallow,ideep, * nsurfs,bmiss,xovertop,xunderbot,ieof) numlevels(1)=nsurfs call grid(rlat,rlon,latt,lonn,grids) c if ( jj .eq. 106017 ) then if ( latt .eq. 119 .and. lonn .eq. 271 ) then write(6,*) 'idense',idenserecord(nprec(ik,m)),ik write(6,*) 'idense',idenserecord(25) write(6,*) 'ip2',ireqpars,ip2(1) write(6,*) 'sec',nopas(m),jj,jjd endif if ( idenserecord(nprec(ik,m)) .lt. 0 ) then ieof=2 return endif c if ( jj .eq. 106017 ) if ( latt .eq. 119 .and. lonn .eq. 271 ) * write(6,*) 'hereout', jj,sntemp(16,ip2(1)), * ip2(1),kdim1,kdimx,xovertop,xunderbot c convert density surfaces which are too heavy or too light c to missing value do 736 nz=kdim1,kdimx do 737 ipi=1,ireqpars 737 if ( sntemp(nz,ip2(ipi)) .le. xovertop +1.E-5 .and. * sntemp(nz,ip2(ipi)) .ge. xunderbot -1.E-5 ) * sntemp(nz,ip2(ipi))=bmiss 736 continue endif c********************************************************* c c Special for nitrate c c********************************************************* if ( ik .eq. 8 .and. * irecord((nopas(m)*isoor)+nprec(ik,m)) .lt. 0 .and. * irecord((nopas(m)*isoor)+nprec(ik2,m)) .ge. 0 ) then do 555 nilev=1,numlevels(isoor) call overmax(maxlevel,numeach(1),nilev,0,ic,np) sntemp(ic,ik+np)=sntemp(ic,ik2+np) 555 continue irecord((nopas(m)*isoor)+nprec(ik,m))= * irecord((nopas(m)*isoor)+nprec(ik2,m)) endif c********************************************************** c c Special adjustment for XBTs c c********************************************************** itsprobe=0 itsprobe2=0 if ( sntemp(numtot2+29,numeach(1)) .gt. bmiss ) then itsprobe=sntemp(numtot2+29,numeach(1)) endif if ( sntemp(numtot2+noffset+5,numeach(1)) .gt. bmiss ) then itsprobe2=sntemp(numtot2+noffset+5,numeach(1)) if ( itsprobe2 .eq. 204. .or. itsprobe2 .eq. 210 .or. * itsprobe2 .eq. 219 .or. itsprobe2 .eq. 213 .or. * itsprobe .eq. 236) then itsprobe=0 itsprobe2=-2 endif endif if ( itsprobe .eq. 2 ) then c remove XBT depths deeper than 700 meters if ( numlevels(1) .gt. levelmaxxbt .and. * itsprobe2 .ne. -2 ) then do 577 nxx=levelmaxxbt,numlevels(1) 577 sntemp(nxx,1)=bmiss numlevels(1)=levelmaxxbt endif if ( ixbt_adjust .eq. 1 ) then call XBTbiaschoose(5,iyear,month,numlevels,isoor,ifix,0,1) c do 678 nxx=1,numlevels(1) c78 write(6,*) jj,dz(nxx),sntemp(nxx,1) c call XBTbiasant(itsprobe,sntemp,maxlevel,maxcalc, c * numlevels(1),iyear,levelmaxbias) endif elseif ( itsprobe .eq. 1 ) then if ( numlevels(1) .gt. levelmaxmbt ) then do 578 nxx=levelmaxmbt+1,numlevels(1) 578 sntemp(nxx,1)=bmiss numlevels(1)=levelmaxmbt endif if ( ixbt_adjust .eq. 1 ) then call MBTbiasant(itsprobe,sntemp,maxlevel,maxcalc, * numlevels(1),iyear,levelmaxbias) endif endif c********************************************************* c c Special binning for high density data c c********************************************************* if ( icruise .eq. 34062 ) write(6,*) 'ibin',ibinyes,ieof if ( ibinyes .eq. 1 ) then call cruisecon(cc,icruise,icrcon) do 735 ncf=1,ncruisebin 735 if ( icrbin(ncf) .eq. icrcon ) ieof=2 if ( ieof .eq. 2 ) return call binprofs(jj,jj1,irecord,ioldrecord,m,bmiss, * maxlevel,maxcalc,numlevels,itwo,depth, * sntemp,isignif,nprec,nopas,maxparm, * icrcon,ibinmax,kdim,snbin,nbins,npos, * grids,idimt,icruisefile,icruisefilep,nbdate, * icruisefile2,isea,nyc,nyc2,maxprobe, * ibinsurf,ik,ireqtots,ip2,ireqpars2,iannual, * itoto,icounter,nmask,ifn,ipar2,ipro2,maxmask) ncruisebin=ncruisebin+1 icrbin(ncruisebin)=icrcon if ( ncruisebin .gt. ncruisebinmax ) then write(6,*) 'Number of binned cruises exceeds maximum', * ncruisebinmax write(6,*) 'Program will be terminated' stop endif c********************************************************* c c Mask out depths as per depth files in maskchoice.d c c Only do this if not all masks have been completely c read in. c c depthmask: c maxlevel - maximum number of levels c jj - present record number c ifn - file identification numbers for masks c nmask(2) - starting mask number for depth masks (minus one) c nmask(3) - ending mask number for depth masks c temp - parameter data c ipar2 - parameter which corresponds to each mask c icounter - counter of depths which are masked c for each parameter c itoto - counter of number of profiles which are masked c for each parameter c ienda - counter of number of end of files reached c ireqtots - number of requested parameters c ip2 - parameter code numbers for requested parameters c numlevels(1) - number of standard levels c bmiss - missing value marker c c*************************************************************** elseif (ienda.lt.nmask(3)-nmask(2)) then call depthmask(maxlevel,jj,ifn,nmask(2),nmask(3),sntemp, * ipar2,icounter,itoto,ienda,ireqtotz,ip2, * numlevels(isoor),bmiss,ireqpars2) endif if ( jj .eq. 1085883 ) then write(6,*) 'zip',sntemp(1,2),ienda,nmask(3),nmask(2),ip2(1) endif c****************************************************************** c c Calculate geopotential thickness if desired c c****************************************************************** if ( igeot .eq. 1 .and. isurf .eq. 0 ) * call geothick(dz,kdim2,sntemp,maxlevel,maxcalc,ik,rlat, * numlevels(1)) c****************************************************************** c c Bin surface data c c****************************************************************** if ( isurf .gt. 0 ) then if ( ju .eq. 7806336 ) write(6,*) 'subsurf', * iyear,iyx0f,iyx1f,imx0,imx1 call subsurf(iyear,iyx0f,iyx1f,imx0,imx1, * 0,0,idimgs,idimge,jdimgs,jdimge, * depth,sntemp,numlevels(0),nlevout,numeach(1), * ireqtots, ip2,latindex,lonindex, * itimeindex,grids,1) numlevels(0)=nlevout call binsurf(depth,sntemp,maxcalc,maxlevel, * numlevels(0),latindex,lonindex,itimeindex, * ik,grids,numeach(1),ibinsurf) kall=numlevels(0) c write(6,*)'sntempout',sntemp(1,2),kall nlvs=1 iyear0=iyear elseif ( ibinyes .eq. 1 ) then kall=nbins nlvs=kdimx c****************************************************************** c c Transfer the data from single precision to double precision c c****************************************************************** else do KI = kdim1,kdimx temp(KI,ik) = sntemp(KI,ik) * 1.D+00 enddo kall=1 c****************************************************************** c c If the number of standard levels in this profile is more than c the maximum number of levels, set it to the maximum number c of levels (either user maximum or system maximum. c c****************************************************************** if ( jj .eq. 1085883 ) write(6,*) 'k0', * kdimx,numlevels(1) if (numlevels(1).gt.kdimx .and. * igeot .ne. -1) numlevels(1)=kdimx if (numlevels(1).gt.kdim .and. igeot .ne. -1) * numlevels(1)=kdim nlvs=numlevels(1) if ( igeot .eq. -1 ) then nlvs=kdimx kdima=numlevels(1)+1 endif if ( jj .eq. 1085883 ) write(6,*) 'k', * kdimx,numlevels(1) endif return end SUBROUTINE PRINTOT(mpass,npro,itot,ctf,jp, * ctype,ntype,itty,ittot,nfreq,ndep,dz,ndeptot, * nbad,nalt,sea,numo,itprint,nquad,jquad,maxfreq, * istantotfile) C C PRINTOT PRINTS ACCOUNTING STATISTICS FOR FINALSD C c********************************************************** c c Passed Variables; c mpass - number of pass to print statistics for c npro - number of probes requested c itot - number of profiles measuring desired parameter c for each probe c ctf - parameter name c npr - probe codes c ctype - probe names c ntype - first probe code c itty - number of profiles with outliers c ittot - total number of outliers c nfreq - number of outliers per profile c ndep - outliers per depth c dz - standard depths c ndeptot - total number of outliers c nbad - number of profiles discarded c nalt - number of altered profiles c sea - time period c numo - number of outliers which make a profile bad c itprint - set to one if statistics totals are to be printed out c jquad - number of iterations (1 for 1-degree, 16 for 1/4-degree) c nquad - present section for quarter-degree (1-16) c istantotfile - FORTRAN file ID number for stantots file c c*********************************************************** c*********************************************************** c c Parameters: c maxprobe - maximum number of probes in system c maxparm - maximum number of measured parameters c kdim - number of standard depth levels for which c statistics are compiled c kdim2 - number of depth levels c maxpas - maximum number of passes through statistical check c c************************************************************ parameter (maxprobe=100,maxparm=100,kdim=102,kdim2=137,maxpas=5) c********************************************************* c c Character Arrays: c ctype - probe names c ctf - parameter name c sea - time period c c********************************************************* character ctype(maxprobe)*11,ctf*15,sea*9 c********************************************************** c c Arrays: c c itot - number of profiles with desired parameter for c each probe c itty - number of profiles with outliers c ittot - total number of outliers c nfreq - number of outliers per profile c ndep - outlier per depth counter c dz - standard depths c ndeptot - total number of outliers c jp - probe types c nalt - number of profiles altered c nbad - number of profiles discarded c c********************************************************* dimension itot(maxprobe,maxpas),itty(maxprobe,maxpas), * ittot(maxprobe,maxpas),nfreq(0:maxfreq+1,maxpas), * ndep(2,kdim,maxpas),dz(kdim2),ndeptot(2,maxpas), * jp(0:maxprobe),nalt(maxprobe,maxpas),nbad(maxprobe,maxpas) if ( mpass .eq. 0 .or. itprint .eq. 0 .or. * nquad .ne. jquad ) return c*********************************************************** c c Write out pass number, number of profiles for this c parameter for this time period c c*********************************************************** write(istantotfile,203) mpass 203 format(/,3x,' Pass #',i1,' Through Standard Deviation Check') if (npro.gt.1) then write(istantotfile,805) itot(npro+1,mpass),sea,ctf write(istantotfile,305) do 786 n=1,npro write(istantotfile,305) itot(n,mpass),ctype(jp(n)) 786 continue else nxw=0 do 88 n=1,npro 88 if (jp(n).eq.ntype) nxw=n write(istantotfile,305) itot(nxw,mpass),sea,ctype(ntype),ctf endif 305 format(/,3x,i7,1x,a9,1x,a4,' ',a11,' Profiles') 805 format(/,3x,i7,' Total ',a9,' ',a11,' Profiles') c************************************************************* c c Write number of profiles containing outliers, c per probe if necessary. c c************************************************************* write(istantotfile,806) itty(npro+1,mpass) 806 format(/,3x,i5,' Profiles Contain Standard Deviation * Outliers') if (npro.gt.1) then do 1980 n=1,npro 1980 write(istantotfile,1810) ctype(jp(n)),itty(n,mpass) endif 1810 format(5x,a4,5x,i8) c************************************************************* c c Write total number of outliers c c************************************************************* write(istantotfile,807) ittot(npro+1,mpass) 807 format(/,i8,' Total Standard Deviation Outliers') if (npro.gt.1) then do 2980 n=1,npro 2980 write(istantotfile,2810) ctype(jp(n)),ittot(n,mpass) endif 2810 format(5x,a4,5x,i8) c************************************************************* c c Write number of discarded profiles c c************************************************************* write(istantotfile,907) nbad(npro+1,mpass),numo 907 format(/,i8,' Discarded Profiles (',i2,' or More Outliers)') if (npro.gt.1) then do 2981 n=1,npro 2981 write(istantotfile,2881) ctype(jp(n)),nbad(n,mpass) endif 2881 format(5x,a4,5x,i8) c************************************************************** c c Write number of altered profiles c c************************************************************** write(istantotfile,909) nalt(npro+1,mpass) 909 format(/,i8,' Altered Profiles ') if (npro.gt.1) then do 2991 n=1,npro 2991 write(istantotfile,2891) ctype(jp(n)),nalt(n,mpass) endif 2891 format(5x,a4,5x,i8) c************************************************************** c c Write out outliers frequency table c c************************************************************** write(istantotfile,340) 340 format(/,' Outliers per Profile Frequency') write(istantotfile,341) 341 format(/,5x,'Outs',5x,'Profiles') write(istantotfile,342) nfreq(0,mpass) 342 format(5x,'None',4x,i7) do 132 i=1,10 write(istantotfile,343) i,nfreq(i,mpass) 132 continue 343 format(7x,i2,4x,i7) write(istantotfile,344) nfreq(11,mpass) 344 format(6x,'>10',4x,i7) c********************************************************** c c Write out errors per level table c c********************************************************** write(istantotfile,605) 605 format(/,10x,' Errors per Depth Level') write(istantotfile,606) 606 format(/,5x,'Num',5x,'Depth',5x,'Outliers',5x,'# of Profiles') do 433 k=1,kdim if (ndep(2,k,mpass).gt.0) then write(istantotfile,607) k,dz(k),ndep(1,k,mpass),ndep(2,k,mpass) endif 433 continue 607 format(6x,i2,5x,f5.0,2(4x,i7)) write(istantotfile,608) ndeptot(1,mpass),ndeptot(2,mpass) 608 format(3x,'Total',10x,2(4x,i8)) return end SUBROUTINE PRINTPRO(ctp,nsea,ctf,nobsdp,rlon,rlat,ik,iyr, * mon,iday,cc,icruz,dz,inume,ity,lon,lat,lon2,lat2,bmiss, * sea,mpass,jj,igeot,ndsurf,istanprofile) C PRINTPRO PRINTS PROFILES FOR FINALSD c******************************************************** c c Passed Variables: c I ctp - dataset name c I nsea - time period (00=annual, 13-16=season, 1-12=month) c I nobsdp - number of levels with measured data c I rlat,rlon - geographic location c I ik - desired measured variable c I iyr,mon,iday - date of measurements c I cc, icruz - cruise information c I dz - standard depth levels c I inume - levels with standard deviation outliers for profile c I ity - number of outliers in profile c I lon,lat - longitude/latitude grid numbers for 1x1 gridbox c I lon2,lat2 - longitude/latitude grid numbers for 5x5 gridbox c I bmiss - missing value marker c I sea - character representation of time period c I mpass - present pass through standard deviation check c I jj - sequential station order c I igeot - 1=calculate geopotential thickness c -1=perform calculations on density surfaces c 0=normal c I ndsurf - index for desired density surface c I istanprofile - FORTRAN file ID number for stanpros file c c************************************************************* c*************************************************************** c c Parameters: c idim - number of longitudinal spacings c jdim - number of latitudinal spacings c kdim - number of depth levels to compute means c kdim2 - number of depth levels in OCL system c idim5 - number of five degree longitudinal spacings c jdim5 - number of five degree latitudinal spacings c maxcalc - maximum number of measured and calculated c parameters c c**************************************************************** parameter (idim=360,jdim=180,kdim=102,kdim2=137,maxparm=100) parameter (idim5=72, jdim5=36, maxcalc=200) c***************************************************************** c c Character Arrays: c ctp - probe name c ctf - 1-6 letter parameter name c cc - country code c sea - time period c c***************************************************************** character ctp*11, ctf*6, cc*2, sea*9 c***************************************************************** c c Arrays: c dz - standard depth levels c inume - depth of each standard deviation outlier c tempav5 - five degree parameter means c tempsq5 - five degree parameter standard deviations c temp - profile parameter data c mask - multiplication factor for each five degree square c at each depth c num5 - number of observations of parameter c c***************************************************************** dimension dz(kdim2),inume(kdim) double precision tempav5 double precision tempsq5 double precision temp common /sdata/ temp(kdim2,maxcalc) common /multy/ mask(idim5,jdim5,kdim) common /five/ tempav5(idim5,jdim5,kdim), * tempsq5(idim5,jdim5,kdim), * num5(idim5,jdim5,kdim) c************************************************************* c c Write out header information for probe, pass number, c time period and probe name. c c************************************************************ write(istanprofile,8303) write(istanprofile,1799) write(istanprofile,1800) rlon,rlat,iyr,mon,iday,cc,icruz,jj 1799 format(/,5x,'Longitude',2x,'Latitude',2x,'Year', * 2x,'Month',2x,'Day',2x,'CC',2x,'Cruise#',2x,'Probe#') 1800 format(6x,f8.2,2x,f8.2,2x,i4,4x,i2,4x,i2,2x,a2,2x,i6,3x,i7/) write(istanprofile,'(a7,i1)') 'Pass # ',mpass write(istanprofile,860) sea,ctp 860 format(/,6x,a9,1x,a11,' Statistics ') write(istanprofile,8803) ctf 8803 format(/,4x,'Num',5x,'Depth',5x,a6,4x,'GM',4x,'StanDev', * 7x,'Mean',3x,'Profiles') c*********************************************************** c c Write out parameter value at each depth, number of c total observations at this depth, multiplication factor, c mean and standard deviaton c c*********************************************************** do 502 j0=1,nobsdp if ( igeot .ne. -1 ) then j=j0 else j=j0+ndsurf-1 endif if (temp(j,ik).gt.bmiss.and. * num5(lon2,lat2,j0).gt.0) then write(istanprofile,803) j0,dz(j0),temp(j,ik), * mask(lon2,lat2,j0),tempsq5(lon2,lat2,j0), * tempav5(lon2,lat2,j0),num5(lon2,lat2,j0) elseif ( temp(j,ik) .gt. bmiss ) then write(istanprofile,903) j0,dz(j0),temp(j,ik), * mask(lon2,lat2,j0),num5(lon2,lat2,j0) elseif ( num5(lon2,lat2,j0) .gt. 0 ) then write(istanprofile,904) j0,dz(j0),mask(lon2,lat2,j0), * tempsq5(lon2,lat2,j0), * tempav5(lon2,lat2,j0), * num5(lon2,lat2,j0) else write(istanprofile,905) j0,dz(j0),mask(lon2,lat2,j0) endif 502 continue 501 continue 803 format (5x,i2,5x,f5.0,5x,f6.2,5x,i1,5x,f6.3,5x,f6.3,5x,i6) 903 format (5x,i2,5x,f5.0,5x,f6.2,5x,i1,5x,' -----',5x,' -----',5x,i6) 904 format (5x,i2,5x,f5.0,5x,' -----',5x,i1,5x,f6.3,5x,f6.3,5x,i6) 905 format (5x,i2,5x,f5.0,5x,' -----',5x,i1,5x,' -----',5x, * ' -----',5x,i6) c************************************************************** c c write out number of errors in file, season of file, c and depth levels of error, for cross reference. c c************************************************************** write(istanprofile,804) ity 804 format(/5x,i3,' Standard Deviation Errors') do 33 i=1,ity if (i.eq.1) then write(istanprofile,853) inume(i) 853 format(5x,' at Depth Level(s) : ',i4) else write(istanprofile,850) inume(i) endif if (i.eq.ity) then write(istanprofile,852) 852 format(/) 850 format(26x,i4) endif 33 continue write(istanprofile,8303) 8303 format(/,79('-')) return end SUBROUTINE PROBESDCLOSE(ifn,maxmask,nmask,ireqtotz,isoor, * icruisefile,icruisefile2,icruisefilep,m,iposdfile, * idensefile,jpotlev,npots,ireqpars,maxparm,maxpot, * ip2,maxlevel,densesurf,maxsurf,ibinyes,igeot) C PROBESDCLOSE CLOSES MASKS AND FILES FOR C PRESENT PROBE c***************************************************** c c Passed Variables: c I ifn - file identification numbers for masks c maxmask - maximun number of masks c I ireqtotz - set to zero for density surfaces, ireqtots otherwise c I isoor - set to 1 for standard levels (0=observed levels) c I icruisefile,icruisefile2,icruisefilep - file ID numbers c for cruise indexes c I m = present dataset code c I iposdfile - file ID numbers for position density surface files c I idensefile - file ID numbers for density surface files c I jpotlev - level of reference for each potential surface c I npots - number of potential surfaces c I ireqpars - number of measured variables to read in c I maxparm - maximum number of measured variables c I maxpot - maximum number of potential surfaces c I ip2 - measured variable codes c I maxlevel = maximum number of depth levels c I densesurf - individual density surfaces c I maxsurf - maximum number of density surfaces c I ibinyes - set to one if data from this dataset should be binned c c*************************************************************** dimension ip2(0:maxlevel),ifn(maxmask) dimension jpotlev(maxpot),densesurf(maxsurf) dimension idensefile(0:maxparm,maxpot),iposdfile(maxpot) dimension nmask(5) c******************************************************************** c c Close all parameter and header files with fileclose c c fileclose: c c ifn - file identification numbers for masks (no masks here) c nmask(5) - number of masks c ireqpars - number of parameters c isoor - standard (1) or observed (0) c 1(0) - second header closing instruction c 0 - no second header c 1 - second header c 1 - file set number c c close date mask c c********************************************************************* if ( igeot .ne. -1 .and. igeot .ne. 2) then call fileclose(ifn,nmask(5),ireqtotz,isoor,1,1) elseif ( igeot .eq. 2 ) then call onefileclose(1) call fileclose(ifn,nmask(5),ireqtotz,isoor,1,1) else call closedensefiles(npots,jpotlev,ireqpars,ip2,maxlevel, * idensefile,maxparm,maxpot,iposdfile,m) call onefileclose(1) endif if ( ibinyes .eq. 1 ) then call onefilecloseft(icruisefile) call onefilecloseft(icruisefilep) call onefilecloseft(icruisefile2) endif return end SUBROUTINE READMIXEDLAYER(ju,sntemp,isign,levels,irecord) C READMIXEDLAYER READS THE MIXED LAYER DEPTH FOR THE GIVEN C UNIQUE STATION NUMBER AND SETS THE NECESSARY VARIABLES parameter (imiss=-1000,bmiss=-1.E10) character*14 firsttime character*80 filename, cmixname save firsttime save imixedfile if ( firsttime .ne. 'readmixedlayer' ) then call mixedlayerchoice(ivar,xmixed,del0,del3, * cmixname,nspmix) cmixname=cmixname(1:nspmix)//'.s'//CHAR(0) call extraname('sys.inf/mld'//CHAR(0),cmixname,filename) imixedfile=-1 call fileopen(filename,'rb+'//CHAR(0),imixedfile) firsttime='readmixedlayer' endif ieof=0 call fileend(imixedfile) call filepos(imixedfile,0,1,irecnum) call posfile(ju,imixedfile) call readmask(mld,imixedfile,ieof) if ( mld .eq. imiss ) then sntemp=bmiss isign=0 irecord=-1 else sntemp=mld call findsignif(isign,0,sntemp) irecord=1 levels=1 endif return end SUBROUTINE RESETMASKS(isea2,iylook,nomal,leof,nmask,nmasko) C RESETMASKS RESETS MASKS TO ORIGINAL C REQUESTED MASKS FROM MASKCHOICE.D c*************************************************** c c Passed Variables: c I isea2 - time period (0 if annual statistics below certain level) c I iylook -look at only certain years (2), months (1), or all values c O nomal - set to one if anomalies are to be calculated, zero c for full variable values c I leof - > 0 if there was data for given time period c nmask - number of each type of mask c c**************************************************** parameter (maxmask=100) dimension ifn(maxmask), ipar2(maxmask), ipro2(maxmask) dimension nmask(5) common /parm2/ ipar2,ifn,ipro2 if ( (isea2 .ne. 0 .or. iylook .ne. 0 .or. nomal .ne. 0) * .and. leof .ne. 0) then do 490 np=1,nmasko 490 call maskerase(1,1,nmask,ifn,ipar2,ipro2) endif return end SUBROUTINE SETFORQUAD(nquad,ngrid,lonquadi,latquadi, * idimgs,jdimgs,idimge,jdimge,idim5gs,jdim5gs,idim5ge, * jdim5ge,ionegrid) C SETFORQUAD SETS FOR EACH SUCCESSIVE SECTION OF C THE QUARTER-DEGREE CHECK/CALCULATION c********************************************************* c c Calculate lonquadi and latquadi. These are the base c numbers for each small segment of the world c c********************************************************* if ( mod(nquad,ngrid).eq.0 ) then lonquadi=ngrid-1 latquadi=nquad/ngrid-1 else lonquadi=mod(nquad,ngrid)-1 latquadi=nquad/ngrid endif c****************************************************** c c set grid boundaries c c****************************************************** if ( ionegrid.lt.1 ) then idimgs=lonquadi*360+1 jdimgs=latquadi*180+1 idimge=idimgs+359 jdimge=jdimgs+179 idim5gs= (idimgs-1)/5 +1 jdim5gs= (jdimgs-1)/5 +1 endif return end SUBROUTINE SETFORSPECPROBE(m,ik,ik2,iyes,ibinyes,icruisefile, * icruisefilep,icruisefile2,ncruisebin,isurf,latindex, * lonindex,itimeindex,isoor,ireqtots,ireqpars,ioldrecord, * ienda,jj,jj1,ireqtotz,nmask,nomal,nyc,nyc2,iend,iuniqn,cyon, * msk) C SET VARIABLES BASED ON GIVEN PROBE c m - probe code c ik - main measured variable c ik2 - secondary measured variable (used for nitrate+nitrite c when calculating means for NO3. c iyes - set to one if probe measured given variable c ibinyes - set to one if data from this probe should be binned c icruisefile,icruisefilep,icruisefile2 - if probe data is binned c by cruise, these files are opened to search by cruise c ncruisebin - number of cruises already binned c initialized here c isurf - set to one for surface data c latindex,lonindex,itimeindex - variable codes for latitude, c longitude,julian year day c isoor - 0 for observed level, 1 for standard level data c ireqtots - total number of variables to examine c ireqpars - number of measured variables to examine c ioldrecord - present file position for variable files c ienda - end of file marker (initialized here) c iuniqn - >0 if probe has unique station numbers c jj - present sequential station number (initialized here) c ireqtotz - set to zero for density surfaces, MLD, ireqtots otherwise c I nyc,nyc2 - beginning, ending anomaly time period c I nmask - number of each type of mask c I nomal - set to one if anomalies are to be calculated, zero c for full parameter values c cyon - 'y' if conventional files, 'n' if not conventionally named c I msk = 0 - simply calculate means c = 1 - calculate means after throwing out c standard deviation outliers c = 2 - calculate means after checking for c standard deviation outliers against c current five-degree statistics c = 3 - same as (2), but standard deviation c outliers masks are automatically created c = 4 - same as (1), but standard deviation c outliers masks are automatically created parameter (maxlevel=6000,maxprobe=100) parameter (maxcalc=200,maxparm=100,maxmask=100) character*1 cyon character*11 ctype(maxprobe) character*80 filename,cfl dimension ip2(0:maxlevel),jp(0:maxprobe) dimension nopas(maxprobe),isur(maxcalc) dimension nprec(maxcalc,maxprobe) dimension ioldrecord(0:2*maxparm) dimension ifn(maxmask), ipar2(maxmask), ipro2(maxmask) dimension nmask(5) dimension nprm(maxprobe,maxparm) common /nopars/ nopas common /extra/ ireqpars2,ireqsecond,ireqbio common /surrogate/ isur common /parmord/ nprec common /eachparm/ numeach(6),numtot2,numtot3,ntax,numsecset common /reqps/ ip2 common /probes/ ctype common /partype/ nprm common /parm2/ ipar2,ifn,ipro2 common /maskfile/ cfl(maxmask) common /reqprobe/ jp common /twoparm/ itwo(0:maxparm) common /geop/ igeot,jlev,dsurf c************************************************************ c c Initialize counters: c c m - probe type number c ioldrecord - position of parameter file pointer c ienda - number of depth masks past end of file c jj - profile number c jj1 - position of header file pointer c c*********************************************************** iyes=0 if ( cyon .eq. 'N' ) iyes=1 do 188 nn=1,nopas(m) 188 if ( ik .eq. nprm(m,nn) ) iyes=1 if ( iyes .ne. 1 ) then do 187 nn=1,nopas(m) 187 if ( isur(ik) .eq. nprm(m,nn) ) iyes=1 endif if ( iyes .ne. 1 ) return c See if this data type needs to be binned c If so, open cruise masks call istobinprof(ibinyes,ctype(m)) ncruisebin=0 icruisefile=-1 icruisefilep=-1 icruisefile2=-1 if ( ibinyes .eq. 1 ) then call clearstring(filename,80) call maskname('cruiseindex'//CHAR(0),0,filename,m) call fileopen(filename,'rb+'//CHAR(0),icruisefile) call maskname('cruisemask.p'//CHAR(0),0,filename,m) call fileopen(filename,'rb+'//CHAR(0),icruisefilep) call maskname('cruisemask.2'//CHAR(0),0,filename,m) call fileopen(filename,'rb+'//CHAR(0),icruisefile2) endif c*********************************************************** c c Check if this is a surface data file c c*********************************************************** call issurf(isurf,m,latindex,lonindex,itimeindex) if ( isurf .eq. 1 ) then isoor=0 ikadd=0 if ( ik2 .gt. 0 ) ikadd=1 ireqtots=4+ikadd ip2(1)=ik if ( ikadd .gt. 0 ) ip2(2)=ik2 ip2(2+ikadd)=latindex ip2(3+ikadd)=lonindex ip2(4+ikadd)=itimeindex else isoor=1 ireqtots=ireqpars+ireqpars2 ip2(1)=ik ip2(2)=ik2 endif do 755 nik=1,ireqtots ioldrecord((nopas(m)*isoor)+nprec(ip2(nik),m)) = 0 755 continue ioldrecord(2*nopas(m)+1)=0 if ( ireqpars2 .gt. 0 ) ioldrecord(2*nopas(m)+3)=0 ienda=0 jj = 0 jj1=1 c set ireqtotz to ireqtots. (This will only change if c data on density surfaces is required or MLD is c requested (igeot=2)) c ireqtotz=ireqtots if ( igeot .eq. 2 ) ireqtotz=0 c******************************************************************* c c Open all masks for probe. c c maskset: c 1 - number of probes to open files for for each mask c m - probe type code c 0 - open old files (one to open new files) c nmask - number of each type of mask c 1. number of inclusive sequential masks to be read c 2. number of inclusive + exclusive sequential to be read c 3. nmask(2) + number of depth exclusion masks to be read c 4. nmask(3) + number of sequential masks to be written c 5. nmask(4) + number of depth exclusion masks to be written c c Common Block used here: reqprobe c Common Block set here: parm2 c c Initialize all mask holder values. A mask holder value c is a C variable which holds the last value read in from c a mask. c c maskinit: c c nmask(5) - total number of requested masks c iend - denotes end of masks c c******************************************************************* call maskset(1,m,0,nmask) if ( nomal .gt. 0 .and. msk .ne. 2 .and. msk .ne. 3) * call anomalymasks(nyc,nyc2,nmask,m,ik,5,msk) if ( nomal .gt. 0 .and. msk .ne. 2 .and. msk .ne. 3 * .and. nyc .eq. nyc2 ) * call anomalymasks(nyc,nyc2,nmask,m,ik,1,msk) call maskinit(nmask(5)+1,iend) call ucheck(m,iuniqn) c******************************************************************* c c If this is a calculated parameter, set mask parameter number to c surrogate parameter number. c c******************************************************************* if ( ik .gt. numeach(1) ) then do 882 n=1,nmask(3) 882 if ( ipar2(n).le.isur(ik)) ipar2(n)=ik endif c**************************************************************** c c open OCL data files c c isoor - observed (0) or standard (1) levels c m - desired probe code c ireqtots - total number of requested parameters c 1 - file set one c 1 - second header opener instruction c -2 - no second header, no header, no depth c -1 - second header, no header, no depth c 0 - no second header, header, depth c 1 - second header, header, depth c 2 - no second header, header, no depth c 3 - no second header, no header, depth c c*************************************************************** if (cyon.ne.'n') call opener(isoor,m,ireqtotz,1,1) return end SUBROUTINE SETFORSPECPROBEDENSE(cyon,igeot,npots,jpotlev, * idensefile,maxpot,iposdfile,ireqtotz,maxparm, * idenserecord,ioldenserecord,m,ireqpars,ip2,maxlevel, * densesurf,maxsurf,jjd) C SETFORSPECPROBEDENSE SETS DENSITY SURFACE C VARIABLES FOR EACH SPECIFIC DATASET c********************************************************** c c Passed Variables: c I cyon - set to 'n' if a non-standard input file name is used c set to 'y' otherwise c I igeot - 1=calculate geopotential thickness c -1=perform calculations on density surfaces c 0=normal c I npots - number of potential surfaces c I jpotlev - level of reference for each potential surface c I maxpot - maximum number of potential surfaces c I iposdfile - file ID numbers for position density surface files c ireqtotz - set to zero for density surfaces, ireqtots otherwise c I maxparm - maximum number of measured variables c I idensefile - file ID numbers for density surface files c I idenserecord,iolddenserecord - holds position for each c density surface file opened (initialized here) c I m - present probe code c I ireqpars - number of measured variables to read in c I ip2 - variable codes c I maxlevel - array elements in ip2 c I densesurf - individual density surfaces c I maxsurf - maximum number of density surfaces c O jjd - present position of density file c c*********************************************************** parameter (maxprobe=100,maxcalc=200) parameter (maxlevloc=6000) character*1 cyon character*3 mtype dimension jpotlev(maxpot),densesurf(maxsurf) dimension idensefile(0:maxparm,maxpot),iposdfile(maxpot) dimension idenserecord(0:maxparm),iolddenserecord(0:maxparm) dimension ip2(0:maxlevel) dimension nprec(maxcalc,maxprobe) common /parmord/ nprec common /parname1/ mtype(maxlevloc) c Shut off ireqtotz. For density surfaces, no measured c variables should be read in if ( igeot .eq. -1 ) ireqtotz=0 c*************************************************************** c c Open density surface files c c*************************************************************** if ( cyon .ne. 'n' .and. cyon .ne. 'N' .and. igeot .eq. -1 ) then inewfile=0 call opendensefiles(npots,jpotlev,ireqpars,ip2,maxlevel, * idensefile,maxparm,maxpot,inewfile,iposdfile,m,1) jjd=1 do 731 nz=1,ireqpars idenserecord(nprec(ip2(nz),m))=-1 iolddenserecord(nprec(ip2(nz),m))=0 731 continue endif return end SUBROUTINE SETGRADIENTX(numlevels, * sntemp,maxlevel,maxcalc,bmiss,ik,iwrite) parameter (kdim=137) character*12 firsttime dimension numlevels(0:1),dzdiff(kdim) dimension sntemp(maxlevel,maxcalc) common /levelsfull/ nfulldeps,dz(kdim) save firsttime save dzdiff if ( firsttime .ne. 'setgradientx' ) then do 60 nnx=2,kdim 60 dzdiff(nnx)=dz(nnx)-dz(nnx-1) firsttime='setgradientx' endif do 75 nnx=numlevels(1),2,-1 if ( sntemp(nnx,ik) .gt. bmiss .and. * sntemp(nnx-1,ik) .gt. bmiss ) then if ( iwrite .eq. 1 ) then write(6,*) 'snt',nnx,sntemp(nnx,ik),sntemp(nnx-1,ik),dzdiff(nnx) endif sntemp(nnx,ik)=(sntemp(nnx-1,ik)-sntemp(nnx,ik))/dzdiff(nnx) if ( iwrite .eq. 1 ) write(6,*) 'sntx',sntemp(nnx,ik) else sntemp(nnx,ik)=bmiss endif 75 continue sntemp(1,ik)=bmiss return end SUBROUTINE SETSDDATE(m,iyx,imx,nmask,nmaskg,nmasko, * ik,iend,isea2,iylook,nomal,jyear,jyearmax,isurf, * nyc,leof,iyx0f,iyx1f,imx0,imx1) C SETSDDATE SETS DATE MASK FOR SPECIFIC DATE PERIOD C AND ADDS THIS MASK TO ALL OTHER USED MASKS c m - probe type code c iyx - year under examination c imx - month under examination c nmask - list of types of mask from maskchoice.d c nmaskg - set to one if geographic mask is present c nmasko - number of sequential inclusive masks, including c geography and date c ik - variable code c iend - initialized mask position pointer c isea2 - time period (0 if annual statistics below certain level) c iylook - look at certain years c nomal - > 0 if looking at anomalies c jyear,jyearmax - first/last possible year c isurf - 1=surface-only data,0=no c nyc - first year period to perform check c leof - set to zero if no data for this time period c I iyx0f,iyx1f,imx0,imx1 - beginning,ending years/months c for subsetting surface-only data parameter (maxmask=100) parameter (ndays=366) character*80 filename dimension ifn(maxmask), ipar2(maxmask), ipro2(maxmask) dimension monthday(13) dimension nmask(5) common /parm2/ ipar2,ifn,ipro2 data monthday/0,31,60,91,121,152,182,213,244,274,305,335,366/ leof=0 iyx0f=0 iyx1f=0 imx0=0 imx1=0 if ( nomal .gt. 0 ) then iyx0f=iyx iyx1f=iyx imx0=imx imx1=imx endif if ( isea2 .ne. 0 .or. iylook .ne. 0 .or. nomal .ne. 0) then idatemask=-1 call tempfile(idatemask) if ( isurf .eq. 0 ) then idatefile=-1 call maskname('dates.l'//CHAR(0),0,filename,m) call fileopen(filename,'rb+'//CHAR(0),idatefile) c write(6,*) 'here' c write(6,*) 'datefile',js,je,iyx,imx,leof call loopread(js, je, idatefile, iyx, imx, leof) c write(6,*) 'loopread',idatefile,leof call onefilecloseft(idatefile) if ( leof .gt. 0 ) then call onefilecloseft(idatemask) leof=-1 return endif c write(6,*) js,je,iyx,imx,leof,idatemask ndatecount=0 do 300 nj=js,je 300 call writetomask(nj,ndatecount,idatemask) call filerewind(idatemask) else if ( imx .eq. 13 ) then mday1=1 mday2=ndays else mday1=monthday(imx) mday2=monthday(imx+1) endif if ( nyc .gt. 1 ) then iyx0=iyx-jyear+1 iyx1=iyx-jyear+1 iyx0f=nyc iyx1f=nyc else iyx0=1 iyx1=jyearmax-jyear+1 iyx0f=0 iyx1f=0 endif call gridset(mday1,mday2,iyx0,iyx1, * ndays,m,1,idatemask,ndatecount) c write(6,*) 'gd',idatemask endif c************************************************************ c c Combine date mask with any other input inclusive masks c c************************************************************ c write(6,*) 'nmasks is ',nmask if ( ndatecount .gt. 0 ) then c write(6,*) 'tempcombine_1',idatemask,ndatecount,nmasko,nmaskg, c * nmask call tempcombine(nmasko,nmaskg,idatemask,ndatecount,nmask) c write(6,*) 'tempcombine',nmask,ifn(2),ipar2(2),ipro2(2), c * nmasko call maskinit(nmasko,iend) c See if there is any data here (leof > 0) iendxx=0 call readmask(jxxx,ifn(1),iendxx) c write(6,*) 'jxxx is',jxxx,iendxx,ifn(1) call fileend(ifn(1)) call filepos(ifn(1),0,1,leof) call filerewind(ifn(1)) c write(6,*) 'fileend',ifn(1),leof if ( leof .eq. 0 ) leof=-2 call maskinit(1,idumend) else call onefilecloseft(idatemask) leof=-1 endif endif return end SUBROUTINE SETSDGEOGRAPHY(nmaskg,ionegrid,jquad, * idimgs,idimge,jdimgs,jdimge,idimt,m,ik,nmask) C SET GEOGRAPHIC REGION (FOR QUARTER-DEGREE OR C SELECTED 5-DEGREE BOX. INSERT REGION INTO NORMAL C MASK SCHEME parameter (maxmask=100) dimension nmask(5) dimension ifn(maxmask),ipro2(maxmask),ipar2(maxmask) common /parm2/ ipar2,ifn,ipro2 c******************************************************************* c c Call onegridset to make a temporary mask of all profiles c within the boundaries of the specified area. c c onegridset: c nmask - number of each type of mask c 1. number of inclusive sequential masks to be read c 2. number of inclusive + exclusive sequential to be read c 3. nmask(2) + number of depth exclusion masks to be read c 4. nmask(3) + number of sequential masks to be written c 5. nmask(4) + number of depth exclusion masks to be written c c idimgs - starting longitude number c idimge - ending longitude number c jdimgs - starting latitude number c jdimge - ending latitude number c c idim - number of longitude spacings c c m - probe type c c******************************************************************* nmaskg=0 igeofile=-1 if (ionegrid .ge. 1 .or. jquad.gt.1 ) then call tempfile(igeofile) call gridset(idimgs,idimge,jdimgs,jdimge,idimt, * m,0,igeofile,ngeocount) c************************************************************ c c Combine geographic mask with any user input inclusive masks c c************************************************************ write(6,*) 'tempcg',nmaskg,igeofile,ngeocount,nmask,ifn(1) call tempcombine(nmaskg,0,igeofile,ngeocount,nmask) write(6,*) 'tempcg',nmaskg,igeofile,ngeocount,nmask, * ifn(1),ifn(2) endif return end SUBROUTINE SETTIMELOOPS(imn1,imn2,isea,ibinsurf,numsuc, * nsub,ihalf,isea2,iannual,nomal,nomaly,nlo,nyc,nyc2, * jyear,jyearmax,iylook,imn1x,imn2x) C SET YEAR AND MONTH TIME LOOPS c IO imn1 - first month(season) to check c IO imn2 - last month(season) to check c I isea - time period of statistical checks: c 00 - annual c 01-12 - months c 13-16 - seasons c 17-40 - half months c IO ibinsurf - surface data binning:-1=no date or geographic binning c 0=no date binning,gridbox c 1=bin by day,gridbox c 2=bin by month,gridbox c 3=bin by season,gridbox c I numsuc - number of successive years to composite c I nsub - which fraction of year (used only if numsuc < 0 ) c O ihalf - for biweekly statistics, ihalf is set to 1 for first half c of month, 2 for second half of month c O imn1x,imn2x - original monthly values. Used for statistics c compilation, as imn, imn2 may change during program run c O isea2 - if annually compiled statistics are used below a certain c level (iannual=1), set iseas to zero, otherwise to isea c I nomal - set to one if anomalies are to be calculated, zero c for full parameter values c I nomaly - read which type of anomaly: annual - 0, seasonal - 1, c monthly - 2 c I iylook -look at only certain years (2), months (1), or all values c IO nyc - starting year for statistics checks c IO nyc2 - ending year of statistics checks character*12 firsttime save firsttime c The first part was moved from inputf subroutine to c keep all manipulations of year and month loops together if ( firsttime .ne. 'settimeloops' ) then firsttime='settimeloops' c*********************************************************** c c set month loop according to time period. c If annual time period set loop ends to 13 to read in c only annual data. c c*********************************************************** if ( isea .gt. 12 .and. isea .lt. 17 ) then imn1 = ( isea - 13 ) * 3 + 1 imn2 = imn1 + 2 elseif ( isea .gt. 0 .and. isea .lt. 13 ) then imn1 = isea imn2 = isea if ( ibinsurf .eq. 3 ) ibinsurf=2 elseif ( isea .gt. 16 ) then imn1=-(((isea-16)/2)+mod(isea,2)) imn2=imn1-mod(isea,2) else imn1 = 13 imn2 = 13 endif if ( numsuc .lt. 0 ) then numsucx = -numsuc imn1= (nsub-1)* numsucx + 1 imn2=nsub*numsucx endif c******************************************************************* c c Set imn1x to beginning month and imn2x to ending month. c These values will be used to compile needed statistics and c data counts when annual statistics are being used to c calculate standard deviation outliers for a lesser period c ( a month or season ) c c******************************************************************* ihalf=0 if ( imn1 .lt. 0 ) then if ( imn1 .ne. imn2 ) then ihalf=1 else ihalf=2 endif imn1=-imn1 imn2=imn1 endif c write(6,*) 'timeloop',imn1,imn2,ihalf isea2=isea imn1x=abs(imn1) imn2x=abs(imn2) if ( iannual.eq.1 .and. nomal .eq. 0) isea2=0 endif c*************************************************************** c c set month loop if anomalies are being calculated. c If anomalies are being calculated, only the months (season) c involved is pertinent. c c If no anomalies are being calculated, user input time c period is used to determine starting and ending months c (in inputf subroutine) c c*************************************************************** if ( nomal .gt. 0 ) then if ( nomaly .eq. 2 ) then c imn1 = (nlo-1) * 3 c imn2 = imn1 + 2 c This is wrong. Fix before doing seasonal anomalies imn1=13 imn2=16 elseif ( nomaly .eq. 3 .and. imn1 .eq. imn2 .and. * numsuc .ge. 0 ) then c imn1 =nlo c imn2 = nlo imn1=1 imn2=12 endif endif c**************************************************************** c c If no specific years are being investigated, set the years c to the largest range possible. c c**************************************************************** if (iylook.eq.0) then nyc=jyear nyc2=jyearmax endif c*************************************************************** c c If annual statistics are being compiled for a partial annual c period, set the months to annual (13) c c*************************************************************** if (iannual.eq.1 .and. nomal .eq. 0 ) then imn1=13 imn2=13 endif c*************************************************************** c c If all data is being looked at set years, to one, c months to one, and profiles counters to maximum and minimum c c*************************************************************** if ( isea2 .eq. 0 .and. iylook .eq. 0 .and. nomal .eq. 0) then nyc = 1 nyc2 = 1 imn1 = 1 imn2 = 1 endif return end SUBROUTINE SORTQMASKS(nquad,jquad,msk,npro, * istanoutfile,jp,maxprobe,fileout,fileoutd) C SORTQMASKS SORTS QUARTER-DEGREE C SDOUT MASKS c*********************************************************** c c Passed Variables: c I nquad - present section for quarter-degree (1-16) c I jquad - number of iterations (1 for 1-degree, 16 for 1/4-degree) c I msk = 0 - simply calculate means c = 1 - calculate means after throwing out c standard deviation outliers c = 2 - calculate means after checking for c standard deviation outliers against c current five-degree statistics c = 3 - same as (2), but standard deviation c outliers masks are automatically created c = 4 - same as (1), but standard deviation c outliers masks are automatically created c I npro - number of datasets requested c I istanoutfile - C file ID number for standard deviation outlier c masks for quarter-degree. This is necessary c because these masks must be kept open through c 16 iterations through database c I jp - requested dataset codes c I maxprobe - maximum number of datasets c O fileout - sdout[].s file name c O fileoutd - sdout[].d file name c c********************************************************* character*80 fileout,fileoutd,filetemp,filetempd dimension jp(0:maxprobe) dimension istanoutfile(maxprobe,2) c***************************************************** c c Sort 1/4 degree sdout masks if necessary c c***************************************************** if ( nquad .eq. jquad .and. jquad .gt. 1. * .and. msk .ge. 3 ) then call clearstring(filetemp,80) filetemp=fileout call addCend(filetemp,80,0) call clearstring(filetempd,80) filetempd=fileoutd call addCend(filetempd,80,0) do 91 nn=1,npro isdsfile=istanoutfile(jp(nn),1) isddfile=istanoutfile(jp(nn),2) itempfile=isdsfile call masksift(-1,isdsfile,itempfile,jp(nn),ik,filetemp) itempfile=isddfile call masksift(-2,isddfile,itempfile,jp(nn),ik,filetempd) call onefileclose(isdsfile) call onefileclose(isddfile) 91 continue endif return end SUBROUTINE SWITCHFIVEDEGREE(lon2,lat2,lon,lat,lon5,lat5) C SWITCHFIVEDEGREE SWITCHES FIVE DEGREE BOX FOR C STATISTICS ONLY FOR ONE PARTICULAR BOX NEAR C THE GULF STREAM C FIVE DEGREE GRIDBOX 59,27 SOUTHEAST CORNER C WILL BE CONSIDERED PART OF GRIDBOX 60,27 C THESE ONE-DEGREE SQUARES: C (295,131), (294,131), (293,131), (292,131), C (291,131), (295,132), (294,132), (293,132), C (295,133), (294,133), (295,134) lon5=lon2 lat5=lat2 if ( lon5 .eq. 59 .and. lat5 .eq. 27 ) then if ( lon .eq. 295 .and. (lat .eq. 131 .or. * lat .eq. 132 .or. lat .eq. 133 .or. * lat .eq. 134) ) then lon5=60 elseif ( lon .eq. 294 .and. (lat .eq. 131 .or. * lat .eq. 132 .or. lat .eq. 133) ) then lon5=60 elseif ( lon .eq. 293 .and. (lat .eq. 131 .or. * lat .eq. 132) ) then lon5=60 elseif ( (lon .eq. 292 .or. lon .eq. 291) * .and. lat .eq. 131 ) then lon5=60 endif endif return end SUBROUTINE TEMPCOMBINE(nmaskx,nmasky,inewfile,ncount,nmask) C TEMPCOMBINE COMBINES A GIVEN INCLUSIVE MASK WITH ALL OTHER C EXISTING MASKS. parameter (maxmask=100) dimension ifnc(maxmask),iprof(maxmask) dimension ifn(maxmask),ipro2(maxmask),ipar2(maxmask) dimension nmask(5) common /parm2/ ipar2,ifn,ipro2 c nmaskx - number of sequential inclusive masks c nmasky - a previous mask has already been combined c with all existing sequential inclusive masks. In this c case the new mask need only be combined with this c previous mask c inewfile - file ID number (C) for mask to be combined c with existing sequential inclusive masks c ncount - count of stations in mask to be combined c nmask - count of each type of existing mask (nmask(1)=sequential c inclusive) nmaskx=nmask(1) if ( nmasky .gt. 0 ) then do 310 nzp=2,nmask(1) 310 call fileend(ifn(nzp)) nmaskx=1 endif do 321 np=1,nmaskx ifnc(1)=inewfile iprof(1)=ncount ifnc(2)=ifn(np) call fileend(ifn(np)) call filepos(ifn(np),0,1,iprof(2)) call filerewind(ifn(np)) irewind=-1 c write(6,*) 'every',ifnc(1),ifnc(2),iprof(1),iprof(2) call everymask(2,0,iprof,ifnc,irewind) call onefilecloseft(inewfile) call fileend(ifn(np)) c write(6,*) 'np',ifn(np),np,ifnc(1),ifnc(2) call masksee(1,nmask,ifn,ipar2,ipro2,ifnc(1), * 0,ipro2(1)) ipx=0 call fileend(ifn(1)) call filepos(ifn(1),0,1,ipx) call filerewind(ifn(1)) c write(6,*) 'ifn',ifn(1),ipx,iprof(1),iprof(2),iprof(3) 321 continue c write(6,*) 'here2',nmask(1),ifn(1),ifn(2) if ( nmask(1) .eq. 0 ) then call masksee(1,nmask,ifn,ipar2,ipro2,inewfile, * 0,ipro2(1)) nmaskx=1 endif return end SUBROUTINE UPDATESDSTATS(iax,ity,ialready,nx,mpass, * npro,numo,isurf,ndepcheck,nalt,nbad,ndep,ndeptot, * itti,nfreq,maxfreq,namis,nam,numlevels,isoor, * igeot,kdim1,kdimx,isea,m,ik,rlon,rlat,iyear,month, * iday,cc,icruise,inume,lon,lat,lon2,lat2,bmiss,season, * jj,ndsurf,npass,nmask,isdsfile,maxpas,kdim,kdim2, * maxprobe,itty,iprint,msk,istanprofile,ju) C UPDATESDSTATS UPDATES ACCOUNTING STATISTICS C FOR PROGRAM c**************************************************** c c Passed Variables: c O itty - number of profiles with outliers (per probe/pass) c O nalt - number of altered profiles (per probe/pass) c altered profile is profile with standard deviation c outliers masked out, but still used c O nbad - number of profiles not used due to outliers (per probe/pass) c O ndep - 1=frequency of error (per depth/pass) c 2=total number of measurements (per depth/pass) c O ndeptot - 1=total number of errors (per pass) c 2=total number of measurements (per pass) c O nfreq - number of profiles by the number of outliers c contained (per pass) c I maxpas - maximum number of passes through standard deviation check c I maxfreq - maximum number of outliers per profile listed in statistics c I maxprobe - maximum number of data sets c I kdim - number of depth levels examined in standard deviation check c I iax - 1=check annual statistics below certain level,0=do not c I ity - number of outliers in profile c O ialready - if 1, record profile specific statistics, c if > 1, statistics are already recorded c I nx - sequential dataset order c I mpass - present pass through standard deviation check c I npro - number of datasets requested c I numo - number of standard deviation outliers which invalidate an c entire profile from use c I isurf - if 1, dataset is surface-only c O ndepcheck - set to one to discard standard deviation outliers c I itti - set to one when first value of variable is found c I namis - number of individual profiles to write to stanpros.d file c IO nam - number of individual profiles presently written to c stanpros.d file c I numlevels - number of levels with measured data in profile c I isoor - set to 1 for standard levels (0=observed levels) c I igeot - 1=calculate geopotential thickness c -1=perform calculations on density surfaces c 0=normal c I kdim1 - first level to perform calculations c I kdimx - deepest level to perform calculations c I isea - time period (00=annual, 13-16=seasons, 1-12=months) c I m - present dataset code c I ik - measured variable code c I rlon,rlat,iyear,month,iday,cc,icruise - information for station c I inume - levels with standard deviation outliers for profile c I lon,lat - longitude/latitude grid numbers for 1x1 gridbox c I lon2,lat2 - longitude/latitude grid numbers for 5x5 gridbox c I bmiss - missing value marker c I season - character representation of time period c I jj - sequential station order c I ndsurf - index for desired density surface c I npass - total number of passes through standard deviation check c I nmask - number of each type of mask c I isdsfile - file ID number for sdout[].s mask c I kdim2 - total number of depth levels c I iprint - set to one to print out stations with outliers c I msk = 0 - simply calculate means c = 1 - calculate means after throwing out c standard deviation outliers c = 2 - calculate means after checking for c standard deviation outliers against c current five-degree statistics c = 3 - same as (2), but standard deviation c outliers masks are automatically created c = 4 - same as (1), but standard deviation c outliers masks are automatically created c I istanprofile - FORTRAN file ID number for stanpros file c I ju - unique station number c c**************************************************************** parameter (maxprobel=100,maxcalc=200,kdiml=102) parameter (idim5=72,jdim5=36,kdim2l=137) parameter (maxlevel=6000) character*2 cc character*9 season character ctype(maxprobel)*11 character citf2(maxlevel)*6 common /parname4/ citf2 common /probes/ ctype dimension nfreq(0:maxfreq+1,maxpas) dimension ndep(2,kdim,maxpas),inume(kdim) dimension itty(maxprobe,maxpas) dimension ndeptot(2,maxpas) dimension dz(kdim2l),nbad(maxprobe,maxpas),nalt(maxprobe,maxpas) dimension nmask(5) dimension numlevels(0:1) double precision tempav5 double precision tempsq5 double precision temp common /sdata/ temp(kdim2l,maxcalc) common /multy/ mask(idim5,jdim5,kdiml) c common /levels/ dz common /levelsfull/ ndepsfull,dz common /five/ tempav5(idim5,jdim5,kdiml), * tempsq5(idim5,jdim5,kdiml), * num5(idim5,jdim5,kdiml) icountsd=0 c*********************************************************** c c If this is the correct time period, add to counters: c c itty - number of profiles with outliers c nalt - number of altered profiles c nbad - number of profiles not used due to outliers c ndep - frequency of error per depth c ndeptot - frequency of errors c nfreq - number of profiles with the number of outliers c contained in this profile c c********************************************************* if (iax.eq.0) then if (ity.gt.0) then ialready=ialready+1 if ( ialready .eq. 1 ) itty(nx,mpass)=itty(nx,mpass)+1 if ( ialready .eq. 1 ) itty(npro+1,mpass)=itty(npro+1,mpass)+1 endif if (ity.lt.numo .or. isurf .eq. 1) then if (ity.gt.0.and.ndepcheck.gt.0 .and. ialready .eq. 1) then nalt(nx,mpass)=nalt(nx,mpass)+1 nalt(npro+1,mpass)=nalt(npro+1,mpass)+1 endif elseif ( isurf .eq. 0 ) then nbad(nx,mpass)=nbad(nx,mpass)+1 nbad(npro+1,mpass)=nbad(npro+1,mpass)+1 endif do 454 k=1,ity if ( isurf .eq. 0 ) then ndep(1,inume(k),mpass)=ndep(1,inume(k),mpass)+1 else ndep(1,1,mpass)=ndep(1,1,mpass) endif ndeptot(1,mpass)=ndeptot(1,mpass)+1 454 continue if ( itti .eq. 1 ) then if (ity.le.maxfreq) then nfreq(ity,mpass)=nfreq(ity,mpass)+1 else nfreq(maxfreq+1,mpass)=nfreq(maxfreq+1,mpass)+1 endif endif c*********************************************************** c c write out first twenty-five outlier profiles c c*********************************************************** if (ity.gt.0.and.nam.lt.namis.and. * mpass.gt.1.and.iprint.gt.0 ) then if ( ialready .eq. 1 ) nam=nam+1 numls=numlevels(isoor) if ( isurf .eq. 1 ) numls=1 if ( igeot .eq. -1 ) numls=kdimx-kdim1+1 call printpro(ctype(m),isea,citf2(ik), * numls,rlon,rlat,ik,iyear, * month,iday,cc,icruise,dz,inume,ity,lon,lat, * lon2,lat2,bmiss,season,mpass,jj,igeot, * ndsurf,istanprofile) endif c*********************************************************** c c record in mask if profile was discarded (if mask is c present) c c writetomask: c c jj - record number c ic - counter of number of profiles written to file c ifn - file identification number for duplicate mask c*********************************************************** if ( mpass .eq. npass .and. * (msk .ge. 3 .or. nmask(4).gt.nmask(3)) * .and.ity.ge.numo ) then call writetomask(ju,icountsd,isdsfile) endif endif return end FUNCTION USEPRESTATS(grids,ik,nomal,iylook) C DECISION ON WHETHER TO USE PRE STATS C PRE-STATS ARE THE STATISTICS (5-GRID BOX) FROM C THE STANDARD DEVIATION CHECK DIRECTLY BEFORE C THE LAST STANDARD DEVIATION CHECK. c I ik - code for requested measured variable c I grids - grid size of statistical check c I nomal - set to one if anomalies are to be calculated, zero c for full variable values c I iylook -look at only certain years (2), months (1), or all values character*80 filename maxparm=100 useprestats=0 c make sure this is 1 degree grid if ( grids .lt. 1.00 - .1 .or. grids .gt. 1.00 + .1 ) return c make sure this is not anomaly if ( nomal .eq. 1 ) return c make sure this is not year specific if ( iylook .ne. 0 ) return c make sure this is a qc'ed variable call clearstring(filename,80) call extraname('sys.inf'//CHAR(0),'infinalsd.i'//CHAR(0), * filename) open(10,file=filename,status='unknown') do 50 n=1,maxparm read(10,*,END=4) iparm if ( iparm .eq. 7 ) goto 50 if ( iparm .eq. ik ) then useprestats=1 goto 4 endif 50 continue 4 continue close(10) return end SUBROUTINE USERSD( nsuct, isea, nomal, nomaly, * iylook, nsucy, numsuc, nyc, nyc2, istc, grids, msk, numo, * ndepcheck, cyon, cyon2, uninfile, unoutfile, statpath, * anlypath, iprint, namis, itprint, uninparm, kdimx, npass, * ionegrid, latgrid, longrid, kdim1, iannual, kdima, nsub, * igeot,ibinsurf,jlev,dsurf ) C USERSD READS IN USER INPUT INFORMATION FROM DEVFILE.D c********************************************************************* c c Passed Variables: c c nsuct - set to zero for one time period, one for all time periods c requested in devshell c isea - time period of statistical checks: c 00 - annual c 01-12 - months c 13-16 - seasons c 17-40 - half months c nomal - set to one if anomalies are to be calculated, zero c for full parameter values c nomaly - read which type of anomaly: annual - 0, seasonal - 1, c monthly - 2 c iylook -look at only certain years (2), months (1), or all values c nsucy - set to one if successive years are desired (using c devshell2) c nyc - starting year for statistics checks c nyc2 - ending year of statistics checks c istc - information on which statistics to save: c 0. one grid box means c 1. one grid box data distributions c 2. one grid box standard deviations c 3. five grid box means c 4. five grid box data distributions c 5. five grid box standard deviations c each array member is set to one if the user wants the c statistic saved, zero otherwise. c grids - grid size of statistical check c msk - set to two for a statistical check against existing five c grid statistics, one for normal statistical checks, and c zero to do simple statistical compilations c numo - number of standard deviation outliers which invalidate a c profile c ndepcheck - set to one to discard standard deviation outliers c in useable profiles c cyon - set to 'n' if a non-standard input file name is used c set to 'y' otherwise c cyon2 - set to 'n' if a non-standard output file name is used c set to 'y' otherwise c unoutfile - non-standard output file name c uninfile - non-standard input file name c statpath - standard statistics output path c anlypath - path for analyzed field c iprint - set to one if the user requests individual profile c printouts, zero otherwise c namis - the number of individual profiles requested c itprint - set to one if statistics totals are to be printed out c uninparm - name of non-standard input file c kdimx - lowest level to perform checks c npass - number of passes through standard deviation check. c One pass is all that is completed for simple c calculations, npass + 1 passes for deviation checks. c ionegrid - set to one if only one five grid box is to be checked, c zero otherwise c latgrid - latitude grid box for one grid check c longrid - longitude grid box for one longitude check c kdim1 - first level to perform checks on c iannual - set to one if annual statistics are to be used below c a certain depth level c kdima - depth level at which annual statistics are to be used c minstan - minimum number of standard deviations c nsub - which fraction of year desired (used only if numsuc < 0) c igeot - set to one if geopotential thickness is to be c calculated c set to 2 to calculate stats for mixed layer depth c set to 3 to calculate a gradient c set to -1 to calculate variable on given density surface c ibinsurf - surface data binning:-1=no date or geographic binning c 0=no date binning,gridbox c 1=bin by day,gridbox c 2=bin by month,gridbox c 3=bin by season,gridbox c jlev - reference level for density surface c dsurf - desired density surface c c******************************************************************* parameter (icmax=10) common /stanmin/ minstan c******************************************************************* c c Character Arrays: c statpath - standard statistics output path c uninparm - name of non-standard input file c anlypath - path for analyzed field c unoutfile - non-standard output file name c uninfile - non-standard input file name c cyon - set to 'n' if non-standard input file name c cyon2 - set to 'n' if non-standard output file name c text - dummy variable c c********************************************************************** character*80 unoutfile, uninfile, statpath character*80 anlypath, uninparm, text character*1 cyon,cyon2 dimension latgrid(2),longrid(2) dimension icomma(icmax,2) call initsdvars( nsuct, isea, nomal, nomaly, * iylook, nsucy, numsuc, nyc, nyc2, istc, grids, msk, numo, * ndepcheck, cyon, cyon2, uninfile, unoutfile, statpath, * anlypath, iprint, namis, itprint, uninparm, kdimx, npass, * ionegrid, latgrid, longrid, kdim1, iannual, kdima, nsub, * igeot,ibinsurf,jlev,dsurf ) idebug = 0 !- 0 = off 1 = on ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c c READ IN INSTRUCTIONS c c c ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c************************************************************** c c Open devfile.d, the instruction file c c************************************************************** open(11,file='./devfile.d',status='old') read(11,'(a1)') text read(11,'(a1)') text read(11,'(a1)') text read(11,'(a1)') text c*************************************************************** c c read in whether to do only one time period or all time periods c 00 through 16 c c*************************************************************** if (idebug .gt. 0) print *,'nsuct' read(11,*) nsuct read(11,'(a1)') text read(11,'(a1)') text read(11,'(a1)') text c*************************************************************** c c read in time period c c************************************************************** if (idebug .gt. 0) print *,'isea' read(11,*) isea read(11,'(a1)') text c************************************************************** c c Annual statistics at low depths (1) c c************************************************************** if (idebug .gt. 0) print *,'iannual' read(11,*) iannual read(11,'(a1)') text c************************************************************** c c Depth level to start using Annual statitstics c c************************************************************** if (idebug .gt. 0) print *,'kdima' read(11,*) kdima read(11,'(a1)') text read(11,'(a1)') text c************************************************************** c c read in whether to do anomalies(1) or full values(0) c c************************************************************** if (idebug .gt. 0) print *,'nomal' read(11,*) nomal read(11,'(a1)') text c*************************************************************** c c read which type of anomaly: annual - 0, seasonal - 1, monthly - 2 c c*************************************************************** if (idebug .gt. 0) print *,'nomaly' read(11,*) nomaly read(11,'(a1)') text read(11,'(a1)') text c*************************************************************** c c Calculate geopotential thickness (1) c Calculate statistics on density surface (-1) c for density surface add two more entries on line c - reference surface (-1=insitu, 0,500,1000, etc..) c - desired density surface (ex. 31.85) c c ex: -1 1000 31.85 c calculate variable on density surface 31.85 relative c to 1000 decibars c c*************************************************************** if (idebug .gt. 0) print *,'igeot' ierr=1 read(11,*) igeot if ( igeot .eq. -1 ) then read(11,*) jlev,dsurf endif read(11,'(a1)') text read(11,'(a1)') text c*************************************************************** c c look at only certain years (2), months (1), or all values c c*************************************************************** if (idebug .gt. 0) print *,'iylook' read(11,*) iylook read(11,'(a1)') text c*************************************************************** c c do successive years (2), months (1), or all values c c*************************************************************** if (idebug .gt. 0) print *,'nsucy' read(11,*) nsucy read(11,'(a1)') text c*************************************************************** c c Number of successive years (months) to composite c If number of successive years is a fraction of a c year, make numsuc negative c c*************************************************************** if (idebug .gt. 0) print *,'xnumsuc' read(11,*) xnumsuc if ( xnumsuc .lt. 1. .and. xnumsuc .gt. 0) then numsuc = -12/( 1./xnumsuc) else numsuc = xnumsuc endif read(11,'(a1)') text c*************************************************************** c c read in starting and ending years to be looked at c if fraction of year is desired (numsuc < 0 ) then c read nsub, which part of year is desired (1=first fraction) c c*************************************************************** if (idebug .gt. 0) print *,'nyc,nyc2' if ( numsuc .lt. 0. ) then read(11,*) nyc,nyc2, nsub else read(11,*) nyc,nyc2 endif read(11,'(a1)') text read(11,'(a1)') text c**************************************************************** c c first level to look at c c**************************************************************** if (idebug .gt. 0) print *,'kdim1' read(11,*) kdim1 read(11,'(a1)') text c**************************************************************** c c lowest level to look at c c**************************************************************** if (idebug .gt. 0) print *,'kdimx' read(11,*) kdimx read(11,'(a1)') text read(11,'(a1)') text c***************************************************************** c c one gridbox (1) or entire world (0) c c***************************************************************** if (idebug .gt. 0) print *,'ionegrid' read(11,*) ionegrid read(11,'(a1)') text c****************************************************************** c c latitude and longitude gridboxes (5X5) c c****************************************************************** if (idebug .gt. 0) print *,'latgrid,longrid' call clearstring(text,80) read(11,'(a)') text call findlastchar(nsp,text,80) call spacefindl(nsp,icmax,ncomma,icomma,text) if ( ncomma. ge. 1 ) * call readintfromchar(latgrid(1),icomma(1,2), * icomma(1,1),icomma(1,1)+icomma(1,2)-1,text) if ( ncomma .ge. 2 ) * call readintfromchar(longrid(1),icomma(2,2), * icomma(2,1),icomma(2,1)+icomma(2,2)-1,text) if ( ncomma .ge. 3 ) * call readintfromchar(latgrid(2),icomma(3,2), * icomma(3,1),icomma(3,1)+icomma(3,2)-1,text) if ( ncomma .ge. 4 ) * call readintfromchar(longrid(2),icomma(4,2), * icomma(4,1),icomma(4,1)+icomma(4,2)-1,text) if ( ionegrid .gt. 0 .and. ncomma .gt. 2 ) ionegrid=2 read(11,'(a1)') text read(11,'(a1)') text c****************************************************************** c c number of passes through deviation check c c****************************************************************** if (idebug .gt. 0) print *,'npass' read(11,*) npass read(11,'(a1)') text read(11,'(a1)') text read(11,'(a1)') text read(11,'(a1)') text read(11,'(a1)') text read(11,'(a1)') text read(11,'(a1)') text read(11,'(a1)') text read(11,'(a1)') text read(11,'(a1)') text read(11,'(a1)') text c************************************************************** c c save statistics or not c c************************************************************** if (idebug .gt. 0) print *,'istc' read(11,*) istc read(11,'(a1)') text read(11,'(a1)') text c****************************************************************** c c read grid spacing c c****************************************************************** if (idebug .gt. 0) print *,'grids' read(11,*) grids read(11,'(a1)') text read(11,'(a1)') text c****************************************************************** c c do checks on 5X5 statistics (2), c discard deviation outliers (1) or simply calculate means (0) c c****************************************************************** if (idebug .gt. 0) print *,'msk' read(11,*) msk read(11,'(a1)') text c****************************************************************** c c number of outliers which invalidate a profile c c****************************************************************** if (idebug .gt. 0) print *,'numo' read(11,*) numo read(11,'(a1)') text c****************************************************************** c c discard outliers in useable profiles (1) or not (0) c c****************************************************************** if (idebug .gt. 0) print *,'ndepcheck' read(11,*) ndepcheck read(11,'(a1)') text read(11,'(a1)') text c****************************************************************** c c minimum number of standard deviatons c c****************************************************************** if (idebug .gt. 0) print *,'minstan' read(11,*) minstan open(12,file='./dumfin',status='unknown') if ( minstan .lt. 0 ) then read(12,*) minstan else write(12,*) minstan endif close(12) read(11,'(a1)') text read(11,'(a1)') text c******************************************************************* c c standard input file (y) or not (n) c c******************************************************************* if (idebug .gt. 0) print *,'cyon' read(11,'(a1)') cyon read(11,'(a1)') text c******************************************************************* c c name of non-standard input file c c******************************************************************* if (idebug .gt. 0) print *,'uninfile' read(11,*) uninfile read(11,*) uninparm read(11,'(a1)') text c******************************************************************* c c standard output file (y) or not (n) c c******************************************************************* if (idebug .gt. 0) print *,'cyon2' read(11,'(a1)') cyon2 read(11,'(a1)') text c******************************************************************* c c name of non-standard output file c c******************************************************************* if (idebug .gt. 0) print *,'unoutfile' read(11,*) unoutfile read(11,'(a1)') text read(11,'(a1)') text c******************************************************************* c c standard statistics output pathname c c******************************************************************* if (idebug .gt. 0) print *,'statpath' read(11,*) statpath read(11,'(a1)') text c******************************************************************* c c standard analyzed data pathname c c******************************************************************* if (idebug .gt. 0) print *,'anlypath' read(11,*) anlypath c*************************************************************** c c Surface binning c c*************************************************************** read(11,'(a1)') text read(11,'(a1)') text read(11,*) ibinsurf c cant do daily yet if ( ibinsurf .eq. 1 ) ibinsurf=0 read(11,'(a1)') text read(11,'(a1)') text read(11,'(a1)') text read(11,'(a1)') text c******************************************************************** c c make individual profile printout (1) or not (0) c c******************************************************************** if (idebug .gt. 0) print *,'iprint' read(11,*) iprint read(11,'(a1)') text c******************************************************************** c c number of individual profiles to print c c******************************************************************** if (idebug .gt. 0) print *,'namis' read(11,*) namis read(11,'(a1)') text c******************************************************************** c c print out totals (1) or not (0) c c******************************************************************** if (idebug .gt. 0) print *,'itprint' read(11,'(i1)') itprint close(11) c convert lat/lon if ( ionegrid .eq. 2 ) then rlat1=latgrid(1) rlon1=longrid(1) rlat2=latgrid(2) rlon2=longrid(2) call grid(rlat1,rlon1,latgrid(1),longrid(1),grids) call grid(rlat2,rlon2,latgrid(2),longrid(2),grids) write(6,*) 'rlat',latgrid,longrid endif return end SUBROUTINE WRITEOUTSDSTATS(igeot,kdim1,ndsurf,kdimx,istc, * idimgs,jdimgs,idimge,jdimge,idimt,jdimt,idim5gs,jdim5gs, * idim5ge,jdim5ge,idimt5,jdimt5,npre) C WRITEOUTSDSTATS WRITES OUT FINAL STATISTICS c**************************************************** c c Passed Variables: c I igeot - 1=calculate geopotential thickness c -1=perform calculations on density surfaces c 0=normal c I kdim1 - first level to perform calculations c I kdimx - deepest level to perform calculations c I istc - information on which statistics to save: c 0. one grid box means c 1. one grid box data distributions c 2. one grid box standard deviations c 3. five grid box means c 4. five grid box data distributions c 5. five grid box standard deviations c each array member is set to file ID number if c statistic are to be saved, zero otherwise. c I idimt,jdimt - ngrid*idim, ngrid*jdim c I idimt5,jdimt5 - ngrid*idim5, ngrid*jdimt5 c I idimgs - starting longitude to perform checks on c I idimge - ending longitude to perform checks on c I jdimgs - starting latitude to perform checks on c I jdimge - ending latitude to perform checks on c I idim5gs - starting five grid longitude to perform checks on c I idim5gs - ending five grid longitude to perform checks on c I jdim5gs - starting five grid latitude to perform checks on c I jdim5ge - ending five grid latitude to perform checks on c I npre - if 0 - final statistics c if 1 - pre stats (statistics from next to last run c through standard deviation check) c c************************************************************ parameter (jdim=180, idim=360, kdim=102) parameter (idim5=72,jdim5=36) dimension istc(0:5,0:1) dimension num5(idim5,jdim5,kdim) dimension num(idim,jdim,kdim) double precision tempsq5(idim5,jdim5,kdim) double precision tempsq(idim,jdim,kdim) double precision tempav(idim,jdim,kdim) double precision tempav5(idim5,jdim5,kdim) common /mean/ tempav common /squared/ tempsq common /dist/ num common/five/ tempav5,tempsq5,num5 c************************************************************* c c Write out data distributions if desired c c************************************************************* if ( igeot .ne. -1 ) then kd1=kdim1 else kd1=kdim1-ndsurf+1 endif if ( istc(0,npre) .gt. 0 ) then call writeleveli(1,1,kd1,idimgs,jdimgs,idimge,jdimge, * idimt,jdimt,idim,jdim,kdim1,kdimx,istc(0,npre),num) endif if ( istc(3,npre) .gt. 0 ) then call writeleveli(1,1,kd1,idim5gs,jdim5gs,idim5ge,jdim5ge, * idimt5,jdimt5,idim5,jdim5,kdim1,kdimx,istc(3,npre),num5) endif do 550 m7 = 1,5 c******************************************* c c Put proper double precision information c into single precision variable: c if m7=1 one grid means c if m7=2 one grid standard deviations c if m7=4 five grid means c if m7=5 five grid standard deviations c c******************************************* if ( m7.ne.3.and. istc(m7,npre).gt.0 ) then c********************************************* c c Save statistics to correct file c c********************************************* if ( m7 .eq. 1 ) then call writeleveld(1,1,kd1,idimgs,jdimgs,idimge,jdimge,idimt, * jdimt,idim,jdim,kdim1,kdimx,istc(m7,npre),tempav) write(6,*) 'm7',kd1,kdim1,kdimx,idimgs,jdimgs, * idimge,jdimge,idimt,jdimt,idim,jdim, * tempav(255,65,1) elseif ( m7 .eq. 2 ) then call writeleveld(1,1,kd1,idimgs,jdimgs,idimge,jdimge,idimt, * jdimt,idim,jdim,kdim1,kdimx,istc(m7,npre),tempsq) elseif ( m7 .eq. 4 ) then call writeleveld(1,1,kd1,idim5gs,jdim5gs,idim5ge,jdim5ge, * idimt5,jdimt5,idim5,jdim5,kdim1,kdimx,istc(m7,npre),tempav5) elseif ( m7 .eq. 5 ) then call writeleveld(1,1,kd1,idim5gs,jdim5gs,idim5ge,jdim5ge, * idimt5,jdimt5,idim5,jdim5,kdim1,kdimx,istc(m7,npre),tempsq5) endif endif 550 continue return end