--- MITgcm/model/src/ini_cg2d.F 1998/05/25 16:58:49 1.6 +++ MITgcm/model/src/ini_cg2d.F 2001/06/07 15:53:21 1.35 @@ -1,8 +1,8 @@ -C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/model/src/ini_cg2d.F,v 1.6 1998/05/25 16:58:49 cnh Exp $ +C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/model/src/ini_cg2d.F,v 1.35 2001/06/07 15:53:21 adcroft Exp $ +C $Name: $ -#include "CPP_EEOPTIONS.h" +#include "CPP_OPTIONS.h" -CStartOfInterface SUBROUTINE INI_CG2D( myThid ) C /==========================================================\ C | SUBROUTINE INI_CG2D | @@ -11,18 +11,23 @@ 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" #include "EEPARAMS.h" #include "PARAMS.h" #include "GRID.h" +c#include "DYNVARS.h" +#include "SURFACE.h" #include "CG2D.h" +#ifdef ALLOW_OBCS +#include "OBCS.h" +#endif C === Routine arguments === C myThid - Thread no. that called this routine. INTEGER myThid -CEndOfInterface C === Local variables === C xG, yG - Global coordinate location. @@ -31,14 +36,52 @@ C bi,bj - Loop counters C faceArea - Temporary used to hold cell face areas. C I,J,K -C myNorm - Work variable used in clculating normalisation factor +C myNorm - Work variable used in calculating normalisation factor +C sumArea - Work variable used to compute the total Domain Area CHARACTER*(MAX_LEN_MBUF) msgBuf INTEGER bi, bj INTEGER I, J, K - real faceArea - _RL myNorm + _RL faceArea + _RS myNorm + _RL sumArea _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 @@ -51,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) @@ -71,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) @@ -100,16 +163,45 @@ 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) +c _EXCH_XY_R4(aW2d, myThid) +c _EXCH_XY_R4(aS2d, myThid) + CALL EXCH_UV_XY_RS(aW2d,aS2d,.FALSE.,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-- Define the solver tolerance in the appropriate Unit : + cg2dNormaliseRHS = cg2dTargetResWunit.LE.0 + IF (cg2dNormaliseRHS) THEN +C- when using a normalisation of RHS, tolerance has no unit => no conversion + cg2dTolerance = cg2dTargetResidual + ELSE +C- compute the total Area of the domain : + sumArea = 0. + DO bj=myByLo(myThid),myByHi(myThid) + DO bi=myBxLo(myThid),myBxHi(myThid) + DO j=1,sNy + DO i=1,sNx + IF (Ro_surf(i,j,bi,bj).GT.R_low(i,j,bi,bj)) THEN + sumArea = sumArea + rA(i,j,bi,bj) + ENDIF + ENDDO + ENDDO + ENDDO + ENDDO +c WRITE(*,*) ' mythid, sumArea = ', mythid, sumArea + _GLOBAL_SUM_R8( sumArea, myThid ) +C- convert Target-Residual (in W unit) to cg2d-solver residual unit [m^2/s^2] + cg2dTolerance = cg2dNorm * cg2dTargetResWunit + & * sumArea / deltaTMom + WRITE(*,'(2A,1P2E22.14)') ' ini_cg2d: ', + & 'sumArea,cg2dTolerance =', sumArea,cg2dTolerance + ENDIF + C-- Initialise preconditioner C Note. 20th May 1998 C I made a weird discovery! In the model paper we argue @@ -131,24 +223,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 @@ -173,25 +265,14 @@ ENDDO C-- Update overlap regions _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 - ENDDO - ENDDO - ENDDO - ENDDO -C-- Update overlap regions - _EXCH_XY_R8(cg2d_x, myThid) - _EXCH_XY_R8(cg2d_b, myThid) - ENDIF +c _EXCH_XY_R4(pW, myThid) +c _EXCH_XY_R4(pS, myThid) + CALL EXCH_UV_XY_RS(pW,pS,.FALSE.,myThid) +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