C $Id: ini_cartesian_grid.F,v 1.1 1998/04/22 19:15:30 cnh Exp $ #include "CPP_EEOPTIONS.h" CStartOfInterface SUBROUTINE INI_CARTESIAN_GRID( myThid ) C /==========================================================\ C | SUBROUTINE INI_CARTESIAN_GRID | C | o Initialise model coordinate system | C |==========================================================| C | These arrays are used throughout the code in evaluating | C | gradients, integrals and spatial avarages. This routine | C | is called separately by each thread and initialise only | C | the region of the domain it is "responsible" for. | C | Notes: | C | Two examples are included. One illustrates the | C | initialisation of a cartesian grid. The other shows the | C | inialisation of a spherical polar grid. Other orthonormal| C | grids can be fitted into this design. In this case | C | custom metric terms also need adding to account for the | C | projections of velocity vectors onto these grids. | C | The structure used here also makes it possible to | C | implement less regular grid mappings. In particular | C | o Schemes which leave out blocks of the domain that are | C | all land could be supported. | C | o Multi-level schemes such as icosohedral or cubic | C | grid projections onto a sphere can also be fitted | C | within the strategy we use. | C | Both of the above also require modifying the support | C | routines that map computational blocks to simulation | C | domain blocks. | C | Under the cartesian grid mode primitive distances in X | C | and Y are in metres. Disktance in Z are in m or Pa | C | depending on the vertical gridding mode. | C \==========================================================/ C === Global variables === #include "SIZE.h" #include "EEPARAMS.h" #include "PARAMS.h" #include "GRID.h" C == Routine arguments == C myThid - Number of this instance of INI_CARTESIAN_GRID INTEGER myThid CEndOfInterface C == Local variables == C xG, yG - Global coordinate location. C zG C xBase - South-west corner location for process. C yBase C zUpper - Work arrays for upper and lower C zLower cell-face heights. C phi - Temporary scalar C xBase - Temporaries for lower corner coordinate C yBase C iG, jG - Global coordinate index. Usually used to hold C the south-west global coordinate of a tile. C bi,bj - Loop counters C zUpper - Temporary arrays holding z coordinates of C zLower upper and lower faces. C I,J,K _RL xG, yG, zG _RL phi _RL zUpper(Nz), zLower(Nz) _RL xBase, yBase INTEGER iG, jG INTEGER bi, bj INTEGER I, J, K C-- Simple example of inialisation on cartesian grid C-- First set coordinates of cell centers C This operation is only performed at start up so for more C complex configurations it is usually OK to pass iG, jG to a custom C function and have it return xG and yG. C Set up my local grid first DO bj = myByLo(myThid), myByHi(myThid) jG = myYGlobalLo + (bj-1)*sNy DO bi = myBxLo(myThid), myBxHi(myThid) iG = myXGlobalLo + (bi-1)*sNx yBase = 0. _d 0 xBase = 0. _d 0 DO i=1,iG-1 xBase = xBase + delX(i) ENDDO DO j=1,jG-1 yBase = yBase + delY(j) ENDDO yG = yBase DO J=1,sNy xG = xBase DO I=1,sNx xc(I,J,bi,bj) = xG + delX(iG+i-1)*0.5 _d 0 yc(I,J,bi,bj) = yG + delY(jG+j-1)*0.5 _d 0 xG = xG + delX(iG+I-1) dxF(I,J,bi,bj) = delX(iG+i-1) dyF(I,J,bi,bj) = delY(jG+j-1) ENDDO yG = yG + delY(jG+J-1) ENDDO ENDDO ENDDO C Now sync. and get edge regions from other threads and/or processes. C Note: We could just set the overlap regions ourselves here but C exchanging edges is safer and is good practice! _EXCH_XY_R4( xc, myThid ) _EXCH_XY_R4( yc, myThid ) _EXCH_XY_R4(dxF, myThid ) _EXCH_XY_R4(dyF, myThid ) C-- Calculate separation between other points C dxG, dyG are separations between cell corners along cell faces. DO bj = myByLo(myThid), myByHi(myThid) DO bi = myBxLo(myThid), myBxHi(myThid) DO J=1,sNy DO I=1,sNx dxG(I,J,bi,bj) = (dxF(I,J,bi,bj)+dxF(I,J-1,bi,bj))*0.5 _d 0 dyG(I,J,bi,bj) = (dyF(I,J,bi,bj)+dyF(I-1,J,bi,bj))*0.5 _d 0 ENDDO ENDDO ENDDO ENDDO _EXCH_XY_R4(dxG, myThid ) _EXCH_XY_R4(dyG, myThid ) C dxV, dyU are separations between velocity points along cell faces. DO bj = myByLo(myThid), myByHi(myThid) DO bi = myBxLo(myThid), myBxHi(myThid) DO J=1,sNy DO I=1,sNx dxV(I,J,bi,bj) = (dxG(I,J,bi,bj)+dxG(I-1,J,bi,bj))*0.5 _d 0 dyU(I,J,bi,bj) = (dyG(I,J,bi,bj)+dyG(I,J-1,bi,bj))*0.5 _d 0 ENDDO ENDDO ENDDO ENDDO _EXCH_XY_R4(dxV, myThid ) _EXCH_XY_R4(dyU, myThid ) C dxC, dyC is separation between cell centers DO bj = myByLo(myThid), myByHi(myThid) DO bi = myBxLo(myThid), myBxHi(myThid) DO J=1,sNy DO I=1,sNx dxC(I,J,bi,bj) = (dxF(I,J,bi,bj)+dxF(I-1,J,bi,bj))*0.5 D0 dyC(I,J,bi,bj) = (dyF(I,J,bi,bj)+dyF(I,J-1,bi,bj))*0.5 D0 ENDDO ENDDO ENDDO ENDDO _EXCH_XY_R4(dxC, myThid ) _EXCH_XY_R4(dyC, myThid ) C Calculate recipricols DO bj = myByLo(myThid), myByHi(myThid) DO bi = myBxLo(myThid), myBxHi(myThid) DO J=1,sNy DO I=1,sNx rDxG(I,J,bi,bj)=1.d0/dxG(I,J,bi,bj) rDyG(I,J,bi,bj)=1.d0/dyG(I,J,bi,bj) rDxC(I,J,bi,bj)=1.d0/dxC(I,J,bi,bj) rDyC(I,J,bi,bj)=1.d0/dyC(I,J,bi,bj) rDxF(I,J,bi,bj)=1.d0/dxF(I,J,bi,bj) rDyF(I,J,bi,bj)=1.d0/dyF(I,J,bi,bj) rDxV(I,J,bi,bj)=1.d0/dxV(I,J,bi,bj) rDyU(I,J,bi,bj)=1.d0/dyU(I,J,bi,bj) ENDDO ENDDO ENDDO ENDDO _EXCH_XY_R4(rDxG, myThid ) _EXCH_XY_R4(rDyG, myThid ) _EXCH_XY_R4(rDxC, myThid ) _EXCH_XY_R4(rDyC, myThid ) _EXCH_XY_R4(rDxF, myThid ) _EXCH_XY_R4(rDyF, myThid ) _EXCH_XY_R4(rDxV, myThid ) _EXCH_XY_R4(rDyU, myThid ) C Calculate vertical face area DO bj = myByLo(myThid), myByHi(myThid) DO bi = myBxLo(myThid), myBxHi(myThid) DO J=1,sNy DO I=1,sNx zA(I,J,bi,bj) = dxF(I,J,bi,bj)*dyF(I,J,bi,bj) ENDDO ENDDO ENDDO ENDDO DO bj = myByLo(myThid), myByHi(myThid) DO bi = myBxLo(myThid), myBxHi(myThid) DO K=1,Nz DO J=1,sNy DO I=1,sNx IF (HFacC(I,J,K,bi,bj) .NE. 0. D0 ) THEN rHFacC(I,J,K,bi,bj) = 1. D0 / HFacC(I,J,K,bi,bj) ELSE rHFacC(I,J,K,bi,bj) = 0. D0 ENDIF IF (HFacW(I,J,K,bi,bj) .NE. 0. D0 ) THEN rHFacW(I,J,K,bi,bj) = 1. D0 / HFacW(I,J,K,bi,bj) maskW(I,J,K,bi,bj) = 1. D0 ELSE rHFacW(I,J,K,bi,bj) = 0. D0 maskW(I,J,K,bi,bj) = 0.0 D0 ENDIF IF (HFacS(I,J,K,bi,bj) .NE. 0. D0 ) THEN rHFacS(I,J,K,bi,bj) = 1. D0 / HFacS(I,J,K,bi,bj) maskS(I,J,K,bi,bj) = 1. D0 ELSE rHFacS(I,J,K,bi,bj) = 0. D0 maskS(I,J,K,bi,bj) = 0. D0 ENDIF ENDDO ENDDO ENDDO ENDDO ENDDO C Now sync. and get/send edge regions that are shared with C other threads. _EXCH_XYZ_R4(rHFacC , myThid ) _EXCH_XYZ_R4(rHFacW , myThid ) _EXCH_XYZ_R4(rHFacS , myThid ) _EXCH_XYZ_R4(maskW , myThid ) _EXCH_XYZ_R4(maskS , myThid ) C RETURN END C $Id: ini_cartesian_grid.F,v 1.1 1998/04/22 19:15:30 cnh Exp $