psmile_write_byid_dble.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_dble
00008 !
00009 ! !INTERFACE:
00010       subroutine psmile_write_byid_dble(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_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 (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       double precision,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 double precision 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_dble.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_dble: 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_dble: 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_dble: 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 
00137       il_varid=id_varid
00138       il_method_id=Fields(id_varid)%method_id
00139       il_grid_id=Methods(il_method_id)%grid_id
00140       il_grid_type=Grids(il_grid_id)%grid_type
00141       il_ndim=Grids(il_grid_id)%n_dim
00142       il_blockid=0
00143       is_block=.false.
00144       dl_ju_day=ju_day
00145       dl_ju_sec=ju_sec
00146 !
00147 !     Time offset with respect to the date of the job start
00148 !
00149       dl_time=86400.d0 *(dl_ju_day - Fields(id_varid)%io_infos%ju_start_day) &
00150              +          (dl_ju_sec - Fields(id_varid)%io_infos%ju_start_sec)
00151 !
00152 !     Convert Julian date to Gregorian date in order to construct a file name
00153 !     if a new file has to be created due to an exceed of the file size
00154 !     threshold.
00155 !
00156       call psmile_ju2date ( cur_date, dl_ju_day, dl_ju_sec )
00157 !   
00158 !     Check file size and open transparently a new file if necessary.
00159 !
00160       call psmile_open_file_byid(id_varid,id_taskid,cur_date,ierror)
00161 !
00162 !     Vector field ?
00163 !
00164       nvcomp=1
00165       vectorfield=0
00166       if(associated(Methods(il_method_id)%vector_pointer)) then
00167         nvcomp=size(Methods(il_method_id)%vector_pointer%array_of_point_ids)
00168         vectorfield=1
00169       endif
00170 
00171 !
00172 !     Set up shapes depending on the grid type.
00173 !
00174       if ( il_grid_type == PRISM_unstructlonlatvrt) then
00175          len = 1
00176          myvarshape(1:2,1)=Fields(id_varid)%var_shape(1:2,1)
00177          myvarshape(1:2,2)=1
00178          myvarshape(1:2,3)=1
00179          mygrdshape(1:2,1)=Grids(il_grid_id)%grid_shape(1:2,1)
00180          mygrdshape(1:2,2)=1
00181          mygrdshape(1:2,3)=1
00182       else if ( il_grid_type == PRISM_unstructlonlat_sigmavrt .or. &
00183             il_grid_type == PRISM_unstructlonlat_regvrt   .or. &
00184             il_grid_type == PRISM_Gaussreduced_regvrt   .or. &
00185             il_grid_type == PRISM_Gaussreduced_sigmavrt   ) then
00186          len = 3
00187          myvarshape(1:2,1)=Fields(id_varid)%var_shape(1:2,1)
00188          myvarshape(1:2,2)=1
00189          myvarshape(1:2,3)=Fields(id_varid)%var_shape(1:2,3)
00190 
00191          mygrdshape(1:2,1)=Grids(il_grid_id)%grid_shape(1:2,1)
00192          mygrdshape(1:2,2)=1
00193          mygrdshape(1:2,3)=Grids(il_grid_id)%grid_shape(1:2,3)
00194       else
00195          len = il_ndim
00196          myvarshape(1:2,1:len)=Fields(id_varid)%var_shape(1:2,1:len)
00197          mygrdshape(1:2,1:len)=Grids(il_grid_id)%grid_shape(1:2,1:len)
00198       endif
00199 !
00200 !     Is it a bundle ?
00201 !
00202        
00203       lenb=Fields(id_varid)%var_shape(2,len+vectorfield+1) &
00204          -Fields(id_varid)%var_shape(1,len+vectorfield+1)+1
00205 
00206       if(lenb.gt.1) then
00207         is_bundle=.true.
00208         myvarshape(1:2,4)=Fields(id_varid)%var_shape(1:2,len+vectorfield+1)
00209         mygrdshape(1:2,4)=Fields(id_varid)%var_shape(1:2,len+vectorfield+1)
00210       else ! No bundle
00211         is_bundle=.false.
00212       endif
00213 !
00214 !     Calculate  offset to address vector components.
00215 !     Components are written as separate NetCDF variables with 3 axis
00216 !     declared (Sophie Valckes Nov. 2003)
00217 !
00218       offset=1
00219       do ii = 1, len
00220          offset = offset * (Fields(id_varid)%var_shape(2,ii) &
00221                 - Fields(id_varid)%var_shape(1,ii)+1)
00222       enddo
00223 !
00224 !    Set the first location in the user provided field
00225 !
00226 
00227       il_loc=1
00228 !
00229 !     For GME like grids with a roundrobin distribution of fields over
00230 !     blocks I match the given varid with the set of related ids.
00231 !     The position of a match will be used as a block id.
00232 !     The file and field  information will be take from the last
00233 !     varid of the set of related ids because each block of the GME grid
00234 !     is partitioned in the same manner.
00235 !
00236       if(associated(Fields(id_varid)%io_infos%related_ids)) then
00237         fullsize=size(Fields(id_varid)%io_infos%related_ids)
00238         if(fullsize.gt.1) then
00239           do il_i=1,fullsize
00240             if(id_varid.eq. Fields(id_varid)%io_infos%related_ids(il_i)) exit
00241           enddo
00242 !
00243 !         The MPP_IO file descriptor is always taken from the first entry
00244 !         of the sets of related ids.
00245 !
00246           il_varid=Fields(id_varid)%io_infos%related_ids(1)
00247 !
00248 !         Take the block number from the matching entry within the set
00249 !         of related ids. 
00250 !
00251           il_blockid=Fields(il_i)%io_infos%block_id
00252           is_block=.true.
00253         endif
00254       endif 
00255 
00256       if(is_bundle) then
00257         offset=offset*lenb
00258         if(nvcomp.gt.1) then
00259           fullsize=offset*nvcomp
00260           allocate(afield(1:fullsize),stat=ierror)
00261           if (ierror > 0) then
00262             ierrp (1) = ierror
00263             ierrp (2) = 1
00264             ierror = PRISM_Error_Alloc
00265             call psmile_error ( ierror, 'afield', &
00266                          ierrp, 2, __FILE__, __LINE__ )
00267             return
00268           endif
00269 
00270 !
00271 !         Bundles appear as the last set of indices
00272 !         In order to write  individual vector components as bundles the
00273 !         vector indices have to be last dimension.
00274 !         Therefore, a transpose of the last two indices has to take place.
00275 
00276           if(len.eq.1) call trs_bundle_vec_1d_dble(dd_a,afield &
00277                                                   ,Fields(id_varid)%var_shape &
00278                                                   ,mygrdshape &
00279                                                   ,ierror)
00280           if(len.eq.2) call trs_bundle_vec_2d_dble(dd_a,afield &
00281                                                   ,Fields(id_varid)%var_shape &
00282                                                   ,mygrdshape &
00283                                                   ,ierror)
00284           if(len.eq.3) call trs_bundle_vec_3d_dble(dd_a,afield &
00285                                                   ,Fields(id_varid)%var_shape &
00286                                                   ,mygrdshape &
00287                                                   ,ierror)
00288         endif
00289       endif
00290 
00291 !
00292 !     Dump it
00293 !
00294       
00295       do ii=1,nvcomp
00296         
00297         if(.not.is_bundle) then
00298         call psmile_write_3d_dble(Fields(il_varid)%io_infos%file_unit &
00299                                  ,Fields(il_varid)%io_infos%p_mpp_io%field(ii) &
00300                                  ,Fields(il_varid)%io_infos%p_mpp_io%domain(1) &
00301                                  ,dd_a(il_loc) &
00302                                  ,myvarshape &
00303                                  ,mygrdshape &
00304                                  ,dl_time,.TRUE.,il_blockid,is_block,ierror)
00305         else if (is_bundle) then
00306           if(vectorfield.eq.0) then
00307             call psmile_write_4d_dble(Fields(il_varid)%io_infos%file_unit &
00308                                  ,Fields(il_varid)%io_infos%p_mpp_io%field(ii) &
00309                                  ,Fields(il_varid)%io_infos%p_mpp_io%domain(1) &
00310                                      ,dd_a(il_loc) &
00311                                      ,myvarshape &
00312                                      ,mygrdshape &
00313                                      ,dl_time,.TRUE.,il_blockid,is_block,ierror)
00314           else
00315             call psmile_write_4d_dble(Fields(il_varid)%io_infos%file_unit &
00316                                  ,Fields(il_varid)%io_infos%p_mpp_io%field(ii) &
00317                                  ,Fields(il_varid)%io_infos%p_mpp_io%domain(1) &
00318                                      ,afield(il_loc) &
00319                                      ,myvarshape &
00320                                      ,mygrdshape &
00321                                      ,dl_time,.TRUE.,il_blockid,is_block,ierror)
00322           endif
00323         endif
00324         il_loc=il_loc+offset
00325       enddo
00326 !
00327 !     Update the file size in kbytes (sizeof(double)/1024)=0.0078125)
00328 !     Note that the file size is here an rough estimate only. 
00329 !
00330       Fields(il_varid)%io_infos%old_filesize= &
00331             Fields(il_varid)%io_infos%current_filesize
00332       Fields(il_varid)%io_infos%current_filesize= &
00333             Fields(il_varid)%io_infos%current_filesize+max(0.0078125 &
00334                                                        *offset*nvcomp,1.)
00335 #ifdef VERBOSE
00336       print*,trim(ch_id),' : psmile_write_byid_dble: end'
00337       call psmile_flushstd
00338 #endif
00339 #endif
00340 
00341       end subroutine psmile_write_byid_dble
00342 
00343       subroutine trs_bundle_vec_1d_dble(user,xio,shape,vshape,ierror)
00344       integer,intent(out)         :: ierror
00345       integer,intent(in)      :: shape(2,3)
00346       integer,intent(in)      :: vshape(2,3)
00347       double precision,intent(in) :: user(shape(1,1):shape(2,1) 
00348                                          ,shape(1,2):shape(2,2) 
00349                                          ,shape(1,3):shape(2,3))
00350       double precision,intent(out)::  xio(shape(1,1):shape(2,1) 
00351                                          ,shape(1,3):shape(2,3) 
00352                                          ,shape(1,2):shape(2,2))
00353 !
00354 !     Local
00355 !
00356       integer                     :: i,j
00357       
00358 !
00359 !     The statement below is semantically the same as the manually coded loops.
00360 !     However, don't know what is faster. For now I trust my coding.
00361 !
00362 !      xio(shape(1,1):shape(2,1)                                           &
00363 !         ,shape(1,3):shape(2,3)                                           &
00364 !         ,shape(1,2):shape(2,2))=reshape(user,(/shape(2,1)-shape(1,1)+1    &
00365 !                                               ,shape(2,3)-shape(1,3)+1    &
00366 !                                               ,shape(2,1)-shape(1,2)+1 /) &
00367 !                                             ,ORDER=(/1,3,2/))
00368       ierror=0
00369 
00370       do j=shape(1,3),shape(2,3)
00371         do i=shape(1,2),shape(2,2)
00372           xio(vshape(1,1):vshape(2,1),j,i)=user(vshape(1,1):vshape(2,1),i,j)
00373         enddo
00374       enddo
00375       
00376       end subroutine trs_bundle_vec_1d_dble
00377 
00378       subroutine trs_bundle_vec_2d_dble(user,xio,shape,vshape,ierror)
00379       integer,intent(out)         :: ierror
00380       integer,intent(in)      :: shape(2,4)
00381       integer,intent(in)      :: vshape(2,4)
00382       double precision,intent(in) :: user(shape(1,1):shape(2,1) 
00383                                          ,shape(1,2):shape(2,2) 
00384                                          ,shape(1,3):shape(2,3) 
00385                                          ,shape(1,4):shape(2,4))
00386       double precision,intent(out)::  xio(shape(1,1):shape(2,1) 
00387                                          ,shape(1,2):shape(2,2) 
00388                                          ,shape(1,4):shape(2,4) 
00389                                          ,shape(1,3):shape(2,3))
00390 !
00391 !     Local
00392 !
00393       integer                     :: i,j
00394       
00395       ierror=0
00396 
00397       do j=shape(1,4),shape(2,4)
00398         do i=shape(1,3),shape(2,3)
00399 
00400           xio(vshape(1,1):vshape(2,1) &
00401              ,vshape(1,2):vshape(2,2),j,i)=user(vshape(1,1):vshape(2,1) &
00402                                                ,vshape(1,2):vshape(2,2), i,j)
00403         enddo
00404       enddo
00405       
00406       end subroutine trs_bundle_vec_2d_dble
00407 
00408       subroutine trs_bundle_vec_3d_dble(user,xio,shape,vshape,ierror)
00409       integer,intent(out)         :: ierror
00410       integer,intent(in)      :: shape(2,5)
00411       integer,intent(in)      :: vshape(2,5)
00412       double precision,intent(in) :: user(shape(1,1):shape(2,1) 
00413                                          ,shape(1,2):shape(2,2) 
00414                                          ,shape(1,3):shape(2,3) 
00415                                          ,shape(1,4):shape(2,4) 
00416                                          ,shape(1,5):shape(2,5))
00417       double precision,intent(out)::  xio(shape(1,1):shape(2,1) 
00418                                          ,shape(1,2):shape(2,2) 
00419                                          ,shape(1,3):shape(2,2) 
00420                                          ,shape(1,5):shape(2,5) 
00421                                          ,shape(1,4):shape(2,4))
00422 !
00423 !     Local
00424 !
00425       integer                     :: i,j
00426       
00427       ierror=0
00428 
00429       do j=shape(1,5),shape(2,5)
00430         do i=shape(1,4),shape(2,4)
00431 
00432           xio(vshape(1,1):vshape(2,1) &
00433              ,vshape(1,2):vshape(2,2) &
00434              ,vshape(1,3):vshape(2,3),j,i)=user(vshape(1,1):vshape(2,1) &
00435                                                ,vshape(1,2):vshape(2,2) &
00436                                                ,vshape(1,3):vshape(2,3) , i,j)
00437         enddo
00438       enddo
00439       
00440       end subroutine trs_bundle_vec_3d_dble