module mod_hycomfile_io ! Module for reading several types of files produced by hycom, ! including restart, archv, nersc_daily and nersc_weekly. ! ! Types: ! hycomfile - contains info on contents of the hycomfile ! ! Subroutines : ! --subroutine initHF(hycomfile, file, filetype, [template]) ! This normally is one of the first steps for setting up the ! HYCOM I/O module. It initializes the hycomfile type by parsing ! the contents of filename, given a hint on what file type this is. ! ! --subroutine HFReadField(hycomfile,field,idm,jdm,cfld,coord,tlevel) ! Use this routine after calling initHF. It reads a 2D field "field" ! with horizontal dimensions (idm,jdm), having the name "cfld", ! at vertical coordinate "coord" and time level "tlevel" ! ! --subroutine HFReadField3D(df,field,idm,jdm,kdm,cfld,tlevel) ! Use this routine after calling initHF. It reads a 3D field "field" ! with 3D dimensions (idm,jdm,kdm), having the name "cfld", ! at vertical coordinate "coord" and time level "tlevel". It ! uses the routine HFReadField ! ! --subroutine HFReadDPField(df,dp,idm,jdm,coord,tlevel) ! Use this routine after calling initHF. Since variabe names are ! different across hycom files custom routines are available for some ! fields. This routine will read the layer thickness variable ! ("dp" in restart files, "thknss" in archive files, etc) of vertical ! coordinaye "coord" at time level "tlevel". ! ! --subroutine HFReaduvtot(df,ut,vt,idm,jdm,vlevel,tlevel) ! Use this routine after calling initHF. Since variabe names are ! different across hycom files custom routines are available for some ! fields. This routine will read the total velocity of layer ! coordinaye "coord" at time level "tlevel". ! ! --subroutine HFReaduvbaro(df,ub,vb,idm,jdm,tlevel) ! Use this routine after calling initHF. Since variabe names are ! different across hycom files custom routines are available for some ! fields. This routine will read the barotropic velocity at ! time level "tlevel". ! ! ! ! Usage example: ! ------------------------------------------------------------------- ! use mod_hycomfile_io ! ... ! type(hycomfile) :: hfile ! real, dimension(:,:), allocatable :: field ! real, dimension(:,:,:), allocatable :: field3D ! ! ! Retrieve grid size first - this uses the other modules in libhycnersc.a ! call xcspmd() ! sets idm, jdm ! call zaiost() ! call get_grid() ! ! ! Reading a 2D field ! fname='TP3restart1990_250_00.a' ! A test file ! ftype=getfiletype(filename) ! 'restart' in this case ! call initHF(hfile, file, filetype) ! kdm=vDim(hfile) ! Retrieves vertical dimension in file ! ! ! Allocate variables to be read ! allocate(field3D(idm,jdm,kdm)) ! allocate(field(idm,jdm)) ! ! ! Read mixed layer thickness from that file ! call HFReadField(hfile,field,idm,jdm'dpmixl ',0,1) ! ! ! Read salinity in layer 1 from the same file ! call HFReadField(hfile,field,idm,jdm'saln ',1,1) ! ! ! Read temperature in all layers from the same file ! call HFReadField(hfile,field3D,idm,jdm,kdm,'temp ',1) ! ! ------------------------------------------------------------------- ! TODO: More info ! TODO: Loosen restriction on character length when calling routines use mod_parameters implicit none type hycomfile integer :: iversn = 0 integer :: iexpt = 0 integer :: nstep = 0 integer :: yrflag = 0 real*8 :: dtime = 0 real*8 :: fyear = 0 integer :: count = 0 integer :: nrec = 0 character(len=80) :: filebase ='', ftype='' character(len=8), pointer :: cfld (:) integer , pointer, dimension(:) :: coord, tlevel ! Time info integer :: iyear = 0 integer :: imonth= 0 integer :: iweek = 0 integer :: iday = 0 ! day starting at 0 in year integer :: ihour = 0 ! day starting at 0 in year real :: ftime ! floating point year integer :: start_iyear integer :: start_iday end type logical, parameter, private :: silent=.true. ! TODO: Document ! TODO: Make sure DP routine returns in pressure coords always private :: readHeader, readVarHeaders, readFieldEntry, writeFieldEntry, & HFwriteHeader contains subroutine initHF(df,filename,ftype,template) implicit none type(hycomfile) , intent(out) :: df character(len=*), intent( in) :: filename,ftype type(hycomfile) , intent( in),optional :: template integer :: ind if (present(template)) then df=template ind=max(index(filename,'.a'),index(filename,'.b')) df%filebase=filename(1:ind-1) df%ftype=ftype return end if ! Get file base ind=max(index(filename,'.a'),index(filename,'.b')) df%filebase=filename(1:ind-1) df%ftype=ftype ! Read header call ReadHeader(df) ! Get var info call readVarHeaders(df) end subroutine !!!!!!!!!!! Header processing - parameters !!!!!!!!!!!!! subroutine readHeader(df) use mod_year_info implicit none type(hycomfile), intent(inout) :: df character(len=80) :: c80 character(len=80) :: ctitle(4) integer :: nop,k,i, ind integer :: iversn, iexpt, yrflag,lidm,ljdm,lkdm, dmonth integer :: syear, sday, dyear, dday, itmp, ios(3), indx nop=777 !TODO : check date setup df%start_iyear=0 df%start_iday =0 open(nop,file=trim(df%filebase)//'.b',status='old') if (trim(df%ftype)=='restart') then read(nop,'(a)') c80 ; ind=index(c80,'='); read(c80(ind+1:),*) df%iexpt, df%iversn, df%yrflag read(nop,'(a)') c80 ; ind=index(c80,'='); read(c80(ind+1:),*) df%nstep, df%dtime if (df%yrflag==3) then call juliantodate(floor(df%dtime),df%iyear,itmp,itmp,1901,1,1) df%iday=floor(df%dtime) - datetojulian(df%iyear,1,1,1901,1,1) ! ordinal day df%ihour=nint((df%dtime - floor(df%dtime))*24.) else ! Due to a timing bug file name is our safest bet is to parse file ! names read(df%filebase(11:14),'(i4.4)',iostat=ios(1)) df%iyear read(df%filebase(16:18),'(i3.3)',iostat=ios(2)) df%iday read(df%filebase(20:21),'(i3.3)',iostat=ios(3)) df%ihour if (any(ios/=0)) then ! Drop it if that failed df%iyear=0 df%iday=0 df%ihour=0 end if end if df%fyear = df%iyear + min((df%iday + df%ihour)/365.,1.) else if (trim(df%ftype)=='nersc_daily') then read(nop,116) ctitle,df%iversn,df%iexpt,df%yrflag, & lidm,ljdm,lkdm,df%start_iyear,df%start_iday , & df%iyear,df%iday,df%count df%ihour=12 ! hardcoded for daily average ! "Floating point year" df%fyear = df%iyear + min((df%iday + df%ihour)/365.,1.) !print * , df%iyear, df%iday, df%ihour, df%fyear else if (trim(df%ftype)=='nersc_weekly') then read(nop,216) ctitle,df%iversn,df%iexpt,df%yrflag, & lidm,ljdm,lkdm,dmonth,dyear , & df%count read(df%filebase(5 :8),'(i4.4)',iostat=ios(1)) df%iyear read(df%filebase(10:12),'(i3.3)',iostat=ios(2)) df%imonth read(df%filebase(14:15),'(i3.3)',iostat=ios(3)) df%iweek if (any(ios/=0)) then ! Drop it if that failed df%iyear=0 df%iday=0 df%ihour=0 end if ! Floating point years - only useful as indicator df%fyear = df%iyear + (df%imonth -1) / 12. + (df%iweek -1. ) /(4.*12.) + 1./96. else if (trim(df%ftype)=='archv') then read(nop,316) ctitle,df%iversn,df%iexpt,df%yrflag, & lidm,ljdm ! Time from file name indx=index(df%filebase,'archv.') read(df%filebase(indx+6 :indx+9 ),'(i4)') df%iyear read(df%filebase(indx+11:indx+13),'(i3)') df%iday read(df%filebase(indx+15:indx+16),'(i2)') df%ihour df%fyear = df%iyear + min((df%iday + df%ihour)/365.,1.) else print *,'(Unknown file type : '//trim(df%ftype)//')' stop end if close (nop) ! nersc_daily format 116 format (a80/a80/a80/a80/ & i5,4x,'''iversn'' = hycom version number x10'/ & i5,4x,'''iexpt '' = experiment number x10'/ & i5,4x,'''yrflag'' = days in year flag'/ & i5,4x,'''idm '' = longitudinal array size'/ & i5,4x,'''jdm '' = latitudinal array size'/ & i5,4x,'''kdm '' = Vertical array size'/ & i5,4x,'''syear '' = Year of integration start '/ & i5,4x,'''sday '' = Day of integration start'/ & i5,4x,'''dyear '' = Year of this dump '/ & i5,4x,'''dday '' = Day of this dump '/ & i5,4x,'''count '' = Ensemble counter ') ! nersc_weekly format 216 format (a80/a80/a80/a80/ & i5,4x,'''iversn'' = hycom version number x10'/ & i5,4x,'''iexpt '' = experiment number x10'/ & i5,4x,'''yrflag'' = days in year flag'/ & i5,4x,'''idm '' = longitudinal array size'/ & i5,4x,'''jdm '' = latitudinal array size'/ & i5,4x,'''kdm '' = Vertical array size'/ & i5,4x,'''month '' = Month of this dump '/ & i5,4x,'''year '' = Year of this dump '/ & i5,4x,'''count '' = Averaging counter '/ & 'field time step model day', & ' k dens min max') ! Archive format 316 format (a80/a80/a80/a80/ & i5,4x,'''iversn'' = hycom version number x10'/ & i5,4x,'''iexpt '' = experiment number x10'/ & i5,4x,'''yrflag'' = days in year flag'/ & i5,4x,'''idm '' = longitudinal array size'/ & i5,4x,'''jdm '' = latitudinal array size'/ & 'field time step model day', & ' k dens min max') end subroutine readHeader subroutine skipHeader(ftype,nop) implicit none character(len=*), intent(in) :: ftype integer , intent(in) :: nop character(len=5) :: char5 integer :: ios if ( trim(ftype)=='nersc_weekly' .or. trim(ftype) == 'nersc_daily' & .or.trim(ftype)=='archv') then ios=0 ; char5='' do while (char5/='field' .and. ios==0) read(nop,'(a5)',iostat=ios) char5 end do else if (trim(ftype)=='restart' ) then read(nop,*) ; read(nop,*) else print *,'unknown file type '//trim(ftype) stop end if end subroutine !!!!!!!!!!! Header processing - fields !!!!!!!!!!!!! subroutine readVarHeaders(df) implicit none type(hycomfile), intent(inout) :: df character(len=5) :: char5 character(len=8) :: char8 integer :: ios, nop, nrec, coord, nstep, indx, irec,tlevel real :: xmin, xmax logical :: ex inquire(exist=ex,file=trim(df%filebase)//'.b') if (.not. ex) then print *,'file does not exist :' print *,trim(df%filebase)//'.b' stop '()' end if open(nop,file=trim(df%filebase)//'.b',status='old') ! Open input file call skipHeader(df%ftype,nop) ! Read until we get the index we want nrec=0 ; ios=0 do while(ios==0) call readFieldEntry(df%ftype,char8,coord,tlevel,xmin,xmax,nop,ios) nrec=nrec+1 !print *,nrec,char8,coord,tlevel,ios end do nrec=nrec-1; df%nrec=nrec rewind(nop) call skipHeader(df%ftype,nop) allocate(df%cfld (nrec)) allocate(df%coord (nrec)) allocate(df%tlevel (nrec)) ios=0 do irec=1,nrec call readFieldEntry(df%ftype,char8,coord,tlevel,xmin,xmax,nop,ios) !print *,irec,char8,coord,tlevel df%cfld (irec)=char8 df%coord (irec)=coord df%tlevel (irec)=tlevel end do close(nop) end subroutine !!!!!!!!!!! Header processing - one variable item !!!!!!!!!!!! ! Read one line of variable info and parse it subroutine readFieldEntry(ftype,cfld,coord,tlevel,xmin,xmax,nop,ios) implicit none character(len=*), intent(in) :: ftype character(len=8), intent(out) :: cfld integer , intent(out) :: coord,tlevel,ios real , intent(out) :: xmin, xmax integer , intent(in) :: nop integer :: nstep real :: rday,dens ! TODO: make sure this works properly - make it more robust if (trim(ftype)=='restart') then read(nop,4100,iostat=ios) cfld,coord,tlevel,xmin,xmax else if (trim(ftype)=="nersc_daily" .or. trim(ftype)=="nersc_weekly") then read(nop,117,iostat=ios) cfld,nstep,rday,coord,dens,xmin,xmax tlevel=1 else if (trim(ftype)=="archv") then read(nop,118,iostat=ios) cfld,nstep,rday,coord,dens,xmin,xmax tlevel=1 else print *,'unknown file type '//trim(ftype) stop end if 4100 format(a,': layer,tlevel,range = ',i3,i3,2x,1p2e16.7) 117 format (a8,' = ',i11,f11.2,i3,f7.3,1p2e16.7) 118 format (a8,' =',i11,f11.3,i3,f7.3,1p2e16.7) end subroutine ! Read one line of variable info and parse it subroutine writeFieldEntry(ftype,cfld,coord,tlevel,xmin,xmax,nop,ios) implicit none character(len=*), intent(in) :: ftype character(len=8), intent(in) :: cfld integer , intent(in) :: coord,tlevel real , intent(in) :: xmin, xmax integer , intent(in) :: nop integer , intent(out) :: ios if (trim(ftype)=='restart') then write(nop,4100,iostat=ios) cfld,coord,tlevel,xmin,xmax else if (trim(ftype)=="nersc_daily" .or. trim(ftype)=="nersc_weekly") then write(nop,117,iostat=ios) cfld,0,0.,coord,0.,xmin,xmax else if (trim(ftype)=="archv") then write(nop,118,iostat=ios) cfld,0,0.,coord,0.,xmin,xmax else print *,'unknown file type '//trim(ftype) stop end if 4100 format(a,': layer,tlevel,range = ',i3,i3,2x,1p2e16.7) 117 format (a8,' = ',i11,f11.2,i3,f7.3,1p2e16.7) 118 format (a8,' =',i11,f11.3,i3,f7.3,1p2e16.7) end subroutine !!!!!!!!!!! Header processing - write !!!!!!!!!!!!! subroutine HFwriteHeader(df,kdm) use mod_xc, only : idm, jdm implicit none type(hycomfile), intent(in) :: df integer , intent(in) :: kdm character(len=80) :: c80, cline character(len=80) :: ctitle(4) integer :: nop,k,i, ind integer :: iversn, iexpt, yrflag,lidm,ljdm,lkdm, dmonth integer :: syear, sday, dyear, dday nop=777 ctitle(1)='Generated by hycave' ctitle(2)='Generated by hycave' ctitle(3)='Generated by hycave' ctitle(4)='Generated by hycave' open(nop,file=trim(df%filebase)//'.b',status='replace') if (trim(df%ftype)=='restart') then write(nop,'(a,3i6)') 'RESTART: iexpt,iversn,yrflag = ', df%iexpt,df%iversn,df%yrflag write(cline,*) df%nstep,df%dtime write(nop,'(a,a)') 'RESTART: nstep,dtime = ',cline(1:len_trim(cline)) close(nop) else if (trim(df%ftype)=='nersc_daily') then write(nop,116) ctitle,df%iversn,df%iexpt,df%yrflag, & idm,jdm,kdm,9999,0 , & 9999,0,df%count do k=1,1000 write(nop,118) k,0 end do write(nop,119) else if (trim(df%ftype)=='nersc_weekly') then write(nop,216) ctitle,df%iversn,df%iexpt,df%yrflag, & idm,jdm,kdm,dmonth,9999, df%count else if (trim(df%ftype)=='archv') then write(nop,316) ctitle,df%iversn,df%iexpt,df%yrflag, & idm,jdm else print *,'(Unknown file type : '//trim(df%ftype)//')' end if close (nop) 116 format (a80/a80/a80/a80/ & i5,4x,'''iversn'' = hycom version number x10'/ & i5,4x,'''iexpt '' = experiment number x10'/ & i5,4x,'''yrflag'' = days in year flag'/ & i5,4x,'''idm '' = longitudinal array size'/ & i5,4x,'''jdm '' = latitudinal array size'/ & i5,4x,'''kdm '' = Vertical array size'/ & i5,4x,'''syear '' = Year of integration start '/ & i5,4x,'''sday '' = Day of integration start'/ & i5,4x,'''dyear '' = Year of this dump '/ & i5,4x,'''dday '' = Day of this dump '/ & i5,4x,'''count '' = Ensemble counter ') 216 format (a80/a80/a80/a80/ & i5,4x,'''iversn'' = hycom version number x10'/ & i5,4x,'''iexpt '' = experiment number x10'/ & i5,4x,'''yrflag'' = days in year flag'/ & i5,4x,'''idm '' = longitudinal array size'/ & i5,4x,'''jdm '' = latitudinal array size'/ & i5,4x,'''kdm '' = Vertical array size'/ & i5,4x,'''month '' = Month of this dump '/ & i5,4x,'''year '' = Year of this dump '/ & i5,4x,'''count '' = Averaging counter '/ & 'field time step model day', & ' k dens min max') 118 format ('member ',i5.5,' = ',i1,' Ensemble member flag') 119 format('field time step model day', & ' k dens min max') ! Archive format 316 format (a80/a80/a80/a80/ & i5,4x,'''iversn'' = hycom version number x10'/ & i5,4x,'''iexpt '' = experiment number x10'/ & i5,4x,'''yrflag'' = days in year flag'/ & i5,4x,'''idm '' = longitudinal array size'/ & i5,4x,'''jdm '' = latitudinal array size'/ & 'field time step model day', & ' k dens min max') end subroutine HFwriteHeader subroutine HFReadField3D(df,field,idm,jdm,kdm,cfld,tlevel) implicit none type(hycomfile), intent(in) :: df integer, intent(in) :: idm,jdm,kdm,tlevel real, intent(out) :: field(idm,jdm,kdm) character(len=8), intent(in) :: cfld integer :: k do k=1,kdm call HFReadField(df,field(1,1,k),idm,jdm,cfld,k,tlevel) end do end subroutine subroutine HFReadField(df,field,idm,jdm,cfld,coord,tlevel) implicit none type(hycomfile), intent(in) :: df integer, intent(in) :: idm,jdm,coord,tlevel real, intent(out) :: field(idm,jdm) character(len=8), intent(in) :: cfld character(len=8) :: cfld2 integer, dimension(idm,jdm) :: mask real, dimension(idm,jdm) :: dp real*4 :: A(idm,jdm), AMN, AMX, spval real :: xmin,dens,xmax,dens0,rday0 integer :: indx,nstep0 ! Some (unfortunate) corrections for ice variables cfld2=cfld if (trim(df%ftype)=="nersc_daily" .and. trim(cfld2)=='hicem') then cfld2='hice' elseif (trim(df%ftype)=="nersc_daily" .and. trim(cfld2)=='ficem') then cfld2='fice' end if call indexFromHeader(df,cfld2,coord,tlevel,indx) if (indx/=-1) then spval=undef call READRAW(A,AMN,AMX,IDM,JDM,.false.,spval,trim(df%filebase)//'.a',indx) field=A else print '(a,i3.3)', 'Could not get field "'//cfld//'" at coordinate',coord field=undef end if where(field > 0.5*huge) field=0. if (trim(df%ftype)=="nersc_daily") then field=field/df%count else if (trim(df%ftype)=="nersc_weekly") then if (is3DVar(df,cfld2,tlevel) .and. .not. isDPVar(df,cfld2)) then call HFReadDPField(df,dp,idm,jdm,coord,tlevel) field=field/max(dp,onecm) else field=field/df%count end if end if end subroutine HFReadField subroutine HFReadDPField(df,dp,idm,jdm,coord,tlevel) implicit none type(hycomfile), intent(in) :: df integer, intent(in) :: idm,jdm,coord,tlevel real, intent(out) :: dp(idm,jdm) real, dimension(idm,jdm) :: lw, up if (trim(df%ftype)=='restart') then call HFReadField(df,dp,idm,jdm,'dp ',coord,tlevel) elseif (trim(df%ftype)=="nersc_daily") then if (coord==1) then up=0. else call HFReadField(df,up,idm,jdm,'pres ',coord-1,tlevel) end if call HFReadField(df,lw,idm,jdm,'pres ',coord,tlevel) dp=lw-up elseif (trim(df%ftype)=="nersc_weekly") then call HFReadField(df,dp,idm,jdm,'pres ',coord,tlevel) elseif (trim(df%ftype)=="archv") then call HFReadField(df,dp,idm,jdm,'thknss ',coord,tlevel) else print *,'unknown file type '//trim(df%ftype) stop end if end subroutine HFReadDPField ! Special routine for reading total velocity subroutine HFReaduvtot(df,ut,vt,idm,jdm,vlevel,tlevel) implicit none type(hycomfile), intent(in) :: df integer, intent(in) :: idm,jdm,tlevel, vlevel real, intent(out) :: ut(idm,jdm) real, intent(out) :: vt(idm,jdm) real, dimension(idm,jdm) :: ub, vb if (trim(df%ftype)=='restart') then call HFReadField(df,ub,idm,jdm,'ubavg ',0,tlevel) call HFReadField(df,vb,idm,jdm,'vbavg ',0,tlevel) call HFReadField(df,ut,idm,jdm,'ut ',vlevel,tlevel) call HFReadField(df,vt,idm,jdm,'vt ',vlevel,tlevel) ut=ut+ub vt=vt+vb elseif (trim(df%ftype)=="nersc_daily") then call HFReadField(df,ut,idm,jdm,'utot ',vlevel,1) call HFReadField(df,vt,idm,jdm,'vtot ',vlevel,1) elseif (trim(df%ftype)=="nersc_weekly") then call HFReadField(df,ut,idm,jdm,'utot ',vlevel,1) call HFReadField(df,vt,idm,jdm,'vtot ',vlevel,1) elseif (trim(df%ftype)=="archv") then call HFReadField(df,ut,idm,jdm,'u-vel. ',vlevel,1) call HFReadField(df,vt,idm,jdm,'v-vel. ',vlevel,1) call HFReadField(df,ub,idm,jdm,'u_btrop ',0,1) call HFReadField(df,vb,idm,jdm,'v_btrop ',0,1) else print *,'unknown file type '//trim(df%ftype) stop end if ut=ut+ub vt=vt+vb end subroutine ! Special routine for reading barotropic velocity subroutine HFReaduvbaro(df,ub,vb,idm,jdm,tlevel) implicit none type(hycomfile), intent(in) :: df integer, intent(in) :: idm,jdm,tlevel real, intent(out) :: ub(idm,jdm) real, intent(out) :: vb(idm,jdm) real, dimension(idm,jdm) :: ut, vt, dp, dpsumu,dpsumv integer :: i,j,k, kdm real :: dpu, dpv kdm=vDim(df) if (trim(df%ftype)=='restart') then call HFReadField(df,ub,idm,jdm,'ubavg ',0,1) call HFReadField(df,vb,idm,jdm,'vbavg ',0,1) elseif (trim(df%ftype)=="nersc_daily") then ub=0. vb=0. dpsumu=0. dpsumv=0. do k=1,kdm call HFReadField(df,ut,idm,jdm,'utot ',k,tlevel) call HFReadField(df,vt,idm,jdm,'vtot ',k,tlevel) call HFReadDPField(df,dp,idm,jdm,k,tlevel) !print *,minval(ut),maxval(ut) ! TODO - need periodic fix do j=2,jdm do i=2,idm dpu=0.5*(dp(i,j)+dp(i-1,j)) dpv=0.5*(dp(i,j)+dp(i,j-1)) ub(i,j)=ub(i,j)+dpu*ut(i,j) vb(i,j)=vb(i,j)+dpv*vt(i,j) dpsumu(i,j)= dpsumu(i,j) + dpu dpsumv(i,j)= dpsumv(i,j) + dpv end do end do end do do j=1,jdm do i=1,idm ub(i,j)=ub(i,j)/max(0.1,dpsumu(i,j)) vb(i,j)=vb(i,j)/max(0.1,dpsumv(i,j)) end do end do !print *,minval(ub),maxval(ub) !print *,minval(vb),maxval(vb) elseif (trim(df%ftype)=="nersc_weekly") then call HFReadField(df,ub,idm,jdm,'ubavg ',0,1) call HFReadField(df,vb,idm,jdm,'vbavg ',0,1) elseif (trim(df%ftype)=="archv") then call HFReadField(df,ub,idm,jdm,'u_btrop ',0,1) call HFReadField(df,vb,idm,jdm,'v_btrop ',0,1) else print *,'unknown file type '//trim(df%ftype) stop end if end subroutine HFReaduvbaro subroutine HFWriteField(df,field,idm,jdm,cfld,coord,tlevel,indx) implicit none type(hycomfile) , intent(in) ::df integer, intent(in) :: idm,jdm,coord,indx,tlevel real, intent(in) :: field(idm,jdm) character(len=8), intent(in) :: cfld integer, dimension(idm,jdm) :: mask real*4 :: A(idm,jdm), AMN, AMX,spval real :: xmin,xmax integer :: nop,ios A=field spval=undef call WRITERAW(A,AMN,AMX,IDM,JDM,.false.,spval,trim(df%filebase)//'.a',indx) xmax=AMX ; xmin=AMN open(nop,file=trim(df%filebase)//'.b',action='write',form='formatted',status='old',position='append',iostat=ios) call writeFieldEntry(df%ftype,cfld,coord,tlevel,xmin,xmax,nop,ios) close(nop) end subroutine HFWriteField subroutine HFHeaderFromIndex(df,indx,ierr,varname,level,timelevel) implicit none type(hycomfile) , intent(in) :: df integer , intent(in) :: indx integer , intent(out) :: ierr character(len=8), intent(out), optional :: varname integer , intent(out), optional :: level integer , intent(out), optional :: timelevel ierr=0 if (indx> df%nrec .or. indx <1) then ierr=-1 else if (present(varname)) varname =df%cfld (indx) if (present(level )) level =df%coord(indx) if (present(timelevel )) timelevel =df%tlevel(indx) end if end subroutine ! Retrieve a-file index for variable subroutine indexFromHeader(df,cfld,coord,tlevel,indx) implicit none type(hycomfile) , intent(in) :: df character(len=8), intent(in) :: cfld integer , intent(in) :: coord,tlevel integer , intent(out) :: indx integer :: irec indx=-1 do irec=1,df%nrec if (trim(cfld)==trim(df%cfld(irec)) .and. coord==df%coord(irec) .and. df%tlevel(irec)==tlevel ) then indx=irec end if end do end subroutine subroutine HFUpdateAverage(dfave,df1) implicit none type(hycomfile), intent(in) :: df1 type(hycomfile), intent(inout) :: dfave integer :: k ! Consistency check and "addition" of restart average files if (df1%iversn/=dfave%iversn) then print *, 'iversn differ ' stop '(df_plus_df)' end if if (df1%iexpt/=dfave%iexpt) then print *, 'iexpt differ ' stop '(df_plus_df)' end if if (df1%yrflag/=dfave%yrflag) then print *, 'yrflag differ ' stop '(df_plus_df)' end if ! Ambigous dfave%nstep=9999 dfave%dtime=9999 !dfave%count=df1%count+dfave%count dfave%count=1 end subroutine logical function is3DVar(df,cfld,timelevel) implicit none type(hycomfile) , intent(in) :: df character(len=*), intent(in) :: cfld integer , intent(in) :: timelevel character(len=8) char8 char8=adjustl(cfld) !print *,count( df%cfld == char8 .and. df%tlevel==timelevel ) is3DVar=count( df%cfld == char8 .and. df%tlevel==timelevel ) > 1 end function logical function isDPVar(df,cfld) implicit none type(hycomfile) , intent(in) :: df character(len=*), intent(in) :: cfld if (trim(df%ftype)=='restart') then isDPVar=trim(cfld)=='dp' else if (trim(df%ftype)=='nersc_daily' .or. trim(df%ftype)=='nersc_weekly') then isDPVar=trim(cfld)=='pres' elseif (trim(df%ftype)=='archv') then isDPVar=trim(cfld)=='tknss' else print *,'Unknown file type '//trim(df%ftype) stop '(mod_hycomfile_io:isDPVar)' end if end function integer function vDim(df) implicit none type(hycomfile) , intent(in) :: df if (trim(df%ftype)=='restart') then vDim=count( df%cfld == 'dp ' .and. df%tlevel==1 ) else if (trim(df%ftype)=='nersc_daily' .or. trim(df%ftype)=='nersc_weekly') then vDim=count( df%cfld == 'pres ' .and. df%tlevel==1 ) elseif (trim(df%ftype)=='archv') then vDim=count( df%cfld == 'thknss ' .and. df%tlevel==1 ) else print *,'Unknown file type '//trim(df%ftype) stop '(mod_hycomfile_io:isDPVar)' end if end function real function fyear(df) implicit none type(hycomfile) , intent(in) :: df fyear=df%fyear end function ! Returns file type based on file name function getfiletype(filename) implicit none character(len=*), intent(in) :: filename character(len=20) :: getfiletype integer :: findhdr,findab,finddaily,findweek,findrst, findarchv ! Check for type ... findhdr =index(filename,'.hdr') findab=max(index(filename,'.a'),index(filename,'.b')) if (.not. findab>0 .and. .not. findhdr>0) then print *,'No .ab or .hdr files' stop '(mod_hycomfile_io:getfiletype)' else ! We have a .ab-file. Now figure out what type. findrst =index(filename,'restart') finddaily=index(filename,'DAILY') findweek =index(filename,'AVE') findhdr =index(filename,'.hdr') findarchv=index(filename,'archv.') if (findrst==4) then getfiletype='restart' elseif (finddaily==4) then getfiletype='nersc_daily' elseif (findweek==4) then getfiletype='nersc_weekly' elseif (findarchv>0) then getfiletype='archv' elseif (findhdr>0) then getfiletype='pak' print *,'pak files no longer supported in this version' stop '(mod_hycomfile_io:getfiletype)' else print *,'Can not deduce file type from file name' stop '(mod_hycomfile_io:getfiletype)' end if end if end function getfiletype subroutine netcdfInfo(vnamein,gridrotate,stdname,units,vname,cellmethod,limits) implicit none character(len=*), intent(in) :: vnamein logical , intent(in) :: gridrotate character(len=*) ,intent(out) :: stdName ! standard_name attrib in netcdf file character(len=*) ,intent(out) :: units ! Units attribute in netcdf file character(len=*) ,intent(out) :: vname ! Name of variable in netcdf file character(len=*) ,intent(out) :: cellmethod ! Name of cell method real ,intent(out) :: limits(2) ! Lower and upper limits ! Standard name and unit lookup table for a given var name units='' cellmethod='area: mean' stdname='' limits=(/0,0/) select case (trim(vnamein)) case ('saln','salin') stdname='sea_water_salinity' ; units='1e-3' ; vname='salinity' limits=(/0,45/) case ('temp') stdname='sea_water_potential_temperature' ; units='K' ; vname='temperature' limits=(/-3,50/) case ('levsaln') stdname='sea_water_salinity' ; units='1e-3' ; vname='levitus_salinity' limits=(/0,45/) case ('levtemp') stdname='sea_water_potential_temperature' ; units='K' ; vname='levitus_temperature' limits=(/-3,50/) case ('ssh','srfhgt') stdname='sea_surface_elevation' ; units='m' ; vname='ssh' limits=(/-4,4/) case ('utot') if (.not.gridrotate) then stdname='eastward_sea_water_velocity' else stdname='x_sea_water_velocity' end if units='m s-1' ; vname='utot' limits=(/-3,3/) case ('vtot') if (.not.gridrotate) then stdname='northward_sea_water_velocity' else stdname='y_sea_water_velocity' end if units='m s-1' ; vname='vtot' limits=(/-3,3/) case ('u','u-vel.') if (.not.gridrotate) then stdname='baroclinic_eastward_sea_water_velocity' else stdname='baroclinic_x_sea_water_velocity' end if units='m s-1' ; vname='utot' limits=(/-3,3/) case ('v','v-vel.') if (.not.gridrotate) then stdname='baroclinic_northward_sea_water_velocity' else stdname='baroclinic_y_sea_water_velocity' end if units='m s-1' ; vname='vtot' limits=(/-3,3/) case ('hice','hicem') stdname='sea_ice_thickness' ; units='m' ; vname='hice' cellmethod='area: mean where sea_ice' limits=(/0,20/) case ('hsnw','hsnwm') stdname='surface_snow_thickness' ; units='m' ; vname='hsnw' cellmethod='area: mean where sea_ice' limits=(/0,3/) case ('fice','ficem') stdname='sea_ice_concentration' ; units='1' ; vname='fice' limits=(/0,1/) case ('ubavg','u_btrop') if (.not.gridrotate) then stdname='barotropic_eastward_sea_water_velocity' else stdname='barotropic_x_sea_water_velocity' end if units='m s-1' ; vname='ub' limits=(/-3,3/) case ('vbavg','v_btrop') if (.not.gridrotate) then stdname='barotropic_northward_sea_water_velocity' else stdname='barotropic_y_sea_water_velocity' end if units='m s-1' ; vname='vb' limits=(/-3,3/) case ('uice') if (.not.gridrotate) then stdname='eastward_sea_ice_velocity' else stdname='sea_ice_x_velocity' end if cellmethod='area: mean where sea_ice' units='m s-1' ; vname='uice' limits=(/-3,3/) case ('vice') if (.not.gridrotate) then stdname='northward_sea_ice_velocity' else stdname='sea_ice_y_velocity' end if cellmethod='area: mean where sea_ice' units='m s-1' ; vname='vice' limits=(/-3,3/) case ('taux') if (.not.gridrotate) then stdname='surface_downward_eastward_stress' else stdname='surface_downward_x_stress' end if units='pascal' ; vname='taux' limits=(/-3,3/) case ('tauy') if (.not.gridrotate) then stdname='surface_downward_north_stress' else stdname='surface_downward_y_stress' end if units='pascal' ; vname='tauy' limits=(/-3,3/) case default vname=trim(vnamein) end select end subroutine !TODO: Fix for different yrflag subroutine forecastDate(hfile,rt) use mod_year_info implicit none type(hycomfile), intent(in) :: hfile type(year_info), intent(out) :: rt integer diy if (trim(hfile%ftype)=='nersc_daily') then call year_day(real(hfile%iday),hfile%iyear,rt,'ecmwf') else if (trim(hfile%ftype)=='nersc_weekly') then ! Only approxmate, day in year diy=3+(hfile%iweek-1)*7 + (hfile%imonth-1)*30 call year_day(real(diy),hfile%iyear,rt,'ecmwf') else if (trim(hfile%ftype)=='restart') then call year_day(real(hfile%iday),hfile%iyear,rt,'ecmwf') else if (trim(hfile%ftype)=='archv') then call year_day(real(hfile%iday),hfile%iyear,rt,'ecmwf') else write(6,*) 'Unknown file type '//trim(hfile%ftype) call exit(1) end if end subroutine !TODO: Fix for different yrflag subroutine startDate(hfile,rt) use mod_year_info implicit none type(hycomfile), intent(in) :: hfile type(year_info), intent(out) :: rt if (trim(hfile%ftype)=='nersc_daily') then call year_day(real(hfile%start_iday),hfile%start_iyear,rt,'ecmwf') else if (trim(hfile%ftype)=='nersc_weekly') then call forecastDate(hfile,rt) else if (trim(hfile%ftype)=='restart') then call forecastDate(hfile,rt) else if (trim(hfile%ftype)=='archv') then call forecastDate(hfile,rt) else write(6,*) 'Unknown file type '//trim(hfile%ftype) call exit(1) end if end subroutine ! Modified from Alan Wallcraft's RAW routine by Knut Liseter @ NERSC ! So far only the "I" in "IO" is present SUBROUTINE READRAW(A,AMN,AMX,IDM,JDM,LSPVAL,SPVAL,CFILE1,K) IMPLICIT NONE ! REAL*4 SPVALH PARAMETER (SPVALH=1e30) ! REAL*4, INTENT(OUT) :: A(IDM,JDM) REAL*4, INTENT(OUT) :: AMN,AMX INTEGER, INTENT(IN) :: IDM,JDM LOGICAL, INTENT(IN) :: LSPVAL REAL*4, INTENT(INOUT) :: SPVAL INTEGER, INTENT(IN) :: K CHARACTER(len=*), INTENT(IN) :: CFILE1 ! REAL*4 :: PADA(4096) ! MOST OF WORK IS DONE HERE. ! INTEGER I,J,IOS,NRECL INTEGER NPAD ! IF(.NOT.LSPVAL) THEN SPVAL = SPVALH ENDIF ! !!! Calculate the number of elements padded!!!!!!!!!!!!!!!!!!!!!!!! NPAD=GET_NPAD(IDM,JDM) INQUIRE( IOLENGTH=NRECL) A,PADA(1:NPAD) ! ! OPEN(UNIT=11, FILE=CFILE1, FORM='UNFORMATTED', STATUS='old', & ACCESS='DIRECT', RECL=NRECL, IOSTAT=IOS) IF (IOS.NE.0) THEN write(6,*) 'Error: can''t open ',CFILE1(1:LEN_TRIM(CFILE1)) write(6,*) 'ios = ',ios write(6,*) 'nrecl = ',nrecl CALL EXIT(3) ENDIF ! READ(11,REC=K,IOSTAT=IOS) A close(11) ! IF (IOS.NE.0) THEN WRITE(6,*) 'can''t read record ',K, & ' from '//CFILE1(1:LEN_TRIM(CFILE1)) CALL EXIT(4) ENDIF ! AMN = SPVALH AMX = -SPVALH DO J= 1,JDM DO I=1,IDM IF (A(I,J).LE.SPVALH) THEN AMN = MIN( AMN, A(I,J) ) AMX = MAX( AMX, A(I,J) ) ELSEIF (LSPVAL) THEN A(I,J) = SPVAL ENDIF END DO END DO ! RETURN END SUBROUTINE ! Modified from Alan Wallcraft's RAW routine by Knut Liseter @ NERSC SUBROUTINE WRITERAW(A,AMN,AMX,IDM,JDM,LSPVAL,SPVAL,CFILE1,K) IMPLICIT NONE REAL*4 SPVALH PARAMETER (SPVALH=1e30) REAL*4, INTENT(INOUT) :: A(IDM,JDM) REAL*4, INTENT(OUT) :: AMN,AMX INTEGER, INTENT(IN) :: IDM,JDM LOGICAL, INTENT(IN) :: LSPVAL REAL*4, INTENT(INOUT) :: SPVAL INTEGER, INTENT(IN) :: K CHARACTER(len=*), INTENT(IN) :: CFILE1 ! REAL*4 :: PADA(4096) ! ! MOST OF WORK IS DONE HERE. ! CHARACTER*18 CASN INTEGER LEN_TRIM INTEGER I,J,IOS,NRECL INTEGER NPAD ! IF(.NOT.LSPVAL) THEN SPVAL = SPVALH ENDIF ! !!! Calculate the number of elements padded!!!!!!!!!!!!!!!!!!!!!!!! NPAD=GET_NPAD(IDM,JDM) PADA=0. INQUIRE( IOLENGTH=NRECL) A,PADA(1:NPAD) OPEN(UNIT=11, FILE=CFILE1, FORM='UNFORMATTED', STATUS='unknown', & ACCESS='DIRECT', RECL=NRECL, IOSTAT=IOS) IF (IOS.NE.0) THEN write(6,*) 'Error: can''t open ',CFILE1(1:LEN_TRIM(CFILE1)) write(6,*) 'ios = ',ios write(6,*) 'nrecl = ',nrecl CALL EXIT(3) ENDIF ! WRITE(11,REC=K,IOSTAT=IOS) A,PADA(1:NPAD) close(11) IF (IOS.NE.0) THEN WRITE(6,*) 'can''t write record ',K, & ' from '//CFILE1(1:LEN_TRIM(CFILE1)) CALL EXIT(4) ENDIF ! AMN = SPVALH AMX = -SPVALH DO J= 1,JDM DO I=1,IDM IF (A(I,J).LE.SPVALH) THEN AMN = MIN( AMN, A(I,J) ) AMX = MAX( AMX, A(I,J) ) ELSEIF (LSPVAL) THEN A(I,J) = SPVAL ENDIF END DO END DO ! RETURN END SUBROUTINE !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! INTEGER FUNCTION GET_NPAD(IDM,JDM) IMPLICIT NONE INTEGER, INTENT(IN) :: IDM,JDM GET_NPAD = 4096 - MOD(IDM*JDM,4096) GET_NPAD = mod(GET_NPAD,4096) END FUNCTION end module