Back to home page

MITgcm
 
 

    

File indexing completed on 2018-03-02 18:45:27 UTC

view on githubraw file Latest commit 95ec9058 on 2009-04-02 17:49:59 UTC
95ec9058d9 Davi*0001 program cvfloat
                0002    !
                0003    !=======================================================================
                0004    !
                0005    !     converts binary float trajectories to netCDF
                0006    !
                0007    !=======================================================
                0008    !
                0009    !     * uses namelist data.float
                0010    !
                0011    !     Arne Biastoch, abiastoch@ucsd.edu, 11/16/2000
                0012    !     Updated 03/16/2009 Christopher Wolfe, clwolfe@ucsd.edu
                0013    !
                0014    !=======================================================================
                0015    !
                0016    use netcdf
                0017    implicit none
                0018    integer, parameter :: STDIN = 5, STDOUT = 6, STDERR = 6
                0019    
                0020    character(80) :: iotext
                0021    character(32) :: stamp
                0022    integer :: ioun,iargc,m,n,ilen,irec,ip,np
                0023    integer :: icount,itotalrecord,imaxrecord
                0024 
                0025    ! tile counters
                0026    integer :: iGmax,jGmax,iG,jG
                0027 
                0028    ! variables for filenames
                0029    integer :: narg, npart
                0030    character(80) :: dataFName
                0031    logical :: exst
                0032 
                0033    ! number of variables per record
                0034    integer, parameter :: IMAX = 13
                0035    !
                0036    !      integer narg
                0037    logical preflag
                0038    !
                0039    ! netCDF ids
                0040    !
                0041    integer :: ncid
                0042    integer :: partdim,Timedim
                0043    integer :: partid, Timeid
                0044    integer :: xpartid, ypartid, zpartid
                0045    integer :: ipartid, jpartid, kpartid
                0046    integer  tempid, saltid, uvelid, vvelid, presid
                0047    !
                0048    ! attribute vectors
                0049    !
                0050    character(24) :: name*24
                0051    character(16) :: unit*16
                0052    ! 
                0053    ! data variables for NetCDF
                0054    !
                0055    double precision :: rcountstart, rcountdelta
                0056    double precision, dimension(:),   allocatable :: pnum,time
                0057    double precision, dimension(:,:), allocatable :: xpart,ypart,zpart, &
                0058                                                    ipart,jpart,kpart, &
                0059                                                    pres,uvel,vvel,temp,salt
                0060    double precision, dimension(:), allocatable ::  tmp
                0061    !
                0062    ! namelist contains
                0063    !
                0064    character(50) :: outname2, fName="float_trajectories", outname
                0065    character(20) :: startDate="01-JAN-1992:12:00:00"
                0066    character(60) :: expnam = "Experiment name not set"
                0067    logical :: usingSphericalPolarGrid=.true.
                0068    
                0069    namelist /dimensions/ expnam, startDate, usingSphericalPolarGrid
                0070    namelist /floats/ fName
                0071 
                0072    ioun=11
                0073    open(ioun,file="data.float",status="old",form="formatted")
                0074    read  (unit=ioun, end=666, nml=dimensions)
                0075    write (STDOUT,dimensions)
                0076    close (ioun)
                0077 666  continue
                0078    open(ioun,file="data.float",status="old",form="formatted")
                0079    read  (unit=ioun, end=999, nml=floats)
                0080    write (STDOUT,floats)
                0081    close (ioun)
                0082 999  continue
                0083 
                0084    !
                0085    !     PRELIMINARY USE:
                0086    !     IF FLOATS SHOULD BE VIEWED DURING A CURRENT MODEL RUN THE FIRST
                0087    !     LINE OF THE FILE MAY NOT BE UPDATED CORRECTLY, I.E. THERE MIGHT
                0088    !     BE MORE TIMES THAN STATED AT THE BEGINNING. BY GIVING A FLAG
                0089    !     ONLY ICOUNT-1 TIMESTEPS ARE USED
                0090    !
                0091    preflag = .false.
                0092    narg=iargc()
                0093    if ( narg > 0 ) preflag = .true.
                0094    !
                0095    ! check existent files
                0096    !
                0097    iGmax=1
                0098    do m=1,100
                0099       write(dataFname(1:80),"(2a,i3.3,a,i3.3,a)") trim(fName),".",iGmax,".",1,".data"
                0100       inquire( file=dataFname, exist=exst )
                0101       if (exst)  iGmax = iGmax + 1
                0102    enddo
                0103    
                0104    jGmax=1
                0105    do m=1,100
                0106       write(dataFname(1:80),"(2a,i3.3,a,i3.3,a)") trim(fName),".",1,".",jGmax,".data"
                0107       inquire( file=dataFname, exist=exst )
                0108       if (exst)  jGmax = jGmax + 1
                0109    enddo
                0110    
                0111    iGmax = iGmax - 1
                0112    jGmax = jGmax - 1
                0113    print*, "There are ",iGmax," x ",jGmax," files"
                0114 
                0115    ! open first file and get dimensions (float number and time)
                0116    ilen=IMAX*8
                0117    allocate (tmp(IMAX))
                0118 
                0119    write(dataFname(1:80),"(2a,a)") trim(fName),".001.001.data"
                0120    open(1,file=dataFname,status="old",form="unformatted",access="direct",recl=ilen)
                0121    
                0122    read(1,rec=1) tmp
                0123 !   print*,"tmp:", tmp
                0124    
                0125    rcountstart = tmp(2)
                0126    rcountdelta = tmp(4)
                0127    icount      = INT(tmp(5))
                0128    npart       = INT(tmp(6))
                0129    close(1)
                0130    
                0131    print*, "npart    = ",npart
                0132    print*, "timesteps= ",icount
                0133 !    print*,"rcountstart=",rcountstart,"rcountdelta",rcountdelta
                0134    if (preflag) then
                0135       icount=icount-1
                0136       print*, "preliminary --> use timesteps= ",icount
                0137    endif
                0138 
                0139    !-----------------------------------------------------------------------
                0140    !     allocate variables
                0141    !-----------------------------------------------------------------------
                0142    !
                0143    allocate (pnum(npart))
                0144    allocate (time(icount))
                0145    allocate (xpart(npart,icount))
                0146    allocate (ypart(npart,icount))
                0147    allocate (zpart(npart,icount))
                0148    allocate (ipart(npart,icount))
                0149    allocate (jpart(npart,icount))
                0150    allocate (kpart(npart,icount))
                0151    allocate (uvel(npart,icount))
                0152    allocate (vvel(npart,icount))
                0153    allocate (temp(npart,icount))
                0154    allocate (salt(npart,icount))
                0155    allocate (pres(npart,icount))
                0156 
                0157    ! initialize arrays
                0158    !
                0159    do m=1,npart
                0160       do n=1,icount
                0161          xpart(m,n) = NF90_FILL_DOUBLE
                0162          ypart(m,n) = NF90_FILL_DOUBLE
                0163          zpart(m,n) = NF90_FILL_DOUBLE
                0164          ipart(m,n) = NF90_FILL_DOUBLE
                0165          jpart(m,n) = NF90_FILL_DOUBLE
                0166          kpart(m,n) = NF90_FILL_DOUBLE
                0167          uvel(m,n) = NF90_FILL_DOUBLE
                0168          vvel(m,n) = NF90_FILL_DOUBLE
                0169          temp(m,n) = NF90_FILL_DOUBLE
                0170          salt(m,n) = NF90_FILL_DOUBLE
                0171          pres(m,n) = NF90_FILL_DOUBLE
                0172       enddo
                0173    enddo
                0174 
                0175    ! generate axes
                0176    !
                0177    time(1)=rcountstart
                0178    do m=2,icount
                0179       time(m) = time(m-1)+rcountdelta
                0180    enddo
                0181    print*, "time: ",time(1:icount)
                0182    
                0183    do ip=1,npart
                0184       pnum(ip) = FLOAT(ip)
                0185    enddo
                0186    !
                0187    !-----------------------------------------------------------------------
                0188    !     open files and read input
                0189    !-----------------------------------------------------------------------
                0190    !
                0191 
                0192    itotalrecord = 0
                0193    
                0194    do iG=1,iGmax
                0195       do jG=1,jGmax
                0196    
                0197          write(dataFname(1:80),"(2a,i3.3,a,i3.3,a)") trim(fName),".",iG,".",jG,".data"
                0198          open(1,file=dataFname,status="old",form="unformatted",access="direct",recl=ilen)
                0199    
                0200          read(1,rec=1) tmp
                0201          imaxrecord = INT(tmp(1))
                0202          print "(1x,2a)","read ",dataFname
                0203          itotalrecord = itotalrecord + imaxrecord
                0204    
                0205          do irec=2,imaxrecord+1
                0206    
                0207             read(1,rec=irec) tmp
                0208             ip = INT(tmp(1))
                0209             if (ip > npart) then
                0210                print*,"ip out of order: ",ip,npart
                0211                stop
                0212             endif
                0213             
                0214             ! Note: If the float initial conditions are also written to the data files the time interval
                0215             ! between record 1 and 2 may be somewhat less than rcountdelta. The +0.9999 deals with
                0216             ! this case without upsetting the rest of the indexing.
                0217             np = INT((tmp(2)-rcountstart)/rcountdelta + 1 + 0.9999)
                0218    
                0219             ! this is only for prelimiray results. Use only icount-1 timesteps
                0220             if (preflag .and. (np > icount .or. np < 1)) cycle
                0221    
                0222             if (usingSphericalPolarGrid) then
                0223                xpart(ip,np)  = tmp(3)
                0224                ypart(ip,np)  = tmp(4)
                0225             else
                0226                xpart(ip,np)  = tmp(3)
                0227                ypart(ip,np)  = tmp(4)
                0228             endif
                0229             zpart(ip,np)  = tmp(5)
                0230             ipart(ip,np)  = tmp(6)
                0231             jpart(ip,np)  = tmp(7)
                0232             kpart(ip,np)  = tmp(8)
                0233             pres(ip,np)   = tmp(9)
                0234             uvel(ip,np)   = tmp(10)
                0235             vvel(ip,np)   = tmp(11)
                0236             temp(ip,np)   = tmp(12)
                0237             salt(ip,np)   = tmp(13)
                0238          enddo
                0239    
                0240          close(1)
                0241       enddo
                0242    enddo
                0243 
                0244    print*,icount," x ",npart," = ",icount*npart," records expected,",& 
                0245       " found ",itotalrecord," float records"
                0246    print*,"==> ",icount*npart-itotalrecord," float records missing"
                0247    !
                0248    !
                0249    !-----------------------------------------------------------------------
                0250    !     define netCDF variables
                0251    !-----------------------------------------------------------------------
                0252    !
                0253    write(STDOUT,*) " Start Converting"
                0254    !
                0255    ! enter define mode
                0256    !
                0257    outname2=trim(fname)//".cdf"
                0258    write (STDOUT,*)" ==>  Writing a trajectories to file ",trim(outname2)
                0259    call nc_check(nf90_create(trim(outname2),NF90_CLOBBER,ncid),__LINE__)
                0260    ! 
                0261    ! define dimensions
                0262    !
                0263    call nc_check(nf90_def_dim(ncid, "Particles", npart, partdim),__LINE__)
                0264    call nc_check(nf90_def_dim(ncid, "Time",      NF90_UNLIMITED, Timedim),__LINE__)
                0265    !
                0266    ! define variables
                0267    !
                0268    call nc_check(nf90_def_var(ncid, "Particles",NF90_DOUBLE, (/partdim/), partid),__LINE__)
                0269    call nc_check(nf90_def_var(ncid, "Time",     NF90_DOUBLE, (/Timedim/), Timeid),__LINE__)
                0270    
                0271    call nc_check(nf90_def_var(ncid, "x",        NF90_DOUBLE, (/partdim,Timedim/), xpartid),__LINE__)
                0272    call nc_check(nf90_def_var(ncid, "y",        NF90_DOUBLE, (/partdim,Timedim/), ypartid),__LINE__)
                0273    call nc_check(nf90_def_var(ncid, "z",        NF90_DOUBLE, (/partdim,Timedim/), zpartid),__LINE__)
                0274    call nc_check(nf90_def_var(ncid, "i",        NF90_DOUBLE, (/partdim,Timedim/), ipartid),__LINE__)
                0275    call nc_check(nf90_def_var(ncid, "j",        NF90_DOUBLE, (/partdim,Timedim/), jpartid),__LINE__)
                0276    call nc_check(nf90_def_var(ncid, "k",        NF90_DOUBLE, (/partdim,Timedim/), kpartid),__LINE__)
                0277    call nc_check(nf90_def_var(ncid, "pressure", NF90_DOUBLE, (/partdim,Timedim/), presid),__LINE__)
                0278    call nc_check(nf90_def_var(ncid, "u",        NF90_DOUBLE, (/partdim,Timedim/), uvelid),__LINE__)
                0279    call nc_check(nf90_def_var(ncid, "v",        NF90_DOUBLE, (/partdim,Timedim/), vvelid),__LINE__)
                0280    call nc_check(nf90_def_var(ncid, "T",        NF90_DOUBLE, (/partdim,Timedim/), tempid),__LINE__)
                0281    call nc_check(nf90_def_var(ncid, "S",        NF90_DOUBLE, (/partdim,Timedim/), saltid),__LINE__)
                0282    !
                0283    !-----------------------------------------------------------------------
                0284    !     assign attributes
                0285    !-----------------------------------------------------------------------
                0286    !
                0287    name = "Particle Number"
                0288    unit = "particle number"
                0289    call nc_check(nf90_put_att(ncid,partid,"long_name",name),__LINE__)
                0290    call nc_check(nf90_put_att(ncid,partid,"units",    unit),__LINE__)
                0291    
                0292    name = "Time"
                0293    unit = "seconds"
                0294    call nc_check(nf90_put_att(ncid,Timeid,"long_name",  name),__LINE__)
                0295    call nc_check(nf90_put_att(ncid,Timeid,"units",      unit),__LINE__)
                0296    call nc_check(nf90_put_att(ncid,Timeid,"time_origin",startDate),__LINE__)
                0297    
                0298    if (usingSphericalPolarGrid) then
                0299       name = "LONGITUDE "
                0300       unit = "degrees_W "
                0301    else
                0302       name = "X_t "
                0303       unit = "meter "
                0304    endif
                0305    call nc_check(nf90_put_att(ncid,xpartid,"long_name",  name),__LINE__)
                0306    call nc_check(nf90_put_att(ncid,xpartid,"units",      unit),__LINE__)
                0307    call nc_check(nf90_put_att(ncid,xpartid,"missing_value",NF90_FILL_DOUBLE),__LINE__)
                0308    call nc_check(nf90_put_att(ncid,xpartid,"_FillValue",   NF90_FILL_DOUBLE),__LINE__)
                0309    
                0310    if (usingSphericalPolarGrid) then
                0311       name = "LATITUDE "
                0312       unit = "degrees_N "
                0313    else
                0314       name = "Y_t "
                0315       unit = "meter "
                0316    endif
                0317    call nc_check(nf90_put_att(ncid,ypartid,"long_name",  name),__LINE__)
                0318    call nc_check(nf90_put_att(ncid,ypartid,"units",      unit),__LINE__)
                0319    call nc_check(nf90_put_att(ncid,ypartid,"missing_value",NF90_FILL_DOUBLE),__LINE__)
                0320    call nc_check(nf90_put_att(ncid,ypartid,"_FillValue",   NF90_FILL_DOUBLE),__LINE__)
                0321    
                0322    name = "DEPTH "
                0323    unit = "meter "
                0324    call nc_check(nf90_put_att(ncid,zpartid,"long_name",  name),__LINE__)
                0325    call nc_check(nf90_put_att(ncid,zpartid,"units",      unit),__LINE__)
                0326    call nc_check(nf90_put_att(ncid,zpartid,"missing_value",NF90_FILL_DOUBLE),__LINE__)
                0327    call nc_check(nf90_put_att(ncid,zpartid,"_FillValue",   NF90_FILL_DOUBLE),__LINE__)
                0328    
                0329    name = "iINDEX "
                0330    unit = "index "
                0331    call nc_check(nf90_put_att(ncid,ipartid,"long_name",  name),__LINE__)
                0332    call nc_check(nf90_put_att(ncid,ipartid,"units",      unit),__LINE__)
                0333    call nc_check(nf90_put_att(ncid,ipartid,"missing_value",NF90_FILL_DOUBLE),__LINE__)
                0334    call nc_check(nf90_put_att(ncid,ipartid,"_FillValue",   NF90_FILL_DOUBLE),__LINE__)
                0335    
                0336    name = "jINDEX "
                0337    unit = "index "
                0338    call nc_check(nf90_put_att(ncid,jpartid,"long_name",  name),__LINE__)
                0339    call nc_check(nf90_put_att(ncid,jpartid,"units",      unit),__LINE__)
                0340    call nc_check(nf90_put_att(ncid,jpartid,"missing_value",NF90_FILL_DOUBLE),__LINE__)
                0341    call nc_check(nf90_put_att(ncid,jpartid,"_FillValue",   NF90_FILL_DOUBLE),__LINE__)
                0342    
                0343    name = "LEVEL "
                0344    unit = "level "
                0345    call nc_check(nf90_put_att(ncid,kpartid,"long_name",  name),__LINE__)
                0346    call nc_check(nf90_put_att(ncid,kpartid,"units",      unit),__LINE__)
                0347    call nc_check(nf90_put_att(ncid,kpartid,"missing_value",NF90_FILL_DOUBLE),__LINE__)
                0348    call nc_check(nf90_put_att(ncid,kpartid,"_FillValue",   NF90_FILL_DOUBLE),__LINE__)
                0349    
                0350    name = "POTENTIAL TEMPERATURE "
                0351    unit = "deg C "
                0352    call nc_check(nf90_put_att(ncid,tempid,"long_name",  name),__LINE__)
                0353    call nc_check(nf90_put_att(ncid,tempid,"units",      unit),__LINE__)
                0354    call nc_check(nf90_put_att(ncid,tempid,"missing_value",NF90_FILL_DOUBLE),__LINE__)
                0355    call nc_check(nf90_put_att(ncid,tempid,"_FillValue",   NF90_FILL_DOUBLE),__LINE__)
                0356    
                0357    name = "SALINITY "
                0358    unit = "PSU "
                0359    call nc_check(nf90_put_att(ncid,saltid,"long_name",  name),__LINE__)
                0360    call nc_check(nf90_put_att(ncid,saltid,"units",      unit),__LINE__)
                0361    call nc_check(nf90_put_att(ncid,saltid,"missing_value",NF90_FILL_DOUBLE),__LINE__)
                0362    call nc_check(nf90_put_att(ncid,saltid,"_FillValue",   NF90_FILL_DOUBLE),__LINE__)
                0363    
                0364    name = "U VELOCITY "
                0365    unit = "m/sec"
                0366    call nc_check(nf90_put_att(ncid,uvelid,"long_name",  name),__LINE__)
                0367    call nc_check(nf90_put_att(ncid,uvelid,"units",      unit),__LINE__)
                0368    call nc_check(nf90_put_att(ncid,uvelid,"missing_value",NF90_FILL_DOUBLE),__LINE__)
                0369    call nc_check(nf90_put_att(ncid,uvelid,"_FillValue",   NF90_FILL_DOUBLE),__LINE__)
                0370    
                0371    name = "V VELOCITY "
                0372    unit = "m/sec"
                0373    call nc_check(nf90_put_att(ncid,vvelid,"long_name",  name),__LINE__)
                0374    call nc_check(nf90_put_att(ncid,vvelid,"units",      unit),__LINE__)
                0375    call nc_check(nf90_put_att(ncid,vvelid,"missing_value",NF90_FILL_DOUBLE),__LINE__)
                0376    call nc_check(nf90_put_att(ncid,vvelid,"_FillValue",   NF90_FILL_DOUBLE),__LINE__)
                0377    
                0378    name = "PRESSURE "
                0379    unit = "N/m^2 "
                0380    call nc_check(nf90_put_att(ncid,presid,"long_name",  name),__LINE__)
                0381    call nc_check(nf90_put_att(ncid,presid,"units",      unit),__LINE__)
                0382    call nc_check(nf90_put_att(ncid,presid,"missing_value",NF90_FILL_DOUBLE),__LINE__)
                0383    call nc_check(nf90_put_att(ncid,presid,"_FillValue",   NF90_FILL_DOUBLE),__LINE__)
                0384    
                0385    
                0386 !   expnam= " "
                0387 !   stamp = " "
                0388    call nc_check(nf90_put_att(ncid,NF90_GLOBAL,"title",  trim(expnam)),__LINE__)
                0389 !   call nc_check(nf90_put_att(ncid,NF90_GLOBAL,"history",stamp),__LINE__)
                0390    !
                0391    !-----------------------------------------------------------------------
                0392    !     leave define mode
                0393    !-----------------------------------------------------------------------
                0394    !
                0395    call nc_check(nf90_enddef(ncid),__LINE__)
                0396    !
                0397    !
                0398    !-----------------------------------------------------------------------
                0399    !     put variables in netCDF file
                0400    !-----------------------------------------------------------------------
                0401    !
                0402    ! store Particles
                0403    call nc_check(nf90_put_var(ncid,partid,pnum),__LINE__)
                0404    !
                0405    ! store Time      
                0406    call nc_check(nf90_put_var(ncid,Timeid,Time),__LINE__)
                0407    !
                0408    ! store values
                0409    call nc_check(nf90_put_var(ncid,xpartid,xpart),__LINE__)
                0410    call nc_check(nf90_put_var(ncid,ypartid,ypart),__LINE__)
                0411    call nc_check(nf90_put_var(ncid,zpartid,zpart),__LINE__)
                0412    call nc_check(nf90_put_var(ncid,ipartid,ipart),__LINE__)
                0413    call nc_check(nf90_put_var(ncid,jpartid,jpart),__LINE__)
                0414    call nc_check(nf90_put_var(ncid,kpartid,kpart),__LINE__)
                0415    call nc_check(nf90_put_var(ncid,tempid, temp),__LINE__)
                0416    call nc_check(nf90_put_var(ncid,saltid, salt),__LINE__)
                0417    call nc_check(nf90_put_var(ncid,uvelid, uvel),__LINE__)
                0418    call nc_check(nf90_put_var(ncid,vvelid, vvel),__LINE__)
                0419    call nc_check(nf90_put_var(ncid,presid, pres),__LINE__)
                0420    
                0421    call nc_check(nf90_close(ncid),__LINE__)
                0422    
                0423    write(STDOUT,*) " End "
                0424       
                0425    !-----------------------------------------------------------------------
                0426    !     internal subroutines
                0427    !-----------------------------------------------------------------------
                0428 contains
                0429    subroutine nc_check(status,lineno)
                0430       integer, intent ( in) :: status
                0431       integer, intent ( in) :: lineno
                0432    
                0433       if(status /= nf90_noerr) then
                0434          print *, "Error at line number ",lineno,":"
                0435          print *, trim(nf90_strerror(status))
                0436          stop 2
                0437       end if
                0438    end subroutine nc_check
                0439 end program cvfloat
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated from https://github.com/MITgcm/MITgcm by the 2.2.1-MITgcm-0.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!