C $Header: /u/gcmpack/MITgcm/pkg/flt/flt_init_varia.F,v 1.13 2017/03/24 23:38:56 jmc Exp $ C $Name: $ #include "FLT_OPTIONS.h" SUBROUTINE FLT_INIT_VARIA ( myThid ) C ================================================================== C SUBROUTINE FLT_INIT_VARIA C ================================================================== C o This routine initializes the start/restart positions. C o Either read initial position from file "flt_file" or C read pickup file. The 2 type of files are similar, except C initial positions are given on grid-coordinate (distance/degree C depending on the grid) whereas in pickup file, positions are C fractional indices along the grid and local to the tile. C For this reason global pickup file is not supported. C Initialisation: C o First it check for global file, and when found, reads the global file C (that has the same format as local files) and sorts those floats C that exist on the specific tile into the local array. C o If no global file is available or in a case of a restart (pickup C file from a previous integration) then read tiled file without C any further check (because they exist on the specific tile). C ================================================================== C !USES: IMPLICIT NONE #include "SIZE.h" #include "EEPARAMS.h" #include "PARAMS.h" #include "FLT_SIZE.h" #include "FLT.h" C == routine arguments == C myThid - thread number for this instance of the routine. INTEGER myThid C == Functions == INTEGER ILNBLNK EXTERNAL _RL FLT_MAP_R2K EXTERNAL C == local variables == INTEGER bi, bj INTEGER ip, iL INTEGER imax PARAMETER(imax=9) _RL tmp(imax) _RS dummyRS(1) _RL ix, jy, kz _RL iLo, iHi, jLo, jHi INTEGER fp, ioUnit CHARACTER*(MAX_LEN_FNAM) fn CHARACTER*(MAX_LEN_MBUF) msgBuf CHARACTER*(10) suff #ifdef DEVEL_FLT_EXCH2 INTEGER iG,jG,i,j #endif C number of active record in the file (might be lower than the C total number of records because the tile could have contained C more floats at an earlier restart INTEGER npart_read _RL npart_dist C == end of interface == C- Tile boundary on index map: iLo = 0.5 _d 0 iHi = 0.5 _d 0 + DFLOAT(sNx) jLo = 0.5 _d 0 jHi = 0.5 _d 0 + DFLOAT(sNy) C- all threads initialise local var: npart_read = 0 npart_dist = 0. _BEGIN_MASTER(myThid) DO bj = 1,nSy DO bi = 1,nSx npart_tile(bi,bj) = 0 ENDDO ENDDO C read floats initial condition from file IF ( nIter0.EQ.FLT_Iter0 ) THEN fn = flt_file fp = readBinaryPrec ELSEIF ( nIter0.GT.FLT_Iter0 ) THEN IF ( pickupSuff .EQ. ' ' ) THEN IF ( rwSuffixType.EQ.0 ) THEN WRITE(suff,'(I10.10)') nIter0 ELSE CALL RW_GET_SUFFIX( suff, startTime, nIter0, myThid ) ENDIF ELSE WRITE(suff,'(A10)') pickupSuff ENDIF WRITE(fn,'(A,A10)') 'pickup_flt.',suff fp = precFloat64 ELSE WRITE(msgBuf,'(2A,I3,A)') 'FLT_INIT_VARIA:', & ' wrong setting of FLT_Iter0 :' CALL PRINT_ERROR( msgBuf, myThid ) WRITE(msgBuf,'(2A,I3,A)') 'FLT_INIT_VARIA:', & ' nIter0 < FLT_Iter0 not supported' CALL PRINT_ERROR( msgBuf, myThid ) STOP 'ABNORMAL END: S/R FLT_INIT_VARIA' ENDIF iL = ILNBLNK(fn) WRITE(msgBuf,'(2A)') & 'FLT_INIT_VARIA: reading Floats from: ', fn(1:iL) CALL PRINT_MESSAGE( msgBuf, standardMessageUnit, & SQUEEZE_RIGHT, myThid ) C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| C-- Initial position: first try to read from a global file. ioUnit = -2 bi = 0 bj = 0 IF ( nIter0.EQ.FLT_Iter0 ) THEN C- read actual number of floats from file CALL MDS_READVEC_LOC( fn, fp, ioUnit, & 'RL', imax, tmp, dummyRS, & bi, bj, 1, myThid ) ENDIF IF ( ioUnit.GT.0 .AND. mapIniPos2Index ) THEN C-- Found a global file WRITE(msgBuf,'(A,2I4,A,1P2E15.8)') & ' bi,bj=', bi, bj, ' , npart,max_npart=', tmp(1), tmp(6) CALL PRINT_MESSAGE( msgBuf, standardMessageUnit, & SQUEEZE_RIGHT, myThid ) npart_read = NINT(tmp(1)) max_npart = tmp(6) DO ip=1,npart_read C- read individual float position from file CALL MDS_READVEC_LOC( fn, fp, ioUnit, & 'RL', imax, tmp, dummyRS, & bi, bj, ip+1, myThid ) DO bj = 1,nSy DO bi = 1,nSx C- For initial condition only, convert coordinates to index map: CALL FLT_MAP_XY2IJLOCAL( ix, jy, I tmp(3), tmp(4),bi,bj,myThid ) kz = FLT_MAP_R2K( tmp(5), bi, bj, myThid ) C- Check if float exists on this tile. If not, try next tile IF ( ix.GE.iLo .AND. ix.LT.iHi .AND. & jy.GE.jLo .AND. jy.LT.jHi ) THEN npart_tile(bi,bj) = npart_tile(bi,bj) + 1 IF ( npart_tile(bi,bj).LE.max_npart_tile ) THEN npart( npart_tile(bi,bj),bi,bj) = tmp(1) tstart(npart_tile(bi,bj),bi,bj) = tmp(2) ipart( npart_tile(bi,bj),bi,bj) = ix jpart( npart_tile(bi,bj),bi,bj) = jy kpart( npart_tile(bi,bj),bi,bj) = kz kfloat(npart_tile(bi,bj),bi,bj) = tmp(6) iup( npart_tile(bi,bj),bi,bj) = tmp(7) itop( npart_tile(bi,bj),bi,bj) = tmp(8) tend( npart_tile(bi,bj),bi,bj) = tmp(9) ENDIF ENDIF C- end bi,bj loops ENDDO ENDDO ENDDO CLOSE( ioUnit ) ELSEIF ( ioUnit.GT.0 ) THEN WRITE(msgBuf,'(2A)') 'FLT_INIT_VARIA:', & ' need mapIniPos2Index=T for global file' CALL PRINT_ERROR( msgBuf , myThid) STOP 'ABNORMAL END: S/R FLT_INIT_VARIA' ELSE C-- then try to read from a tiled file: DO bj = 1,nSy DO bi = 1,nSx #ifdef DEVEL_FLT_EXCH2 iG=bi+(myXGlobalLo-1)/sNx jG=bj+(myYGlobalLo-1)/sNy WRITE(msgBuf,'(A,I4,A,I4)') 'iG,jG=',iG,',',jG CALL PRINT_MESSAGE( msgBuf, standardMessageUnit, & SQUEEZE_RIGHT, myThid ) c npart_tile(bi,bj) = 6*6 max_npart = 6*6*10 npart_read = 1 c DO j=1,6 DO i=1,6 ip=i+(j-1)*6 npart(ip,bi,bj) = (iG-1)*6*6+ip tstart(ip,bi,bj) = -1. _d 0 ipart(ip,bi,bj) = 1+5*(i-1) jpart(ip,bi,bj) = 1+5*(j-1) kpart(ip,bi,bj) = 1. _d 0 kfloat(ip,bi,bj) = 0. _d 0 iup( ip,bi,bj) = 0. _d 0 itop( ip,bi,bj) = 0. _d 0 tend( ip,bi,bj) = -1. _d 0 ENDDO ENDDO #else ioUnit = -1 C- read actual number floats from file CALL MDS_READVEC_LOC( fn, fp, ioUnit, & 'RL', imax, tmp, dummyRS, & bi, bj, 1, myThid ) WRITE(msgBuf,'(A,2I4,A,1P2E15.8)') & ' bi,bj=', bi, bj, ' , npart,max_npart=', tmp(1), tmp(6) CALL PRINT_MESSAGE( msgBuf, standardMessageUnit, & SQUEEZE_RIGHT, myThid ) npart_tile(bi,bj) = NINT(tmp(1)) max_npart = tmp(6) npart_read = MIN( npart_tile(bi,bj), max_npart_tile ) DO ip=1,npart_read C- read individual float position from file CALL MDS_READVEC_LOC( fn, fp, ioUnit, & 'RL', imax, tmp, dummyRS, & bi, bj, ip+1, myThid ) IF ( nIter0.EQ.FLT_Iter0 .AND. mapIniPos2Index ) THEN C-- For initial condition only, convert coordinates to index map: CALL FLT_MAP_XY2IJLOCAL( ix, jy, I tmp(3), tmp(4),bi,bj,myThid ) kz = FLT_MAP_R2K( tmp(5), bi, bj, myThid ) ELSE ix = tmp(3) jy = tmp(4) kz = tmp(5) ENDIF C not a global file: assume that all particles from this tiled-file C belong to this current tile (=> do not no check) npart(ip,bi,bj) = tmp(1) tstart(ip,bi,bj) = tmp(2) ipart(ip,bi,bj) = ix jpart(ip,bi,bj) = jy kpart(ip,bi,bj) = kz kfloat(ip,bi,bj) = tmp(6) iup( ip,bi,bj) = tmp(7) itop( ip,bi,bj) = tmp(8) tend( ip,bi,bj) = tmp(9) ENDDO CLOSE( ioUnit ) #endif C- end bi,bj loops ENDDO ENDDO C-- end global-file / tiled-file separated treatment ENDIF DO bj = 1,nSy DO bi = 1,nSx npart_dist = npart_dist + DBLE(npart_tile(bi,bj)) IF ( npart_tile(bi,bj).GT.max_npart_tile ) THEN WRITE(msgBuf,'(2A,2I4,2(A,I8))') 'FLT_INIT_VARIA:', & ' bi,bj=', bi, bj, & ' npart_tile=', npart_tile(bi,bj), & ' > max_npart_tile=', max_npart_tile CALL PRINT_ERROR( msgBuf , myThid) STOP 'ABNORMAL END: S/R FLT_INIT_VARIA' ENDIF ENDDO ENDDO _END_MASTER( myThid ) _BARRIER _GLOBAL_SUM_RL( npart_dist, myThid ) _BEGIN_MASTER( myThid ) WRITE(msgBuf,'(A,2(A,I9))') 'FLT_INIT_VARIA:', & ' max npart=', NINT(max_npart), & ' , sum npart_tile=', NINT(npart_dist) CALL PRINT_MESSAGE( msgBuf, standardMessageUnit, & SQUEEZE_RIGHT, myThid ) WRITE(msgBuf,'(A)') ' ' CALL PRINT_MESSAGE( msgBuf, standardMessageUnit, & SQUEEZE_RIGHT, myThid ) _END_MASTER( myThid ) C-- Initial call just to check which variables to write IF ( flt_int_prof.NE.0. ) & CALL FLT_UP( startTime, nIter0, myThid ) IF ( flt_int_traj.NE.0. ) & CALL FLT_TRAJ( startTime, nIter0, myThid ) RETURN END