psmile_write_byid_int.F90

Go to the documentation of this file.
00001 !-----------------------------------------------------------------------
00002 ! Copyright 2007-2010, SGI Germany, Munich, Germany.
00003 ! All rights reserved. Use is subject to OASIS4 license terms.
00004 !-----------------------------------------------------------------------
00005 !BOP
00006 !
00007 ! !ROUTINE:  psmile_write_byid_int
00008 !
00009 ! !INTERFACE:
00010       subroutine psmile_write_byid_int(id_varid,id_taskid  &
00011                                       ,dd_a,ju_day,ju_sec,ierror)
00012 !
00013 ! !USES:
00014 !
00015       use prism_constants
00016       use psmile ! , dummy_interface=> psmile_write_byid_int
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       integer,target, Intent (In) :: dd_a(*)
00026       double precision, Intent (In) :: ju_day, ju_sec
00027 !
00028 ! !OUTPUT PARAMETERS:
00029 !
00030       integer,intent(out):: ierror
00031 !
00032 ! !LOCAL VARIABLES:
00033 !
00034       integer                           :: ierrp(2)
00035       integer               :: il_grid_id,il_grid_type
00036       integer               :: il_method_id,il_varid
00037       integer                           :: il_blockid
00038       integer                           :: iicomp_id
00039       logical                           :: is_block
00040       integer               :: nvcomp,il_ndim,offset,vectorfield
00041       integer                           :: fullsize
00042       integer               :: il_i,lenb,len,ii,il_loc
00043       integer               :: myvarshape(2,ndim_3d+2)
00044       integer               :: mygrdshape(2,ndim_3d+2)
00045       integer,allocatable   :: afield(:)
00046       double precision          :: dl_time,dl_ju_day,dl_ju_sec
00047       logical               :: is_bundle
00048       Type(PRISM_Time_struct)           :: cur_date
00049 !
00050 ! !DESCRIPTION:
00051 !
00052 !  Routine writes a integer field variable.
00053 !  It can handle scalar fields and vector fields as well as bundles of
00054 !  the two types of fields.
00055 !  Last but not least it writes subblocks of the above in case of a
00056 !  GME grid.
00057 !  Time stamps are written relative to the job start date.
00058 !  A check of the file size is performed calling psmile_open_file_byid.
00059 !  In case of bundles of vector fields transpose routines are called
00060 !  which transpose the indices of bundles with vector components.
00061 
00062 !
00063 ! !REVISION HISTORY:
00064 !
00065 !   Date      Programmer   Description
00066 ! ----------  ----------   -----------
00067 !   21.12.03  R. Vogeslang created
00068 !   23.3.04     R. Vogelsang Added the taskid to the interface and related code.
00069 !   02.10.05  R. Vogelsang Gaussian reduced grids added
00070 !EOP
00071 !----------------------------------------------------------------------
00072 !
00073       character(len=len_cvs_string),save :: mcvs = 
00074 '$Id: psmile_write_byid_int.F90 2899 2011-01-20 10:11:51Z hanke $'
00075 
00076 !
00077       ierror=0
00078 #ifdef __PSMILE_WITH_IO
00079 #ifdef VERBOSE
00080       print*,trim(ch_id),' : psmile_write_byid_int: start'
00081       call psmile_flushstd
00082 #endif
00083 !
00084 !     Return from this routine if the field is not subject to file I/O
00085 !     or the I/O action is not MPP_OVERWR (write).
00086 !
00087       if (.not.associated(Fields(id_varid)%io_chan_infos)) then
00088 #ifdef VERBOSE
00089          print*,trim(ch_id),' : psmile_write_byid_int: end'
00090          call psmile_flushstd
00091 #endif
00092          return
00093       endif
00094 !
00095 !     IO infos are pointing to io_chan_infos(taskid)
00096 !
00097       if(id_taskid .le. size(Fields(id_varid)%io_task_lookup)) then
00098 
00099         il_i=Fields(id_varid)%io_task_lookup(id_taskid)
00100 
00101       else
00102 
00103         ierror = PRISM_Error_IO_Meta
00104         ierrp (1) = id_taskid
00105         call psmile_error ( ierror &
00106                           , 'Task id out of range! ', &
00107                             ierrp, 1, __FILE__, __LINE__ )
00108         return
00109 
00110       endif
00111 
00112       if (il_i.gt.0) then
00113 
00114          Fields(id_varid)%io_infos =>  Fields(id_varid)%io_chan_infos(il_i)
00115 
00116       else
00117         ierror = PRISM_Error_IO_Meta
00118         ierrp (1) = id_taskid
00119         call psmile_error ( ierror &
00120                           , 'Negative task id! ', &
00121                             ierrp, 1, __FILE__, __LINE__ )
00122         return
00123       endif
00124 
00125       if (Fields(id_varid)%io_infos%action .ne. MPP_OVERWR ) then
00126 #ifdef VERBOSE
00127          print*,trim(ch_id),' : psmile_write_byid_int: end'
00128          call psmile_flushstd
00129 #endif
00130          return
00131       endif
00132 !
00133 !     Set the pelist for the current active component
00134 !
00135       iicomp_id=Fields(id_varid)%comp_id
00136       call mpp_set_current_pelist(IO_Comps_infos(iicomp_id)%pelist)
00137 
00138 
00139       il_varid=id_varid
00140       il_method_id=Fields(id_varid)%method_id
00141       il_grid_id=Methods(il_method_id)%grid_id
00142       il_grid_type=Grids(il_grid_id)%grid_type
00143       il_ndim=Grids(il_grid_id)%n_dim
00144       il_blockid=0
00145       is_block=.false.
00146       dl_ju_day=ju_day
00147       dl_ju_sec=ju_sec
00148 !
00149 !     Time offset with respect to the date of the job start
00150 !
00151       dl_time=86400.d0 *(dl_ju_day - Fields(id_varid)%io_infos%ju_start_day) &
00152              +          (dl_ju_sec - Fields(id_varid)%io_infos%ju_start_sec)
00153 !
00154 !     Convert Julian date to Gregorian date in order to construct a file name
00155 !     if a new file has to be created due to an exceed of the file size
00156 !     threshold.
00157 !
00158       call psmile_ju2date ( cur_date, dl_ju_day, dl_ju_sec )
00159 !   
00160 !     Check file size and open transparently a new file if necessary.
00161 !
00162       call psmile_open_file_byid(id_varid,id_taskid,cur_date,ierror)
00163 !
00164 !     Vector field ?
00165 !
00166       nvcomp=1
00167       vectorfield=0
00168       if(associated(Methods(il_method_id)%vector_pointer)) then
00169         nvcomp=size(Methods(il_method_id)%vector_pointer%array_of_point_ids)
00170         vectorfield=1
00171       endif
00172 
00173 !
00174 !     Set up shapes depending on the grid type.
00175 !
00176       if ( il_grid_type == PRISM_unstructlonlatvrt) then
00177          len = 1
00178          myvarshape(1:2,1)=Fields(id_varid)%var_shape(1:2,1)
00179          myvarshape(1:2,2)=1
00180          myvarshape(1:2,3)=1
00181          mygrdshape(1:2,1)=Grids(il_grid_id)%grid_shape(1:2,1)
00182          mygrdshape(1:2,2)=1
00183          mygrdshape(1:2,3)=1
00184       else if ( il_grid_type == PRISM_unstructlonlat_sigmavrt .or. &
00185             il_grid_type == PRISM_unstructlonlat_regvrt   .or. &
00186             il_grid_type == PRISM_Gaussreduced_regvrt   .or. &
00187             il_grid_type == PRISM_Gaussreduced_sigmavrt   ) then
00188 
00189          len = 3
00190          myvarshape(1:2,1)=Fields(id_varid)%var_shape(1:2,1)
00191          myvarshape(1:2,2)=1
00192          myvarshape(1:2,3)=Fields(id_varid)%var_shape(1:2,3)
00193 
00194          mygrdshape(1:2,1)=Grids(il_grid_id)%grid_shape(1:2,1)
00195          mygrdshape(1:2,2)=1
00196          mygrdshape(1:2,3)=Grids(il_grid_id)%grid_shape(1:2,3)
00197       else
00198          len = il_ndim
00199          myvarshape(1:2,1:len)=Fields(id_varid)%var_shape(1:2,1:len)
00200          mygrdshape(1:2,1:len)=Grids(il_grid_id)%grid_shape(1:2,1:len)
00201       endif
00202 !
00203 !     Is it a bundle ?
00204 !
00205        
00206       lenb=Fields(id_varid)%var_shape(2,len+vectorfield+1) &
00207          -Fields(id_varid)%var_shape(1,len+vectorfield+1)+1
00208 
00209       if(lenb.gt.1) then
00210         is_bundle=.true.
00211         myvarshape(1:2,4)=Fields(id_varid)%var_shape(1:2,len+vectorfield+1)
00212         mygrdshape(1:2,4)=Fields(id_varid)%var_shape(1:2,len+vectorfield+1)
00213       else ! No bundle
00214         is_bundle=.false.
00215       endif
00216 !
00217 !     Calculate  offset to address vector components.
00218 !     Components are written as separate NetCDF variables with 3 axis
00219 !     declared (Sophie Valckes Nov. 2003)
00220 !
00221       offset=1
00222       do ii = 1, len
00223          offset = offset * (Fields(id_varid)%var_shape(2,ii) &
00224                 - Fields(id_varid)%var_shape(1,ii)+1)
00225       enddo
00226 !
00227 !    Set the first location in the user provided field
00228 !
00229 
00230       il_loc=1
00231 !
00232 !     For GME like grids with a roundrobin distribution of fields over
00233 !     blocks I match the given varid with the set of related ids.
00234 !     The position of a match will be used as a block id.
00235 !     The file and field  information will be take from the last
00236 !     varid of the set of related ids because each block of the GME grid
00237 !     is partitioned in the same manner.
00238 !
00239       if(associated(Fields(id_varid)%io_infos%related_ids)) then
00240         fullsize=size(Fields(id_varid)%io_infos%related_ids)
00241         if(fullsize.gt.1) then
00242           do il_i=1,fullsize
00243             if(id_varid.eq. Fields(id_varid)%io_infos%related_ids(il_i)) exit
00244           enddo
00245 !
00246 !         The MPP_IO file descriptor is always taken from the first entry
00247 !         of the sets of related ids.
00248 !
00249           il_varid=Fields(id_varid)%io_infos%related_ids(1)
00250 !
00251 !         Take the block number from the matching entry within the set
00252 !         of related ids. 
00253 !
00254           il_blockid=Fields(il_i)%io_infos%block_id
00255           is_block=.true.
00256         endif
00257       endif 
00258 
00259       if(is_bundle) then
00260         offset=offset*lenb
00261         if(nvcomp.gt.1) then
00262           fullsize=offset*nvcomp
00263           allocate(afield(1:fullsize),stat=ierror)
00264           if (ierror > 0) then
00265             ierrp (1) = ierror
00266             ierrp (2) = 1
00267             ierror = PRISM_Error_Alloc
00268             call psmile_error ( ierror, 'afield', &
00269                          ierrp, 2, __FILE__, __LINE__ )
00270             return
00271           endif
00272 
00273 !
00274 !         Bundles appear as the last set of indices
00275 !         In order to write  individual vector components as bundles the
00276 !         vector indices have to be last dimension.
00277 !         Therefore, a transpose of the last two indices has to take place.
00278 
00279           if(len.eq.1) call trs_bundle_vec_1d_integer(dd_a,afield &
00280                                                   ,Fields(id_varid)%var_shape &
00281                                                   ,mygrdshape &
00282                                                   ,ierror)
00283           if(len.eq.2) call trs_bundle_vec_2d_integer(dd_a,afield &
00284                                                   ,Fields(id_varid)%var_shape &
00285                                                   ,mygrdshape &
00286                                                   ,ierror)
00287           if(len.eq.3) call trs_bundle_vec_3d_integer(dd_a,afield &
00288                                                   ,Fields(id_varid)%var_shape &
00289                                                   ,mygrdshape &
00290                                                   ,ierror)
00291         endif
00292       endif
00293 
00294 !
00295 !     Dump it
00296 !
00297       
00298       do ii=1,nvcomp
00299         
00300         if(.not.is_bundle) then
00301         call psmile_write_3d_int(Fields(il_varid)%io_infos%file_unit &
00302                                  ,Fields(il_varid)%io_infos%p_mpp_io%field(ii) &
00303                                  ,Fields(il_varid)%io_infos%p_mpp_io%domain(1) &
00304                                  ,dd_a(il_loc) &
00305                                  ,myvarshape &
00306                                  ,mygrdshape &
00307                                  ,dl_time,.TRUE.,il_blockid,is_block,ierror)
00308         else if (is_bundle) then
00309           if(vectorfield.eq.0) then
00310             call psmile_write_4d_int(Fields(il_varid)%io_infos%file_unit &
00311                                  ,Fields(il_varid)%io_infos%p_mpp_io%field(ii) &
00312                                  ,Fields(il_varid)%io_infos%p_mpp_io%domain(1) &
00313                                      ,dd_a(il_loc) &
00314                                      ,myvarshape &
00315                                      ,mygrdshape &
00316                                      ,dl_time,.TRUE.,il_blockid,is_block,ierror)
00317           else
00318             call psmile_write_4d_int(Fields(il_varid)%io_infos%file_unit &
00319                                  ,Fields(il_varid)%io_infos%p_mpp_io%field(ii) &
00320                                  ,Fields(il_varid)%io_infos%p_mpp_io%domain(1) &
00321                                      ,afield(il_loc) &
00322                                      ,myvarshape &
00323                                      ,mygrdshape &
00324                                      ,dl_time,.TRUE.,il_blockid,is_block,ierror)
00325           endif
00326         endif
00327         il_loc=il_loc+offset
00328       enddo
00329 !
00330 !     Update the file size in kbytes (sizeof(integer)/1024)=0.00390625)
00331 !     Note that the file size is here an rough estimate only. 
00332 !
00333       Fields(il_varid)%io_infos%old_filesize= &
00334             Fields(il_varid)%io_infos%current_filesize
00335       Fields(il_varid)%io_infos%current_filesize= &
00336             Fields(il_varid)%io_infos%current_filesize+max(0.00390625 &
00337                                                        *offset*nvcomp,1.)
00338 #ifdef VERBOSE
00339       print*,trim(ch_id),' : psmile_write_byid_int: end'
00340       call psmile_flushstd
00341 #endif
00342 #endif
00343 
00344       end subroutine psmile_write_byid_int
00345 
00346       subroutine trs_bundle_vec_1d_integer(user,xio,shape,vshape,ierror)
00347       integer,intent(out)         :: ierror
00348       integer,intent(in)      :: shape(2,3)
00349       integer,intent(in)      :: vshape(2,3)
00350       integer,intent(in) :: user(shape(1,1):shape(2,1) 
00351                                          ,shape(1,2):shape(2,2) 
00352                                          ,shape(1,3):shape(2,3))
00353       integer,intent(out)::  xio(shape(1,1):shape(2,1) 
00354                                          ,shape(1,3):shape(2,3) 
00355                                          ,shape(1,2):shape(2,2))
00356 !
00357 !     Local
00358 !
00359       integer                     :: i,j
00360       
00361 !
00362 !     The statement below is semantically the same as the manually coded loops.
00363 !     However, don't know what is faster. For now I trust my coding.
00364 !
00365 !      xio(shape(1,1):shape(2,1)                                           &
00366 !         ,shape(1,3):shape(2,3)                                           &
00367 !         ,shape(1,2):shape(2,2))=reshape(user,(/shape(2,1)-shape(1,1)+1    &
00368 !                                               ,shape(2,3)-shape(1,3)+1    &
00369 !                                               ,shape(2,1)-shape(1,2)+1 /) &
00370 !                                             ,ORDER=(/1,3,2/))
00371       ierror=0
00372 
00373       do j=shape(1,3),shape(2,3)
00374         do i=shape(1,2),shape(2,2)
00375           xio(vshape(1,1):vshape(2,1),j,i)=user(vshape(1,1):vshape(2,1),i,j)
00376         enddo
00377       enddo
00378       
00379       end subroutine trs_bundle_vec_1d_integer
00380 
00381       subroutine trs_bundle_vec_2d_integer(user,xio,shape,vshape,ierror)
00382       integer,intent(out)         :: ierror
00383       integer,intent(in)      :: shape(2,4)
00384       integer,intent(in)      :: vshape(2,4)
00385       integer,intent(in) :: user(shape(1,1):shape(2,1) 
00386                                          ,shape(1,2):shape(2,2) 
00387                                          ,shape(1,3):shape(2,3) 
00388                                          ,shape(1,4):shape(2,4))
00389       integer,intent(out)::  xio(shape(1,1):shape(2,1) 
00390                                          ,shape(1,2):shape(2,2) 
00391                                          ,shape(1,4):shape(2,4) 
00392                                          ,shape(1,3):shape(2,3))
00393 !
00394 !     Local
00395 !
00396       integer                     :: i,j
00397       
00398       ierror=0
00399 
00400       do j=shape(1,4),shape(2,4)
00401         do i=shape(1,3),shape(2,3)
00402 
00403           xio(vshape(1,1):vshape(2,1) &
00404              ,vshape(1,2):vshape(2,2),j,i)=user(vshape(1,1):vshape(2,1) &
00405                                                ,vshape(1,2):vshape(2,2), i,j)
00406         enddo
00407       enddo
00408       
00409       end subroutine trs_bundle_vec_2d_integer
00410 
00411       subroutine trs_bundle_vec_3d_integer(user,xio,shape,vshape,ierror)
00412       integer,intent(out)         :: ierror
00413       integer,intent(in)      :: shape(2,5)
00414       integer,intent(in)      :: vshape(2,5)
00415       integer,intent(in) :: user(shape(1,1):shape(2,1) 
00416                                          ,shape(1,2):shape(2,2) 
00417                                          ,shape(1,3):shape(2,3) 
00418                                          ,shape(1,4):shape(2,4) 
00419                                          ,shape(1,5):shape(2,5))
00420       integer,intent(out)::  xio(shape(1,1):shape(2,1) 
00421                                          ,shape(1,2):shape(2,2) 
00422                                          ,shape(1,3):shape(2,2) 
00423                                          ,shape(1,5):shape(2,5) 
00424                                          ,shape(1,4):shape(2,4))
00425 !
00426 !     Local
00427 !
00428       integer                     :: i,j
00429       
00430       ierror=0
00431 
00432       do j=shape(1,5),shape(2,5)
00433         do i=shape(1,4),shape(2,4)
00434 
00435           xio(vshape(1,1):shape(2,1) &
00436              ,vshape(1,2):shape(2,2) &
00437              ,vshape(1,3):shape(2,3),j,i)=user(vshape(1,1):vshape(2,1) &
00438                                               ,vshape(1,2):shape(2,2) &
00439                                               ,vshape(1,3):shape(2,3) , i,j)
00440         enddo
00441       enddo
00442       
00443       end subroutine trs_bundle_vec_3d_integer