Blame | Last modification | View Log | Download | RSS feed
PROGRAM differencec **************************************************************************c * Get the difference between two runs *c * Michael Sprenger / Autumn 2006 *c **************************************************************************implicit nonec --------------------------------------------------------------------------c Declaration of subroutine parametersc --------------------------------------------------------------------------c Physical and numerical parametersreal epsparameter (eps=0.01)integer nmaxparameter (nmax=300*300*200)c Variables for input file 1character*80 i1_filenamereal i1_varmin(4),i1_varmax(4),i1_stag(4)integer i1_vardim(4)real i1_mdvinteger i1_ndiminteger i1_nx,i1_ny,i1_nzreal i1_timeinteger i1_nvarscharacter*80 i1_vnam(100)real i1_field(nmax)c Variables for input file 2character*80 i2_filenamereal i2_varmin(4),i2_varmax(4),i2_stag(4)integer i2_vardim(4)real i2_mdvinteger i2_ndiminteger i2_nx,i2_ny,i2_nzreal i2_timeinteger i2_nvarscharacter*80 i2_vnam(100)real i2_field(nmax)c Variables for output filecharacter*80 o_filenamereal o_varmin(4),o_varmax(4),o_stag(4)integer o_vardim(4)real o_mdvinteger o_ndiminteger o_nx,o_ny,o_nzreal o_timereal o_field(nmax)c Auxiliary variablesinteger i,j,kinteger cdfid,ierrinteger isokcharacter*80 cfn,varnamec --------------------------------------------------------------------------------c Inputc --------------------------------------------------------------------------------print*,'********************************************************'print*,'* DIFFERENCE *'print*,'********************************************************'c Read parameter fileopen(10,file='fort.10')read(10,*) i1_filenameread(10,*) i2_filenameread(10,*) o_filenameclose(10)print*print*,trim(i1_filename)print*,trim(i2_filename)print*,trim(o_filename)print*c Get a list of all variables on the two input filescall cdfopn(i1_filename,cdfid,ierr)if (ierr.ne.0) goto 998call getvars(cdfid,i1_nvars,i1_vnam,ierr)if (ierr.ne.0) goto 998call getdef(cdfid,i1_vnam(i1_nvars),i1_ndim,i1_mdv,i1_vardim,> i1_varmin,i1_varmax,i1_stag,ierr)call clscdf(cdfid,ierr)print*,(trim(i1_vnam(i))//' ',i=1,i1_nvars)call cdfopn(i2_filename,cdfid,ierr)if (ierr.ne.0) goto 997call getvars(cdfid,i2_nvars,i2_vnam,ierr)if (ierr.ne.0) goto 997call clscdf(cdfid,ierr)print*,(trim(i2_vnam(i))//' ',i=1,i2_nvars)print*c Create the output fileo_ndim=i1_ndimdo i=1,i1_ndimo_varmin(i)=i1_varmin(i)o_varmax(i)=i1_varmax(i)enddocfn=trim(o_filename)//'_cst'call crecdf(trim(o_filename),cdfid,o_varmin,o_varmax,> o_ndim,trim(cfn),ierr)if (ierr.ne.0) goto 996call clscdf(cdfid,ierr)c --------------------------------------------------------------------------------c Loop through all variablesc --------------------------------------------------------------------------------do i=1,i1_nvarsc Check wether the variable is available on both filesvarname=i1_vnam(i)if (varname.eq.'time') goto 100isok=0call check_varok (isok,varname,i2_vnam,i2_nvars)if (isok.eq.0) goto 100c Read first input filecall cdfopn(i1_filename,cdfid,ierr)if (ierr.ne.0) goto 998call getdef(cdfid,varname,i1_ndim,i1_mdv,i1_vardim,> i1_varmin,i1_varmax,i1_stag,ierr)if (ierr.ne.0) goto 998call getdat(cdfid,varname,0.,0,i1_field,ierr)if (ierr.ne.0) goto 998call clscdf(cdfid,ierr)c Read second input filecall cdfopn(i2_filename,cdfid,ierr)if (ierr.ne.0) goto 998call getdef(cdfid,varname,i2_ndim,i2_mdv,i2_vardim,> i2_varmin,i2_varmax,i2_stag,ierr)if (ierr.ne.0) goto 998call getdat(cdfid,varname,0.,0,i2_field,ierr)if (ierr.ne.0) goto 998call clscdf(cdfid,ierr)c Consistency checkif (i1_ndim.ne.i2_ndim) thenprint*,'Inconsistent input files... Stop',i1_ndim,i2_ndimstopendifdo j=1,3if ( (i1_vardim(j).ne.i2_vardim(j)).or.> (abs(i1_varmin(j)-i2_varmin(j)).gt.eps).or.> (abs(i1_varmax(j)-i2_varmax(j)).gt.eps).or.> (abs(i1_stag(j)-i2_stag(j)).gt.eps)) thenprint*,'Inconsistent input files... Stop'print*,j,i1_varmin(j),i2_varmin(j)print*,j,i1_varmax(j),i2_varmax(j)print*,j,i1_vardim(j),i2_vardim(j)print*,j,i1_stag(j), i2_stag(j)stopendifenddoc Get the differencedo j=1,i1_vardim(1)*i1_vardim(2)*i1_vardim(3)if ( (abs(i1_field(j)-i1_mdv).gt.eps).and.> (abs(i2_field(j)-i2_mdv).gt.eps) ) theno_field(j)=i1_field(j)-i2_field(j)elseo_field(j)=i1_mdvendifenddoc Write to output fileo_ndim=i1_ndimo_mdv =i1_mdvdo j=1,i1_ndimo_vardim(j)=i1_vardim(j)o_varmin(j)=i1_varmin(j)o_varmax(j)=i1_varmax(j)o_stag(j) =i1_stag(j)enddocall cdfwopn(o_filename,cdfid,ierr)if (ierr.ne.0) goto 996call putdef(cdfid,varname,o_ndim,o_mdv,o_vardim,> o_varmin,o_varmax,o_stag,ierr)if (ierr.ne.0) goto 996call putdat(cdfid,varname,0.,0,o_field,ierr)if (ierr.ne.0) goto 996print*,'W ',trim(varname),' ',trim(o_filename)call clscdf(cdfid,ierr)if (ierr.ne.0) goto 996c Next100 continueenddoc -----------------------------------------------------------------c Exception handling and format specsc -----------------------------------------------------------------stop998 print*,'Problems with input file 1 ',trim(i1_filename)stop997 print*,'Problems with input file 2 ',trim(i2_filename)stop996 print*,'Problems with output file ',trim(o_filename)stopendc --------------------------------------------------------------------------------c Check whether variable is found on netcdf filec --------------------------------------------------------------------------------subroutine check_varok (isok,varname,varlist,nvars)c Check whether the variable <varname> is in the list <varlist(nvars)>. If this isC the case, <isok> is incremented by 1. Otherwise <isok> keeps its value.implicit nonec Declaraion of subroutine parametersinteger isokinteger nvarscharacter*80 varnamecharacter*80 varlist(nvars)c Auxiliary variablesinteger ic Maindo i=1,nvarsif (trim(varname).eq.trim(varlist(i))) isok=isok+1enddoend