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