psmile_read_byid_real.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_real
00008 !
00009 ! !INTERFACE:
00010       subroutine psmile_read_byid_real(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_real
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       real,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       real,allocatable                  :: afield(:)
00049       real                              :: 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 real 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. Vogeslang 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_rundef.
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_real.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_real: 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_real: 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 #endif
00176       dl_time=86400.d0 *(dl_ju_day - Fields(id_varid)%io_infos%ju_start_day) &
00177              +          (dl_ju_sec - Fields(id_varid)%io_infos%ju_start_sec)
00178 #ifdef VERBOSE
00179       print*,trim(ch_id),' : PSMILe_read_byid_real: called: ju_day= ' &
00180             ,ju_day ,' ju_sec= ',ju_sec,' ju_dayb=',ju_dayb,' ju_secb=' &
00181             ,ju_secb,' offset=',dl_time
00182       call psmile_flushstd
00183 #endif
00184 !
00185 !     Time offset of the date bounds with regards to the job start date
00186 !     read from the currentfile
00187 !
00188       dl_timeb(1:2)=86400.d0 *(ju_dayb(1:2) &
00189                    - Fields(id_varid)%io_infos%ju_start_day) &
00190                    +          (ju_secb(1:2)  &
00191                    - Fields(id_varid)%io_infos%ju_start_sec)
00192 !
00193 !     Search the I/O cache of time stamps for the closed intervall containing
00194 !     the current time offset.
00195 !     
00196       if(associated(Fields(id_varid)%io_infos%p_cache%time_stamps)) then
00197         pl_times=>Fields(id_varid)%io_infos%p_cache%time_stamps
00198         pl_cache=>Fields(id_varid)%io_infos%p_cache
00199       else
00200         ierror=PRISM_Error_IO_Read
00201         call psmile_error(ierror, &
00202               'IO cache for time stamps  not allocated!'  &
00203                         ,ierrp,0, __FILE__, __LINE__)
00204       endif
00205 
00206       lmatch=.FALSE.
00207       llast=>pl_cache%llast
00208       ju_sec_last=>pl_cache%ju_sec_last
00209       ju_day_last=>pl_cache%ju_day_last
00210       time_last=>pl_cache%time_last
00211 !
00212       do while(.not.lmatch)
00213         if(size(pl_times)-1.eq.1) then
00214 !
00215 !       Special treatment if the p_cache contains only 1 time stamp
00216 !
00217           il_i=1
00218           if(abs(dl_time-pl_times(il_i)).lt.1.d-8) then
00219             lmatch=.TRUE.
00220             il_times(1)=il_i
00221             ibl=1
00222             ibu=1
00223             pl_cache%ilast=il_i
00224           else
00225             lower_in=.false.
00226             upper_in=.false.
00227             if(pl_times(il_i).le.dl_timeb(1).and. &
00228                pl_times(il_i).ge.dl_timeb(2)) lower_in=.true.
00229             if(lower_in) then
00230               lmatch=.TRUE.
00231               il_times(1)=il_i
00232               ibl=1
00233               ibu=1
00234               pl_cache%ilast=il_i
00235             else
00236               lmatch=.FALSE.
00237             endif
00238           endif
00239 !
00240         else
00241 !
00242 !       More than 1 time stamp within pcache. Find matching closed intervall.
00243 !
00244           if(llast) then
00245 !
00246 !           From a preceeding file we kept the last time stamp
00247 !
00248             pl_times(0)=86400.d0 *(ju_day_last  &
00249                    - Fields(id_varid)%io_infos%ju_start_day) &
00250                    +          (ju_sec_last  &
00251                    - Fields(id_varid)%io_infos%ju_start_sec)
00252 
00253           else
00254 !
00255 !           Just to skip the time stamp 0 the very first time if no last
00256 !           cached time stamp is present.
00257 !
00258             pl_times(0)=huge(aone)
00259           endif
00260 
00261 !rv          do il_i=pl_cache%ilast,size(pl_times)-1
00262 !
00263 !         Loop over ubound(time stamps)-1 not  ubound(time stamps)
00264 !
00265           do il_i=pl_cache%ilast,size(pl_times)-2
00266             if ((dl_time.ge.pl_times(il_i)) .and.  &
00267                 (dl_time.le.pl_times(il_i+1)) )exit
00268           enddo
00269 !
00270 !         We don't need the stored data of the preceeding file.
00271 !
00272           if(il_i.gt.1) then
00273             llast=.false.
00274             if(associated(pl_cache%data_real)) &
00275               deallocate(pl_cache%data_real,stat=ierror)
00276             if (ierror > 0) then
00277               ierrp (1) = ierror
00278                ierrp (2) = 1
00279                ierror = PRISM_Error_Alloc
00280                call psmile_error ( ierror, 'deallocate pl_cache%data_real', &
00281                             ierrp, 2, __FILE__, __LINE__ )
00282                return
00283             endif
00284           endif
00285 !
00286 !       Now check the found boundaries of the closed intervall 
00287 !       against the date bounds
00288 !
00289           lower_in=.false.
00290           upper_in=.false.
00291           w=-1.
00292           if(pl_times(il_i).ge.dl_timeb(1).and. &
00293              pl_times(il_i).le.dl_timeb(2)) lower_in=.true.
00294           if(il_i.lt.size(pl_times)-1) then
00295             if(pl_times(il_i+1).ge.dl_timeb(1).and. &
00296                pl_times(il_i+1).le.dl_timeb(2)) upper_in=.true.
00297           endif
00298 !
00299           if(il_i.lt.size(pl_times)-1) then
00300             lmatch=.TRUE.
00301 !
00302 !           1st case: None of the intervall boundaries are inside the
00303 !                     date bounds. Return nearest time step.
00304 !
00305             if((.not.lower_in).and.(.not.upper_in)) then
00306 !
00307 !             Which intervall boundary is nearest?
00308 !
00309               if((dl_time-pl_times(il_i)).lt.(pl_times(il_i+1)-dl_time)) then
00310                 il_times(1)=il_i
00311               else
00312                 il_times(1)=il_i+1
00313               endif
00314               ibl=1
00315               ibu=1
00316 !
00317 !           2nd case: Both intervall boundaries are inside the
00318 !                     date bounds. 
00319 !                     Fields(id_varid)%io_infos%ilag_mode = PSMILe_time_nnghbr:
00320 !                     Return nearest neighbour
00321 !                     Fields(id_varid)%io_infos%ilag_mode = PSMILe_time_linear:
00322 !                     Calculate a weight to return
00323 !                     a weighted mean.
00324 !                     DEFAULT:
00325 !                     Calculate a weight to return
00326 !                     a weighted mean.
00327 !
00328             elseif(lower_in.and.upper_in) then
00329 
00330               Select Case(Fields(id_varid)%io_infos%ilag_mode)
00331               Case (PSMILe_time_nnghbr)
00332                 if((dl_time-pl_times(il_i)).lt.(pl_times(il_i+1)-dl_time)) then
00333                   il_times(1)=il_i
00334                 else
00335                   il_times(1)=il_i+1
00336                 endif
00337                 ibl=1
00338                 ibu=1
00339                 timeop=.true.
00340               Case DEFAULT
00341                 w=pl_times(il_i+1)-pl_times(il_i)
00342                 w=(dl_time-pl_times(il_i))/w
00343                 il_times(1)=il_i
00344                 il_times(2)=il_i+1
00345                 ibl=1
00346                 ibu=2
00347                 timeop=.true.
00348               End Select
00349 !
00350 !           3rd case: Only one of the intervall boundaries is inside the
00351 !                     date bounds. Return what was found.
00352 !
00353             elseif(lower_in) then
00354               il_times(1)=il_i
00355               ibl=1
00356               ibu=1
00357             elseif(upper_in) then
00358               il_times(1)=il_i+1
00359               ibl=1
00360               ibu=1
00361             endif
00362             pl_cache%ilast=il_i
00363 #ifdef VERBOSE
00364       print*,trim(ch_id),' : PSMILe_read_byid_real: Match: ju_day= ',dl_ju_day &
00365             ,' ju_sec= ',dl_ju_sec,' offsets of date bounds=',dl_timeb(1:2) &
00366             ,' time_levels= ',il_times(ibl:ibu),pl_times(il_times(ibl)) &
00367                              ,pl_times(il_times(ibu)),' weight(upper)=',w
00368       call psmile_flushstd
00369 #endif
00370           else
00371             
00372 !
00373 !        We did not find a closed intervall in the current file.
00374 !        Cache all time information and trigger the caching of the
00375 !        data as well.
00376 !
00377             lmatch=.FALSE.
00378             llast=.TRUE.
00379             il_times(1)=il_i
00380             ibl=1
00381             ibu=1
00382             time_last=pl_times(il_i)
00383             ju_sec_last=Fields(id_varid)%io_infos%ju_start_sec+time_last
00384             days=dble(int(ju_sec_last/86400.d0))
00385             ju_sec_last=ju_sec_last-days*86400.d0
00386             ju_day_last=Fields(id_varid)%io_infos%ju_start_day + days
00387 #ifdef VERBOSE
00388       print*,trim(ch_id),' : PSMILe_read_byid_real: No match (EOF): ju_day= ' &
00389             ,dl_ju_day &
00390             ,' ju_sec= ',dl_ju_sec,' offsets of date bounds=',dl_timeb(1:2) &
00391             ,' time_levels= ',il_times(ibl:ibu)
00392       call psmile_flushstd
00393 #endif
00394           endif
00395 !
00396         endif
00397 !
00398 !
00399 !     Vector field ?
00400 !
00401       nvcomp=1
00402       vectorfield=0
00403       if(associated(Methods(il_method_id)%vector_pointer)) then
00404         nvcomp=size(Methods(il_method_id)%vector_pointer%array_of_point_ids)
00405         vectorfield=1
00406       endif
00407 
00408 !
00409 !     Set up shapes depending on the grid type.
00410 !
00411       if ( il_grid_type == PRISM_unstructlonlatvrt) then
00412          len = 1
00413          myvarshape(1:2,1)=Fields(id_varid)%var_shape(1:2,1)
00414          myvarshape(1:2,2)=1
00415          myvarshape(1:2,3)=1
00416          mygrdshape(1:2,1)=Grids(il_grid_id)%grid_shape(1:2,1)
00417          mygrdshape(1:2,2)=1
00418          mygrdshape(1:2,3)=1
00419       else if ( il_grid_type == PRISM_unstructlonlat_sigmavrt .or. &
00420             il_grid_type == PRISM_unstructlonlat_regvrt   .or. &
00421             il_grid_type == PRISM_Gaussreduced_regvrt   .or. &
00422             il_grid_type == PRISM_Gaussreduced_sigmavrt   ) then
00423          len = 2
00424          myvarshape(1:2,1)=Fields(id_varid)%var_shape(1:2,1)
00425          myvarshape(1:2,2)=1
00426          myvarshape(1:2,3)=Fields(id_varid)%var_shape(1:2,3)
00427 
00428          mygrdshape(1:2,1)=Grids(il_grid_id)%grid_shape(1:2,1)
00429          mygrdshape(1:2,2)=1
00430          mygrdshape(1:2,3)=Grids(il_grid_id)%grid_shape(1:2,3)
00431       else
00432          len = il_ndim
00433          myvarshape(1:2,1:len)=Fields(id_varid)%var_shape(1:2,1:len)
00434          mygrdshape(1:2,1:len)=Grids(il_grid_id)%grid_shape(1:2,1:len)
00435       endif
00436 !
00437 !     Is it a bundle ?
00438 !
00439        
00440       lenb=Fields(id_varid)%var_shape(2,len+vectorfield+1) &
00441          -Fields(id_varid)%var_shape(1,len+vectorfield+1)+1
00442 
00443       if(lenb.gt.1) then
00444         is_bundle=.true.
00445         myvarshape(1:2,4)=Fields(id_varid)%var_shape(1:2,len+vectorfield+1)
00446         mygrdshape(1:2,4)=Fields(id_varid)%var_shape(1:2,len+vectorfield+1)
00447       else ! No bundle
00448         is_bundle=.false.
00449       endif
00450 !
00451 !     Calculate  offset to address vector components.
00452 !     Components are written as separate NetCDF variables with 3 axis
00453 !     declared (Sophie Valckes Nov. 2003)
00454 !
00455       offset=1
00456       do ii = 1, len
00457          offset = offset * (Fields(id_varid)%var_shape(2,ii) &
00458                 - Fields(id_varid)%var_shape(1,ii)+1)
00459       enddo
00460 !
00461 !     For GME like grids with a roundrobin distribution of fields over
00462 !     blocks I match the given varid with the set of related ids.
00463 !     The position of a match will be used as a block id.
00464 !     The file and field  information will be take from the last
00465 !     varid of the set of related ids because each block of the GME grid
00466 !     is partitioned in the same manner.
00467 !
00468       if(associated(Fields(id_varid)%io_infos%related_ids)) then
00469         fullsize=size(Fields(id_varid)%io_infos%related_ids)
00470         if(fullsize.gt.1) then
00471           do il_i=1,fullsize
00472             if(id_varid.eq. Fields(id_varid)%io_infos%related_ids(il_i)) exit
00473           enddo
00474 !
00475 !         The MPP_IO file descriptor is always taken from the first entry
00476 !         of the sets of related ids.
00477 !
00478           il_varid=Fields(id_varid)%io_infos%related_ids(1)
00479 !
00480 !         Take the block number from the matching entry within the set
00481 !         of related ids. 
00482 !
00483           il_blockid=Fields(il_i)%io_infos%block_id
00484           is_block=.true.
00485         endif
00486       endif 
00487 
00488       if(is_bundle) then
00489         offset=offset*lenb
00490         if(nvcomp.gt.1) then
00491           fullsize=offset*nvcomp
00492           allocate(afield(1:fullsize),stat=ierror)
00493           if (ierror > 0) then
00494             ierrp (1) = ierror
00495             ierrp (2) = 1
00496             ierror = PRISM_Error_Alloc
00497             call psmile_error ( ierror, 'afield', &
00498                          ierrp, 2, __FILE__, __LINE__ )
00499             return
00500           endif
00501 
00502         endif
00503       endif
00504 
00505 !
00506 !     Read it
00507 !
00508       
00509       do jj=ibl,ibu
00510       il_time=il_times(jj)
00511 #ifdef VERBOSE
00512       print*,trim(ch_id),' : PSMILe_read_byid_real: ju_day= ',dl_ju_day &
00513             ,' ju_sec= ',dl_ju_sec,' time_level= ',il_time
00514       call psmile_flushstd
00515 #endif
00516       
00517 !
00518 !    Set the first location in the user provided field dd_a.
00519 !
00520 
00521       il_loc=1
00522       
00523       do il_i=1,nvcomp
00524         if(il_time.eq.0.and.llast) then
00525 !
00526 !        We had a match with the last time_stamp of the preceeding file.
00527 !        Return what was cached.
00528 !
00529          dd_a(il_loc:il_loc+offset-1)=pl_cache &
00530                                      %data_real(il_loc:il_loc+offset-1)
00531         else
00532         
00533         ii=Fields(il_varid)%io_infos%p_mpp_io%findex(il_i)
00534 
00535         if(.not.is_bundle) then
00536         call psmile_read_3d_real(Fields(il_varid)%io_infos%file_unit &
00537                                  ,Fields(il_varid)%io_infos%p_mpp_io%field(ii) &
00538                                  ,Fields(il_varid)%io_infos%p_mpp_io%domain(1) &
00539                                  ,dd_a(il_loc) &
00540                                  ,myvarshape &
00541                                  ,mygrdshape &
00542                                  ,il_time,.TRUE.,il_blockid,is_block &
00543                                  ,use_domain,ierror)
00544         else if (is_bundle) then
00545           if(vectorfield.eq.0) then
00546             call psmile_read_4d_real(Fields(il_varid)%io_infos%file_unit &
00547                                  ,Fields(il_varid)%io_infos%p_mpp_io%field(ii) &
00548                                  ,Fields(il_varid)%io_infos%p_mpp_io%domain(1) &
00549                                      ,dd_a(il_loc) &
00550                                      ,myvarshape &
00551                                      ,mygrdshape &
00552                                      ,il_time,.TRUE.,il_blockid,is_block &
00553                                      ,use_domain,ierror)
00554           else
00555             call psmile_read_4d_real(Fields(il_varid)%io_infos%file_unit &
00556                                  ,Fields(il_varid)%io_infos%p_mpp_io%field(ii) &
00557                                  ,Fields(il_varid)%io_infos%p_mpp_io%domain(1) &
00558                                      ,afield(il_loc) &
00559                                      ,myvarshape &
00560                                      ,mygrdshape &
00561                                      ,il_time,.TRUE.,il_blockid,is_block &
00562                                      ,use_domain,ierror)
00563           endif
00564         endif
00565         endif
00566         il_loc=il_loc+offset
00567       enddo
00568 
00569       if(is_bundle) then
00570         if(nvcomp.gt.1) then
00571 !
00572 !         Individual vector components as bundles were read by looping over
00573 !         vector indices as the last dimension.
00574 !         However, in PSMILE the bundles have to appear as the last 
00575 !         set of indices.
00576 !         Therefore, a transpose of the last two indices has to take place.
00577 !
00578 
00579           if(len.eq.1) call trs_vec_bundle_1d_real(afield,dd_a &
00580                                                   ,Fields(id_varid)%var_shape &
00581                                                   ,mygrdshape &
00582                                                   ,ierror)
00583           if(len.eq.2) call trs_vec_bundle_2d_real(afield,dd_a &
00584                                                   ,Fields(id_varid)%var_shape &
00585                                                   ,mygrdshape &
00586                                                   ,ierror)
00587           if(len.eq.3) call trs_vec_bundle_3d_real(afield,dd_a &
00588                                                   ,Fields(id_varid)%var_shape &
00589                                                   ,mygrdshape &
00590                                                   ,ierror)
00591         endif
00592       endif
00593 !
00594 !
00595 !     In case that the end of the time stamps in one file was hit or
00596 !     or two time stamps fall into the date bounds we have to cache one
00597 !     set of data.
00598 !
00599         if((jj.eq.1.and.ibu.eq.2).or.llast) then
00600 #ifdef DEBUG
00601       print*,trim(ch_id),' : PSMILe_read_byid_real:',nvcomp*offset &
00602             ,' reals must be cached . Allocating pcache!'
00603       call psmile_flushstd
00604 #endif
00605 !
00606           if(associated(pl_cache%data_real)) then
00607              if(size(pl_cache%data_real).lt.nvcomp*offset) &
00608                 deallocate(pl_cache%data_real,stat=ierror)
00609              if (ierror > 0) then
00610                ierrp (1) = ierror
00611                ierrp (2) = 1
00612                ierror = PRISM_Error_Alloc
00613                call psmile_error ( ierror, 'deallocate pl_cache%data_real', &
00614                              ierrp, 2, __FILE__, __LINE__ )
00615                return
00616              endif
00617           endif
00618 
00619 !
00620           if(.not.associated(pl_cache%data_real)) &
00621              allocate(pl_cache%data_real(1:nvcomp*offset),stat=ierror)
00622           if (ierror > 0) then
00623               ierrp (1) = ierror
00624               ierrp (2) = 1
00625               ierror = PRISM_Error_Alloc
00626               call psmile_error ( ierror, 'allocate pl_cache%data_real', &
00627                             ierrp, 2, __FILE__, __LINE__ )
00628               return
00629           endif
00630 !
00631           pl_cache%data_real(1:nvcomp*offset)=dd_a(1:nvcomp*offset)
00632 #ifdef DEBUG
00633       print*,trim(ch_id),' : PSMILe_read_byid_real:',nvcomp*offset &
00634             ,' reals cached'
00635       call psmile_flushstd
00636 #endif
00637         endif
00638 !
00639       enddo  ! End of loop over matching time stamps
00640 !
00641 !     Compute the waited mean of the time stamps i and i+1
00642 !
00643       if(ibu.eq.2.and.lmatch) then
00644         dd_a(1:nvcomp*offset)=(aone-w)*pl_cache%data_real(1:nvcomp*offset) &
00645                              +w*dd_a(1:nvcomp*offset)
00646       endif
00647       
00648 
00649       if (.not.lmatch) then
00650 !
00651 !         Convert Julian date to Gregorian date in order to construct 
00652 !         a file name if a search for a file containing a new, 
00653 !         advanced job start date has
00654 !         to be performed in order to match the  current date.
00655 !
00656           call psmile_ju2date ( cur_date, dl_ju_day, dl_ju_sec )
00657 !   
00658 !         Check file size and open transparently a new file if necessary.
00659 !
00660           call psmile_open_file_byid(id_varid,id_taskid,cur_date,ierror)
00661       endif
00662       enddo   !End of do while
00663 !
00664 !     Save some memory!
00665 !
00666       if(.not.llast) then
00667         if(associated(pl_cache%data_real)) &
00668            deallocate(pl_cache%data_real,stat=ierror)
00669         if (ierror > 0) then
00670           ierrp (1) = ierror
00671           ierrp (2) = 1
00672           ierror = PRISM_Error_Alloc
00673           call psmile_error ( ierror, 'deallocate pl_cache%data_real', &
00674                               ierrp, 2, __FILE__, __LINE__ )
00675           return
00676         endif
00677       endif
00678       
00679 
00680 #ifdef VERBOSE
00681       print*,trim(ch_id),' : PSMILe_read_byid_real: end'
00682       call psmile_flushstd
00683 #endif
00684 #endif
00685 
00686       end subroutine psmile_read_byid_real
00687 
00688       subroutine trs_vec_bundle_1d_real(user,xio,shape,vshape,ierror)
00689       integer,intent(out)         :: ierror
00690       integer,intent(in)      :: shape(2,3)
00691       integer,intent(in)      :: vshape(2,3)
00692       real,intent(in)             :: user(shape(1,1):shape(2,1) 
00693                                          ,shape(1,3):shape(2,3) 
00694                                          ,shape(1,2):shape(2,2))
00695       real,intent(out)            ::  xio(shape(1,1):shape(2,1) 
00696                                          ,shape(1,2):shape(2,2) 
00697                                          ,shape(1,3):shape(2,3))
00698 !
00699 !     Local
00700 !
00701       integer                     :: i,j
00702       
00703       ierror=0
00704 
00705       do j=vshape(1,3),vshape(2,3)
00706         do i=vshape(1,2),vshape(2,2)
00707           xio(vshape(1,1):vshape(2,1),i,j)=user(vshape(1,1):vshape(2,1),j,i)
00708         enddo
00709       enddo
00710       
00711       end subroutine trs_vec_bundle_1d_real
00712 
00713       subroutine trs_vec_bundle_2d_real(user,xio,shape,vshape,ierror)
00714       integer,intent(out)         :: ierror
00715       integer,intent(in)      :: shape(2,4)
00716       integer,intent(in)      :: vshape(2,4)
00717       real,intent(in)             :: user(shape(1,1):shape(2,1) 
00718                                          ,shape(1,2):shape(2,2) 
00719                                          ,shape(1,4):shape(2,4) 
00720                                          ,shape(1,3):shape(2,3))
00721       real,intent(out)            ::  xio(shape(1,1):shape(2,1) 
00722                                          ,shape(1,2):shape(2,2) 
00723                                          ,shape(1,3):shape(2,3) 
00724                                          ,shape(1,4):shape(2,4))
00725 !
00726 !     Local
00727 !
00728       integer                     :: i,j
00729       
00730       ierror=0
00731 
00732       do j=vshape(1,4),vshape(2,4)
00733         do i=vshape(1,3),vshape(2,3)
00734 
00735           xio(vshape(1,1):vshape(2,1) &
00736              ,vshape(1,2):vshape(2,2),i,j)=user(vshape(1,1):vshape(2,1) &
00737                                                ,vshape(1,2):vshape(2,2), j,i)
00738         enddo
00739       enddo
00740       
00741       end subroutine trs_vec_bundle_2d_real
00742 
00743       subroutine trs_vec_bundle_3d_real(user,xio,shape,vshape,ierror)
00744       integer,intent(out)         :: ierror
00745       integer,intent(in)      :: shape(2,5)
00746       integer,intent(in)      :: vshape(2,5)
00747       real,intent(in)             :: user(shape(1,1):shape(2,1) 
00748                                          ,shape(1,2):shape(2,2) 
00749                                          ,shape(1,3):shape(2,3) 
00750                                          ,shape(1,5):shape(2,5) 
00751                                          ,shape(1,4):shape(2,4))
00752       real,intent(out)            ::  xio(shape(1,1):shape(2,1) 
00753                                          ,shape(1,2):shape(2,2) 
00754                                          ,shape(1,3):shape(2,2) 
00755                                          ,shape(1,4):shape(2,4) 
00756                                          ,shape(1,5):shape(2,5))
00757 !
00758 !     Local
00759 !
00760       integer                     :: i,j
00761       
00762       ierror=0
00763 
00764       do j=vshape(1,5),vshape(2,5)
00765         do i=vshape(1,4),vshape(2,4)
00766 
00767           xio(vshape(1,1):vshape(2,1) &
00768              ,vshape(1,2):vshape(2,2) &
00769              ,vshape(1,3):vshape(2,3),i,j)=user(vshape(1,1):vshape(2,1) &
00770                                                ,vshape(1,2):vshape(2,2) &
00771                                                ,vshape(1,3):vshape(2,3) , j,i)
00772         enddo
00773       enddo
00774       
00775       end subroutine trs_vec_bundle_3d_real