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	C $Header: /u/gcmpack/MITgcm/model/src/load_ref_files.F,v 1.3 2014/05/06 15:47:05 jmc Exp $
2	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	#include "GRID.h"
26
27	C !INPUT/OUTPUT PARAMETERS:
28	C == Routine arguments ==
29	C myThid :: my Thread Id number
30	INTEGER myThid
31
32	C !FUNCTIONS:
33	INTEGER ILNBLNK
34	EXTERNAL ILNBLNK
35
36	C !LOCAL VARIABLES:
37	C == Local variables ==
38	C k :: loop index
39	C msgBuf :: Informational/error message buffer
40	_RL tracerDefault
41	INTEGER k, kLen
42	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	IF ( thetaConst.NE.UNSET_RL ) tracerDefault = thetaConst
53	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	WRITE(msgBuf,'(2A,I4,A)') 'S/R LOAD_REF_FILES: ',
62	& 'Cannot set both tRef(k=', k, ') and tRefFile'
63	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	CALL READ_GLVEC_RL( tRefFile, ' ', tRef, Nr, 1, myThid )
72	WRITE(msgBuf,'(3A)') 'S/R LOAD_REF_FILES: ',
73	& 'tRef loaded from file: ', tRefFile(1:kLen)
74	CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
75	& SQUEEZE_RIGHT, myThid )
76	ENDIF
77	IF ( thetaConst.EQ.UNSET_RL ) thetaConst = tRef(1)
78
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	WRITE(msgBuf,'(2A,I4,A)') 'S/R LOAD_REF_FILES: ',
93	& 'Cannot set both sRef(k=', k, ') and sRefFile'
94	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	CALL READ_GLVEC_RL( sRefFile, ' ', sRef, Nr, 1, myThid )
103	WRITE(msgBuf,'(3A)') 'S/R LOAD_REF_FILES: ',
104	& 'sRef loaded from file: ', sRefFile(1:kLen)
105	CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
106	& SQUEEZE_RIGHT, myThid )
107	ENDIF
108
109	C-- Set reference Density
110	IF ( rhoRefFile .NE. ' ' ) THEN
111	kLen = ILNBLNK(rhoRefFile)
112	C- read from file:
113	CALL READ_GLVEC_RL( rhoRefFile, ' ', rho1Ref, Nr, 1, myThid )
114	WRITE(msgBuf,'(3A)') 'S/R LOAD_REF_FILES: ',
115	& 'rho1Ref loaded from file: ', rhoRefFile(1:kLen)
116	CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
117	& 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	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