MODULE CloudTab ! This module produces cloud properties from tables of gridded data. ! The user provides a netCDF file that contains the data. ! The netCDF file must contain these anciliary data: ! - the nodes of the grid for each dimension of the data ! - strings that define the type of interpolation (e.g. linear) to use ! for each dimension (it is assumed that the table was built to achieve ! required accuracy with a certain type of interpolation) ! - a flag for each grid that indicates whether the grid is uniformly ! spaced, which speeds computation ! The data may represent scattering and absorption coefficients, asymmetry ! factor, etc. There may be more than one type of data in the file. ! Each type of data has a distinct netCDF variable name. ! The data must be on a grid of 4 dimensions, where the dimensions ! respectively cover: primary size parameter (may be effective diameter, ! mass-diameter, etc.), secondary size parameter (may be variance or width of ! size ditribution), temperature (K), and spectral point (may be wavenumber, ! frequency, etc.). ! With the exception of spectral point, the dimension may be 1, whereby the ! property is assumed not to vary in that dimension. ! The subroutine getCloudTab provides cloud property data for all ! spectral points interpolated to the specified ! values of size parameters and temperature. use ncdf_module implicit none private public:: cloudTabInit,getCloudTab,cloudTabDestroy,ReadMie ! Global private data integer, parameter :: mxiv=12 ! length of variable types integer, parameter :: mxty=15 ! max # number of data types character(len=mxiv), dimension(mxty) :: dataType integer, parameter :: mxim=12 ! length of interp types integer, parameter :: mxit=2 ! number of interp types character(len=mxim), dimension(mxit) :: validInterpTypes=(/'none ','linear'/) character(len=mxim), dimension(mxty) :: tySzP,tySzS,tyT integer, dimension(mxty) :: unigSzP,unigSzS,unigT real, dimension(mxty) :: stepSzP,stepSzS,stepT logical, dimension(mxty) :: initialized=.false. logical :: first_init=.true. logical :: outBndsExit=.false. ! To allow cloud data tabulation to use different grid for each ! data type (see EmisTab): ! * move myspecPtsGrd into structure Variables ! * add function to interpolate from spectral point grid to input ! spectral position, such that parent program does not need to ! know spectral point grid, the grid is ! * remove nSpcPts and spectPtsGrd from cloudTabInit calling arguments real, dimension(:), allocatable :: myspectPtsGrd integer :: nSPsave=0 type Variables real, dimension(:,:,:,:), pointer :: datArr real, dimension(:,:,:), pointer :: tmp3 real, dimension(:,:), pointer :: tmp2 real, dimension(:), pointer :: grSzP,grSzS,grT end type Variables type(Variables), dimension(mxty) :: allVar CONTAINS SUBROUTINE cloudTabInit(fname,dataType_in,outBndsExit_in,nSpcPts,spectPtsGrd) ! Initialize access to tabulated cloud data. Data table is assumed ! to represent a spectral point set and the user may verify that ! the set is as expected by returning spectPtsGrd and nSpcPts. character(len=*), intent(in) :: fname character(len=*), dimension(:), intent(in) :: dataType_in logical, intent(in), optional :: outBndsExit_in integer, intent(inout), optional :: nSpcPts real, dimension(:), pointer, optional :: spectPtsGrd !-- Local variables integer :: indx integer :: ixkabs,ixkscat,ixgasym integer :: mynSpcPts integer, dimension(6) :: mydimid integer :: tmpvarid integer :: ncid,ncStatus,istatus integer :: neSzP,neSzS,neT character(len=60) :: htName if (size(dataType_in) > mxty) then print *,'err[CloudTab::cloudTabInit]: ' print *,' Trying to load more variables than the global table can hold' print *,' size(dataType_in) > mxty:',size(dataType_in),mxty call errorHalt(1) endif !-- Read file call openNcdfFile(ncid, fname) do indx=1,size(dataType_in) call cloudTabDestroy(indx) ! Start fresh with memory allocation !--The "indx" in the statement below indicates that !--the field is the indx-th variable in the file ncStatus = nf_inq_varid(ncid,trim(dataType_in(indx)),tmpvarid) if (ncStatus /= nf_noerr) call handle_nc_err(ncStatus,'cloudTabInit', & 'File does not contain dataType '//trim(dataType_in(indx))) dataType(indx)=dataType_in(indx) ! Hold in global data for reference ! Dimensions depend ultimately on dataType; file could contain > 1 datatype ncStatus = nf_inq_vardimid(ncid,tmpvarid,mydimid) if (ncStatus /= nf_noerr) call handle_nc_err(ncStatus,'cloudTabInit', & 'Failed to get dimension IDs') ncStatus = nf_inq_dimlen(ncid,mydimid(1),neSzP) if (ncStatus /= nf_noerr) call handle_nc_err(ncStatus,'cloudTabInit', & 'Failed to get dimension of sizP') ncStatus = nf_inq_dimlen(ncid,mydimid(2),neSzS) if (ncStatus /= nf_noerr) call handle_nc_err(ncStatus,'cloudTabInit', & 'Failed to get dimension of sizS') ncStatus = nf_inq_dimlen(ncid,mydimid(3),neT) if (ncStatus /= nf_noerr) call handle_nc_err(ncStatus,'cloudTabInit', & 'Failed to get dimension of T') ncStatus = nf_inq_dimlen(ncid,mydimid(4),mynSpcPts) if (ncStatus /= nf_noerr) call handle_nc_err(ncStatus,'cloudTabInit', & 'Failed to get nSpcPts') if (present(nSpcPts)) nSpcPts=mynSpcPts allocate(allVar(indx)%grSzP(neSzP), & allVar(indx)%grSzS(neSzS), & allVar(indx)%grT(neT), & allVar(indx)%datArr(neSzP,neSzS,neT,mynSpcPts)) allocate(allVar(indx)%tmp2(neT,mynSpcPts), & allVar(indx)%tmp3(neSzS,neT,mynSpcPts), stat=istatus) if (istatus /=0 ) then print *,'err[CloudTab::cloudTabInit]: ', & 'Not enough memory for tmp2,tmp3; istatus:',istatus call errorHalt(1) endif ncStatus = nf_get_var_real(ncid,tmpvarid,allVar(indx)%datArr) if (ncStatus /= nf_noerr) then print *,'err[CloudTab::cloudTabInit]: '// & '(get_var_real:'//trim(dataType(indx))// & ') '//trim(nf_strerror(ncStatus)) call errorHalt(1) endif ixkabs =index(dataType_in(indx),'_kabs') ixkscat=index(dataType_in(indx),'_kscat') ixgasym=index(dataType_in(indx),'_gasym') htName=dataType_in(indx) if (ixkabs /= 0) htName=dataType_in(indx)(1:ixkabs-1) if (ixkscat /= 0) htName=dataType_in(indx)(1:ixkscat-1) if (ixgasym /= 0) htName=dataType_in(indx)(1:ixgasym-1) call readNcdfAttr(ncid,attrName='grid'//trim(htName)//'SizP', & attr=allVar(indx)%grSzP) call readNcdfAttr(ncid,attrName='grid'//trim(htName)//'SizS', & attr=allVar(indx)%grSzS) call readNcdfAttr(ncid,attrName='grid'//trim(htName)//'Temp', & attr=allVar(indx)%grT) if (.not. allocated(mySpectPtsGrd)) then allocate(mySpectPtsGrd(mynSpcPts)) call readNcdfAttr(ncid,attrName='spectralPoints',attr=mySpectPtsGrd) nSPsave=mynSpcPts else if (mynSpcPts /= nSPsave) then print *,'err[CloudTab::cloudTabInit]: '// & 'Inconsistent number of spectral points' print *,'mynSpcPts /= nSPsave:',mynSpcPts,nSPsave call errorHalt(1) endif endif call readNcdfAttr(ncid,attrName='interp'//trim(htName)//'SizP', & attr=tySzP(indx)) call readNcdfAttr(ncid,attrName='interp'//trim(htName)//'SizS', & attr=tySzS(indx)) call readNcdfAttr(ncid,attrName='interp'//trim(htName)//'Temp', & attr=tyT(indx)) call interpTypeCheck('SzP',neSzP,tySzP(indx)) call interpTypeCheck('SzS',neSzS,tySzS(indx)) call interpTypeCheck('Temp',neT, tyT(indx)) call readNcdfAttr(ncid,attrName='unifGrid'//trim(htName)//'SizP', & attr=unigSzP(indx)) call readNcdfAttr(ncid,attrName='unifGrid'//trim(htName)//'SizS', & attr=unigSzS(indx)) call readNcdfAttr(ncid,attrName='unifGrid'//trim(htName)//'Temp', & attr=unigT(indx)) stepSzP(indx)=stepUniform('SzP',unigSzP(indx),allVar(indx)%grSzP) stepSzS(indx)=stepUniform('SzS',unigSzS(indx),allVar(indx)%grSzS) stepT(indx) =stepUniform('Temp',unigT(indx),allVar(indx)%grT) if (present(spectPtsGrd)) then if (associated(spectPtsGrd)) deallocate(spectPtsGrd) allocate (spectPtsGrd(mynSpcPts)) spectPtsGrd=mySpectPtsGrd endif initialized(indx)=.true. enddo call closeNcdfFile(ncid) if(present(outBndsExit_in))outBndsExit=outBndsExit_in return END SUBROUTINE cloudTabInit !-------------------------------------------------------------------------- SUBROUTINE handle_nc_err(ncStatus,subproName,msg) integer, intent(in) :: ncStatus character (len=*), intent(in) :: subproName character (len=*), intent(in), optional :: msg character (len=256) :: lmsg if (present(msg)) then if (len(msg) > 256) print*,'message is truncated' lmsg=trim(msg) else lmsg=':::' endif print *, 'err[ncdf_module::handle_nc_err(',trim(subproName),')]: ',& ' ('//trim(msg)//') '//trim(nf_strerror(ncStatus)) call errorHalt(1) END SUBROUTINE handle_nc_err !-------------------------------------------------------------------------- SUBROUTINE interpTypeCheck(name,nelem,type) character(len=*), intent(in) :: name integer, intent(in) :: nelem character(len=*), intent(inout) :: type integer :: itype character(len=mxim) :: tmptype integer :: i itype=0 if (nelem == 1) then type='none' ! Override interpolation type when only one node print *,'msg[CloudTab::interpTypeCheck]: ' print *,name,' grid has only 1 node; will be assumed constant' endif ! Get integer index of interpolation type do i=1,size(validInterpTypes) tmptype=validInterpTypes(i) if (trim(adjustl(type)) == trim(adjustl(tmptype))) then itype=i exit endif enddo if (itype == 0) then print *,'err[CloudTab::interpTypeCheck]: ' print *,'For '//trim(name)//', interp type "'//trim(adjustl(type))// & '" not in valid types:' print *,validInterpTypes call errorHalt(1) endif return END SUBROUTINE interpTypeCheck !-------------------------------------------------------------------------- FUNCTION stepUniform(name,unig,grid) character(len=*), intent(in) :: name integer, intent(in) :: unig real, dimension(:), intent(in) :: grid real :: stepUniform real, parameter :: deplim=0.001 integer :: i,ngrid real :: step,dep step=0. if (unig /= 1) return ngrid=size(grid) if (ngrid == 1) return step=(grid(ngrid)-grid(1))/float(ngrid-1) do i=2,size(grid) dep=abs(grid(i)-grid(1)-step*float(i-1))/step if (dep > deplim) then print*,grid print *,'err[CloudTab::stepUniform]: ', & name,' grid claimed to be uniform but was not; i,dep:',i,dep call errorHalt(1) endif enddo stepUniform=step return END FUNCTION stepUniform !-------------------------------------------------------------------------- SUBROUTINE getCloudTab(indx,sizP,sizS,temp, & prop,dPropdSzP,dPropdSzS,dPropdTemp) integer, intent(in) :: indx real, intent(in) :: sizP,sizS,temp real, dimension(:), intent(inout) :: prop real, dimension(:), intent(inout), optional :: dPropdSzP,dPropdSzS,dPropdTemp integer :: iSzP1,iSzP2,iSzS1,iSzS2,iT1,iT2 integer :: iSzS,iT if (indx > mxty) then print *,'err[CloudTab::getCloudTab]: ', & 'Trying to access field beyond the rank of the global table' call errorHalt(1) endif if (.not. initialized(indx)) then print *,'err[CloudTab::getCloudTab]: ', & 'Must call cloudTabInit before first call to getCloudTab' call errorHalt(1) endif if (size(prop) /= size(allVar(indx)%datArr,dim=4)) then print *,'err[CloudTab::getCloudTab]: ', & 'vector length mismatch; size(prop) /= nSpcPts:', & size(prop) /= size(allVar(indx)%datArr,dim=4) call errorHalt(1) endif ! Find bounding grid indices call boundGrid('SzP',allVar(indx)%grSzP,unigSzP(indx),stepSzP(indx), & sizP,iSzP1,iSzP2) call boundGrid('SzS',allVar(indx)%grSzS,unigSzS(indx),stepSzS(indx), & sizS,iSzS1,iSzS2) call boundGrid('Temp',allVar(indx)%grT, unigT(indx), stepT(indx) , & temp, iT1, iT2) do iT=iT1,iT2 do iSzS=iSzS1,iSzS2 call interpolate(tySzP(indx),allVar(indx)%grSzP,iSzP1,iSzP2, & sizP,allVar(indx)%datArr(:,iSzS,iT,:),allVar(indx)%tmp3(iSzS,iT,:)) enddo call interpolate(tySzS(indx),allVar(indx)%grSzS,iSzS1,iSzS2, & sizS,allVar(indx)%tmp3(:,iT,:),allVar(indx)%tmp2(iT,:)) enddo call interpolate(tyT(indx),allVar(indx)%grT,iT1,iT2,temp,allVar(indx)%tmp2(:,:),prop) return END SUBROUTINE getCloudTab !-------------------------------------------------------------------------- SUBROUTINE boundGrid(name,grid,unig,step,val,i1,i2) ! Find the indices of the grid that bound the given value character(len=*), intent(in) :: name real, dimension(:), intent(in) :: grid integer, intent(in) :: unig real, intent(in) :: step real, intent(in) :: val integer, intent(inout) :: i1 integer, intent(inout) :: i2 integer :: i,ngrid ngrid=size(grid) ! Grid has size 1 if (ngrid == 1) then i1=1 i2=1 ! Grid is uniformly spaced; saves time without loop over if block elseif (unig == 1) then i1=floor((val-grid(1))/step)+1 i2=i1+1 if (i1 < 1) then if (outBndsExit) then print *,'err[CloudTab::boundGrid]: ', & name,' value',val,' is out of uniform range:',grid(1),grid(ngrid) call errorHalt(1) elseif (val < 0.) then print *,'err[CloudTab::boundGrid]: ', & name,' value',val,' is unphysical negative; min grid=',grid(1) call errorHalt(1) else i1=1 i2=2 endif endif if (i2 > ngrid) then if (val <= grid(ngrid)) then ! Accommodate roundoff i2=ngrid i1=i2-1 elseif (outBndsExit) then print *,'err[CloudTab::boundGrid]: ', & name,' value',val,' is out of uniform range:',grid(1),grid(ngrid) call errorHalt(1) else i2=ngrid i1=i2-1 endif endif else ! Grid is not uniformly spaced if ((val < grid(1) .or. val > grid(ngrid)) .and. outBndsExit) then print *,'err[CloudTab::boundGrid]: ', & name,' value',val,' is out of range:',grid(1),grid(ngrid) call errorHalt(1) elseif (val < 0.) then print *,'err[CloudTab::boundGrid]: ', & name,' value',val,' is unphysical negative; min grid=',grid(1) call errorHalt(1) endif do i=2,ngrid if (grid(i) >= val) exit enddo i2=min(i,ngrid) i1=i2-1 endif return END SUBROUTINE boundGrid !-------------------------------------------------------------------------- SUBROUTINE interpolate(type,grid,i1,i2,xin,arr,arrIntrp,arrDer) character(len=*), intent(in) :: type real, dimension(:), intent(in) :: grid integer, intent(in) :: i1 integer, intent(in) :: i2 real, intent(in) :: xin real, dimension(:,:), intent(in) :: arr real, dimension(:), intent(inout) :: arrIntrp real, dimension(:), intent(inout), optional :: arrDer select case(type) case ('none') arrIntrp=arr(i1,:) if(present(arrDer))arrDer=0. case ('linear') arrIntrp=(arr(i1,:)*(grid(i2)-xin)+arr(i2,:)*(xin-grid(i1)))/ & (grid(i2)-grid(i1)) if(present(arrDer))arrDer =(arr(i2,:)-arr(i1,:))/(grid(i2)-grid(i1)) case default print *,'err[CloudTab::interpolate]: ', & 'Unsupported interpolation method: ',type call errorHalt(1) end select return END SUBROUTINE interpolate !-------------------------------------------------------------------------- SUBROUTINE cloudTabDestroy(indx) integer, intent(in) :: indx integer :: ix ! Deallocate memory that can be allocated by cloudTabInit, except ! for optional arguments of cloudTabInit that are handled there if(first_init)then do ix=1,mxty nullify(allVar(ix)%grSzP,allVar(ix)%grSzS,allVar(ix)%grT) nullify(allVar(ix)%datArr) nullify(allVar(ix)%tmp2,allVar(ix)%tmp3) enddo first_init=.false. endif if(indx == 0)then ! destroy all arrays do ix=1,mxty if(associated(allVar(ix)%grSzP))deallocate(allVar(ix)%grSzP) if(associated(allVar(ix)%grSzS))deallocate(allVar(ix)%grSzS) if(associated(allVar(ix)%grT))deallocate(allVar(ix)%grT) if(associated(allVar(ix)%datArr))deallocate(allVar(ix)%datArr) if(associated(allVar(ix)%tmp2))deallocate(allVar(ix)%tmp2) if(associated(allVar(ix)%tmp3))deallocate(allVar(ix)%tmp3) enddo if (allocated(mySpectPtsGrd)) deallocate(mySpectPtsGrd) else ! destroy arrays pertaining to a particular file if (indx > mxty) then print *,'err[CloudTab::cloudTabDestroy]: ', & 'Trying to destroy field beyond the rank of the global table' call errorHalt(1) endif if(associated(allVar(indx)%grSzP))deallocate(allVar(indx)%grSzP,allVar(indx)%grSzS, & allVar(indx)%grT) if(associated(allVar(indx)%datArr))deallocate(allVar(indx)%datArr) if(associated(allVar(indx)%tmp2))deallocate(allVar(indx)%tmp2,allVar(indx)%tmp3) endif return END SUBROUTINE cloudTabDestroy SUBROUTINE ReadMie(FrqMie,QMie) real, dimension(:), pointer :: FrqMie real, dimension(:,:), pointer :: QMie !---Local variables character(len=100) :: filename real, dimension(:), allocatable :: FrqMie_loc real, dimension(:,:), allocatable :: QMie_loc filename='/storm/cris/DataBases/MIE/cris.new.mier01v10.wat' filename=trim(filename) print*,filename stop return END SUBROUTINE ReadMie END MODULE CloudTab