psmile_read_byid_dble.F90

Go to the documentation of this file.
00001 !-----------------------------------------------------------------------
00002 ! Copyright 2006-2010, SGI Germany, Munich, Germany.
00003 ! All rights reserved. Use is subject to OASIS4 license terms.
00004 !-----------------------------------------------------------------------
00005 !BOP
00006 !
00007 ! !ROUTINE:  psmile_read_byid_dble
00008 !
00009 ! !INTERFACE:
00010       subroutine psmile_read_byid_dble(id_varid,id_taskid,dd_a,ju_day,ju_sec &
00011                                       ,ju_dayb,ju_secb,timeop,ierror)
00012 !
00013 ! !USES:
00014 !
00015       use prism_constants
00016       use psmile, dummy_interface=> psmile_read_byid_dble
00017       use prism_calendar
00018 !
00019       implicit none
00020 !
00021 ! !INPUT PARAMETERS:
00022 !
00023       integer,intent(in)::id_varid
00024       integer,intent(in)::id_taskid
00025       double precision,target, Intent (inout) :: dd_a(*)
00026       double precision, Intent (In) :: ju_day, ju_sec
00027       double precision, Intent (In) :: ju_dayb(2), ju_secb(2)
00028 !
00029 ! !OUTPUT PARAMETERS:
00030 !
00031       logical,intent(out):: timeop
00032       integer,intent(out):: ierror
00033 !
00034 ! !LOCAL VARIABLES:
00035 !
00036       integer                           :: ierrp(2)
00037       integer               :: il_grid_id,il_grid_type
00038       integer               :: il_method_id,il_varid
00039       integer                           :: il_blockid
00040       integer                           :: iicomp_id
00041       logical                           :: is_block
00042       integer               :: nvcomp,il_ndim,offset,vectorfield
00043       integer                           :: fullsize
00044       integer               :: il_i,lenb,len,ii,jj,il_loc 
00045                                           ,il_time,il_times(2),ibl,ibu
00046       integer               :: myvarshape(2,ndim_3d+2)
00047       integer               :: mygrdshape(2,ndim_3d+2)
00048       double precision,allocatable  :: afield(:)
00049       double precision                  :: tol,fill
00050       double precision          :: dl_time,dl_ju_day,dl_ju_sec
00051       double precision          :: dl_timeb(2),w,days,sec
00052       double precision,pointer          :: ju_sec_last,ju_day_last 
00053                                          , time_last
00054       double precision,pointer      :: pl_times(:)
00055       double precision              :: aone = 1.d0
00056       logical               :: is_bundle,lmatch,use_domain
00057       logical                           :: lower_in,upper_in
00058       logical,pointer                   :: llast
00059 #ifdef __PSMILE_WITH_IO
00060       Type(IO_cache),pointer            :: pl_cache
00061 #endif
00062       Type(PRISM_Time_struct)           :: cur_date
00063 !
00064 ! !DESCRIPTION:
00065 !  Routine reads a dble field variable.
00066 !  It can handle scalar fields and vector fields as well as bundles of
00067 !  the two types of fields.
00068 !  Last but not least it reads subblocks of the above in case of a
00069 !  GME grid.
00070 !  If the current time stamp can not be read by matching the time offset
00071 !  with respect to the current job start date a new file will be opened
00072 !  transparently if the user has given a base name for the input file.
00073 !  The search is based on the extension of the file basename by the job 
00074 !  start dates defining the origin to the time axes contained in those
00075 !  files under consideration.  
00076 !  For a completely specified filename no search is performed.
00077 !  An unsuccessful match of a
00078 !  time stamp leads to a termination of the model component code.
00079 !  
00080 !  The routine honors date bounds for matching time stamps (time lags).
00081 !
00082 !  In case of bundles of vector fields transpose routines are called
00083 !  which transpose the indices of vector components with bundles.
00084 !
00085 ! !REVISION HISTORY:
00086 !
00087 !   Date      Programmer   Description
00088 ! ----------  ----------   -----------
00089 !   21.12.03  R. Vogeslang created
00090 !   29.01.04  R. Vogelsang Time interpolation using date bounds.
00091 !   23.03.04  R. Vogelsang Added the taskid to the interface and related code.
00092 !   12.10.04  R. Vogelsang Values = fill value are set to psmile_dundef.
00093 !   2.10.05   R. Vogelsang Added Gaussian reduced grids
00094 !EOP
00095 !----------------------------------------------------------------------
00096 !
00097       character(len=len_cvs_string),save :: mcvs = 
00098 '$Id: psmile_read_byid_dble.F90 2920 2011-01-27 10:55:28Z coquart $'
00099 
00100 !
00101       ierror=0
00102       timeop=.false.
00103 #ifdef __PSMILE_WITH_IO
00104 #ifdef VERBOSE
00105       print*,trim(ch_id),' : PSMILe_read_byid_dble: start'
00106       call psmile_flushstd
00107 #endif
00108 !
00109 !     Return from this routine if the field is not subject to file I/O
00110 !     or the I/O action is not MPP_RDONLY (read).
00111 !
00112       if (.not.associated(Fields(id_varid)%io_chan_infos)) return
00113 !
00114 !     IO infos are pointing to io_chan_infos(taskid)
00115 !
00116       if(id_taskid .le. size(Fields(id_varid)%io_task_lookup)) then
00117 
00118         il_i=Fields(id_varid)%io_task_lookup(id_taskid)
00119 
00120       else
00121 
00122         ierror = PRISM_Error_IO_Meta
00123         ierrp (1) = id_taskid
00124         call psmile_error ( ierror &
00125                           , 'Task id out of range! ', &
00126                             ierrp, 1, __FILE__, __LINE__ )
00127         return
00128 
00129       endif
00130 
00131       if (il_i.gt.0) then
00132 
00133          Fields(id_varid)%io_infos =>  Fields(id_varid)%io_chan_infos(il_i)
00134 
00135       else
00136         ierror = PRISM_Error_IO_Meta
00137         ierrp (1) = id_taskid
00138         call psmile_error ( ierror &
00139                           , 'Negative task id! ', &
00140                             ierrp, 1, __FILE__, __LINE__ )
00141         return
00142       endif
00143 
00144       if (Fields(id_varid)%io_infos%action .ne. MPP_RDONLY ) return
00145 !
00146 !     Set the pelist for the current active component
00147 !
00148       iicomp_id=Fields(id_varid)%comp_id
00149       call mpp_set_current_pelist(IO_Comps_infos(iicomp_id)%pelist)
00150 
00151 
00152 
00153       il_varid=id_varid
00154       il_method_id=Fields(id_varid)%method_id
00155       il_grid_id=Methods(il_method_id)%grid_id
00156       il_grid_type=Grids(il_grid_id)%grid_type
00157       il_ndim=Grids(il_grid_id)%n_dim
00158       il_blockid=0
00159       is_block=.false.
00160 
00161       use_domain=.true.
00162       if(Fields(id_varid)%io_infos%threading .eq. MPP_MULTI .and. &
00163          Fields(id_varid)%io_infos%fileset .eq. MPP_MULTI )use_domain=.FALSE.
00164 
00165       dl_ju_day=ju_day
00166       dl_ju_sec=ju_sec
00167 !
00168 !     Time offset with respect to the date of the job start contained
00169 !     in the current file.
00170 !
00171 #ifdef VERBOSE
00172       print*,trim(ch_id),' : PSMILe_read_byid_dble: called: ju_start_day= ' &
00173             ,Fields(id_varid)%io_infos%ju_start_day &
00174             ,' ju_start_sec= ',Fields(id_varid)%io_infos%ju_start_sec
00175       print*,trim(ch_id),' : PSMILe_read_byid_dble: id_varid=',id_varid &
00176                         ,' id_taskid=',id_taskid,' il_i= ', il_i
00177       call psmile_flushstd
00178 #endif
00179       dl_time=86400.d0 *(dl_ju_day - Fields(id_varid)%io_infos%ju_start_day) &
00180              +          (dl_ju_sec - Fields(id_varid)%io_infos%ju_start_sec)
00181 #ifdef VERBOSE
00182       print*,trim(ch_id),' : PSMILe_read_byid_dble: called: ju_day= ' &
00183             ,ju_day ,' ju_sec= ',ju_sec,' ju_dayb=',ju_dayb,' ju_secb=' &
00184             ,ju_secb,' offset=',dl_time
00185       call psmile_flushstd
00186 #endif
00187 !
00188 !     Time offset of the date bounds with regards to the job start date
00189 !     read from the currentfile
00190 !
00191       dl_timeb(1:2)=86400.d0 *(ju_dayb(1:2) &
00192                    - Fields(id_varid)%io_infos%ju_start_day) &
00193                    +          (ju_secb(1:2)  &
00194                    - Fields(id_varid)%io_infos%ju_start_sec)
00195 #ifdef VERBOSE
00196       print*,trim(ch_id),' : PSMILe_read_byid_dble: called' &
00197             ,' ju_dayb=',ju_dayb,' ju_secb=' &
00198             ,ju_secb,' offsetsb=',dl_timeb
00199       call psmile_flushstd
00200 #endif
00201 !
00202 !     Search the I/O cache of time stamps for the closed intervall containing
00203 !     the current time offset.
00204 !     
00205       if(associated(Fields(id_varid)%io_infos%p_cache)) then
00206         pl_cache=>Fields(id_varid)%io_infos%p_cache
00207       else
00208         ierror=PRISM_Error_IO_Read
00209         call psmile_error(ierror, &
00210               'IO cache not allocated!'  &
00211                         ,ierrp,0, __FILE__, __LINE__)
00212       endif
00213       if(associated(Fields(id_varid)%io_infos%p_cache%time_stamps)) then
00214         pl_times=>Fields(id_varid)%io_infos%p_cache%time_stamps
00215         pl_cache=>Fields(id_varid)%io_infos%p_cache
00216       else
00217         ierror=PRISM_Error_IO_Read
00218         call psmile_error(ierror, &
00219               'IO cache for time stamps  not allocated!'  &
00220                         ,ierrp,0, __FILE__, __LINE__)
00221       endif
00222 
00223       lmatch=.FALSE.
00224       llast=>pl_cache%llast
00225       ju_sec_last=>pl_cache%ju_sec_last
00226       ju_day_last=>pl_cache%ju_day_last
00227       time_last=>pl_cache%time_last
00228 !
00229       do while(.not.lmatch)
00230         if(size(pl_times)-1.eq.1) then
00231 !
00232 !       Special treatment if the p_cache contains only 1 time stamp
00233 !
00234           il_i=1
00235           if(abs(dl_time-pl_times(il_i)).lt.1.d-8) then
00236             lmatch=.TRUE.
00237             il_times(1)=il_i
00238             ibl=1
00239             ibu=1
00240             pl_cache%ilast=il_i
00241           else
00242             lower_in=.false.
00243             upper_in=.false.
00244             if(pl_times(il_i).le.dl_timeb(1).and. &
00245                pl_times(il_i).ge.dl_timeb(2)) lower_in=.true.
00246             if(lower_in) then
00247               lmatch=.TRUE.
00248               il_times(1)=il_i
00249               ibl=1
00250               ibu=1
00251               pl_cache%ilast=il_i
00252             else
00253               lmatch=.FALSE.
00254             endif
00255           endif
00256 !
00257         else
00258 !
00259 !       More than 1 time stamp within pcache. Find matching closed intervall.
00260 !
00261           if(llast) then
00262 !
00263 !           From a preceeding file we kept the last time stamp
00264 !
00265             pl_times(0)=86400.d0 *(ju_day_last  &
00266                    - Fields(id_varid)%io_infos%ju_start_day) &
00267                    +          (ju_sec_last  &
00268                    - Fields(id_varid)%io_infos%ju_start_sec)
00269 
00270           else
00271 !
00272 !           Just to skip the time stamp 0 the very first time if no last
00273 !           cached time stamp is present.
00274 !
00275             pl_times(0)=huge(aone)
00276           endif
00277 
00278 !rv          do il_i=pl_cache%ilast,size(pl_times)-1
00279 !
00280 !         Loop over ubound(time stamps)-1 not  ubound(time stamps)
00281 !
00282           do il_i=pl_cache%ilast,size(pl_times)-2
00283             if ((dl_time.ge.pl_times(il_i)) .and.  &
00284                 (dl_time.le.pl_times(il_i+1)) )exit
00285           enddo
00286 !
00287 !         We don't need the stored data of the preceeding file.
00288 !
00289           if(il_i.gt.1) then
00290             llast=.false.
00291             if(associated(pl_cache%data_dble)) &
00292               deallocate(pl_cache%data_dble,stat=ierror)
00293             if (ierror > 0) then
00294               ierrp (1) = ierror
00295                ierrp (2) = 1
00296                ierror = PRISM_Error_Alloc
00297                call psmile_error ( ierror, 'deallocate pl_cache%data_dble', &
00298                             ierrp, 2, __FILE__, __LINE__ )
00299                return
00300             endif
00301           endif
00302 !
00303 !       Now check the found boundaries of the closed intervall 
00304 !       against the date bounds
00305 !
00306           lower_in=.false.
00307           upper_in=.false.
00308           w=-1.
00309           if(pl_times(il_i).ge.dl_timeb(1).and. &
00310              pl_times(il_i).le.dl_timeb(2)) lower_in=.true.
00311           if(il_i.lt.size(pl_times)-1) then
00312             if(pl_times(il_i+1).ge.dl_timeb(1).and. &
00313                pl_times(il_i+1).le.dl_timeb(2)) upper_in=.true.
00314           endif
00315 !
00316           if(il_i.lt.size(pl_times)-1) then
00317             lmatch=.TRUE.
00318 !
00319 !           1st case: None of the intervall boundaries are inside the
00320 !                     date bounds. Return nearest time step.
00321 !
00322             if((.not.lower_in).and.(.not.upper_in)) then
00323 !
00324 !             Which intervall boundary is nearest?
00325 !
00326               if((dl_time-pl_times(il_i)).lt.(pl_times(il_i+1)-dl_time)) then
00327                 il_times(1)=il_i
00328               else
00329                 il_times(1)=il_i+1
00330               endif
00331               ibl=1
00332               ibu=1
00333 !
00334 !           2nd case: Both intervall boundaries are inside the
00335 !                     date bounds. 
00336 !                     Fields(id_varid)%io_infos%ilag_mode = PSMILe_time_nnghbr:
00337 !                     Return nearest neighbour
00338 !                     Fields(id_varid)%io_infos%ilag_mode = PSMILe_time_linear:
00339 !                     Calculate a weight to return
00340 !                     a weighted mean.
00341 !                     DEFAULT:
00342 !                     Calculate a weight to return
00343 !                     a weighted mean.
00344 !
00345             elseif(lower_in.and.upper_in) then
00346 
00347               Select Case(Fields(id_varid)%io_infos%ilag_mode)
00348               Case (PSMILe_time_nnghbr)
00349                 if((dl_time-pl_times(il_i)).lt.(pl_times(il_i+1)-dl_time)) then
00350                   il_times(1)=il_i
00351                 else
00352                   il_times(1)=il_i+1
00353                 endif
00354                 ibl=1
00355                 ibu=1
00356                 timeop=.true.
00357               Case DEFAULT
00358                 w=pl_times(il_i+1)-pl_times(il_i)
00359                 w=(dl_time-pl_times(il_i))/w
00360                 il_times(1)=il_i
00361                 il_times(2)=il_i+1
00362                 ibl=1
00363                 ibu=2
00364                 timeop=.true.
00365               End Select
00366 !
00367 !           3rd case: Only one of the intervall boundaries is inside the
00368 !                     date bounds. Return what was found.
00369 !
00370             elseif(lower_in) then
00371               il_times(1)=il_i
00372               ibl=1
00373               ibu=1
00374             elseif(upper_in) then
00375               il_times(1)=il_i+1
00376               ibl=1
00377               ibu=1
00378             endif
00379             pl_cache%ilast=il_i
00380 #ifdef VERBOSE
00381       print*,trim(ch_id),' : PSMILe_read_byid_dble: Match: ju_day= ',dl_ju_day &
00382             ,' ju_sec= ',dl_ju_sec,' offsets of date bounds=',dl_timeb(1:2) &
00383             ,' time_levels= ',il_times(ibl:ibu),pl_times(il_times(ibl)) &
00384                              ,pl_times(il_times(ibu)),' weight(upper)=',w &
00385                              ,' il_i=',il_i,'ibl,ibu',ibl,ibu,lower_in,upper_in
00386       call psmile_flushstd
00387 #endif
00388           else
00389             
00390 !
00391 !        We did not find a closed intervall in the current file.
00392 !        Cache all time information and trigger the caching of the
00393 !        data as well.
00394 !
00395             lmatch=.FALSE.
00396             llast=.TRUE.
00397             il_times(1)=il_i
00398             ibl=1
00399             ibu=1
00400             time_last=pl_times(il_i)
00401             ju_sec_last=Fields(id_varid)%io_infos%ju_start_sec+time_last
00402             days=dble(int(ju_sec_last/86400.d0))
00403             ju_sec_last=ju_sec_last-days*86400.d0
00404             ju_day_last=Fields(id_varid)%io_infos%ju_start_day + days
00405 #ifdef VERBOSE
00406       print*,trim(ch_id),' : PSMILe_read_byid_dble: No match (EOF): ju_day= ' &
00407             ,dl_ju_day &
00408             ,' ju_sec= ',dl_ju_sec,' offsets of date bounds=',dl_timeb(1:2) &
00409             ,' time_levels= ',il_times(ibl:ibu)
00410       call psmile_flushstd
00411 #endif
00412           endif
00413 !
00414         endif
00415 !
00416 !
00417 !     Vector field ?
00418 !
00419       nvcomp=1
00420       vectorfield=0
00421       if(associated(Methods(il_method_id)%vector_pointer)) then
00422         nvcomp=size(Methods(il_method_id)%vector_pointer%array_of_point_ids)
00423         vectorfield=1
00424       endif
00425 
00426 !
00427 !     Set up shapes depending on the grid type.
00428 !
00429       if ( il_grid_type == PRISM_unstructlonlatvrt) then
00430          len = 1
00431          myvarshape(1:2,1)=Fields(id_varid)%var_shape(1:2,1)
00432          myvarshape(1:2,2)=1
00433          myvarshape(1:2,3)=1
00434          mygrdshape(1:2,1)=Grids(il_grid_id)%grid_shape(1:2,1)
00435          mygrdshape(1:2,2)=1
00436          mygrdshape(1:2,3)=1
00437       else if ( il_grid_type == PRISM_unstructlonlat_sigmavrt .or. &
00438             il_grid_type == PRISM_unstructlonlat_regvrt   .or. &
00439             il_grid_type == PRISM_Gaussreduced_regvrt   .or. &
00440             il_grid_type == PRISM_Gaussreduced_sigmavrt   ) then
00441          len = 2
00442          myvarshape(1:2,1)=Fields(id_varid)%var_shape(1:2,1)
00443          myvarshape(1:2,2)=1
00444          myvarshape(1:2,3)=Fields(id_varid)%var_shape(1:2,3)
00445 
00446          mygrdshape(1:2,1)=Grids(il_grid_id)%grid_shape(1:2,1)
00447          mygrdshape(1:2,2)=1
00448          mygrdshape(1:2,3)=Grids(il_grid_id)%grid_shape(1:2,3)
00449       else
00450          len = il_ndim
00451          myvarshape(1:2,1:len)=Fields(id_varid)%var_shape(1:2,1:len)
00452          mygrdshape(1:2,1:len)=Grids(il_grid_id)%grid_shape(1:2,1:len)
00453       endif
00454 !
00455 !     Is it a bundle ?
00456 !
00457        
00458       lenb=Fields(id_varid)%var_shape(2,len+vectorfield+1) &
00459          -Fields(id_varid)%var_shape(1,len+vectorfield+1)+1
00460 
00461       if(lenb.gt.1) then
00462         is_bundle=.true.
00463         myvarshape(1:2,4)=Fields(id_varid)%var_shape(1:2,len+vectorfield+1)
00464         mygrdshape(1:2,4)=Fields(id_varid)%var_shape(1:2,len+vectorfield+1)
00465       else ! No bundle
00466         is_bundle=.false.
00467       endif
00468 !
00469 !     Calculate  offset to address vector components.
00470 !     Components are written as separate NetCDF variables with 3 axis
00471 !     declared (Sophie Valckes Nov. 2003)
00472 !
00473       offset=1
00474       do ii = 1, len
00475          offset = offset * (Fields(id_varid)%var_shape(2,ii) &
00476                 - Fields(id_varid)%var_shape(1,ii)+1)
00477       enddo
00478 !
00479 !     For GME like grids with a roundrobin distribution of fields over
00480 !     blocks I match the given varid with the set of related ids.
00481 !     The position of a match will be used as a block id.
00482 !     The file and field  information will be take from the last
00483 !     varid of the set of related ids because each block of the GME grid
00484 !     is partitioned in the same manner.
00485 !
00486       if(associated(Fields(id_varid)%io_infos%related_ids)) then
00487         fullsize=size(Fields(id_varid)%io_infos%related_ids)
00488         if(fullsize.gt.1) then
00489           do il_i=1,fullsize
00490             if(id_varid.eq. Fields(id_varid)%io_infos%related_ids(il_i)) exit
00491           enddo
00492 !
00493 !         The MPP_IO file descriptor is always taken from the first entry
00494 !         of the sets of related ids.
00495 !
00496           il_varid=Fields(id_varid)%io_infos%related_ids(1)
00497 !
00498 !         Take the block number from the matching entry within the set
00499 !         of related ids. 
00500 !
00501           il_blockid=Fields(il_i)%io_infos%block_id
00502           is_block=.true.
00503         endif
00504       endif 
00505 
00506       if(is_bundle) then
00507         offset=offset*lenb
00508         if(nvcomp.gt.1) then
00509           fullsize=offset*nvcomp
00510           allocate(afield(1:fullsize),stat=ierror)
00511           if (ierror > 0) then
00512             ierrp (1) = ierror
00513             ierrp (2) = 1
00514             ierror = PRISM_Error_Alloc
00515             call psmile_error ( ierror, 'afield', &
00516                          ierrp, 2, __FILE__, __LINE__ )
00517             return
00518           endif
00519 
00520         endif
00521       endif
00522 
00523 !
00524 !     Read it
00525 !
00526       
00527       do jj=ibl,ibu
00528       il_time=il_times(jj)
00529 #ifdef VERBOSE
00530       print*,trim(ch_id),' : PSMILe_read_byid_dble: ju_day= ',dl_ju_day &
00531             ,' ju_sec= ',dl_ju_sec,' time_level= ',il_time
00532       call psmile_flushstd
00533 #endif
00534       
00535 !
00536 !    Set the first location in the user provided field dd_a.
00537 !
00538 
00539       il_loc=1
00540       
00541       do il_i=1,nvcomp
00542         if(il_time.eq.0.and.llast) then
00543 !
00544 !        We had a match with the last time_stamp of the preceeding file.
00545 !        Return what was cached.
00546 !
00547          dd_a(il_loc:il_loc+offset-1)=pl_cache &
00548                                      %data_dble(il_loc:il_loc+offset-1)
00549         else
00550         
00551         ii=Fields(il_varid)%io_infos%p_mpp_io%findex(il_i)
00552 
00553         if(.not.is_bundle) then
00554         call psmile_read_3d_dble(Fields(il_varid)%io_infos%file_unit &
00555                                  ,Fields(il_varid)%io_infos%p_mpp_io%field(ii) &
00556                                  ,Fields(il_varid)%io_infos%p_mpp_io%domain(1) &
00557                                  ,dd_a(il_loc) &
00558                                  ,myvarshape &
00559                                  ,mygrdshape &
00560                                  ,il_time,.TRUE.,il_blockid,is_block &
00561                                  ,use_domain,ierror)
00562         else if (is_bundle) then
00563           if(vectorfield.eq.0) then
00564             call psmile_read_4d_dble(Fields(il_varid)%io_infos%file_unit &
00565                                  ,Fields(il_varid)%io_infos%p_mpp_io%field(ii) &
00566                                  ,Fields(il_varid)%io_infos%p_mpp_io%domain(1) &
00567                                      ,dd_a(il_loc) &
00568                                      ,myvarshape &
00569                                      ,mygrdshape &
00570                                      ,il_time,.TRUE.,il_blockid,is_block &
00571                                      ,use_domain,ierror)
00572           else
00573             call psmile_read_4d_dble(Fields(il_varid)%io_infos%file_unit &
00574                                  ,Fields(il_varid)%io_infos%p_mpp_io%field(ii) &
00575                                  ,Fields(il_varid)%io_infos%p_mpp_io%domain(1) &
00576                                      ,afield(il_loc) &
00577                                      ,myvarshape &
00578                                      ,mygrdshape &
00579                                      ,il_time,.TRUE.,il_blockid,is_block &
00580                                      ,use_domain,ierror)
00581           endif
00582         endif
00583         endif
00584         il_loc=il_loc+offset
00585       enddo
00586 
00587       if(is_bundle) then
00588         if(nvcomp.gt.1) then
00589 !
00590 !         Individual vector components as bundles were read by looping over
00591 !         vector indices as the last dimension.
00592 !         However, in PSMILE the bundles have to appear as the last 
00593 !         set of indices.
00594 !         Therefore, a transpose of the last two indices has to take place.
00595 !
00596 
00597           if(len.eq.1) call trs_vec_bundle_1d_dble(afield,dd_a &
00598                                                   ,Fields(id_varid)%var_shape &
00599                                                   ,mygrdshape &
00600                                                   ,ierror)
00601           if(len.eq.2) call trs_vec_bundle_2d_dble(afield,dd_a &
00602                                                   ,Fields(id_varid)%var_shape &
00603                                                   ,mygrdshape &
00604                                                   ,ierror)
00605           if(len.eq.3) call trs_vec_bundle_3d_dble(afield,dd_a &
00606                                                   ,Fields(id_varid)%var_shape &
00607                                                   ,mygrdshape &
00608                                                   ,ierror)
00609         endif
00610       endif
00611 !
00612 !
00613 !     In case that the end of the time stamps in one file was hit or
00614 !     or two time stamps fall into the date bounds we have to cache one
00615 !     set of data.
00616 !
00617         if((jj.eq.1.and.ibu.eq.2).or.llast) then
00618 #ifdef DEBUG
00619       print*,trim(ch_id),' : PSMILe_read_byid_dble:',nvcomp*offset &
00620             ,' doubles must be cached . Allocating pcache!'
00621       call psmile_flushstd
00622 #endif
00623 
00624 !
00625           if(associated(pl_cache%data_dble)) then
00626              if(size(pl_cache%data_dble).lt.nvcomp*offset) &
00627                 deallocate(pl_cache%data_dble,stat=ierror)
00628              if (ierror > 0) then
00629                ierrp (1) = ierror
00630                ierrp (2) = 1
00631                ierror = PRISM_Error_Alloc
00632                call psmile_error ( ierror, 'deallocate pl_cache%data_dble', &
00633                              ierrp, 2, __FILE__, __LINE__ )
00634                return
00635              endif
00636           endif
00637 
00638 !
00639           if(.not.associated(pl_cache%data_dble)) &
00640              allocate(pl_cache%data_dble(1:nvcomp*offset),stat=ierror)
00641           if (ierror > 0) then
00642               ierrp (1) = ierror
00643               ierrp (2) = 1
00644               ierror = PRISM_Error_Alloc
00645               call psmile_error ( ierror, 'allocate pl_cache%data_dble', &
00646                             ierrp, 2, __FILE__, __LINE__ )
00647               return
00648           endif
00649 !
00650           pl_cache%data_dble(1:nvcomp*offset)=dd_a(1:nvcomp*offset)
00651 #ifdef DEBUG
00652       print*,trim(ch_id),' : PSMILe_read_byid_dble:',nvcomp*offset &
00653             ,' dbles cached'
00654       call psmile_flushstd
00655 #endif
00656         endif
00657 !
00658       enddo  ! End of loop over matching time stamps
00659 !
00660 !     Compute the waited mean of the time stamps i and i+1
00661 !
00662       if(ibu.eq.2.and.lmatch) then
00663         dd_a(1:nvcomp*offset)=(aone-w)*pl_cache%data_dble(1:nvcomp*offset) &
00664                              +w*dd_a(1:nvcomp*offset)
00665       endif
00666 
00667       if (.not.lmatch) then
00668 !
00669 !         Convert Julian date to Gregorian date in order to construct 
00670 !         a file name if a search for a file containing a new, 
00671 !         advanced job start date has
00672 !         to be performed in order to match the  current date.
00673 !
00674           call psmile_ju2date ( cur_date, dl_ju_day, dl_ju_sec )
00675 !   
00676 !         Check file size and open transparently a new file if necessary.
00677 !
00678           call psmile_open_file_byid(id_varid,id_taskid,cur_date,ierror)
00679       endif
00680       enddo   !End of do while
00681 !
00682 !     Save some memory!
00683 !
00684       if(.not.llast) then
00685         if(associated(pl_cache%data_dble)) &
00686            deallocate(pl_cache%data_dble,stat=ierror)
00687         if (ierror > 0) then
00688           ierrp (1) = ierror
00689           ierrp (2) = 1
00690           ierror = PRISM_Error_Alloc
00691           call psmile_error ( ierror, 'deallocate pl_cache%data_dble', &
00692                               ierrp, 2, __FILE__, __LINE__ )
00693           return
00694         endif
00695       endif
00696       
00697 
00698 #ifdef VERBOSE
00699       print*,trim(ch_id),' : PSMILe_read_byid_dble: end'
00700       call psmile_flushstd
00701 #endif
00702 #endif
00703 
00704       end subroutine psmile_read_byid_dble
00705 
00706       subroutine trs_vec_bundle_1d_dble(user,xio,shape,vshape,ierror)
00707       integer,intent(out)         :: ierror
00708       integer,intent(in)      :: shape(2,3)
00709       integer,intent(in)      :: vshape(2,3)
00710       double precision,intent(in) :: user(shape(1,1):shape(2,1) 
00711                                          ,shape(1,3):shape(2,3) 
00712                                          ,shape(1,2):shape(2,2))
00713       double precision,intent(out)::  xio(shape(1,1):shape(2,1) 
00714                                          ,shape(1,2):shape(2,2) 
00715                                          ,shape(1,3):shape(2,3))
00716 !
00717 !     Local
00718 !
00719       integer                     :: i,j
00720       
00721       ierror=0
00722 
00723       do j=vshape(1,3),vshape(2,3)
00724         do i=vshape(1,2),vshape(2,2)
00725           xio(vshape(1,1):vshape(2,1),i,j)=user(vshape(1,1):vshape(2,1),j,i)
00726         enddo
00727       enddo
00728       
00729       end subroutine trs_vec_bundle_1d_dble
00730 
00731       subroutine trs_vec_bundle_2d_dble(user,xio,shape,vshape,ierror)
00732       integer,intent(out)         :: ierror
00733       integer,intent(in)      :: shape(2,4)
00734       integer,intent(in)      :: vshape(2,4)
00735       double precision,intent(in) :: user(shape(1,1):shape(2,1) 
00736                                          ,shape(1,2):shape(2,2) 
00737                                          ,shape(1,4):shape(2,4) 
00738                                          ,shape(1,3):shape(2,3))
00739       double precision,intent(out)::  xio(shape(1,1):shape(2,1) 
00740                                          ,shape(1,2):shape(2,2) 
00741                                          ,shape(1,3):shape(2,3) 
00742                                          ,shape(1,4):shape(2,4))
00743 !
00744 !     Local
00745 !
00746       integer                     :: i,j
00747       
00748       ierror=0
00749 
00750       do j=vshape(1,4),vshape(2,4)
00751         do i=vshape(1,3),vshape(2,3)
00752 
00753           xio(vshape(1,1):vshape(2,1) &
00754              ,vshape(1,2):vshape(2,2),i,j)=user(vshape(1,1):vshape(2,1) &
00755                                                ,vshape(1,2):vshape(2,2), j,i)
00756         enddo
00757       enddo
00758       
00759       end subroutine trs_vec_bundle_2d_dble
00760 
00761       subroutine trs_vec_bundle_3d_dble(user,xio,shape,vshape,ierror)
00762       integer,intent(out)         :: ierror
00763       integer,intent(in)      :: shape(2,5)
00764       integer,intent(in)      :: vshape(2,5)
00765       double precision,intent(in) :: user(shape(1,1):shape(2,1) 
00766                                          ,shape(1,2):shape(2,2) 
00767                                          ,shape(1,3):shape(2,3) 
00768                                          ,shape(1,5):shape(2,5) 
00769                                          ,shape(1,4):shape(2,4))
00770       double precision,intent(out)::  xio(shape(1,1):shape(2,1) 
00771                                          ,shape(1,2):shape(2,2) 
00772                                          ,shape(1,3):shape(2,2) 
00773                                          ,shape(1,4):shape(2,4) 
00774                                          ,shape(1,5):shape(2,5))
00775 !
00776 !     Local
00777 !
00778       integer                     :: i,j
00779       
00780       ierror=0
00781 
00782       do j=vshape(1,5),vshape(2,5)
00783         do i=vshape(1,4),vshape(2,4)
00784 
00785           xio(vshape(1,1):vshape(2,1) &
00786              ,vshape(1,2):vshape(2,2) &
00787              ,vshape(1,3):vshape(2,3),i,j)=user(vshape(1,1):vshape(2,1) &
00788                                                ,vshape(1,2):vshape(2,2) &
00789                                                ,vshape(1,3):vshape(2,3) , j,i)
00790         enddo
00791       enddo
00792       
00793       end subroutine trs_vec_bundle_3d_dble