model/src/load_ref_files.F

C $Header: /u/gcmpack/MITgcm/model/src/load_ref_files.F,v 1.3 2014/05/06 15:47:05 jmc Exp $
C $Name:  $

c #include "PACKAGES_CONFIG.h"
#include "CPP_OPTIONS.h"

CBOP
C     !ROUTINE: LOAD_REF_FILES
C     !INTERFACE:
      SUBROUTINE LOAD_REF_FILES( myThid )
C     !DESCRIPTION: \bv
C     *==========================================================*
C     | SUBROUTINE LOAD_REF_FILES
C     | o Read reference vertical profile from files
C     |   (Pot.Temp., Salinity/Specif.Humid., density ... )
C     *==========================================================*
C     \ev

C     !USES:
      IMPLICIT NONE
C     === Global variables ===
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"

C     !INPUT/OUTPUT PARAMETERS:
C     == Routine arguments ==
C     myThid     :: my Thread Id number
      INTEGER myThid

C     !FUNCTIONS:
      INTEGER  ILNBLNK
      EXTERNAL ILNBLNK

C     !LOCAL VARIABLES:
C     == Local variables ==
C     k          :: loop index
C     msgBuf     :: Informational/error message buffer
      _RL    tracerDefault
      INTEGER  k, kLen
      CHARACTER*(MAX_LEN_MBUF) msgBuf
CEOP

      _BEGIN_MASTER( myThid )

C--   Set reference Potential Temperature
      IF ( tRefFile .EQ. ' ' ) THEN
C-    set default vertical profile for temperature: tRef
        tracerDefault = 20.
        IF ( fluidIsAir ) tracerDefault = 300.
        IF ( thetaConst.NE.UNSET_RL ) tracerDefault = thetaConst
        DO k=1,Nr
          IF (tRef(k).EQ.UNSET_RL) tRef(k) = tracerDefault
          tracerDefault = tRef(k)
        ENDDO
      ELSE
C-    check for multiple definitions:
        DO k=1,Nr
         IF (tRef(k).NE.UNSET_RL) THEN
          WRITE(msgBuf,'(2A,I4,A)') 'S/R LOAD_REF_FILES: ',
     &      'Cannot set both tRef(k=', k, ') and tRefFile'
          CALL PRINT_ERROR( msgBuf, myThid )
          STOP 'ABNORMAL END: S/R INI_PARMS'
         ENDIF
        ENDDO
      ENDIF
C-    read from file:
      IF ( tRefFile .NE. ' ' ) THEN
        kLen = ILNBLNK(tRefFile)
        CALL READ_GLVEC_RL( tRefFile, ' ', tRef, Nr, 1, myThid )
        WRITE(msgBuf,'(3A)') 'S/R LOAD_REF_FILES: ',
     &    'tRef loaded from file: ', tRefFile(1:kLen)
        CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
     &                      SQUEEZE_RIGHT, myThid )
      ENDIF
      IF ( thetaConst.EQ.UNSET_RL ) thetaConst = tRef(1)

C--   Set reference Salinity/Specific Humidity
      IF ( sRefFile .EQ. ' ' ) THEN
C-    set default vertical profile for salinity/water-vapour: sRef
        tracerDefault = 30.
        IF ( fluidIsAir ) tracerDefault = 0.
        DO k=1,Nr
          IF (sRef(k).EQ.UNSET_RL) sRef(k) = tracerDefault
          tracerDefault = sRef(k)
        ENDDO
      ELSE
C-    check for multiple definitions:
        DO k=1,Nr
         IF (sRef(k).NE.UNSET_RL) THEN
          WRITE(msgBuf,'(2A,I4,A)') 'S/R LOAD_REF_FILES: ',
     &      'Cannot set both sRef(k=', k, ') and sRefFile'
          CALL PRINT_ERROR( msgBuf, myThid )
          STOP 'ABNORMAL END: S/R INI_PARMS'
         ENDIF
        ENDDO
      ENDIF
C-    read from file:
      IF ( sRefFile .NE. ' ' ) THEN
        kLen = ILNBLNK(sRefFile)
        CALL READ_GLVEC_RL( sRefFile, ' ', sRef, Nr, 1, myThid )
        WRITE(msgBuf,'(3A)') 'S/R LOAD_REF_FILES: ',
     &    'sRef loaded from file: ', sRefFile(1:kLen)
        CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
     &                      SQUEEZE_RIGHT, myThid )
      ENDIF

C--   Set reference Density
      IF ( rhoRefFile .NE. ' ' ) THEN
        kLen = ILNBLNK(rhoRefFile)
C-    read from file:
        CALL READ_GLVEC_RL( rhoRefFile, ' ', rho1Ref, Nr, 1, myThid )
        WRITE(msgBuf,'(3A)') 'S/R LOAD_REF_FILES: ',
     &    'rho1Ref loaded from file: ', rhoRefFile(1:kLen)
        CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
     &                      SQUEEZE_RIGHT, myThid )
      ENDIF

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

C--   Set gravity vertical profile
      IF ( gravityFile .NE. ' ' ) THEN
        kLen = ILNBLNK(gravityFile)
C-    read from file and store, for now, in gravFacC:
        CALL READ_GLVEC_RL( gravityFile, ' ', gravFacC, Nr, 1, myThid )
        WRITE(msgBuf,'(3A)') 'S/R LOAD_REF_FILES: ',
     &    'gravity vert. prof. loaded from file: ', gravityFile(1:kLen)
        CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
     &                      SQUEEZE_RIGHT, myThid )

C-    Set gravity vertical profile factor: assume surface-interface
C     density-factor to be 1 (grav. acceleration @ rF(1) = gravity)
C- Note: done here (instead of in set_ref_state.F) since gravity fac 
C       might be needed to initialise EOS coeffs (in ini_eos.F)
        gravFacF(1)   = 1. _d 0
C       gravFac(k) = gravity ratio between layer k and top interface
        DO k=1,Nr
          gravFacC(k) = gravFacC(k)*recip_gravity
        ENDDO
        DO k=2,Nr
C       since rkSign=-1, recip_drC(k) = 1./(rC(k-1)-rC(k))
          gravFacF(k) = ( gravFacC(k-1)*(rF(k)-rC(k))
     &                  + gravFacC(k)*(rC(k-1)-rF(k)) )*recip_drC(k)
        ENDDO
C       extrapolate down to the bottom:
        gravFacF(Nr+1) = ( gravFacC(Nr)*(rF(Nr+1)-rF(Nr))
     &                   + gravFacF(Nr)*(rC(Nr)-rF(Nr+1))
     &                   ) / (rC(Nr)-rF(Nr))
C-      set reciprocal gravity-factor:
        DO k=1,Nr
          recip_gravFacC(k) = 1. _d 0/gravFacC(k)
        ENDDO
        DO k=1,Nr+1
          recip_gravFacF(k) = 1. _d 0/gravFacF(k)
        ENDDO

      ELSE
C-    Initialise gravity vertical profile factor:
        DO k=1,Nr
          gravFacC(k) = 1. _d 0
          recip_gravFacC(k) = 1. _d 0
        ENDDO
        DO k=1,Nr+1
          gravFacF(k) = 1. _d 0
          recip_gravFacF(k) = 1. _d 0
        ENDDO
      ENDIF

      _END_MASTER(myThid)
C--   Everyone else must wait for the parameters to be loaded
      _BARRIER

      RETURN
      END
1	jmc	1.4	C $Header: /u/gcmpack/MITgcm/model/src/load_ref_files.F,v 1.3 2014/05/06 15:47:05 jmc Exp $
2	jmc	1.1	C $Name: $
3
4			c #include "PACKAGES_CONFIG.h"
5			#include "CPP_OPTIONS.h"
6
7			CBOP
8			C !ROUTINE: LOAD_REF_FILES
9			C !INTERFACE:
10			SUBROUTINE LOAD_REF_FILES( myThid )
11			C !DESCRIPTION: \bv
12			C ==========================================================
13			C \| SUBROUTINE LOAD_REF_FILES
14			C \| o Read reference vertical profile from files
15			C \| (Pot.Temp., Salinity/Specif.Humid., density ... )
16			C ==========================================================
17			C \ev
18
19			C !USES:
20			IMPLICIT NONE
21			C === Global variables ===
22			#include "SIZE.h"
23			#include "EEPARAMS.h"
24			#include "PARAMS.h"
25	jmc	1.4	#include "GRID.h"
26	jmc	1.1
27			C !INPUT/OUTPUT PARAMETERS:
28			C == Routine arguments ==
29			C myThid :: my Thread Id number
30			INTEGER myThid
31
32	jmc	1.2	C !FUNCTIONS:
33			INTEGER ILNBLNK
34			EXTERNAL ILNBLNK
35
36	jmc	1.1	C !LOCAL VARIABLES:
37			C == Local variables ==
38			C k :: loop index
39			C msgBuf :: Informational/error message buffer
40			_RL tracerDefault
41	jmc	1.2	INTEGER k, kLen
42	jmc	1.1	CHARACTER*(MAX_LEN_MBUF) msgBuf
43			CEOP
44
45			_BEGIN_MASTER( myThid )
46
47			C-- Set reference Potential Temperature
48			IF ( tRefFile .EQ. ' ' ) THEN
49			C- set default vertical profile for temperature: tRef
50			tracerDefault = 20.
51			IF ( fluidIsAir ) tracerDefault = 300.
52	jmc	1.3	IF ( thetaConst.NE.UNSET_RL ) tracerDefault = thetaConst
53	jmc	1.1	DO k=1,Nr
54			IF (tRef(k).EQ.UNSET_RL) tRef(k) = tracerDefault
55			tracerDefault = tRef(k)
56			ENDDO
57			ELSE
58			C- check for multiple definitions:
59			DO k=1,Nr
60			IF (tRef(k).NE.UNSET_RL) THEN
61	jmc	1.4	WRITE(msgBuf,'(2A,I4,A)') 'S/R LOAD_REF_FILES: ',
62			& 'Cannot set both tRef(k=', k, ') and tRefFile'
63	jmc	1.1	CALL PRINT_ERROR( msgBuf, myThid )
64			STOP 'ABNORMAL END: S/R INI_PARMS'
65			ENDIF
66			ENDDO
67			ENDIF
68			C- read from file:
69			IF ( tRefFile .NE. ' ' ) THEN
70			kLen = ILNBLNK(tRefFile)
71	jmc	1.2	CALL READ_GLVEC_RL( tRefFile, ' ', tRef, Nr, 1, myThid )
72	jmc	1.4	WRITE(msgBuf,'(3A)') 'S/R LOAD_REF_FILES: ',
73			& 'tRef loaded from file: ', tRefFile(1:kLen)
74	jmc	1.1	CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
75	jmc	1.4	& SQUEEZE_RIGHT, myThid )
76	jmc	1.1	ENDIF
77	jmc	1.3	IF ( thetaConst.EQ.UNSET_RL ) thetaConst = tRef(1)
78	jmc	1.1
79			C-- Set reference Salinity/Specific Humidity
80			IF ( sRefFile .EQ. ' ' ) THEN
81			C- set default vertical profile for salinity/water-vapour: sRef
82			tracerDefault = 30.
83			IF ( fluidIsAir ) tracerDefault = 0.
84			DO k=1,Nr
85			IF (sRef(k).EQ.UNSET_RL) sRef(k) = tracerDefault
86			tracerDefault = sRef(k)
87			ENDDO
88			ELSE
89			C- check for multiple definitions:
90			DO k=1,Nr
91			IF (sRef(k).NE.UNSET_RL) THEN
92	jmc	1.4	WRITE(msgBuf,'(2A,I4,A)') 'S/R LOAD_REF_FILES: ',
93			& 'Cannot set both sRef(k=', k, ') and sRefFile'
94	jmc	1.1	CALL PRINT_ERROR( msgBuf, myThid )
95			STOP 'ABNORMAL END: S/R INI_PARMS'
96			ENDIF
97			ENDDO
98			ENDIF
99			C- read from file:
100			IF ( sRefFile .NE. ' ' ) THEN
101			kLen = ILNBLNK(sRefFile)
102	jmc	1.2	CALL READ_GLVEC_RL( sRefFile, ' ', sRef, Nr, 1, myThid )
103	jmc	1.4	WRITE(msgBuf,'(3A)') 'S/R LOAD_REF_FILES: ',
104			& 'sRef loaded from file: ', sRefFile(1:kLen)
105	jmc	1.1	CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
106	jmc	1.4	& SQUEEZE_RIGHT, myThid )
107	jmc	1.1	ENDIF
108
109			C-- Set reference Density
110			IF ( rhoRefFile .NE. ' ' ) THEN
111	jmc	1.2	kLen = ILNBLNK(rhoRefFile)
112	jmc	1.1	C- read from file:
113	jmc	1.2	CALL READ_GLVEC_RL( rhoRefFile, ' ', rho1Ref, Nr, 1, myThid )
114	jmc	1.4	WRITE(msgBuf,'(3A)') 'S/R LOAD_REF_FILES: ',
115			& 'rho1Ref loaded from file: ', rhoRefFile(1:kLen)
116	jmc	1.1	CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
117	jmc	1.4	& SQUEEZE_RIGHT, myThid )
118			ENDIF
119
120			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
121
122			C-- Set gravity vertical profile
123			IF ( gravityFile .NE. ' ' ) THEN
124			kLen = ILNBLNK(gravityFile)
125			C- read from file and store, for now, in gravFacC:
126			CALL READ_GLVEC_RL( gravityFile, ' ', gravFacC, Nr, 1, myThid )
127			WRITE(msgBuf,'(3A)') 'S/R LOAD_REF_FILES: ',
128			& 'gravity vert. prof. loaded from file: ', gravityFile(1:kLen)
129			CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
130			& SQUEEZE_RIGHT, myThid )
131
132			C- Set gravity vertical profile factor: assume surface-interface
133			C density-factor to be 1 (grav. acceleration @ rF(1) = gravity)
134			C- Note: done here (instead of in set_ref_state.F) since gravity fac
135			C might be needed to initialise EOS coeffs (in ini_eos.F)
136			gravFacF(1) = 1. _d 0
137			C gravFac(k) = gravity ratio between layer k and top interface
138			DO k=1,Nr
139			gravFacC(k) = gravFacC(k)*recip_gravity
140			ENDDO
141			DO k=2,Nr
142			C since rkSign=-1, recip_drC(k) = 1./(rC(k-1)-rC(k))
143			gravFacF(k) = ( gravFacC(k-1)*(rF(k)-rC(k))
144			& + gravFacC(k)(rC(k-1)-rF(k)) )recip_drC(k)
145			ENDDO
146			C extrapolate down to the bottom:
147			gravFacF(Nr+1) = ( gravFacC(Nr)*(rF(Nr+1)-rF(Nr))
148			& + gravFacF(Nr)*(rC(Nr)-rF(Nr+1))
149			& ) / (rC(Nr)-rF(Nr))
150			C- set reciprocal gravity-factor:
151			DO k=1,Nr
152			recip_gravFacC(k) = 1. _d 0/gravFacC(k)
153			ENDDO
154			DO k=1,Nr+1
155			recip_gravFacF(k) = 1. _d 0/gravFacF(k)
156			ENDDO
157
158			ELSE
159			C- Initialise gravity vertical profile factor:
160			DO k=1,Nr
161			gravFacC(k) = 1. _d 0
162			recip_gravFacC(k) = 1. _d 0
163			ENDDO
164			DO k=1,Nr+1
165			gravFacF(k) = 1. _d 0
166			recip_gravFacF(k) = 1. _d 0
167			ENDDO
168	jmc	1.1	ENDIF
169
170			_END_MASTER(myThid)
171			C-- Everyone else must wait for the parameters to be loaded
172			_BARRIER
173
174			RETURN
175			END