--- MITgcm/model/src/ini_curvilinear_grid.F 2004/06/25 02:49:49 1.14 +++ MITgcm/model/src/ini_curvilinear_grid.F 2005/09/10 18:30:06 1.22 @@ -1,4 +1,4 @@ -C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/model/src/ini_curvilinear_grid.F,v 1.14 2004/06/25 02:49:49 dimitri Exp $ +C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/model/src/ini_curvilinear_grid.F,v 1.22 2005/09/10 18:30:06 edhill Exp $ C $Name: $ #include "PACKAGES_CONFIG.h" @@ -31,26 +31,39 @@ #include "W2_EXCH2_TOPOLOGY.h" #include "W2_EXCH2_PARAMS.h" #endif +#ifdef ALLOW_MNC +#include "MNC_PARAMS.h" +#endif #ifndef ALLOW_EXCH2 C- note: default is to use "new" grid files (OLD_GRID_IO undef) with EXCH2 C but can still use (on 1 cpu) OLD_GRID_IO and EXCH2 independently +#ifdef ALLOW_MDSIO #define OLD_GRID_IO +#endif #endif /* ALLOW_EXCH2 */ C !INPUT/OUTPUT PARAMETERS: C == Routine arguments == -C myThid - Number of this instance of INI_CARTESIAN_GRID +C myThid - Number of this instance of INI_CURVILINEAR_GRID INTEGER myThid C !LOCAL VARIABLES: C == Local variables == - INTEGER bi,bj, myTile, myiter + INTEGER bi,bj, myIter INTEGER I,J - CHARACTER*(15) fName + CHARACTER*(MAX_LEN_FNAM) fName +#ifdef ALLOW_MNC + CHARACTER*(80) mncFn +#endif +#ifdef ALLOW_EXCH2 + _RL buf(sNx*nSx*nPx+1) + INTEGER myTile +#else _RL buf(sNx+1,sNy+1) - INTEGER iG, iL - CHARACTER*(MAX_LEN_MBUF) msgBuf +#endif + INTEGER iG, iL, iLen + CHARACTER*(MAX_LEN_MBUF) msgBuf, tmpBuf INTEGER ILNBLNK EXTERNAL ILNBLNK CEOP @@ -79,6 +92,8 @@ RAS(i,j,bi,bj)=0. tanPhiAtU(i,j,bi,bj)=0. tanPhiAtV(i,j,bi,bj)=0. + angleCosC(i,j,bi,bj)=1. + angleSinC(i,j,bi,bj)=0. cosFacU(J,bi,bj)=1. cosFacV(J,bi,bj)=1. sqcosFacU(J,bi,bj)=1. @@ -89,6 +104,61 @@ ENDDO ENDDO + +#ifdef ALLOW_MNC + IF (useMNC .AND. readgrid_mnc) THEN + + _BEGIN_MASTER(myThid) + DO i = 1,80 + mncFn(i:i) = ' ' + ENDDO + write(mncFn,'(a)') 'mitgrid' + DO i = 1,MAX_LEN_MBUF + msgBuf(i:i) = ' ' + ENDDO + WRITE(msgBuf,'(2A)') msgBuf,' ; Reading grid info using MNC' + CALL PRINT_MESSAGE( msgBuf, standardMessageUnit, + & SQUEEZE_RIGHT , myThid) + CALL MNC_FILE_CLOSE_ALL_MATCHING(mncFn, myThid) + CALL MNC_CW_SET_UDIM(mncFn, 1, myThid) + CALL MNC_CW_SET_CITER(mncFn, 2, -1, -1, -1, myThid) + CALL MNC_CW_RS_R('R',mncFn,0,0,'XC', XC, myThid) + CALL MNC_CW_RS_R('R',mncFn,0,0,'XG', XG, myThid) + CALL MNC_CW_RS_R('R',mncFn,0,0,'YC', YC, myThid) + CALL MNC_CW_RS_R('R',mncFn,0,0,'YG', YG, myThid) + CALL MNC_CW_RS_R('R',mncFn,0,0,'dxC',DXC, myThid) + CALL MNC_CW_RS_R('R',mncFn,0,0,'dyC',DYC, myThid) + CALL MNC_CW_RS_R('R',mncFn,0,0,'dxF',DXF, myThid) + CALL MNC_CW_RS_R('R',mncFn,0,0,'dyF',DYF, myThid) + CALL MNC_CW_RS_R('R',mncFn,0,0,'dxG',DXG, myThid) + CALL MNC_CW_RS_R('R',mncFn,0,0,'dyG',DYG, myThid) + CALL MNC_CW_RS_R('R',mncFn,0,0,'dxV',DXV, myThid) + CALL MNC_CW_RS_R('R',mncFn,0,0,'dyU',DYU, myThid) + CALL MNC_CW_RS_R('R',mncFn,0,0,'rA', RA, myThid) + CALL MNC_CW_RS_R('R',mncFn,0,0,'rAz',RAZ, myThid) + CALL MNC_CW_RS_R('R',mncFn,0,0,'rAw',RAW, myThid) + CALL MNC_CW_RS_R('R',mncFn,0,0,'rAs',RAS, myThid) + + _END_MASTER(myThid) + + CALL EXCH_XY_RS(XC,myThid) + CALL EXCH_XY_RS(YC,myThid) +#ifdef HRCUBE + CALL EXCH_XY_RS(DXF,myThid) + CALL EXCH_XY_RS(DYF,myThid) +#endif + CALL EXCH_XY_RS(RA,myThid ) + CALL EXCH_Z_XY_RS(XG,myThid) + CALL EXCH_Z_XY_RS(YG,myThid) + CALL EXCH_Z_XY_RS(RAZ,myThid) + CALL EXCH_UV_XY_RS(DXC,DYC,.FALSE.,myThid) + CALL EXCH_UV_XY_RS(RAW,RAS,.FALSE.,myThid) + CALL EXCH_UV_XY_RS(DYG,DXG,.FALSE.,myThid) + CALL EXCH_UV_AGRID_XY_RS(angleSinC,angleCosC,.TRUE.,myThid) + + ELSE +#endif + C Here we make no assumptions about grid symmetry and simply C read the raw grid data from files @@ -242,6 +312,7 @@ cs- end block ENDIF CALL EXCH_UV_XY_RS(DYG,DXG,.FALSE.,myThid) + CALL EXCH_UV_AGRID_XY_RS(angleSinC,angleCosC,.TRUE.,myThid) c write(10) XC c write(10) YC @@ -268,75 +339,92 @@ DO bj = 1,nSy DO bi = 1,nSx iG=bi+(myXGlobalLo-1)/sNx - WRITE(fName(1:15),'("tile",I3.3,".mitgrid")') iG - WRITE(msgBuf,'(A,I4)') 'tile:',iG + WRITE(tmpBuf,'(A,I4)') 'tile:',iG #ifdef ALLOW_EXCH2 - myTile = W2_myTileList(bi) - write(fName(1:15),'("tile",I3.3,".mitgrid")') - & exch2_myface(myTile) - WRITE(msgBuf,'(A,I4)') 'tile:',myTile + myTile = W2_myTileList(bi) + WRITE(tmpBuf,'(A,I4)') 'tile:',myTile + iG = exch2_myface(myTile) #endif - iL = ILNBLNK(msgBuf) - WRITE(msgBuf,'(3A)') msgBuf(1:iL), - & ' ; Read from file ',fName(1:15) + iLen = ILNBLNK(horizGridFile) + IF ( iLen .EQ. 0 ) THEN + WRITE(fName,'("tile",I3.3,".mitgrid")') iG + ELSE + WRITE(fName,'(2A,I3.3,A)') horizGridFile(1:iLen), + & '.face',iG,'.bin' + ENDIF + iLen = ILNBLNK(fName) + iL = ILNBLNK(tmpBuf) + WRITE(msgBuf,'(3A)') tmpBuf(1:iL), + & ' ; Read from file ',fName(1:iLen) CALL PRINT_MESSAGE( msgBuf, standardMessageUnit, & SQUEEZE_RIGHT , myThid) WRITE(msgBuf,'(A)') ' =>' CALL READSYMTILE_RS(fName,1,XC,bi,bj,buf,myThid) iL = ILNBLNK(msgBuf) - WRITE(msgBuf,'(A,1X,A)') msgBuf(1:iL),'XC' + WRITE(tmpBuf,'(A,1X,A)') msgBuf(1:iL),'XC' CALL READSYMTILE_RS(fName,2,YC,bi,bj,buf,myThid) - iL = ILNBLNK(msgBuf) - WRITE(msgBuf,'(A,1X,A)') msgBuf(1:iL),'YC' + iL = ILNBLNK(tmpBuf) + WRITE(msgBuf,'(A,1X,A)') tmpBuf(1:iL),'YC' CALL READSYMTILE_RS(fName,3,DXF,bi,bj,buf,myThid) iL = ILNBLNK(msgBuf) - WRITE(msgBuf,'(A,1X,A)') msgBuf(1:iL),'DXF' + WRITE(tmpBuf,'(A,1X,A)') msgBuf(1:iL),'DXF' CALL READSYMTILE_RS(fName,4,DYF,bi,bj,buf,myThid) - iL = ILNBLNK(msgBuf) - WRITE(msgBuf,'(A,1X,A)') msgBuf(1:iL),'DYF' + iL = ILNBLNK(tmpBuf) + WRITE(msgBuf,'(A,1X,A)') tmpBuf(1:iL),'DYF' CALL READSYMTILE_RS(fName,5,RA,bi,bj,buf,myThid) iL = ILNBLNK(msgBuf) - WRITE(msgBuf,'(A,1X,A)') msgBuf(1:iL),'RA' + WRITE(tmpBuf,'(A,1X,A)') msgBuf(1:iL),'RA' CALL READSYMTILE_RS(fName,6,XG,bi,bj,buf,myThid) - iL = ILNBLNK(msgBuf) - WRITE(msgBuf,'(A,1X,A)') msgBuf(1:iL),'XG' + iL = ILNBLNK(tmpBuf) + WRITE(msgBuf,'(A,1X,A)') tmpBuf(1:iL),'XG' CALL READSYMTILE_RS(fName,7,YG,bi,bj,buf,myThid) iL = ILNBLNK(msgBuf) - WRITE(msgBuf,'(A,1X,A)') msgBuf(1:iL),'YG' + WRITE(tmpBuf,'(A,1X,A)') msgBuf(1:iL),'YG' CALL READSYMTILE_RS(fName,8,DXV,bi,bj,buf,myThid) - iL = ILNBLNK(msgBuf) - WRITE(msgBuf,'(A,1X,A)') msgBuf(1:iL),'DXV' + iL = ILNBLNK(tmpBuf) + WRITE(msgBuf,'(A,1X,A)') tmpBuf(1:iL),'DXV' CALL READSYMTILE_RS(fName,9,DYU,bi,bj,buf,myThid) iL = ILNBLNK(msgBuf) - WRITE(msgBuf,'(A,1X,A)') msgBuf(1:iL),'DYU' + WRITE(tmpBuf,'(A,1X,A)') msgBuf(1:iL),'DYU' CALL READSYMTILE_RS(fName,10,RAZ,bi,bj,buf,myThid) - iL = ILNBLNK(msgBuf) - WRITE(msgBuf,'(A,1X,A)') msgBuf(1:iL),'RAZ' + iL = ILNBLNK(tmpBuf) + WRITE(msgBuf,'(A,1X,A)') tmpBuf(1:iL),'RAZ' CALL READSYMTILE_RS(fName,11,DXC,bi,bj,buf,myThid) iL = ILNBLNK(msgBuf) - WRITE(msgBuf,'(A,1X,A)') msgBuf(1:iL),'DXC' + WRITE(tmpBuf,'(A,1X,A)') msgBuf(1:iL),'DXC' CALL READSYMTILE_RS(fName,12,DYC,bi,bj,buf,myThid) - iL = ILNBLNK(msgBuf) - WRITE(msgBuf,'(A,1X,A)') msgBuf(1:iL),'DYC' + iL = ILNBLNK(tmpBuf) + WRITE(msgBuf,'(A,1X,A)') tmpBuf(1:iL),'DYC' CALL READSYMTILE_RS(fName,13,RAW,bi,bj,buf,myThid) iL = ILNBLNK(msgBuf) - WRITE(msgBuf,'(A,1X,A)') msgBuf(1:iL),'RAW' + WRITE(tmpBuf,'(A,1X,A)') msgBuf(1:iL),'RAW' CALL READSYMTILE_RS(fName,14,RAS,bi,bj,buf,myThid) - iL = ILNBLNK(msgBuf) - WRITE(msgBuf,'(A,1X,A)') msgBuf(1:iL),'RAS' + iL = ILNBLNK(tmpBuf) + WRITE(msgBuf,'(A,1X,A)') tmpBuf(1:iL),'RAS' CALL READSYMTILE_RS(fName,15,DXG,bi,bj,buf,myThid) iL = ILNBLNK(msgBuf) - WRITE(msgBuf,'(A,1X,A)') msgBuf(1:iL),'DXG' + WRITE(tmpBuf,'(A,1X,A)') msgBuf(1:iL),'DXG' CALL READSYMTILE_RS(fName,16,DYG,bi,bj,buf,myThid) - iL = ILNBLNK(msgBuf) - WRITE(msgBuf,'(A,1X,A)') msgBuf(1:iL),'DYG' + iL = ILNBLNK(tmpBuf) + WRITE(msgBuf,'(A,1X,A)') tmpBuf(1:iL),'DYG' + + iLen = ILNBLNK(horizGridFile) + IF ( iLen.GT.0 ) THEN + CALL READSYMTILE_RS(fName,17,angleCosC,bi,bj,buf,myThid) + iL = ILNBLNK(msgBuf) + WRITE(tmpBuf,'(A,1X,A)') msgBuf(1:iL),'AngleCS' + CALL READSYMTILE_RS(fName,18,angleSinC,bi,bj,buf,myThid) + iL = ILNBLNK(tmpBuf) + WRITE(msgBuf,'(A,1X,A)') tmpBuf(1:iL),'AngleSN' + ENDIF CALL PRINT_MESSAGE( msgBuf, standardMessageUnit, & SQUEEZE_RIGHT , myThid) ENDDO ENDDO + _END_MASTER(myThid) CALL EXCH_XY_RS(XC,myThid) @@ -347,29 +435,46 @@ CALL EXCH_XY_RS(DYF,myThid) #endif CALL EXCH_XY_RS(RA,myThid ) -#ifndef ALLOW_EXCH2 CALL EXCH_Z_XY_RS(XG,myThid) CALL EXCH_Z_XY_RS(YG,myThid) C !!! _EXCH_ZUV_XY_R4(DXV, DYU, unSigned, myThid) c CALL EXCH_Z_XY_RS(DXV,myThid) c CALL EXCH_Z_XY_RS(DYU,myThid) CALL EXCH_Z_XY_RS(RAZ,myThid) -#endif /* ALLOW_EXCH2 */ CALL EXCH_UV_XY_RS(DXC,DYC,.FALSE.,myThid) CALL EXCH_UV_XY_RS(RAW,RAS,.FALSE.,myThid) CALL EXCH_UV_XY_RS(DYG,DXG,.FALSE.,myThid) + CALL EXCH_UV_AGRID_XY_RS(angleSinC,angleCosC,.TRUE.,myThid) #endif /* OLD_GRID_IO */ -c CALL WRITE_FULLARRAY_RL('DXV',DXV,1,0,myThid) -c CALL WRITE_FULLARRAY_RL('DYU',DYU,1,0,myThid) -c CALL WRITE_FULLARRAY_RL('RAZ',RAZ,1,0,myThid) -c CALL WRITE_FULLARRAY_RL('XG',XG,1,0,myThid) -c CALL WRITE_FULLARRAY_RL('YG',YG,1,0,myThid) +#ifdef ALLOW_MNC + ENDIF +#endif /* ALLOW_MNC */ + +c CALL WRITE_FULLARRAY_RL('DXV',DXV,1,0,0,0,myThid) +c CALL WRITE_FULLARRAY_RL('DYU',DYU,1,0,0,0,myThid) +c CALL WRITE_FULLARRAY_RL('RAZ',RAZ,1,0,0,0,myThid) +c CALL WRITE_FULLARRAY_RL('XG',XG,1,0,0,0,myThid) +c CALL WRITE_FULLARRAY_RL('YG',YG,1,0,0,0,myThid) + +C-- Require that 0 <= longitude < 360 if using exf package +#ifdef ALLOW_EXF + DO bj = 1,nSy + DO bi = 1,nSx + DO J=1-Oly,sNy+Oly + DO I=1-Olx,sNx+Olx + IF (XC(i,j,bi,bj).lt.0.) XC(i,j,bi,bj)=XC(i,j,bi,bj)+360. + IF (XG(i,j,bi,bj).lt.0.) XG(i,j,bi,bj)=XG(i,j,bi,bj)+360. + ENDDO + ENDDO + ENDDO + ENDDO +#endif /* ALLOW_EXF */ C-- Now let's look at all these beasts IF ( debugLevel .GE. debLevB ) THEN - myiter = 1 + myIter = 1 CALL PLOT_FIELD_XYRL( XC , 'Current XC ' , & myIter, myThid ) CALL PLOT_FIELD_XYRL( YC , 'Current YC ' , @@ -404,6 +509,10 @@ & myIter, myThid ) CALL PLOT_FIELD_XYRL( DYG , 'Current DYG ' , & myIter, myThid ) + CALL PLOT_FIELD_XYRL(angleCosC, 'Current AngleCS ' , + & myIter, myThid ) + CALL PLOT_FIELD_XYRL(angleSinC, 'Current AngleSN ' , + & myIter, myThid ) ENDIF RETURN @@ -438,29 +547,34 @@ #endif /* ALLOW_EXCH2 */ C == Local variables == - INTEGER I,J,dUnit + INTEGER I,J,dUnit, iLen INTEGER length_of_rec INTEGER MDS_RECLEN - INTEGER TN, DNX, DNY, TBX, TBY, TNX, TNY, II, iBase +#ifdef ALLOW_EXCH2 + INTEGER TN, dNx, dNy, TBX, TBY, TNX, TNY, II, iBase +#endif + INTEGER ILNBLNK + EXTERNAL ILNBLNK + iLen = ILNBLNK(fName) #ifdef ALLOW_EXCH2 C Figure out offset of tile within face TN = W2_myTileList(bi) - DNX = exch2_mydnx(TN) - DNY = exch2_mydny(TN) + dNx = exch2_mydnx(TN) + dNy = exch2_mydny(TN) TBX = exch2_tbasex(TN) TBY = exch2_tbasey(TN) TNX = exch2_tnx(TN) TNY = exch2_tny(TN) - CALL MDSFINDUNIT( dUnit, mythid ) + CALL MDSFINDUNIT( dUnit, myThid ) length_of_rec=MDS_RECLEN( 64, (dNx+1), myThid ) - OPEN( dUnit, file=fName, status='old', - & access='direct', recl=length_of_rec ) + OPEN( dUnit, file=fName(1:iLen), status='old', + & access='direct', recl=length_of_rec ) J=0 iBase=(irec-1)*(dny+1) - DO I=1+TBY,SNY+1+TBY - READ(dUnit,rec=I+iBase)(buf(ii),ii=1,dnx+1) + DO I=1+TBY,sNy+1+TBY + READ(dUnit,rec=I+iBase)(buf(ii),ii=1,dNx+1) #ifdef _BYTESWAPIO #ifdef REAL4_IS_SLOW CALL MDS_BYTESWAPR8((dNx+1), buf) @@ -477,10 +591,10 @@ #else /* ALLOW_EXCH2 */ - CALL MDSFINDUNIT( dUnit, mythid ) + CALL MDSFINDUNIT( dUnit, myThid ) length_of_rec=MDS_RECLEN( 64, (sNx+1)*(sNy+1), myThid ) - OPEN( dUnit, file=fName, status='old', - & access='direct', recl=length_of_rec ) + OPEN( dUnit, file=fName(1:iLen), status='old', + & access='direct', recl=length_of_rec ) READ(dUnit,rec=irec) buf CLOSE( dUnit )