fizhi-cs-aqualev20/code/ini_salt.F

7ae8fb32b5 Andr*0001 #include "CPP_OPTIONS.h"
                0002
                0003 CBOP
                0004 C     !ROUTINE: INI_SALT
                0005 C     !INTERFACE:
                0006       SUBROUTINE INI_SALT ( myThid )
                0007 C     !DESCRIPTION: \bv
                0008 C     *==========================================================*
                0009 C     | SUBROUTINE INI_SALT
                0010 C     | o Set model initial salinity field.
                0011 C     *==========================================================*
                0012 C     | There are several options for setting the initial
                0013 C     | temperature file
                0014 C     |  1. Inline code
                0015 C     |  2. Vertical profile ( uniform S in X and Y )
                0016 C     |  3. Three-dimensional data from a file. For example from
                0017 C     |     Levitus or from a checkpoint file from a previous
                0018 C     |     integration.
                0019 C     | In addition to setting the salinity field we also
                0020 C     | set the initial salinity tendency term here.
                0021 C     *==========================================================*
                0022 C     \ev
                0023
                0024 C     !USES:
                0025       IMPLICIT NONE
                0026 C     === Global variables ===
                0027 #include "SIZE.h"
                0028 #include "EEPARAMS.h"
                0029 #include "PARAMS.h"
                0030 #include "GRID.h"
                0031 #include "DYNVARS.h"
                0032
                0033 C     !INPUT/OUTPUT PARAMETERS:
                0034 C     == Routine arguments ==
                0035 C     myThid -  Number of this instance of INI_SALT
                0036       INTEGER myThid
                0037
                0038 C     !LOCAL VARIABLES:
                0039 C     == Local variables ==
                0040 C     bi,bj  - Loop counters
                0041 C     I,J,K
                0042       INTEGER bi, bj
                0043       INTEGER  I,  J, K, localWarnings
                0044       CHARACTER*(MAX_LEN_MBUF) msgBuf
                0045 CEOP
                0046
                0047        _RL pedyn(Nr+1), pdyn(Nr), pkappa(Nr)
                0048        integer Lbotij
                0049        _RL getcon, kappa, dum, pinmb
                0050        _RL temperature(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
                0051        _RL rhum(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
                0052        _RL qstar(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
                0053
                0054 C Need to convert from input RH to q -- Need pressure, temperature
                0055 C Build pressures - we know there is no topography and eta = 0
                0056        pedyn(1) = 100000.
                0057        do K = 2,Nr+1
                0058         pedyn(K) = pedyn(K-1) - drF(K-1)
                0059        enddo
                0060        do K = 1,Nr
                0061         pdyn(K)=(pedyn(K+1)+pedyn(K))/2.
                0062        enddo
                0063        kappa = getcon('KAPPA')
                0064        do K = 1,Nr
                0065         pkappa(K)=(pdyn(K)/100000.)**kappa
                0066        enddo
                0067 C Now convert from theta to Temperature
                0068       DO bj = myByLo(myThid), myByHi(myThid)
                0069        DO bi = myBxLo(myThid), myBxHi(myThid)
                0070         DO K=1,Nr
                0071          DO J=1,sNy
                0072           DO I=1,sNx
                0073            temperature(I,J,K,bi,bj) = theta(I,J,K,bi,bj) * pkappa(K)
                0074           ENDDO
                0075          ENDDO
                0076         ENDDO
                0077        ENDDO
                0078       ENDDO
                0079
                0080 C--   Initialise salinity field to the vertical reference profile
                0081       DO bj = myByLo(myThid), myByHi(myThid)
                0082        DO bi = myBxLo(myThid), myBxHi(myThid)
                0083         DO K=1,Nr
                0084          DO J=1-Oly,sNy+Oly
                0085           DO I=1-Olx,sNx+Olx
                0086            salt(I,J,K,bi,bj) = sRef(K)
                0087           ENDDO
                0088          ENDDO
                0089         ENDDO
                0090        ENDDO
                0091       ENDDO
                0092
                0093       IF ( hydrogSaltFile .NE. ' ' ) THEN
                0094        _BEGIN_MASTER( myThid )
                0095        CALL READ_FLD_XYZ_RL( hydrogSaltFile, ' ', rhum, 0, myThid )
                0096        _END_MASTER(myThid)
                0097
                0098 C Now convert from rh to q
                0099       DO bj = myByLo(myThid), myByHi(myThid)
                0100        DO bi = myBxLo(myThid), myBxHi(myThid)
                0101         DO K=1,Nr-1
                0102          DO J=1,sNy
                0103           DO I=1,sNx
                0104            pinmb = pdyn(K)/100.
                0105            call qsat(temperature(i,j,k,bi,bj),pinmb,qstar(i,j,k,bi,bj),
                0106      .               dum,.false.)
                0107            salt(I,J,K,bi,bj) = rhum(i,j,k,bi,bj) * qstar(i,j,k,bi,bj)
                0108           ENDDO
                0109          ENDDO
                0110         ENDDO
                0111        ENDDO
                0112       ENDDO
                0113
12ffad7671 Jean*0114        _EXCH_XYZ_RL(salt  , myThid )
7ae8fb32b5 Andr*0115
                0116       ENDIF
                0117
                0118 C     Apply mask and test consistancy
                0119       localWarnings=0
                0120       DO bj = myByLo(myThid), myByHi(myThid)
                0121        DO bi = myBxLo(myThid), myBxHi(myThid)
                0122         DO K=1,Nr
                0123          DO J=1,sNy
                0124           DO I=1,sNx
                0125            IF (hFacC(I,J,K,bi,bj).EQ.0) salt(I,J,K,bi,bj) = 0.
                0126            IF (hFacC(I,J,K,bi,bj).NE.0.AND.salt(I,J,K,bi,bj).EQ.0.
                0127      &      .AND. sRef(k).NE.0.) THEN
                0128              localWarnings=localWarnings+1
                0129             ENDIF
                0130           ENDDO
                0131          ENDDO
                0132         ENDDO
                0133        ENDDO
                0134       ENDDO
                0135       IF (localWarnings.NE.0) THEN
                0136        WRITE(msgBuf,'(A,A)')
                0137      &  'S/R INI_SALT: salt = 0 identically. If this is intentional',
                0138      &  'you will need to edit ini_salt.F to avoid this safety check'
                0139        CALL PRINT_ERROR( msgBuf , myThid)
                0140        STOP 'ABNORMAL END: S/R INI_SALT'
                0141       ENDIF
                0142
                0143       CALL PLOT_FIELD_XYZRL( salt, 'Initial Salinity' , Nr, 1, myThid )
                0144
                0145       RETURN
                0146       END