module m_read_coriolis private :: errcond, getdimms,getvarstat contains !Reads Coriolis data from a NetCDF file subroutine read_coriolis(filename,var,pres,longitude,latitude,juld,floatid, & errorflag,l_invalid,varname,undef) use netcdf implicit none ! Dummy variables character(len=*), intent(in) :: filename,varname real, pointer, dimension(:) :: longitude,latitude,juld real, pointer, dimension(:,:) :: var,pres logical, pointer, dimension(:,:) :: errorflag character(len=8),dimension(:) , pointer :: floatid logical ,dimension(:) , pointer :: l_invalid real, intent(in) :: undef character(len=16) :: tmp character(len=80), dimension(:), allocatable :: depthname,varname2 ! Local variables character(len=1), allocatable, dimension(:,:,:) :: & sparam character(len=1), allocatable, dimension(:,:) :: & qcvar, qcpres,qcdeph,dstate,dmode character(len=1), allocatable, dimension(:) :: & qcpos,qcjuld logical , allocatable, dimension(:) :: & l_deph, l_adj integer, dimension(nf90_max_dims) :: vdims ! Variable dimensions - not kept integer :: vndim ! number of variable dims - not kept integer :: dstatvarid, dmodevarid,sparamvarid ! Dimension id - not kept integer :: dimid ! Dimension lengths integer :: profdimlen,zlevdimlen,s8dimlen,s4dimlen, nparamdimlen, s16dimlen ! Variable ids - Kept ! integer :: & lonvarid,latvarid, presvarid, qcvarid, pltfvarid, & qcpresvarid, varid,juldvarid,qcposvarid,qcjuldvarid,dephvarid,qcdephvarid integer :: ncid integer :: i,stat,j,k,k2 logical :: ex real :: updiff,lwdiff,diff real :: add_offset,scale_factor, fillvalue logical :: lpres, ldeph,l_present integer :: ncresult real :: fillvalue_pres,fillvalue_var ! Check if variable type is correct if (varname/='PSAL' .and. varname/='TEMP') then print *,'Variable name '//varname//'is unknown' print *,'valid names are PSAL or TEMP' stop '(read_coriolis)' end if ! Check for existence of netcdf file inquire(file=trim(filename),exist=ex) if (.not.ex) then print *,'Can not find file '//trim(filename) stop '(read_coriolis)' end if print *,'reading Coriolis fields from '//trim(filename) stat = nf90_open(trim(filename),nf90_nowrite,ncid) if (stat /= nf90_noerr) then call ncerr( stat) end if ! Check if sought variable is present l_present=probe_var(ncid,varname) !print *,'reading Coriolis fields from '//trim(filename) ! Dimension IDs and lengths (assuming the names are right) call getdimms(ncid,'N_PROF' ,dimid,profdimlen ) call getdimms(ncid,'N_LEVELS',dimid,zlevdimlen ) call getdimms(ncid,'N_PARAM' ,dimid,nparamdimlen) call getdimms(ncid,'STRING16',dimid,s16dimlen ) call getdimms(ncid,'STRING8' ,dimid,s8dimlen ) call getdimms(ncid,'STRING4' ,dimid,s4dimlen ) print *,'# Profiles , # levels ',profdimlen, zlevdimlen !print *,profdimid, zlevdimid ! Better safe than sorry ... if (profdimlen<1 .or. zlevdimlen<1) then print *,'Dimensions are zero for this file ... ' print *,'Exiting..' stop '(read_coriolis)' end if ! Get data state indicator, mode and station parameters call getvarstat(ncid,'DATA_STATE_INDICATOR',dstatvarid ,vndim ,vdims) call getvarstat(ncid,'DATA_MODE' ,dmodevarid ,vndim ,vdims) call getvarstat(ncid,'STATION_PARAMETERS' ,sparamvarid,vndim ,vdims) allocate(sparam (s16dimlen,nparamdimlen,profdimlen)) allocate(dstate (s4dimlen,profdimlen)) allocate(dmode (1 ,profdimlen)) allocate(l_deph (profdimlen)) allocate(l_adj (profdimlen)) allocate(l_invalid(profdimlen)) l_deph=.false. l_adj =.false. l_invalid =.false. do j=1,profdimlen ! Data state indicator call ncerr( nf90_get_var(ncid, dstatvarid, dstate(:,j), & start=(/1,j/),count=(/s4dimlen,1/))) !print *,dstate(:,j) ! Data mode (R)ealtime, (D)elayed or (A)djusted call ncerr( nf90_get_var(ncid, dmodevarid, dmode (:,j), & start=(/1,j/),count=(/1,1/))) !print *,dmode(1,j) ! Station parameters call ncerr( nf90_get_var(ncid, sparamvarid, sparam (:,:,j), & start=(/1,1,j/),count=(/s16dimlen,nparamdimlen,1/))) !print *,sparam(:,:,j) do k=1,nparamdimlen ! Check what variable to use for depth ! KAL - TODO - sparam is a 3D - var => l_deph ! should be as well tmp='' do k2=1,16 tmp(k2:k2)=sparam(k2,k,j) end do !print *,sparam(:,k,j) if (trim(tmp)=='DEPH') then l_deph(j)=.true. elseif (trim(tmp)=='PRES') then l_deph(j)=.false. end if !! Check that variable we want is present !tmp='' !do k2=1,16 ! tmp(k2:k2)=sparam(k2,k,j) !end do !if (trim(tmp)==trim(varname)) then ! l_present(j)=.true. !end if end do ! Check if adjusted variables should be used tmp='' do k2=1,4 tmp(k2:k2)=dstate(k2,j) end do if (dmode(1,j)(1:1)=='A') then l_adj(j)=.true. elseif (dmode(1,j)(1:1)=='R') then if ( trim(tmp) == "2C" .or. trim(tmp)=="2B+") then l_adj(j)=.true. end if end if ! Check if profile is invalid if ( trim(tmp) == "0A") then l_invalid(j)=.true. end if print *,tmp,l_deph(j),l_adj(j),l_invalid(j) end do ! platform, lon, lat and julian day call getvarstat(ncid,'PLATFORM_NUMBER',pltfvarid,vndim ,vdims) call getvarstat(ncid,'LONGITUDE' ,lonvarid ,vndim ,vdims) call getvarstat(ncid,'LATITUDE' ,latvarid ,vndim ,vdims) call getvarstat(ncid,'JULD' ,juldvarid,vndim ,vdims) allocate(floatid (profdimlen)) allocate(longitude(profdimlen)) allocate(latitude (profdimlen)) allocate(juld (profdimlen)) call ncerr( nf90_get_var(ncid, juldvarid, juld ) ) call ncerr( nf90_get_var(ncid, lonvarid, longitude) ) call ncerr( nf90_get_var(ncid, latvarid, latitude ) ) do j=1,profdimlen call ncerr( nf90_get_var(ncid, pltfvarid ,floatid(j)(:), & start=(/1,j/),count=(/s8dimlen,1/)) ) !print *,floatid(j) end do ! QC flags for juld and position call getvarstat(ncid,'POSITION_QC',qcposvarid ,vndim ,vdims) call getvarstat(ncid,'JULD_QC' ,qcjuldvarid,vndim ,vdims) allocate(qcpos (profdimlen)) allocate(qcjuld (profdimlen)) call ncerr( nf90_get_var(ncid, qcposvarid , qcpos, & start=(/1/),count=(/profdimlen/)) ) call ncerr( nf90_get_var(ncid, qcjuldvarid , qcjuld, & start=(/1/),count=(/profdimlen/)) ) print *,'qcpos :',qcpos print *,'qcjuld:',qcjuld ! Adjust invalid flag based on position and juld quality flags do j=1,profdimlen if(qcpos (j) /= '1' .and. qcpos(j) /='0') then l_invalid(j)=.true. end if if(qcjuld(j) /= '1' .and. qcjuld(j) /='0') then l_invalid(j)=.true. end if end do print *,'hei1' ! Set appropriate depth variable name ! KAL - TODO : If l_deph is 3d, then depthname must be allocate(depthname(profdimlen)) allocate(varname2 (profdimlen)) do j=1,profdimlen if (l_deph(j) ) then depthname(j)='DEPH' else depthname(j)='PRES' end if if (l_adj(j)) then depthname(j)=depthname(j)//'_ADJUSTED' varname2 (j)=trim(varname)//'_ADJUSTED' else varname2 (j)=trim(varname) end if print *,j,trim(depthname(j)) ! Security check - we assume that either DEPTH or PRES is present.... ! If this becomes a problem we must gracefully handle it.... !if (j>1) then ! if (depthname(j)(1:4)/=depthname(1)(1:4)) then ! print *,'Error - both PRES and DEPH are used in netcdf file!' ! stop ! end if !end if ! KAL - the routine should be able to handle this end do print *,'hei2' ! Read profile depths and values allocate(var (zlevdimlen,profdimlen)) allocate(pres (zlevdimlen,profdimlen)) allocate(qcvar (zlevdimlen,profdimlen)) allocate(qcpres(zlevdimlen,profdimlen)) var=undef do j=1,profdimlen print *,j,trim(varname2(j)),' ',trim(depthname(j)) ! hmm.. overkill for var stats but... call getvarstat(ncid,trim(depthname(j)) ,presvarid ,vndim,vdims ) call getvarstat(ncid,trim(depthname(j))//'_QC',qcpresvarid ,vndim,vdims ) if (l_present) then call getvarstat(ncid,trim(varname2 (j)) ,varid ,vndim,vdims ) call getvarstat(ncid,trim(varname2 (j))//'_QC',qcvarid,vndim,vdims ) end if call ncerr( nf90_get_var(ncid,presvarid , pres(:,j), & start=(/1,j/),count=(/zlevdimlen,1/))) if (l_present) then call ncerr( nf90_get_var(ncid, varid , var(:,j), & start=(/1,j/),count=(/zlevdimlen,1/))) ! get fillvalue attribute of pressure and variable name ncresult=nf90_get_att(ncid, varid , '_FillValue', fillvalue_var) where (var (:,j)==fillvalue_var ) var (:,j)=undef end if ! set to undef (input) where variable and pressure equals fillvalues ncresult=nf90_get_att(ncid, presvarid, '_FillValue', fillvalue_pres) where (pres(:,j)==fillvalue_pres) pres(:,j)=undef print *,var(:,j) print *,'hei0' call ncerr( nf90_get_var(ncid, qcpresvarid, qcpres(:,j), & start=(/1,j/),count=(/zlevdimlen,1/))) print *,'hei1' if (l_present) then call ncerr( nf90_get_var(ncid, qcvarid , qcvar(:,j), & start=(/1,j/),count=(/zlevdimlen,1/))) end if print *,'hei2' print *,qcvar(:,j) end do print *,'hei' ! Go through data and discard poor observations (errorflag=true) allocate(errorflag(zlevdimlen,profdimlen)) errorflag=.false. do j=1,profdimlen do i=1,zlevdimlen if (l_present) then if(qcvar (i,j) /= '1' .and. qcvar(i,j) /='0') errorflag(i,j)=.true. end if if(qcpres(i,j) /= '1' .and. qcpres(i,j) /='0') errorflag(i,j)=.true. end do end do do j=1,profdimlen if (l_invalid(j)) then errorflag(:,j)=.true. end if end do end subroutine read_coriolis !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !The following are auxillary routines ! Give netcdf error message subroutine ncerr(error) use netcdf implicit none integer, intent(in) :: error if (error/=nf90_noerr) then print *,'read_coriolis: '//nf90_strerror(error) stop end if end subroutine ncerr subroutine getdimms(ncid,dname,dmid,dmlen) use netcdf implicit none integer, intent(in) :: ncid character(len=*), intent(in) :: dname integer, intent(out) :: dmid,dmlen call ncerr( nf90_inq_dimid(ncid,dname,dmid) ) call ncerr( nf90_inquire_dimension(ncid,dmid,len=dmlen) ) end subroutine logical function probe_var(ncid,varname) use netcdf implicit none character(len=*), intent(in) :: varname integer , intent(in) :: ncid integer :: test, varid test=nf90_inq_varid(ncid,varname,varid) probe_var = test==nf90_noerr end function probe_var subroutine getvarstat(ncid,varname,varid,var_ndim,var_dimids) use netcdf implicit none integer, intent(in) :: ncid character(len=*), intent(in) :: varname integer, intent(out) :: varid,var_ndim integer, dimension(nf90_max_dims),intent(out) :: var_dimids call ncerr( nf90_inq_varid(ncid,varname,varid) ) call ncerr( nf90_inquire_variable(ncid,varid,ndims=var_ndim,dimids=var_dimids) ) end subroutine getvarstat end module m_read_coriolis