Back to home page

MITgcm
 
 

    

File indexing completed on 2022-01-06 06:12:04 UTC

view on githubraw file Latest commit 9f5240b5 on 2022-01-05 15:24:45 UTC
b6b8988e60 Jean*0001 #include "PACKAGES_CONFIG.h"
60c223928f Mart*0002 #include "CPP_OPTIONS.h"
                0003 
                0004 CBOP
                0005 C     !ROUTINE: INI_PRESSURE
                0006 C     !INTERFACE:
7e1bd93d4f Jean*0007       SUBROUTINE INI_PRESSURE( myThid )
60c223928f Mart*0008 C     !DESCRIPTION: \bv
                0009 C     *==========================================================*
                0010 C     | SUBROUTINE INI_PRESSURE
d0a4e8c871 Jean*0011 C     | o initialise the pressure field consistently with
60c223928f Mart*0012 C     |   temperature and salinity
d0a4e8c871 Jean*0013 C     |   - needs to be called after ini_theta, ini_salt, and
60c223928f Mart*0014 C     |     ini_psurf
                0015 C     |   - does not include surface pressure loading, because
d0a4e8c871 Jean*0016 C     |     that is only available until after
7e1bd93d4f Jean*0017 C     |     CALL packages_init_variables
60c223928f Mart*0018 C     *==========================================================*
                0019 C     \ev
                0020 
                0021 C     !USES:
7e1bd93d4f Jean*0022       IMPLICIT NONE
60c223928f Mart*0023 C     == Global variables ==
                0024 #include "SIZE.h"
                0025 #include "EEPARAMS.h"
                0026 #include "PARAMS.h"
                0027 #include "EOS.h"
                0028 #include "GRID.h"
                0029 #include "DYNVARS.h"
                0030 
                0031 C     !INPUT/OUTPUT PARAMETERS:
9f5240b52a Jean*0032 C     myThid     :: my Thread Id number
60c223928f Mart*0033       INTEGER myThid
                0034 
                0035 C     !LOCAL VARIABLES:
7e1bd93d4f Jean*0036 C     dPhiHydX,Y :: Gradient (X & Y directions) of Hyd. Potential
d0a4e8c871 Jean*0037 C     bi,bj      :: tile indices
                0038 C     i,j,k      :: Loop counters
60c223928f Mart*0039       INTEGER bi, bj
d0a4e8c871 Jean*0040       INTEGER  i,  j, k
9f5240b52a Jean*0041 #ifndef ALLOW_AUTODIFF
bbff648e65 Jean*0042       INTEGER  iMin, iMax, jMin, jMax, npiter
9bc836469a Jean*0043       _RL PhiHydF (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0044       _RL PhiHydC (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0045       _RL dPhiHydX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0046       _RL dPhiHydY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0047       _RL oldPhi  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
12a36520f9 Jean*0048       _RL count, rmspp, rmsppold
d0a4e8c871 Jean*0049       _RL sumTile, rmsTile
9f5240b52a Jean*0050 #endif
60c223928f Mart*0051 
                0052       CHARACTER*(MAX_LEN_MBUF) msgBuf
e305438401 Mart*0053 CEOP
12a36520f9 Jean*0054 
10022ac442 Jean*0055       _BEGIN_MASTER( myThid )
7e1bd93d4f Jean*0056       WRITE(msgBuf,'(a)')
60c223928f Mart*0057      &     'Start initial hydrostatic pressure computation'
12a36520f9 Jean*0058       CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
9bc836469a Jean*0059      &                    SQUEEZE_RIGHT, myThid )
10022ac442 Jean*0060       _END_MASTER( myThid )
                0061 
7e1bd93d4f Jean*0062       DO bj = myByLo(myThid), myByHi(myThid)
                0063        DO bi = myBxLo(myThid), myBxHi(myThid)
                0064         DO k = 1,Nr
                0065          DO j=1-OLy,sNy+OLy
                0066           DO i=1-OLx,sNx+OLx
bbff648e65 Jean*0067            totPhiHyd(i,j,k,bi,bj) = 0. _d 0
7e1bd93d4f Jean*0068           ENDDO
                0069          ENDDO
                0070         ENDDO
                0071        ENDDO
                0072       ENDDO
60c223928f Mart*0073 
901f12b7bc Jean*0074       IF ( storePhiHyd4Phys ) THEN
60c223928f Mart*0075 
b6b8988e60 Jean*0076 #ifndef ALLOW_AUTODIFF
d74f26ce7e Patr*0077 cph-- deal with this iterative loop for AD once it will
                0078 cph-- really be needed;
                0079 cph-- would need storing of totPhiHyd for each npiter
                0080 
9f5240b52a Jean*0081        iMin = 1-OLx
                0082        iMax = sNx+OLx
                0083        jMin = 1-OLy
                0084        jMax = sNy+OLy
                0085 
60c223928f Mart*0086        rmspp    = 1. _d 0
                0087        rmsppold = 0. _d 0
                0088        npiter = 0
                0089 
                0090 C     iterate pressure p = integral of (g*rho(p)*dz)
bbff648e65 Jean*0091        DO npiter= 1, 15
9bc836469a Jean*0092         IF ( rmspp.GT.zeroRL ) THEN
                0093          rmsppold = rmspp
                0094          rmspp    = 0. _d 0
                0095          count    = 0. _d 0
                0096          DO bj = myByLo(myThid), myByHi(myThid)
                0097           DO bi = myBxLo(myThid), myBxHi(myThid)
                0098            rmsTile = 0. _d 0
                0099            sumTile = 0. _d 0
                0100            DO j=1-OLy,sNy+OLy
                0101             DO i=1-OLx,sNx+OLx
                0102               phiHydF(i,j) = 0. _d 0
                0103             ENDDO
bbff648e65 Jean*0104            ENDDO
9bc836469a Jean*0105            DO k = 1, Nr
d0a4e8c871 Jean*0106 C     for each level save old pressure and compute new pressure
9bc836469a Jean*0107             DO j=jMin,jMax
                0108              DO i=iMin,iMax
                0109               oldPhi(i,j) = totPhiHyd(i,j,k,bi,bj)
                0110              ENDDO
7e1bd93d4f Jean*0111             ENDDO
9bc836469a Jean*0112             CALL CALC_PHI_HYD(
                0113      I           bi, bj, iMin, iMax, jMin, jMax, k,
                0114      U           phiHydF,
                0115      O           phiHydC, dPhiHydX, dPhiHydY,
                0116      I           startTime, -1, myThid )
60c223928f Mart*0117 C     compute convergence criterion
9bc836469a Jean*0118             DO j=1,sNy
                0119              DO i=1,sNx
                0120               rmsTile = rmsTile
bbff648e65 Jean*0121      &            + (totPhiHyd(i,j,k,bi,bj)-oldPhi(i,j))**2
d0a4e8c871 Jean*0122      &             * maskC(i,j,k,bi,bj)
9bc836469a Jean*0123               sumTile = sumTile + maskC(i,j,k,bi,bj)
                0124              ENDDO
7e1bd93d4f Jean*0125             ENDDO
bbff648e65 Jean*0126 C --      end k loop
9bc836469a Jean*0127            ENDDO
                0128            rmspp = rmspp + rmsTile
                0129            count = count + sumTile
7e1bd93d4f Jean*0130           ENDDO
                0131          ENDDO
60c223928f Mart*0132 Cml        WRITE(msgBuf,'(I10.10)') npiter
                0133 Cml        CALL WRITE_FLD_XYZ_RL( 'POLD.',msgBuf,pold,npiter,myThid)
                0134 Cml        CALL WRITE_FLD_XYZ_RL( 'PINI.',msgBuf,pressure,npiter,myThid)
9bc836469a Jean*0135          _GLOBAL_SUM_RL( rmspp, myThid )
                0136          _GLOBAL_SUM_RL( count, myThid )
                0137          IF ( count .EQ. 0. ) THEN
60c223928f Mart*0138            rmspp = 0. _d 0
9bc836469a Jean*0139          ELSE
d0a4e8c871 Jean*0140            rmspp = SQRT(rmspp/count)
9bc836469a Jean*0141          ENDIF
                0142          _BEGIN_MASTER( myThid )
                0143          WRITE(msgBuf,'(A,I4,A,1PE20.12)')
                0144      &        'Iteration', npiter, ', RMS-difference =', rmspp
                0145          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0146      &                       SQUEEZE_RIGHT, myThid )
                0147          _END_MASTER( myThid )
60c223928f Mart*0148 
9bc836469a Jean*0149         ENDIF
bbff648e65 Jean*0150 C --   end npiter loop
7e1bd93d4f Jean*0151        ENDDO
9bc836469a Jean*0152 
30bf074e57 Jean*0153 C     print some diagnostics
10022ac442 Jean*0154        _BEGIN_MASTER( myThid )
d0a4e8c871 Jean*0155        IF ( rmspp .NE. 0. ) THEN
7e1bd93d4f Jean*0156         IF ( rmspp .EQ. rmsppold ) THEN
bbff648e65 Jean*0157          WRITE(msgBuf,'(A)')
12a36520f9 Jean*0158      &      'Initial hydrostatic pressure did not converge perfectly,'
                0159          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
9bc836469a Jean*0160      &                       SQUEEZE_RIGHT, myThid )
bbff648e65 Jean*0161          WRITE(msgBuf,'(A)')
60c223928f Mart*0162      &      'but the RMS-difference is constant, indicating that the'
12a36520f9 Jean*0163          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
9bc836469a Jean*0164      &                       SQUEEZE_RIGHT, myThid )
bbff648e65 Jean*0165          WRITE(msgBuf,'(A)')
12a36520f9 Jean*0166      &      'algorithm converged within machine precision.'
                0167          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
9bc836469a Jean*0168      &                       SQUEEZE_RIGHT, myThid )
7e1bd93d4f Jean*0169         ELSE
bbff648e65 Jean*0170          WRITE(msgBuf,'(A,I2,A)')
12a36520f9 Jean*0171      &       'Initial hydrostatic pressure did not converge after ',
60c223928f Mart*0172      &          npiter-1, ' steps'
                0173          CALL PRINT_ERROR( msgBuf, myThid )
                0174          STOP 'ABNORMAL END: S/R INI_PRESSURE'
7e1bd93d4f Jean*0175         ENDIF
                0176        ENDIF
bbff648e65 Jean*0177        WRITE(msgBuf,'(A)')
12a36520f9 Jean*0178      &     'Initial hydrostatic pressure converged.'
                0179        CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
9bc836469a Jean*0180      &                     SQUEEZE_RIGHT, myThid )
10022ac442 Jean*0181        _END_MASTER( myThid )
60c223928f Mart*0182 
b6b8988e60 Jean*0183 #endif /* ALLOW_AUTODIFF */
d74f26ce7e Patr*0184 
901f12b7bc Jean*0185 c-- else of if storePhiHyd4Phys
7e1bd93d4f Jean*0186       ELSE
                0187 C     print a message and DO nothing
10022ac442 Jean*0188        _BEGIN_MASTER( myThid )
bbff648e65 Jean*0189        WRITE(msgBuf,'(A,A)')
60c223928f Mart*0190      &        'Pressure is predetermined for buoyancyRelation ',
                0191      &        buoyancyRelation(1:11)
12a36520f9 Jean*0192        CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
9bc836469a Jean*0193      &                     SQUEEZE_RIGHT, myThid )
10022ac442 Jean*0194        _END_MASTER( myThid )
60c223928f Mart*0195 
7e1bd93d4f Jean*0196       ENDIF
60c223928f Mart*0197 
10022ac442 Jean*0198       _BEGIN_MASTER( myThid )
bbff648e65 Jean*0199       WRITE(msgBuf,'(A)') ' '
12a36520f9 Jean*0200       CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
9bc836469a Jean*0201      &                    SQUEEZE_RIGHT, myThid )
10022ac442 Jean*0202       _END_MASTER( myThid )
60c223928f Mart*0203 
12a36520f9 Jean*0204       RETURN
7e1bd93d4f Jean*0205       END
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated from https://github.com/MITgcm/MITgcm by the 2.2.1-MITgcm-0.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!