doc-doxygen/html/mpp__read__2_ddecomp_8h_source.html

00001 #ifdef use_CRI_pointers
00002 #undef use_local_CRI_pointers
00003 #endif
00004     subroutine MPP_READ_2DDECOMP_1D_(unit,field,domain,data,tindex,blockid )
00005       integer, intent(in) :: unit
00006       type(fieldtype), intent(in) :: field
00007       type(domain2D), intent(in) :: domain
00008       MPP_TYPE_, intent(inout) :: data(:)
00009       integer, intent(in), optional :: tindex
00010       integer, intent(in), optional :: blockid
00011       integer::il_xbegin,il_xend,il_ybegin,il_yend
00012       MPP_TYPE_,allocatable :: data3D(:,:,:)
00013
00014       call mpp_get_compute_domain(domain &
00015                                 ,xbegin=il_xbegin,xend=il_xend &
00016                                 ,ybegin=il_ybegin,yend=il_yend)
00017
00018       allocate(data3D(1:il_xend-il_xbegin+1,1:il_yend-il_ybegin+1,1:1))
00019       data3D = RESHAPE( data, SHAPE(data3D) )
00020       if(PRESENT(blockid)) then
00021         call mpp_read( unit,field,domain,data3D,tindex,blockid)
00022       else
00023         call mpp_read( unit, field, domain, data3D, tindex )
00024       endif
00025       data   = RESHAPE( data3D, SHAPE(data) )
00026       deallocate(data3D)
00027
00028       return
00029     end subroutine MPP_READ_2DDECOMP_1D_
00030     subroutine MPP_READ_2DDECOMP_2D_( unit, field, domain,data,tindex,blockid )
00031       integer, intent(in) :: unit
00032       type(fieldtype), intent(in) :: field
00033       type(domain2D), intent(in) :: domain
00034       MPP_TYPE_, intent(inout) :: data(:,:)
00035       integer, intent(in), optional :: tindex
00036       integer, intent(in), optional :: blockid
00037       MPP_TYPE_ :: data3D(size(data,1),size(data,2),1)
00038 #ifdef use_local_CRI_pointers
00039       pointer( ptr, data3D )
00040       ptr = LOC(data)
00041       if(PRESENT(blockid)) then
00042       call mpp_read( unit, field, domain, data3D, tindex, blockid )
00043       else
00044       call mpp_read( unit, field, domain, data3D, tindex )
00045       endif
00046 #else
00047       data3D = RESHAPE( data, SHAPE(data3D) )
00048       if(PRESENT(blockid)) then
00049       call mpp_read( unit, field, domain, data3D, tindex, blockid )
00050       else
00051       call mpp_read( unit, field, domain, data3D, tindex )
00052       endif
00053       data   = RESHAPE( data3D, SHAPE(data) )
00054 #endif
00055       return
00056     end subroutine MPP_READ_2DDECOMP_2D_
00057
00058     subroutine MPP_READ_2DDECOMP_3D_(unit,field,domain,data,tindex,blockid )
00059 !mpp_read reads <data> which has the domain decomposition <domain>
00060       integer, intent(in) :: unit
00061       type(fieldtype), intent(in) :: field
00062       type(domain2D), intent(in) :: domain
00063       MPP_TYPE_, intent(inout) :: data(:,:,:)
00064       integer, intent(in), optional :: tindex
00065       integer, intent(in), optional :: blockid
00066       MPP_TYPE_, allocatable :: cdata(:,:,:)
00067       MPP_TYPE_, allocatable :: gdata(:)
00068       integer :: len, lenx,leny,lenz,i,j,k,n
00069 !NEW: data may be on compute OR data domain
00070       logical :: data_has_halos, halos_are_global, x_is_global, y_is_global
00071       integer :: is, ie, js, je, isd, ied, jsd, jed, isg, ieg, jsg, jeg, ioff, joff
00072
00073       if (.NOT. present(tindex) .AND. mpp_file(unit)%time_level .ne. -1) &
00074       call mpp_error(FATAL, 'MPP_READ: need to specify a time level for data with time axis')
00075
00076       if( .NOT.module_is_initialized )call mpp_error( FATAL, 'MPP_READ: must first call mpp_io_init.' )
00077       if( .NOT.mpp_file(unit)%opened )call mpp_error( FATAL, 'MPP_READ: invalid unit number.' )
00078
00079       call mpp_get_compute_domain( domain, is,  ie,  js,  je  )
00080       call mpp_get_data_domain   ( domain, isd, ied, jsd, jed, x_is_global=x_is_global, y_is_global=y_is_global )
00081       call mpp_get_global_domain ( domain, isg, ieg, jsg, jeg )
00082       if( debug ) &
00083       write(stdout(),*) 'MPP_READ_2DDECOMP_3D_:is,  ie,  js,  je, isd, ied, jsd, jed',is,  ie,  js,  je, isd, ied, jsd, jed
00084       if( debug ) &
00085       write(stdout(),*) 'MPP_READ_2DDECOMP_3D_:isg, ieg, jsg, jeg',isg, ieg, jsg, jeg
00086       if( size(data,1).EQ.ie-is+1 .AND. size(data,2).EQ.je-js+1 )then
00087           data_has_halos = .FALSE.
00088       else if( size(data,1).EQ.ied-isd+1 .AND. size(data,2).EQ.jed-jsd+1 )then
00089           data_has_halos = .TRUE.
00090       else
00091           call mpp_error( FATAL, 'MPP_READ: data must be either on compute domain or data domain.' )
00092       end if
00093       halos_are_global = x_is_global .AND. y_is_global
00094       if( debug ) &
00095       write(stdout(),*) 'MPP_READ_2DDECOMP_3D_: halos_are_global', halos_are_global
00096       if( npes.GT.1 .AND. mpp_file(unit)%threading.EQ.MPP_SINGLE )then
00097       if( debug ) &
00098       write(stdout(),*) 'MPP_READ_2DDECOMP_3D_: threading.EQ.MPP_SINGLE'
00099           if( halos_are_global )then !you can read directly into data array
00100               if(PRESENT(blockid)) then
00101               if( pe.EQ.mpp_root_pe() )call read_record_b(unit,field,size(data) &
00102                                            ,data,tindex,block_id=blockid )
00103               else
00104               if(pe.EQ.mpp_root_pe())call read_record(unit,field,size(data),data,tindex )
00105               endif
00106           else
00107               lenx=ieg-isg+1
00108               leny=jeg-jsg+1
00109               lenz=size(data,3)
00110               len=lenx*leny*lenz
00111               write(stdout(),*)'MPP_READ_2DDECOMP_3D: gdata ', lenx,leny,lenz
00112               allocate(gdata(len))
00113 ! read field on pe 0 and pass to all pes
00114               if(PRESENT(blockid)) then
00115               if(pe.EQ.mpp_root_pe()) call read_record_b(unit,field,len,gdata,tindex,block_id=blockid)
00116               else
00117               if( pe.EQ.mpp_root_pe() ) call read_record( unit, field, len, gdata, tindex )
00118               endif
00119 ! broadcasting global array, this can be expensive!
00120               call mpp_transmit( gdata, len, ALL_PES, gdata, len, 0 )
00121               ioff = is; joff = js
00122               if( data_has_halos )then
00123                   ioff = isd; joff = jsd
00124               end if
00125               do k=1,size(data,3)
00126                  do j=js,je
00127                     do i=is,ie
00128                        n=(i-isg+1) + (j-jsg)*lenx + (k-1)*lenx*leny
00129                        data(i-ioff+1,j-joff+1,k)=gdata(n)
00130                     enddo
00131                  enddo
00132               enddo
00133               deallocate(gdata)
00134           end if
00135       else if( data_has_halos )then
00136 ! for uniprocessor or multithreaded read
00137 ! read compute domain as contiguous data
00138
00139           allocate( cdata(is:ie,js:je,size(data,3)) )
00140           if(PRESENT(blockid)) then
00141           call read_record_b(unit,field,size(cdata),cdata,tindex,domain,block_id=blockid)
00142           else
00143           call read_record(unit,field,size(cdata),cdata,tindex,domain)
00144           endif
00145
00146           data(is-isd+1:ie-isd+1,js-jsd+1:je-jsd+1,:) = cdata(:,:,:)
00147           deallocate(cdata)
00148       else
00149           if(PRESENT(blockid)) then
00150           call read_record_b(unit,field,size(data),data,tindex,domain,block_id=blockid)
00151           else
00152           call read_record(unit,field,size(data),data,tindex,domain)
00153           endif
00154       end if
00155
00156       return
00157     end subroutine MPP_READ_2DDECOMP_3D_
00158
00159     subroutine MPP_READ_2DDECOMP_4D_(unit,field,domain,data,tindex,blockid )
00160 !mpp_read reads <data> which has the domain decomposition <domain>
00161       integer, intent(in) :: unit
00162       type(fieldtype), intent(in) :: field
00163       type(domain2D), intent(in) :: domain
00164       MPP_TYPE_, intent(inout) :: data(:,:,:,:)
00165       integer, intent(in), optional :: tindex
00166       integer, intent(in), optional :: blockid
00167       MPP_TYPE_, allocatable :: cdata(:,:,:,:)
00168       MPP_TYPE_, allocatable :: gdata(:)
00169       integer :: len, lenx,leny,lenz,i,j,k,n,l
00170 !NEW: data may be on compute OR data domain
00171       logical :: data_has_halos, halos_are_global, x_is_global, y_is_global
00172       integer :: is, ie, js, je, isd, ied, jsd, jed, isg, ieg, jsg, jeg, ioff, joff
00173
00174       if (.NOT. present(tindex) .AND. mpp_file(unit)%time_level .ne. -1) &
00175       call mpp_error(FATAL, 'MPP_READ: need to specify a time level for data with time axis')
00176
00177       if( .NOT.module_is_initialized )call mpp_error( FATAL, 'MPP_READ: must first call mpp_io_init.' )
00178       if( .NOT.mpp_file(unit)%opened )call mpp_error( FATAL, 'MPP_READ: invalid unit number.' )
00179
00180       call mpp_get_compute_domain( domain, is,  ie,  js,  je  )
00181       call mpp_get_data_domain   ( domain, isd, ied, jsd, jed, x_is_global=x_is_global, y_is_global=y_is_global )
00182       call mpp_get_global_domain ( domain, isg, ieg, jsg, jeg )
00183       if( size(data,1).EQ.ie-is+1 .AND. size(data,2).EQ.je-js+1 )then
00184           data_has_halos = .FALSE.
00185       else if( size(data,1).EQ.ied-isd+1 .AND. size(data,2).EQ.jed-jsd+1 )then
00186           data_has_halos = .TRUE.
00187       else
00188           call mpp_error( FATAL, 'MPP_READ: data must be either on compute domain or data domain.' )
00189       end if
00190       halos_are_global = x_is_global .AND. y_is_global
00191       if( npes.GT.1 .AND. mpp_file(unit)%threading.EQ.MPP_SINGLE )then
00192           if( halos_are_global )then !you can read directly into data array
00193               if(PRESENT(blockid)) then
00194               if( pe.EQ.mpp_root_pe() )call read_record_b(unit,field,size(data) &
00195                                            ,data,tindex,block_id=blockid )
00196               else
00197               if(pe.EQ.mpp_root_pe())call read_record(unit,field,size(data),data,tindex )
00198               endif
00199           else
00200               lenx=ieg-isg+1
00201               leny=jeg-jsg+1
00202               lenz=size(data,3)
00203               len=lenx*leny*lenz*size(data,4)
00204               allocate(gdata(len))
00205 ! read field on pe 0 and pass to all pes
00206               if(PRESENT(blockid)) then
00207               if(pe.EQ.mpp_root_pe()) call read_record_b(unit,field,len,gdata,tindex,block_id=blockid)
00208               else
00209               if( pe.EQ.mpp_root_pe() ) call read_record( unit, field, len, gdata, tindex )
00210               endif
00211 ! broadcasting global array, this can be expensive!
00212               call mpp_transmit( gdata, len, ALL_PES, gdata, len, 0 )
00213               ioff = is; joff = js
00214               if( data_has_halos )then
00215                   ioff = isd; joff = jsd
00216               end if
00217               do l=1,size(data,4)
00218                 do k=1,size(data,3)
00219                   do j=js,je
00220                     do i=is,ie
00221                        n=(i-isg+1)+(j-jsg)*lenx+(k-1)*lenx*leny &
00222                         +(l-1)*lenx*leny*lenz
00223                        data(i-ioff+1,j-joff+1,k,l)=gdata(n)
00224                     enddo
00225                   enddo
00226                 enddo
00227               enddo
00228               deallocate(gdata)
00229           end if
00230       else if( data_has_halos )then
00231 ! for uniprocessor or multithreaded read
00232 ! read compute domain as contiguous data
00233
00234           allocate( cdata(is:ie,js:je,size(data,3),size(data,4)) )
00235           if(PRESENT(blockid)) then
00236           call read_record_b(unit,field,size(cdata),cdata,tindex,domain,block_id=blockid)
00237           else
00238           call read_record(unit,field,size(cdata),cdata,tindex,domain)
00239           endif
00240
00241           data(is-isd+1:ie-isd+1,js-jsd+1:je-jsd+1,:,:) = cdata(:,:,:,:)
00242           deallocate(cdata)
00243       else
00244           if(PRESENT(blockid)) then
00245           call read_record_b(unit,field,size(data),data,tindex,domain,block_id=blockid)
00246           else
00247           call read_record(unit,field,size(data),data,tindex,domain)
00248           endif
00249       end if
00250
00251       return
00252     end subroutine MPP_READ_2DDECOMP_4D_
00253 #define use_local_CRI_pointers