model/src/diags_sound_speed.F

fb2c52cfd4 aver*0001 #include "PACKAGES_CONFIG.h"
                0002 #include "CPP_OPTIONS.h"
                0003
                0004 CBOP
                0005 C     !ROUTINE: DIAGS_SOUND_SPEED
                0006 C     !INTERFACE:
                0007       SUBROUTINE DIAGS_SOUND_SPEED(
                0008      I                       myThid )
                0009 C     !DESCRIPTION: \bv
                0010 C     *==========================================================*
                0011 C     | S/R DIAGS_SOUND_SPEED
                0012 C     | o Diagnose speed of sound in seawater
9590d44ddf aver*0013 C     |   from the algorithm by Del Grosso (1974).
fb2c52cfd4 aver*0014 C     |   This is NOT the sound-speed that can be derived from
                0015 C     |   the equation of state (EOS). It is independent of
                0016 C     |   the model setup specific EOS.
                0017 C     |
                0018 C     | o Reference:
                0019 C     | V. A. Del Grosso, "New equation for the speed of sound in
                0020 C     | natural waters (with comparisons to other equations),"
                0021 C     | J. Acoust. Soc. Am. 56, 1084-1091 (1974).
                0022 C     *==========================================================*
                0023 C     \ev
                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 "SURFACE.h"
                0032 #include "DYNVARS.h"
                0033 #include "EOS.h"
                0034
                0035 C     !INPUT/OUTPUT PARAMETERS:
                0036 C     == Routine Arguments ==
                0037 C     myThid :: Thread number for this instance of the routine.
                0038       INTEGER myThid
                0039
                0040 #ifdef INCLUDE_SOUNDSPEED_CALC_CODE
                0041 C     !FUNCTIONS:
                0042       _RL SW_TEMP
                0043       EXTERNAL SW_TEMP
                0044 #ifdef ALLOW_DIAGNOSTICS
                0045       LOGICAL  DIAGNOSTICS_IS_ON
                0046       EXTERNAL DIAGNOSTICS_IS_ON
                0047 #endif /* ALLOW_DIAGNOSTICS */
                0048
                0049 C     !LOCAL VARIABLES:
                0050 C     == Local variables ==
                0051       INTEGER bi,bj
                0052       INTEGER i,j,k
                0053       _RL c0,ct,cs,cp,cstp
                0054       _RL pres,sal,temp
                0055       LOGICAL calc_soundSpeed
                0056 CEOP
                0057
                0058       calc_soundSpeed = .FALSE.
                0059
                0060 C--   switch on this flag if Sound-Speed needed in any cost fct
                0061 c#ifdef ALLOW_CSOUND_COST
                0062 c     calc_soundSpeed = .TRUE.
                0063 c#endif
                0064
                0065 #ifdef ALLOW_DIAGNOSTICS
                0066       IF ( useDiagnostics ) THEN
                0067         calc_soundSpeed = calc_soundSpeed
                0068      &               .OR. DIAGNOSTICS_IS_ON( 'CSound  ', myThid )
                0069       ENDIF
                0070 #endif /* ALLOW_DIAGNOSTICS */
                0071
                0072       IF ( calc_soundSpeed ) THEN
                0073        c0   = 1402.392 _d 0
                0074        ct   = 0. _d 0
                0075        cs   = 0. _d 0
                0076        cp   = 0. _d 0
                0077        cstp = 0. _d 0
                0078        pres = 0. _d 0
                0079        sal  = 0. _d 0
                0080        temp = 0. _d 0
                0081        DO bj=myByLo(myThid),myByHi(myThid)
                0082         DO bi=myBxLo(myThid),myBxHi(myThid)
                0083          DO k=1,Nr
                0084           DO j=1,sNy
                0085            DO i=1,sNx
                0086             cSound(i,j,k,bi,bj) = 0.0 _d 0
                0087            ENDDO
                0088           ENDDO
                0089          ENDDO
                0090          DO k=1,Nr
                0091           DO j=1,sNy
                0092            DO i=1,sNx
                0093             IF ( maskC(i,j,k,bi,bj).EQ.oneRS ) THEN
                0094 C pressure in dbar (for SW_TEMP)
                0095              pres = rhoConst*( totPhiHyd(i,j,k,bi,bj)
                0096      &               - rC(k)*gravity )*SItodBar
                0097              temp = SW_TEMP( SALT(i,j,k,bi,bj),
                0098      &                 THETA(i,j,k,bi,bj), pres, zeroRL )
                0099              sal  = SALT(i,j,k,bi,bj)
9590d44ddf aver*0100 C convert pressure to kg/cm^2 for Del Grosso algorithm
fb2c52cfd4 aver*0101              pres = pres/gravity
                0102              ct   = ( 5.01109398873 _d 0 - ( 0.550946843172 _d -1
                0103      &              - 0.221535969240 _d -3 * temp ) * temp ) * temp
                0104              cs   = ( 1.32952290781 _d 0 + 0.128955756844 _d -3 * sal
                0105      &              ) * sal
                0106              cp   = ( 0.156059257041 _d 0 + ( 0.244998688441 _d -4
                0107      &              - 0.883392332513 _d -8 * pres ) * pres ) * pres
                0108              cstp = ( - 0.127562783426 _d -1
                0109      &                + 0.968403156410 _d -4 *temp ) * temp * sal
                0110      &            + ((  0.635191613389 _d -2
                0111      &                +(0.265484716608 _d -7 *temp
                0112      &                - 0.159349479045 _d -5
                0113      &                + 0.522116437235 _d -9 *pres ) * pres
                0114      &                - 0.438031096213 _d -6 * temp * temp
                0115      &                +(0.485639620015 _d -5 * sal
                0116      &                - 0.340597039004 _d -3 ) *sal ) * temp
                0117      &                - 0.161674495909 _d -8 * sal * sal * pres
                0118      &              ) * pres
                0119              cSound(i,j,k,bi,bj) = c0+ct+cs+cp+cstp
                0120             ENDIF
                0121            ENDDO
                0122           ENDDO
                0123          ENDDO
                0124         ENDDO
                0125        ENDDO
                0126 #ifdef ALLOW_DIAGNOSTICS
                0127        IF ( useDiagnostics ) THEN
                0128         CALL DIAGNOSTICS_FILL(cSound,'CSound  ',0,Nr,0,1,1,myThid)
                0129        ENDIF
                0130 #endif /* ALLOW_DIAGNOSTICS */
                0131
                0132 C-    end-if calc_soundSpeed
                0133       ENDIF
                0134
                0135 #endif /* INCLUDE_SOUNDSPEED_CALC_CODE */
                0136
                0137       RETURN
                0138       END