C $Header: /u/gcmpack/MITgcm/model/src/ini_cg3d.F,v 1.29 2016/05/04 22:09:54 jmc Exp $ C $Name: $ #include "PACKAGES_CONFIG.h" #include "CPP_OPTIONS.h" CBOP C !ROUTINE: INI_CG3D C !INTERFACE: SUBROUTINE INI_CG3D( myThid ) C !DESCRIPTION: \bv C *==========================================================* C | SUBROUTINE INI_CG3D C | o Initialise 3d conjugate gradient solver operators. C *==========================================================* C | These arrays are purely a function of the basin geom. C | We set then here once and them use then repeatedly. C *==========================================================* C \ev C !USES: IMPLICIT NONE C === Global variables === #include "SIZE.h" #include "EEPARAMS.h" #include "PARAMS.h" #include "GRID.h" #include "SURFACE.h" #include "CG3D.h" #include "SOLVE_FOR_PRESSURE3D.h" C !INPUT/OUTPUT PARAMETERS: C === Routine arguments === C myThid :: My Thread Id number INTEGER myThid #ifdef ALLOW_NONHYDROSTATIC C !LOCAL VARIABLES: C === Local variables === C bi,bj :: tile indices C i,j,k :: Loop counters C faceArea :: Temporary used to hold cell face areas. C myNorm :: Work variable used in clculating normalisation factor CHARACTER*(MAX_LEN_MBUF) msgBuf INTEGER bi, bj INTEGER i, j, k, ks _RL faceArea _RS myNorm _RL theRecip_Dr _RL aU, aL, aW, aE, aN, aS _RL tmpFac, nh_Fac, igwFac _RL locGamma CEOP CcnhDebugStarts c _RL phi(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) CcnhDebugEnds C-- Initialise to zero over the full range of indices DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) DO k=1,Nr C- From common bloc CG3D_R: DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx aW3d(i,j,k,bi,bj) = 0. aS3d(i,j,k,bi,bj) = 0. aV3d(i,j,k,bi,bj) = 0. aC3d(i,j,k,bi,bj) = 0. zMC (i,j,k,bi,bj) = 0. zML (i,j,k,bi,bj) = 0. zMU (i,j,k,bi,bj) = 0. ENDDO ENDDO C- From common bloc CG3D_WK_R: DO j=0,sNy+1 DO i=0,sNx+1 cg3d_q(i,j,k,bi,bj) = 0. cg3d_r(i,j,k,bi,bj) = 0. cg3d_s(i,j,k,bi,bj) = 0. ENDDO ENDDO C- From common bloc SFP3D_COMMON_R: DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx cg3d_b(i,j,k,bi,bj) = 0. ENDDO ENDDO ENDDO ENDDO ENDDO nh_Fac = 0. igwFac = 0. IF ( nonHydrostatic & .AND. nh_Am2.NE.0. ) nh_Fac = 1. _d 0 / nh_Am2 IF ( implicitIntGravWave ) igwFac = 1. _d 0 IF ( use3Dsolver ) THEN C-- Initialise laplace operator C aW3d: Ax/dX C aS3d: Ay/dY C aV3d: Ar/dR myNorm = 0. _d 0 DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) DO k=1,Nr DO j=1,sNy DO i=1,sNx+1 faceArea = _dyG(i,j,bi,bj)*drF(k) & *_hFacW(i,j,k,bi,bj) #ifdef ALLOW_OBCS & *maskInC(i,j,bi,bj)*maskInC(i-1,j,bi,bj) #endif aW3d(i,j,k,bi,bj) = faceArea*recip_dxC(i,j,bi,bj) & *implicitNHPress*implicDiv2DFlow myNorm = MAX(ABS(aW3d(i,j,k,bi,bj)),myNorm) ENDDO ENDDO C deep-model: *deepFacC (faceArea), /deepFacC (recip_dx,y): => no net effect DO j=1,sNy+1 DO i=1,sNx faceArea = _dxG(i,j,bi,bj)*drF(K) & *_hFacS(i,j,k,bi,bj) #ifdef ALLOW_OBCS & *maskInC(i,j,bi,bj)*maskInC(i,j-1,bi,bj) #endif aS3d(i,j,k,bi,bj) = faceArea*recip_dyC(i,j,bi,bj) & *implicitNHPress*implicDiv2DFlow myNorm = MAX(ABS(aS3d(i,j,k,bi,bj)),myNorm) ENDDO ENDDO ENDDO DO k=1,1 DO j=1,sNy DO i=1,sNx aV3d(i,j,k,bi,bj) = 0. ENDDO ENDDO ENDDO DO k=2,Nr tmpFac = nh_Fac*rVel2wUnit(k)*rVel2wUnit(k) & + igwFac*dBdrRef(k)*deltaTMom*dTtracerLev(k) IF (tmpFac.GT.0. ) tmpFac = 1. _d 0 / tmpFac DO j=1,sNy DO i=1,sNx faceArea = _rA(i,j,bi,bj)*maskC(i,j, k ,bi,bj) & *maskC(i,j,k-1,bi,bj) & *deepFac2F(k) #ifdef ALLOW_OBCS & *maskInC(i,j,bi,bj) #endif theRecip_Dr = recip_drC(k) c theRecip_Dr = caja & drF(k )*_hFacC(i,j,k ,bi,bj)*0.5 caja & +drF(k-1)*_hFacC(i,j,k-1,bi,bj)*0.5 c IF ( theRecip_Dr .NE. 0. ) c & theRecip_Dr = 1. _d 0/theRecip_Dr aV3d(i,j,k,bi,bj) = faceArea*theRecip_Dr*tmpFac & *implicitNHPress*implicDiv2DFlow myNorm = MAX(ABS(aV3d(i,j,k,bi,bj)),myNorm) ENDDO ENDDO ENDDO ENDDO ENDDO _GLOBAL_MAX_RS( myNorm, myThid ) IF ( myNorm .NE. 0. _d 0 ) THEN myNorm = 1. _d 0/myNorm ELSE myNorm = 1. _d 0 ENDIF _BEGIN_MASTER( myThid ) C-- set global parameter in common block: cg3dNorm = myNorm WRITE(msgBuf,'(2A,1PE23.16)') 'INI_CG3D: ', & 'CG3D normalisation factor = ', cg3dNorm CALL PRINT_MESSAGE( msgBuf, standardMessageUnit, & SQUEEZE_RIGHT, myThid ) WRITE(msgBuf,*) ' ' CALL PRINT_MESSAGE( msgBuf, standardMessageUnit, & SQUEEZE_RIGHT, myThid ) _END_MASTER( myThid ) DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) C- Set solver main diagonal term DO k=1,Nr DO j=1,sNy DO i=1,sNx aW = aW3d( i, j, k, bi,bj) aE = aW3d(i+1,j, k, bi,bj) aN = aS3d( i,j+1,k, bi,bj) aS = aS3d( i, j, k, bi,bj) aU = aV3d( i, j, k, bi,bj) IF ( k .NE. Nr ) THEN aL = aV3d(i, j,k+1,bi,bj) ELSE aL = 0. ENDIF aC3d(i,j,k,bi,bj) = -aW-aE-aN-aS-aU-aL ENDDO ENDDO ENDDO C- Add free-surface source term IF ( selectNHfreeSurf.GE.1 ) THEN DO j=1,sNy DO i=1,sNx locGamma = drC(1)*recip_Bo(i,j,bi,bj) & /( deltaTMom*deltaTFreeSurf & *implicitNHPress*implicDiv2DFlow ) ks = 1 c ks = kSurfC(i,j,bi,bj) c IF ( ks.LE.Nr ) THEN aC3d(i,j,ks,bi,bj) = aC3d(i,j,ks,bi,bj) & - freeSurfFac*recip_Bo(i,j,bi,bj) & *rA(i,j,bi,bj)*deepFac2F(ks)/deltaTMom/deltaTFreeSurf & / (1. _d 0 + locGamma ) c ENDIF ENDDO ENDDO ELSE DO j=1,sNy DO i=1,sNx ks = kSurfC(i,j,bi,bj) IF ( ks.LE.Nr ) THEN aC3d(i,j,ks,bi,bj) = aC3d(i,j,ks,bi,bj) & - freeSurfFac*recip_Bo(i,j,bi,bj) & *rA(i,j,bi,bj)*deepFac2F(ks)/deltaTMom/deltaTFreeSurf ENDIF ENDDO ENDDO ENDIF C- Matrix solver normalisation DO k=1,Nr DO j=1,sNy DO i=1,sNx aW3d(i,j,k,bi,bj) = aW3d(i,j,k,bi,bj)*myNorm aS3d(i,j,k,bi,bj) = aS3d(i,j,k,bi,bj)*myNorm aV3d(i,j,k,bi,bj) = aV3d(i,j,k,bi,bj)*myNorm aC3d(i,j,k,bi,bj) = aC3d(i,j,k,bi,bj)*myNorm ENDDO ENDDO ENDDO ENDDO ENDDO C-- Update overlap regions CALL EXCH_UV_XYZ_RS(aW3d,aS3d,.FALSE.,myThid) _EXCH_XYZ_RS(aV3d, myThid) _EXCH_XYZ_RS(aC3d, myThid) CcnhDebugStarts C CALL PLOT_FIELD_XYZRS( aW3d, 'AW3D INI_CG3D.1' , Nr, 1, myThid ) C CALL PLOT_FIELD_XYZRS( aS3d, 'AS3D INI_CG3D.1' , Nr, 1, myThid ) CcnhDebugEnds C-- Initialise preconditioner C For now PC is just the identity. Change to C be LU factorization of d2/dz2 later. Note C check for consistency with S/R CG3D before C assuming zML is lower and zMU is upper! DO bj=myByLo(myThid),myByHi(myThid) DO bi=myBxLo(myThid),myBxHi(myThid) DO k=1,Nr DO j=1,sNy DO i=1,sNx IF ( aC3d(i,j,k,bi,bj) .NE. 0. ) THEN zMC(i,j,k,bi,bj) = aC3d(i,j,k,bi,bj) zML(i,j,k,bi,bj) = aV3d(i,j,k,bi,bj) IF ( k.NE.Nr ) THEN zMU(i,j,k,bi,bj)= aV3d(i,j,k+1,bi,bj) ELSE zMU(i,j,k,bi,bj)= 0. ENDIF CcnhDebugStarts C zMC(i,j,k,bi,bj) = 1. C zMU(i,j,k,bi,bj) = 0. C zML(i,j,k,bi,bj) = 0. CcnhDebugEnds ELSE zMC(i,j,k,bi,bj) = 1. _d 0 zMU(i,j,k,bi,bj) = 0. zML(i,j,k,bi,bj) = 0. ENDIF ENDDO ENDDO ENDDO k = 1 DO j=1,sNy DO i=1,sNx zMC(i,j,k,bi,bj) = 1. _d 0 / zMC(i,j,k,bi,bj) zMU(i,j,k,bi,bj) = zMU(i,j,k,bi,bj)*zMC(i,j,k,bi,bj) ENDDO ENDDO DO k=2,Nr DO j=1,sNy DO i=1,sNx zMC(i,j,k,bi,bj) = 1. _d 0 / & (zMC(i,j,k,bi,bj)-zML(i,j,k,bi,bj)*zMU(i,j,k-1,bi,bj)) zMU(i,j,k,bi,bj) = zMU(i,j,k,bi,bj)*zMC(i,j,k,bi,bj) ENDDO ENDDO ENDDO DO k=1,Nr DO j=1,sNy DO i=1,sNx IF ( aC3d(i,j,k,bi,bj) .EQ. 0. ) THEN zMC(i,j,k,bi,bj) = 1. zML(i,j,k,bi,bj) = 0. zMU(i,j,k,bi,bj) = 0. CcnhDebugStarts C ELSE C zMC(i,j,k,bi,bj) = 1. C zML(i,j,k,bi,bj) = 0. C zMU(i,j,k,bi,bj) = 0. CcnhDEbugEnds ENDIF ENDDO ENDDO ENDDO ENDDO ENDDO C-- Update overlap regions _EXCH_XYZ_RS(zMC, myThid) _EXCH_XYZ_RS(zML, myThid) _EXCH_XYZ_RS(zMU, myThid) IF ( debugLevel .GE. debLevC ) THEN CALL WRITE_FLD_XYZ_RS( 'zMC',' ',zMC, 0, myThid ) CALL WRITE_FLD_XYZ_RS( 'zML',' ',zML, 0, myThid ) CALL WRITE_FLD_XYZ_RS( 'zMU',' ',zMU, 0, myThid ) ENDIF CcnhDebugStarts c DO k=1,Nr c DO j=1-OLy,sNy+OLy c DO i=1-OLx,sNx+OLx c phi(i,j,1,1) = zMc(i,j,k,1,1) c ENDDO c ENDDO C CALL PLOT_FIELD_XYRS( phi, 'zMC INI_CG3D.1' , 1, myThid ) c ENDDO C CALL PLOT_FIELD_XYRS( zMU, 'zMU INI_CG3D.1' , Nr, 1, myThid ) C CALL PLOT_FIELD_XYRS( zML, 'zML INI_CG3D.1' , Nr, 1, myThid ) CcnhDebugEnds C-- end if (use3Dsolver) ENDIF #endif /* ALLOW_NONHYDROSTATIC */ RETURN END