/[MITgcm]/MITgcm/pkg/exf/exf_interp.F
ViewVC logotype

Diff of /MITgcm/pkg/exf/exf_interp.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.5 by dimitri, Wed Dec 10 19:37:25 2003 UTC revision 1.10 by heimbach, Mon May 2 15:31:46 2005 UTC
# Line 58  C Line 58  C
58    
59  #include "SIZE.h"  #include "SIZE.h"
60  #include "EEPARAMS.h"  #include "EEPARAMS.h"
61  #ifdef ALLOW_USE_MPI  #include "PARAMS.h"
 # include "EESUPPORT.h"  
 # include "PARAMS.h"  
 #endif /* ALLOW_USE_MPI */  
62    
63  C subroutine variables  C subroutine variables
64        character*(*) infile        character*(*) infile
# Line 74  C subroutine variables Line 71  C subroutine variables
71        integer       method, mythid        integer       method, mythid
72    
73  C local variables  C local variables
       integer  ierr  
74        real*8   ne_fac,nw_fac,se_fac,sw_fac        real*8   ne_fac,nw_fac,se_fac,sw_fac
75        integer  e_ind(snx,sny),w_ind(snx,sny)        integer  e_ind(snx,sny),w_ind(snx,sny)
76        integer  n_ind(snx,sny),s_ind(snx,sny)        integer  n_ind(snx,sny),s_ind(snx,sny)
# Line 84  C local variables Line 80  C local variables
80        real*4   arrayin(-1:nx_in+2 ,      -1:ny_in+2)        real*4   arrayin(-1:nx_in+2 ,      -1:ny_in+2)
81        real*8   x_in   (-1:nx_in+2), y_in(-1:ny_in+2)        real*8   x_in   (-1:nx_in+2), y_in(-1:ny_in+2)
82        integer  i, j, k, l, js, bi, bj, sp, interp_unit        integer  i, j, k, l, js, bi, bj, sp, interp_unit
       real*4   global(nx_in,ny_in)  
83    
84        _BEGIN_MASTER( myThid )         call exf_interp_read(
85         I   infile,
86  C check input arguments       I   filePrec,
87         if ( .NOT. (filePrec .EQ. 32) )       O   arrayin,
88       &     stop 'stop in exf_interp.F: value of filePrec not allowed'       I   irecord, xG, yG,
89         I   lon_0, lon_inc,
90  C read in input data       I   lat_0, lat_inc,
91  #ifdef ALLOW_USE_MPI       I   nx_in, ny_in, method, mythid)
        if (useSingleCPUIO) then  
   
 C master thread of process 0, only, opens a global file  
         IF( mpiMyId .EQ. 0 ) THEN  
          call mdsfindunit( interp_unit, mythid)  
          open(interp_unit,file=infile,status='old',access='direct',  
      &        recl=nx_in*ny_in*4)  
          read(interp_unit,rec=irecord)  
      &        ((global(i,j),i=1,nx_in),j=1,ny_in)  
          close(interp_unit)  
         ENDIF  
   
 C broadcast to all processes  
         call MPI_BCAST(global,nx_in*ny_in,MPI_REAL,  
      &       0,MPI_COMM_MODEL,ierr)  
         do j=1,ny_in  
          do i=1,nx_in  
           arrayin(i,j)=global(i,j)  
          enddo  
         enddo  
92    
93         else        _BEGIN_MASTER( myThid )
 #endif /* ALLOW_USE_MPI */  
   
         call mdsfindunit( interp_unit, mythid)  
         open(interp_unit,file=infile,status='old',access='direct',  
      &       recl=nx_in*ny_in*4)  
         read(interp_unit,rec=irecord)  
      &       ((arrayin(i,j),i=1,nx_in),j=1,ny_in)  
         close(interp_unit)  
   
 #ifdef ALLOW_USE_MPI  
        endif  
 #endif /* ALLOW_USE_MPI */  
   
 #ifdef _BYTESWAPIO  
        call MDS_BYTESWAPR4((nx_in+4)*(ny_in+4), arrayin )  
 #endif /* _BYTESWAPIO */  
94    
95  C setup input grid  C setup input grid
96         do i=-1,nx_in+2         do i=-1,nx_in+2
# Line 166  C enlarge boundary Line 125  C enlarge boundary
125         do bi = mybxlo(mythid), mybxhi(mythid)         do bi = mybxlo(mythid), mybxhi(mythid)
126    
127  C check validity of input/output coordinates  C check validity of input/output coordinates
128          if ( xG(1,1  ,bi,bj) .le. x_in(0)       .or.  #ifdef ALLOW_DEBUG
129       &       xG(snx,1,bi,bj) .ge. x_in(nx_in+1) .or.          if ( debugLevel .ge. debLevB ) then
130       &       yG(1,1  ,bi,bj) .lt. y_in(1)       .or.           do i=1,snx
131       &       yG(1,sny,bi,bj) .gt. y_in(ny_in) ) then            do j=1,sny
132             print*,'ERROR in S/R EXF_INTERP:'             if ( xG(i,j,bi,bj) .lt. x_in(0)         .or.
133             print*,'   input grid must encompass output grid.'       &          xG(i,j,bi,bj) .ge. x_in(nx_in+1)   .or.
134             STOP   ' ABNORMAL END: S/R EXF_INTERP'       &          yG(i,j,bi,bj) .lt. y_in(0)         .or.
135         &          yG(i,j,bi,bj) .ge. y_in(ny_in+1) ) then
136                  print*,'ERROR in S/R EXF_INTERP:'
137                  print*,'   input grid must encompass output grid.'
138                  print*,'i,j,bi,bj'      ,i,j,bi,bj
139                  print*,'xG,yG'          ,xG(i,j,bi,bj),yG(i,j,bi,bj)
140                  print*,'nx_in,ny_in'    ,nx_in        ,ny_in
141                  print*,'x_in(0,nx_in+1)',x_in(0)      ,x_in(nx_in+1)
142                  print*,'y_in(0,ny_in+1)',y_in(0)      ,y_in(ny_in+1)
143                  STOP   ' ABNORMAL END: S/R EXF_INTERP'
144                 endif
145               enddo
146             enddo
147          endif          endif
148    #endif /* ALLOW_DEBUG */
149    
150  C compute interpolation indices  C compute interpolation indices
151          do i=1,snx          do i=1,snx
# Line 184  C compute interpolation indices Line 156  C compute interpolation indices
156             w_ind(i,j) = int((xG(i,j,bi,bj)-x_in(1))/lon_inc)             w_ind(i,j) = int((xG(i,j,bi,bj)-x_in(1))/lon_inc)
157            endif            endif
158            e_ind(i,j) = w_ind(i,j) + 1            e_ind(i,j) = w_ind(i,j) + 1
159            js = ny_in/2            js = ny_in*.5
160            do while (yG(i,j,bi,bj) .lt. y_in(js))            do while (yG(i,j,bi,bj) .lt. y_in(js))
161             js = (js + 1)/2             js = (js - 1)*.5
162            enddo            enddo
163            do while (yG(i,j,bi,bj) .ge. y_in(js+1))            do while (yG(i,j,bi,bj) .ge. y_in(js+1))
164             js = js + 1             js = js + 1
Legend:
Removed from v.1.5  
changed lines
  Added in v.1.10
  ViewVC Help
Powered by ViewVC 1.1.22