model/src/load_ref_files.F

b9306711eb Jean*0001 c #include "PACKAGES_CONFIG.h"
                0002 #include "CPP_OPTIONS.h"
                0003
                0004 CBOP
d87b427b2d Jean*0005 C     !ROUTINE: LOAD_REF_FILES
b9306711eb Jean*0006 C     !INTERFACE:
d87b427b2d Jean*0007       SUBROUTINE LOAD_REF_FILES( myThid )
b9306711eb Jean*0008 C     !DESCRIPTION: \bv
                0009 C     *==========================================================*
d87b427b2d Jean*0010 C     | SUBROUTINE LOAD_REF_FILES
                0011 C     | o Read reference vertical profile from files
                0012 C     |   (Pot.Temp., Salinity/Specif.Humid., density ... )
b9306711eb Jean*0013 C     *==========================================================*
                0014 C     \ev
                0015
                0016 C     !USES:
                0017       IMPLICIT NONE
                0018 C     === Global variables ===
                0019 #include "SIZE.h"
                0020 #include "EEPARAMS.h"
                0021 #include "PARAMS.h"
6ef71429db Jean*0022 #include "GRID.h"
b9306711eb Jean*0023
                0024 C     !INPUT/OUTPUT PARAMETERS:
                0025 C     == Routine arguments ==
                0026 C     myThid     :: my Thread Id number
                0027       INTEGER myThid
                0028
95393692ef Jean*0029 C     !FUNCTIONS:
                0030       INTEGER  ILNBLNK
                0031       EXTERNAL ILNBLNK
                0032
b9306711eb Jean*0033 C     !LOCAL VARIABLES:
                0034 C     == Local variables ==
                0035 C     k          :: loop index
                0036 C     msgBuf     :: Informational/error message buffer
                0037       _RL    tracerDefault
95393692ef Jean*0038       INTEGER  k, kLen
b9306711eb Jean*0039       CHARACTER*(MAX_LEN_MBUF) msgBuf
                0040 CEOP
                0041
                0042       _BEGIN_MASTER( myThid )
                0043
                0044 C--   Set reference Potential Temperature
                0045       IF ( tRefFile .EQ. ' ' ) THEN
                0046 C-    set default vertical profile for temperature: tRef
                0047         tracerDefault = 20.
                0048         IF ( fluidIsAir ) tracerDefault = 300.
aad87ebb4c Jean*0049         IF ( thetaConst.NE.UNSET_RL ) tracerDefault = thetaConst
b9306711eb Jean*0050         DO k=1,Nr
                0051           IF (tRef(k).EQ.UNSET_RL) tRef(k) = tracerDefault
                0052           tracerDefault = tRef(k)
                0053         ENDDO
                0054       ELSE
                0055 C-    check for multiple definitions:
                0056         DO k=1,Nr
                0057          IF (tRef(k).NE.UNSET_RL) THEN
6ef71429db Jean*0058           WRITE(msgBuf,'(2A,I4,A)') 'S/R LOAD_REF_FILES: ',
                0059      &      'Cannot set both tRef(k=', k, ') and tRefFile'
b9306711eb Jean*0060           CALL PRINT_ERROR( msgBuf, myThid )
                0061           STOP 'ABNORMAL END: S/R INI_PARMS'
                0062          ENDIF
                0063         ENDDO
                0064       ENDIF
                0065 C-    read from file:
                0066       IF ( tRefFile .NE. ' ' ) THEN
                0067         kLen = ILNBLNK(tRefFile)
95393692ef Jean*0068         CALL READ_GLVEC_RL( tRefFile, ' ', tRef, Nr, 1, myThid )
6ef71429db Jean*0069         WRITE(msgBuf,'(3A)') 'S/R LOAD_REF_FILES: ',
                0070      &    'tRef loaded from file: ', tRefFile(1:kLen)
b9306711eb Jean*0071         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
6ef71429db Jean*0072      &                      SQUEEZE_RIGHT, myThid )
b9306711eb Jean*0073       ENDIF
aad87ebb4c Jean*0074       IF ( thetaConst.EQ.UNSET_RL ) thetaConst = tRef(1)
b9306711eb Jean*0075
                0076 C--   Set reference Salinity/Specific Humidity
                0077       IF ( sRefFile .EQ. ' ' ) THEN
                0078 C-    set default vertical profile for salinity/water-vapour: sRef
                0079         tracerDefault = 30.
                0080         IF ( fluidIsAir ) tracerDefault = 0.
                0081         DO k=1,Nr
                0082           IF (sRef(k).EQ.UNSET_RL) sRef(k) = tracerDefault
                0083           tracerDefault = sRef(k)
                0084         ENDDO
                0085       ELSE
                0086 C-    check for multiple definitions:
                0087         DO k=1,Nr
                0088          IF (sRef(k).NE.UNSET_RL) THEN
6ef71429db Jean*0089           WRITE(msgBuf,'(2A,I4,A)') 'S/R LOAD_REF_FILES: ',
                0090      &      'Cannot set both sRef(k=', k, ') and sRefFile'
b9306711eb Jean*0091           CALL PRINT_ERROR( msgBuf, myThid )
                0092           STOP 'ABNORMAL END: S/R INI_PARMS'
                0093          ENDIF
                0094         ENDDO
                0095       ENDIF
                0096 C-    read from file:
                0097       IF ( sRefFile .NE. ' ' ) THEN
                0098         kLen = ILNBLNK(sRefFile)
95393692ef Jean*0099         CALL READ_GLVEC_RL( sRefFile, ' ', sRef, Nr, 1, myThid )
6ef71429db Jean*0100         WRITE(msgBuf,'(3A)') 'S/R LOAD_REF_FILES: ',
                0101      &    'sRef loaded from file: ', sRefFile(1:kLen)
b9306711eb Jean*0102         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
6ef71429db Jean*0103      &                      SQUEEZE_RIGHT, myThid )
b9306711eb Jean*0104       ENDIF
                0105
                0106 C--   Set reference Density
                0107       IF ( rhoRefFile .NE. ' ' ) THEN
                0108         kLen = ILNBLNK(rhoRefFile)
95393692ef Jean*0109 C-    read from file:
                0110         CALL READ_GLVEC_RL( rhoRefFile, ' ', rho1Ref, Nr, 1, myThid )
6ef71429db Jean*0111         WRITE(msgBuf,'(3A)') 'S/R LOAD_REF_FILES: ',
                0112      &    'rho1Ref loaded from file: ', rhoRefFile(1:kLen)
b9306711eb Jean*0113         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
6ef71429db Jean*0114      &                      SQUEEZE_RIGHT, myThid )
                0115       ENDIF
                0116
                0117 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0118
                0119 C--   Set gravity vertical profile
                0120       IF ( gravityFile .NE. ' ' ) THEN
                0121         kLen = ILNBLNK(gravityFile)
                0122 C-    read from file and store, for now, in gravFacC:
                0123         CALL READ_GLVEC_RL( gravityFile, ' ', gravFacC, Nr, 1, myThid )
                0124         WRITE(msgBuf,'(3A)') 'S/R LOAD_REF_FILES: ',
                0125      &    'gravity vert. prof. loaded from file: ', gravityFile(1:kLen)
                0126         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0127      &                      SQUEEZE_RIGHT, myThid )
                0128
                0129 C-    Set gravity vertical profile factor: assume surface-interface
                0130 C     density-factor to be 1 (grav. acceleration @ rF(1) = gravity)
                0131 C- Note: done here (instead of in set_ref_state.F) since gravity fac
                0132 C       might be needed to initialise EOS coeffs (in ini_eos.F)
                0133         gravFacF(1)   = 1. _d 0
                0134 C       gravFac(k) = gravity ratio between layer k and top interface
                0135         DO k=1,Nr
                0136           gravFacC(k) = gravFacC(k)*recip_gravity
                0137         ENDDO
                0138         DO k=2,Nr
                0139 C       since rkSign=-1, recip_drC(k) = 1./(rC(k-1)-rC(k))
                0140           gravFacF(k) = ( gravFacC(k-1)*(rF(k)-rC(k))
                0141      &                  + gravFacC(k)*(rC(k-1)-rF(k)) )*recip_drC(k)
                0142         ENDDO
                0143 C       extrapolate down to the bottom:
                0144         gravFacF(Nr+1) = ( gravFacC(Nr)*(rF(Nr+1)-rF(Nr))
                0145      &                   + gravFacF(Nr)*(rC(Nr)-rF(Nr+1))
                0146      &                   ) / (rC(Nr)-rF(Nr))
                0147 C-      set reciprocal gravity-factor:
                0148         DO k=1,Nr
                0149           recip_gravFacC(k) = 1. _d 0/gravFacC(k)
                0150         ENDDO
                0151         DO k=1,Nr+1
                0152           recip_gravFacF(k) = 1. _d 0/gravFacF(k)
                0153         ENDDO
                0154
                0155       ELSE
                0156 C-    Initialise gravity vertical profile factor:
                0157         DO k=1,Nr
                0158           gravFacC(k) = 1. _d 0
                0159           recip_gravFacC(k) = 1. _d 0
                0160         ENDDO
                0161         DO k=1,Nr+1
                0162           gravFacF(k) = 1. _d 0
                0163           recip_gravFacF(k) = 1. _d 0
                0164         ENDDO
b9306711eb Jean*0165       ENDIF
                0166
                0167       _END_MASTER(myThid)
                0168 C--   Everyone else must wait for the parameters to be loaded
                0169       _BARRIER
                0170
                0171       RETURN
                0172       END