C $Header: /u/gcmpack/MITgcm/model/src/ini_pressure.F,v 1.16 2016/02/15 18:00:26 jmc Exp $ C $Name: $ #include "PACKAGES_CONFIG.h" #include "CPP_OPTIONS.h" CBOP C !ROUTINE: INI_PRESSURE C !INTERFACE: SUBROUTINE INI_PRESSURE( myThid ) C !DESCRIPTION: \bv C *==========================================================* C | SUBROUTINE INI_PRESSURE C | o initialise the pressure field consistently with C | temperature and salinity C | - needs to be called after ini_theta, ini_salt, and C | ini_psurf C | - does not include surface pressure loading, because C | that is only available until after C | CALL packages_init_variables C *==========================================================* C \ev C !USES: IMPLICIT NONE C == Global variables == #include "SIZE.h" #include "EEPARAMS.h" #include "PARAMS.h" #include "EOS.h" #include "GRID.h" #include "DYNVARS.h" C !INPUT/OUTPUT PARAMETERS: C == Routine arguments == C myThid - Number of this instance of INI_PRESSURE INTEGER myThid C !LOCAL VARIABLES: C == Local variables == C dPhiHydX,Y :: Gradient (X & Y directions) of Hyd. Potential C bi,bj :: tile indices C i,j,k :: Loop counters INTEGER bi, bj INTEGER i, j, k INTEGER iMin, iMax, jMin, jMax, npiter _RL PhiHydF (1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL PhiHydC (1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL dPhiHydX(1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL dPhiHydY(1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL oldPhi (1-OLx:sNx+OLx,1-OLy:sNy+OLy) _RL count, rmspp, rmsppold _RL sumTile, rmsTile CHARACTER*(MAX_LEN_MBUF) msgBuf CEOP iMin = 1-OLx iMax = sNx+OLx jMin = 1-OLy jMax = sNy+OLy _BEGIN_MASTER( myThid ) WRITE(msgBuf,'(a)') & 'Start initial hydrostatic pressure computation' CALL PRINT_MESSAGE( msgBuf, standardMessageUnit, & SQUEEZE_RIGHT, myThid ) _END_MASTER( myThid ) DO bj = myByLo(myThid), myByHi(myThid) DO bi = myBxLo(myThid), myBxHi(myThid) DO k = 1,Nr DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx totPhiHyd(i,j,k,bi,bj) = 0. _d 0 ENDDO ENDDO ENDDO ENDDO ENDDO IF ( storePhiHyd4Phys ) THEN #ifndef ALLOW_AUTODIFF cph-- deal with this iterative loop for AD once it will cph-- really be needed; cph-- would need storing of totPhiHyd for each npiter rmspp = 1. _d 0 rmsppold = 0. _d 0 npiter = 0 C iterate pressure p = integral of (g*rho(p)*dz) DO npiter= 1, 15 IF ( rmspp.GT.zeroRL ) THEN rmsppold = rmspp rmspp = 0. _d 0 count = 0. _d 0 DO bj = myByLo(myThid), myByHi(myThid) DO bi = myBxLo(myThid), myBxHi(myThid) rmsTile = 0. _d 0 sumTile = 0. _d 0 DO j=1-OLy,sNy+OLy DO i=1-OLx,sNx+OLx phiHydF(i,j) = 0. _d 0 ENDDO ENDDO DO k = 1, Nr C for each level save old pressure and compute new pressure DO j=jMin,jMax DO i=iMin,iMax oldPhi(i,j) = totPhiHyd(i,j,k,bi,bj) ENDDO ENDDO CALL CALC_PHI_HYD( I bi, bj, iMin, iMax, jMin, jMax, k, I theta, salt, U phiHydF, O phiHydC, dPhiHydX, dPhiHydY, I startTime, -1, myThid ) C compute convergence criterion DO j=1,sNy DO i=1,sNx rmsTile = rmsTile & + (totPhiHyd(i,j,k,bi,bj)-oldPhi(i,j))**2 & * maskC(i,j,k,bi,bj) sumTile = sumTile + maskC(i,j,k,bi,bj) ENDDO ENDDO C -- end k loop ENDDO rmspp = rmspp + rmsTile count = count + sumTile ENDDO ENDDO Cml WRITE(msgBuf,'(I10.10)') npiter Cml CALL WRITE_FLD_XYZ_RL( 'POLD.',msgBuf,pold,npiter,myThid) Cml CALL WRITE_FLD_XYZ_RL( 'PINI.',msgBuf,pressure,npiter,myThid) _GLOBAL_SUM_RL( rmspp, myThid ) _GLOBAL_SUM_RL( count, myThid ) IF ( count .EQ. 0. ) THEN rmspp = 0. _d 0 ELSE rmspp = SQRT(rmspp/count) ENDIF _BEGIN_MASTER( myThid ) WRITE(msgBuf,'(A,I4,A,1PE20.12)') & 'Iteration', npiter, ', RMS-difference =', rmspp CALL PRINT_MESSAGE( msgBuf, standardMessageUnit, & SQUEEZE_RIGHT, myThid ) _END_MASTER( myThid ) ENDIF C -- end npiter loop ENDDO C print some diagnostics _BEGIN_MASTER( myThid ) IF ( rmspp .NE. 0. ) THEN IF ( rmspp .EQ. rmsppold ) THEN WRITE(msgBuf,'(A)') & 'Initial hydrostatic pressure did not converge perfectly,' CALL PRINT_MESSAGE( msgBuf, standardMessageUnit, & SQUEEZE_RIGHT, myThid ) WRITE(msgBuf,'(A)') & 'but the RMS-difference is constant, indicating that the' CALL PRINT_MESSAGE( msgBuf, standardMessageUnit, & SQUEEZE_RIGHT, myThid ) WRITE(msgBuf,'(A)') & 'algorithm converged within machine precision.' CALL PRINT_MESSAGE( msgBuf, standardMessageUnit, & SQUEEZE_RIGHT, myThid ) ELSE WRITE(msgBuf,'(A,I2,A)') & 'Initial hydrostatic pressure did not converge after ', & npiter-1, ' steps' CALL PRINT_ERROR( msgBuf, myThid ) STOP 'ABNORMAL END: S/R INI_PRESSURE' ENDIF ENDIF WRITE(msgBuf,'(A)') & 'Initial hydrostatic pressure converged.' CALL PRINT_MESSAGE( msgBuf, standardMessageUnit, & SQUEEZE_RIGHT, myThid ) _END_MASTER( myThid ) #endif /* ALLOW_AUTODIFF */ c-- else of if storePhiHyd4Phys ELSE C print a message and DO nothing _BEGIN_MASTER( myThid ) WRITE(msgBuf,'(A,A)') & 'Pressure is predetermined for buoyancyRelation ', & buoyancyRelation(1:11) CALL PRINT_MESSAGE( msgBuf, standardMessageUnit, & SQUEEZE_RIGHT, myThid ) _END_MASTER( myThid ) ENDIF _BEGIN_MASTER( myThid ) WRITE(msgBuf,'(A)') ' ' CALL PRINT_MESSAGE( msgBuf, standardMessageUnit, & SQUEEZE_RIGHT, myThid ) _END_MASTER( myThid ) RETURN END