--- MITgcm/model/src/ini_cg2d.F 1998/06/08 21:43:01 1.12 +++ MITgcm/model/src/ini_cg2d.F 2001/03/09 20:42:56 1.33 @@ -1,6 +1,7 @@ -C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/model/src/ini_cg2d.F,v 1.12 1998/06/08 21:43:01 cnh Exp $ +C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/model/src/ini_cg2d.F,v 1.33 2001/03/09 20:42:56 jmc Exp $ +C $Name: $ -#include "CPP_EEOPTIONS.h" +#include "CPP_OPTIONS.h" CStartOfInterface SUBROUTINE INI_CG2D( myThid ) @@ -11,6 +12,7 @@ C | These arrays are purely a function of the basin geom. | C | We set then here once and them use then repeatedly. | C \==========================================================/ + IMPLICIT NONE C === Global variables === #include "SIZE.h" @@ -18,7 +20,11 @@ #include "PARAMS.h" #include "GRID.h" #include "DYNVARS.h" -#include "CG2D.h" +#include "SURFACE.h" +#include "CG2D_INTERNAL.h" +#ifdef ALLOW_OBCS +#include "OBCS.h" +#endif C === Routine arguments === C myThid - Thread no. that called this routine. @@ -37,9 +43,45 @@ INTEGER bi, bj INTEGER I, J, K _RL faceArea - _RL myNorm + _RS myNorm _RL aC, aCw, aCs +C-- Initialise -Boyancy at surface level : Bo_surf +C Bo_surf is defined as d/dr(Phi_surf) and set to g/rtoz (linear free surface) +C with rtoz = conversion factor from r-unit to z-unit (=horiVertRatio) +C an acurate formulation will include P_surf and T,S_surf effects on rho_surf: +C z-ocean (rtoz=1) : Bo_surf = - Boyancy = gravity * rho_surf/rho_0 +C p-atmos (rtoz=rho_c*g) : Bo_surf = (1/rho)_surf +C-- + IF ( buoyancyRelation .EQ. 'ATMOSPHERIC' ) THEN + DO bj=myByLo(myThid),myByHi(myThid) + DO bi=myBxLo(myThid),myBxHi(myThid) + DO J=1-Oly,sNy+Oly + DO I=1-Olx,sNx+Olx + Bo_surf(I,J,bi,bj) = recip_rhoConst + recip_Bo(I,J,bi,bj) = rhoConst + ENDDO + ENDDO + ENDDO + ENDDO + ELSE +C- gBaro = gravity (except for External mode test with reduced gravity) + DO bj=myByLo(myThid),myByHi(myThid) + DO bi=myBxLo(myThid),myBxHi(myThid) + DO J=1-Oly,sNy+Oly + DO I=1-Olx,sNx+Olx + Bo_surf(I,J,bi,bj) = gBaro + recip_Bo(I,J,bi,bj) = 1. _d 0 / gBaro + ENDDO + ENDDO + ENDDO + ENDDO + ENDIF + +C-- Update overlap regions + _EXCH_XY_R8(Bo_surf, myThid) + _EXCH_XY_R8(recip_Bo, myThid) + C-- Initialise laplace operator C aW2d: integral in Z Ax/dX C aS2d: integral in Z Ay/dY @@ -52,18 +94,38 @@ aS2d(I,J,bi,bj) = 0. _d 0 ENDDO ENDDO - DO K=1,Nz + DO K=1,Nr DO J=1,sNy DO I=1,sNx - faceArea = _dyG(I,J,bi,bj)*dzF(K)*_hFacW(I,J,K,bi,bj) + faceArea = _dyG(I,J,bi,bj)*drF(K) + & *_hFacW(I,J,K,bi,bj) aW2d(I,J,bi,bj) = aW2d(I,J,bi,bj) + - & gBaro*faceArea*_rdxC(I,J,bi,bj) - faceArea = _dxG(I,J,bi,bj)*dzF(K)*_hFacS(I,J,K,bi,bj) + & implicSurfPress*implicDiv2DFlow + & *faceArea*recip_dxC(I,J,bi,bj) + faceArea = _dxG(I,J,bi,bj)*drF(K) + & *_hFacS(I,J,K,bi,bj) aS2d(I,J,bi,bj) = aS2d(I,J,bi,bj) + - & gBaro*faceArea*_rdyC(I,J,bi,bj) + & implicSurfPress*implicDiv2DFlow + & *faceArea*recip_dyC(I,J,bi,bj) ENDDO ENDDO ENDDO +#ifdef ALLOW_OBCS + IF (useOBCS) THEN + DO I=1,sNx + IF (OB_Jn(I,bi,bj).NE.0) aS2d(I,OB_Jn(I,bi,bj),bi,bj)=0. + IF (OB_Jn(I,bi,bj).NE.0) aS2d(I,OB_Jn(I,bi,bj)+1,bi,bj)=0. + IF (OB_Js(I,bi,bj).NE.0) aS2d(I,OB_Js(I,bi,bj)+1,bi,bj)=0. + IF (OB_Js(I,bi,bj).NE.0) aS2d(I,OB_Js(I,bi,bj),bi,bj)=0. + ENDDO + DO J=1,sNy + IF (OB_Ie(J,bi,bj).NE.0) aW2d(OB_Ie(J,bi,bj),J,bi,bj)=0. + IF (OB_Ie(J,bi,bj).NE.0) aW2d(OB_Ie(J,bi,bj)+1,J,bi,bj)=0. + IF (OB_Iw(J,bi,bj).NE.0) aW2d(OB_Iw(J,bi,bj)+1,J,bi,bj)=0. + IF (OB_Iw(J,bi,bj).NE.0) aW2d(OB_Iw(J,bi,bj),J,bi,bj)=0. + ENDDO + ENDIF +#endif DO J=1,sNy DO I=1,sNx myNorm = MAX(ABS(aW2d(I,J,bi,bj)),myNorm) @@ -72,17 +134,17 @@ ENDDO ENDDO ENDDO - cg2dNbuf(1,myThid) = myNorm - _GLOBAL_MAX_R8( cg2dNbuf, myNorm, myThid ) - IF ( cg2dNbuf(1,1) .NE. 0. _d 0 ) THEN - myNorm = 1. _d 0/cg2dNbuf(1,1) + _GLOBAL_MAX_R4( myNorm, myThid ) + IF ( myNorm .NE. 0. _d 0 ) THEN + myNorm = 1. _d 0/myNorm ELSE myNorm = 1. _d 0 ENDIF cg2dNorm = myNorm _BEGIN_MASTER( myThid ) CcnhDebugStarts - WRITE(msgBuf,'(A,F)') '// CG2D normalisation factor = ', cg2dNorm + WRITE(msgBuf,'(A,E40.25)') '// CG2D normalisation factor = ', + & cg2dNorm CALL PRINT_MESSAGE( msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1) WRITE(msgBuf,*) ' ' CALL PRINT_MESSAGE( msgBuf,standardMessageUnit,SQUEEZE_RIGHT,1) @@ -101,14 +163,14 @@ C-- Update overlap regions CcnhDebugStarts -C CALL PLOT_FIELD_XYR8( aW2d, 'AW2D INI_CG2D.1' , 1, myThid ) -C CALL PLOT_FIELD_XYR8( aS2d, 'AS2D INI_CG2D.1' , 1, myThid ) +C CALL PLOT_FIELD_XYRS( aW2d, 'AW2D INI_CG2D.1' , 1, myThid ) +C CALL PLOT_FIELD_XYRS( aS2d, 'AS2D INI_CG2D.1' , 1, myThid ) CcnhDebugEnds _EXCH_XY_R4(aW2d, myThid) _EXCH_XY_R4(aS2d, myThid) CcnhDebugStarts -C CALL PLOT_FIELD_XYR8( aW2d, 'AW2D INI_CG2D.2' , 1, myThid ) -C CALL PLOT_FIELD_XYR8( aS2d, 'AS2D INI_CG2D.2' , 1, myThid ) +C CALL PLOT_FIELD_XYRS( aW2d, 'AW2D INI_CG2D.2' , 1, myThid ) +C CALL PLOT_FIELD_XYRS( aS2d, 'AS2D INI_CG2D.2' , 1, myThid ) CcnhDebugEnds C-- Initialise preconditioner @@ -132,24 +194,24 @@ pC(I,J,bi,bj) = 1. _d 0 aC = -( & aW2d(I,J,bi,bj) + aW2d(I+1,J ,bi,bj) - & +aS2d(I,J,bi,bj) + aS2D(I ,J+1,bi,bj) - & +freeSurfFac*myNorm* - & zA(I,J,bi,bj)/deltaTMom/deltaTMom + & +aS2d(I,J,bi,bj) + aS2d(I ,J+1,bi,bj) + & +freeSurfFac*myNorm*recip_Bo(I,J,bi,bj)* + & rA(I,J,bi,bj)/deltaTMom/deltaTMom & ) aCs = -( & aW2d(I,J-1,bi,bj) + aW2d(I+1,J-1,bi,bj) & +aS2d(I,J-1,bi,bj) + aS2d(I ,J ,bi,bj) - & +freeSurfFac*myNorm* - & zA(I,J-1,bi,bj)/deltaTMom/deltaTMom + & +freeSurfFac*myNorm*recip_Bo(I,J-1,bi,bj)* + & rA(I,J-1,bi,bj)/deltaTMom/deltaTMom & ) aCw = -( & aW2d(I-1,J,bi,bj) + aW2d(I ,J ,bi,bj) & +aS2d(I-1,J,bi,bj) + aS2d(I-1,J+1,bi,bj) - & +freeSurfFac*myNorm* - & zA(I-1,J,bi,bj)/deltaTMom/deltaTMom + & +freeSurfFac*myNorm*recip_Bo(I-1,J,bi,bj)* + & rA(I-1,J,bi,bj)/deltaTMom/deltaTMom & ) IF ( aC .EQ. 0. ) THEN - pC(I,J,bi,bj) = 0. _d 0 + pC(I,J,bi,bj) = 1. _d 0 ELSE pC(I,J,bi,bj) = 1. _d 0 / aC ENDIF @@ -176,29 +238,11 @@ _EXCH_XY_R4(pC, myThid) _EXCH_XY_R4(pW, myThid) _EXCH_XY_R4(pS, myThid) - -C-- Set default values for initial guess and RHS - IF ( startTime .EQ. 0 ) THEN - DO bj=myByLo(myThid),myByHi(myThid) - DO bi=myBxLo(myThid),myBxHi(myThid) - DO J=1,sNy - DO I=1,sNx - cg2d_x(I,J,bi,bj) = 0. _d 0 - cg2d_b(I,J,bi,bj) = 0. _d 0 -#ifdef ALLOW_CD - cg2d_xNM1(I,J,bi,bj) = 0. _d 0 -#endif - ENDDO - ENDDO - ENDDO - ENDDO -C-- Update overlap regions - _EXCH_XY_R8(cg2d_x, myThid) - _EXCH_XY_R8(cg2d_b, myThid) -#ifdef ALLOW_CD - _EXCH_XY_R8(cg2d_xNM1, myThid) -#endif - ENDIF +CcnhDebugStarts +C CALL PLOT_FIELD_XYRS( pC, 'pC INI_CG2D.2' , 1, myThid ) +C CALL PLOT_FIELD_XYRS( pW, 'pW INI_CG2D.2' , 1, myThid ) +C CALL PLOT_FIELD_XYRS( pS, 'pS INI_CG2D.2' , 1, myThid ) +CcnhDebugEnds RETURN END