high_res_cube/rbcs/rbcs_fields_load.F

C $Header: /u/gcmpack/MITgcm_contrib/high_res_cube/rbcs/rbcs_fields_load.F,v 1.2 2010/10/14 22:32:16 dimitri Exp $
C $Name:  $

#include "RBCS_OPTIONS.h"

C     !ROUTINE: RBCS_FIELDS_LOAD
C     !INTERFACE:
      SUBROUTINE RBCS_FIELDS_LOAD( myTime, myIter, myThid )
C     *==========================================================*
C     | SUBROUTINE RBCS_FIELDS_LOAD
C     | o Control reading of fields from external source.
C     *==========================================================*
C     | Offline External source field loading routine.
C     | This routine is called every time we want to
C     | load a a set of external fields. The routine decides
C     | which fields to load and then reads them in.
C     | This routine needs to be customised for particular
C     | experiments.
C     *==========================================================*

C     !USES:
      IMPLICIT NONE
C     === Global variables ===
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "FFIELDS.h"
#include "GRID.h"
#include "DYNVARS.h"
#ifdef ALLOW_PTRACERS
#include "PTRACERS_SIZE.h"
#include "PTRACERS_PARAMS.h"
#endif
#include "RBCS.h"

C     !INPUT/OUTPUT PARAMETERS:
C     === Routine arguments ===
C     myTime :: Simulation time
C     myIter :: Simulation timestep number
C     myThid :: Thread no. that called this routine.
      _RL     myTime
      INTEGER myIter
      INTEGER myThid

C     !FUNCTIONS:
      INTEGER  IFNBLNK, ILNBLNK
      EXTERNAL IFNBLNK, ILNBLNK

C     !LOCAL VARIABLES:
C     === Local arrays ===
C     [01]      :: End points for interpolation
C     Above use static heap storage to allow exchange.
C     aWght, bWght :: Interpolation weights

      INTEGER bi,bj,i,j,k,intime0,intime1
      _RL aWght,bWght,rdt
      INTEGER nForcingPeriods,Imytm,Ifprd,Ifcyc,Iftm
#ifdef ALLOW_PTRACERS
      INTEGER iTracer
#endif
      INTEGER IL, inIter
      CHARACTER*(MAX_LEN_FNAM) fullName

C     deltaTclock with which the rbcs input files were generated
C     right now hardcoded for output from CS510 integration
      _RL         rbcsTClock
      PARAMETER ( rbcsTclock = 1200. )

#ifdef ALLOW_RBCS
      CALL TIMER_START('RBCS_FIELDS_LOAD      [I/O]', myThid)

C First call requires that we initialize everything to zero for safety
      IF ( myIter .EQ. nIter0 ) THEN
       DO bj = myByLo(myThid), myByHi(myThid)
        DO bi = myBxLo(myThid), myBxHi(myThid)
         DO k=1,Nr
          DO j=1-Oly,sNy+Oly
           DO i=1-Olx,sNx+Olx
             rbct0(i,j,k,bi,bj)=0. _d 0
             rbcs0(i,j,k,bi,bj)=0. _d 0
             rbct1(i,j,k,bi,bj)=0. _d 0
             rbcs1(i,j,k,bi,bj)=0. _d 0
           ENDDO
          ENDDO
         ENDDO
        ENDDO
       ENDDO
#ifdef ALLOW_PTRACERS
       DO iTracer = 1, PTRACERS_numInUse
        DO bj = myByLo(myThid), myByHi(myThid)
         DO bi = myBxLo(myThid), myBxHi(myThid)
          DO k=1,Nr
           DO j=1-Oly,sNy+Oly
            DO i=1-Olx,sNx+Olx
             rbcptr0(i,j,k,bi,bj,iTracer)=0. _d 0
             rbcptr1(i,j,k,bi,bj,iTracer)=0. _d 0
            ENDDO
           ENDDO
          ENDDO
         ENDDO
        ENDDO
       ENDDO
#endif
      ENDIF

C Now calculate whether it is time to update the forcing arrays
      IF (rbcsForcingPeriod.GT.0. _d 0) THEN
       rdt = 1. _d 0 / deltaTclock
       nForcingPeriods = NINT(rbcsForcingCycle/rbcsForcingPeriod)
       Imytm = NINT( (myTime-rbcsIniter*deltaTclock)*rdt )
       Ifprd = NINT(rbcsForcingPeriod*rdt)
       Ifcyc = NINT(rbcsForcingCycle*rdt)
       Iftm  = Imytm-Ifprd/2
       IF (rbcsForcingCycle.GT.0. _d 0)
     &    Iftm  = MOD(Iftm+Ifcyc, Ifcyc)

       intime0 = INT( (Iftm+Ifprd)/Ifprd )
       IF (rbcsForcingCycle.GT.0. _d 0) THEN
        intime1 = 1 + MOD( intime0,nForcingPeriods )
       ELSE
        intime1 = 1 + intime0
       ENDIF
c      aWght = DFLOAT( Iftm-Ifprd*(intime0 - 1) ) / DFLOAT( Ifprd )
       aWght = FLOAT( Iftm-Ifprd*(intime0 - 1) )
       bWght = FLOAT( Ifprd )
       aWght =  aWght / bWght
       bWght = 1. _d 0 - aWght

      ELSE
       intime1 = 1
       intime0 = 1
       Iftm  = 1
       Ifprd = 0
       aWght = .5 _d 0
       bWght = .5 _d 0
      ENDIF

      IF (
     &  Iftm-Ifprd*(intime0-1) .EQ. 0
     &  .OR. myIter .EQ. nIter0
     & ) THEN

C      If the above condition is met then we need to read in
C      data for the period ahead and the period behind myTime.
       _BEGIN_MASTER(myThid)
        WRITE(standardMessageUnit,'(A,2I5,I10,1P1E20.12)')
     &        'S/R RBCS_FIELDS_LOAD: Reading new data:',
     &         intime0, intime1, myIter, myTime
       _END_MASTER(myThid)

       IF ( useRBCtemp .AND. relaxTFile .NE. ' '  ) THEN
        IL=ILNBLNK( relaxTFile )
C       assume files carry iter at end of period in name
        inIter = rbcsIniter + intime0*rbcsForcingPeriod/rbcsTclock
        WRITE(fullName,'(2a,i10.10)') relaxTFile(1:IL),'.',inIter
        CALL READ_REC_XYZ_RS(fullName, rbct0, 1, myIter,myThid)
C
        inIter = rbcsIniter + intime1*rbcsForcingPeriod/rbcsTclock
        WRITE(fullName,'(2a,i10.10)') relaxTFile(1:IL),'.',inIter
        CALL READ_REC_XYZ_RS(fullName, rbct1, 1, myIter,myThid)
       ENDIF
       IF ( useRBCsalt .AND. relaxSFile .NE. ' '  ) THEN
        IL=ILNBLNK( relaxSFile )
        inIter = rbcsIniter + intime0*rbcsForcingPeriod/rbcsTclock
        WRITE(fullName,'(2a,i10.10)') relaxSFile(1:IL),'.',inIter
        CALL READ_REC_XYZ_RS(fullName, rbcs0, 1, myIter,myThid)
        inIter = rbcsIniter + intime1*rbcsForcingPeriod/rbcsTclock
        WRITE(fullName,'(2a,i10.10)') relaxSFile(1:IL),'.',inIter
        CALL READ_REC_XYZ_RS(fullName, rbcs1, 1, myIter,myThid)
       ENDIF

#ifdef ALLOW_PTRACERS
       IF ( usePTRACERS ) THEN
        DO iTracer = 1, PTRACERS_numInUse
         IF ( useRBCptrnum(iTracer) .AND.
     &        relaxPtracerFile(iTracer).NE. ' ' ) THEN
           IL=ILNBLNK( relaxPtracerFile(iTracer) )
           inIter=rbcsIniter+intime0*rbcsForcingPeriod/rbcsTclock
           WRITE(fullName,'(2a,i10.10)') relaxPtracerFile(iTracer)(1:IL)
     &                                   ,'.',inIter
           CALL READ_REC_XYZ_RS( fullName,
     &                           rbcptr0(1-Olx,1-Oly,1,1,1,iTracer),
     &                           1, myIter, myThid )
           inIter=rbcsIniter+intime1*rbcsForcingPeriod/rbcsTclock
           WRITE(fullName,'(2a,i10.10)') relaxPtracerFile(iTracer)(1:IL)
     &                                   ,'.',inIter
           CALL READ_REC_XYZ_RS( fullName,
     &                           rbcptr1(1-Olx,1-Oly,1,1,1,iTracer),
     &                           1, myIter, myThid )
         ENDIF
        ENDDO
       ENDIF
#endif

       IF ( useRBCtemp .AND. relaxTFile .NE. ' '  ) THEN
         CALL EXCH_XYZ_RS( rbct0 , myThid )
         CALL EXCH_XYZ_RS( rbct1 , myThid )
       ENDIF
       IF ( useRBCsalt .AND. relaxSFile .NE. ' '  ) THEN
         CALL EXCH_XYZ_RS( rbcs0 , myThid )
         CALL EXCH_XYZ_RS( rbcs1 , myThid )
       ENDIF
#ifdef ALLOW_PTRACERS
       IF (usePTRACERS) THEN
        DO iTracer = 1, PTRACERS_numInUse
         IF ( useRBCptrnum(iTracer) ) THEN
          CALL EXCH_XYZ_RS( rbcptr0(1-Olx,1-Oly,1,1,1,iTracer),myThid )
          CALL EXCH_XYZ_RS( rbcptr1(1-Olx,1-Oly,1,1,1,iTracer),myThid )
         ENDIF
        ENDDO
       ENDIF
#endif /* ALLOW_PTRACERS */

      ENDIF

C--   Interpolate
      DO bj = myByLo(myThid), myByHi(myThid)
       DO bi = myBxLo(myThid), myBxHi(myThid)
        DO k=1,Nr
         DO j=1-Oly,sNy+Oly
          DO i=1-Olx,sNx+Olx
           RBCtemp(i,j,k,bi,bj) = bWght*rbct0(i,j,k,bi,bj)
     &                           +aWght*rbct1(i,j,k,bi,bj)
           RBCsalt(i,j,k,bi,bj) = bWght*rbcs0(i,j,k,bi,bj)
     &                           +aWght*rbcs1(i,j,k,bi,bj)
          ENDDO
         ENDDO
        ENDDO
       ENDDO
      ENDDO

#ifdef ALLOW_PTRACERS
      IF ( usePTRACERS ) THEN
       DO iTracer = 1, PTRACERS_numInUse
        IF (useRBCptrnum(iTracer)) THEN
         DO bj = myByLo(myThid), myByHi(myThid)
          DO bi = myBxLo(myThid), myBxHi(myThid)
           DO k=1,Nr
            DO j=1-Oly,sNy+Oly
             DO i=1-Olx,sNx+Olx
               RBC_ptracers(i,j,k,bi,bj,iTracer) =
     &                      bWght*rbcptr0(i,j,k,bi,bj,iTracer)
     &                     +aWght*rbcptr1(i,j,k,bi,bj,iTracer)
             ENDDO
            ENDDO
           ENDDO
          ENDDO
         ENDDO
        ENDIF
       ENDDO
      ENDIF
#endif /* ALLOW_PTRACERS */

      CALL TIMER_STOP ('RBCS_FIELDS_LOAD      [I/O]', myThid)

#endif /* ALLOW_RBCS */

      RETURN
      END
1	dimitri	1.3	C $Header: /u/gcmpack/MITgcm_contrib/high_res_cube/rbcs/rbcs_fields_load.F,v 1.2 2010/10/14 22:32:16 dimitri Exp $
2	dimitri	1.2	C $Name: $
3	dimitri	1.1
4	dimitri	1.2	#include "RBCS_OPTIONS.h"
5	dimitri	1.1
6			C !ROUTINE: RBCS_FIELDS_LOAD
7			C !INTERFACE:
8			SUBROUTINE RBCS_FIELDS_LOAD( myTime, myIter, myThid )
9			C ==========================================================
10			C \| SUBROUTINE RBCS_FIELDS_LOAD
11			C \| o Control reading of fields from external source.
12			C ==========================================================
13			C \| Offline External source field loading routine.
14			C \| This routine is called every time we want to
15			C \| load a a set of external fields. The routine decides
16			C \| which fields to load and then reads them in.
17			C \| This routine needs to be customised for particular
18			C \| experiments.
19			C ==========================================================
20
21			C !USES:
22			IMPLICIT NONE
23			C === Global variables ===
24			#include "SIZE.h"
25			#include "EEPARAMS.h"
26			#include "PARAMS.h"
27			#include "FFIELDS.h"
28			#include "GRID.h"
29			#include "DYNVARS.h"
30			#ifdef ALLOW_PTRACERS
31			#include "PTRACERS_SIZE.h"
32			#include "PTRACERS_PARAMS.h"
33			#endif
34			#include "RBCS.h"
35
36			C !INPUT/OUTPUT PARAMETERS:
37			C === Routine arguments ===
38			C myTime :: Simulation time
39			C myIter :: Simulation timestep number
40			C myThid :: Thread no. that called this routine.
41			_RL myTime
42			INTEGER myIter
43			INTEGER myThid
44
45	dimitri	1.2	C !FUNCTIONS:
46			INTEGER IFNBLNK, ILNBLNK
47			EXTERNAL IFNBLNK, ILNBLNK
48
49	dimitri	1.1	C !LOCAL VARIABLES:
50			C === Local arrays ===
51			C [01] :: End points for interpolation
52			C Above use static heap storage to allow exchange.
53			C aWght, bWght :: Interpolation weights
54
55			INTEGER bi,bj,i,j,k,intime0,intime1
56			_RL aWght,bWght,rdt
57			INTEGER nForcingPeriods,Imytm,Ifprd,Ifcyc,Iftm
58	dimitri	1.2	#ifdef ALLOW_PTRACERS
59	dimitri	1.1	INTEGER iTracer
60	dimitri	1.2	#endif
61	dimitri	1.1	INTEGER IL, inIter
62			CHARACTER*(MAX_LEN_FNAM) fullName
63
64	dimitri	1.3	C deltaTclock with which the rbcs input files were generated
65			C right now hardcoded for output from CS510 integration
66			_RL rbcsTClock
67			PARAMETER ( rbcsTclock = 1200. )
68
69	dimitri	1.1	#ifdef ALLOW_RBCS
70			CALL TIMER_START('RBCS_FIELDS_LOAD [I/O]', myThid)
71
72			C First call requires that we initialize everything to zero for safety
73			IF ( myIter .EQ. nIter0 ) THEN
74			DO bj = myByLo(myThid), myByHi(myThid)
75	dimitri	1.2	DO bi = myBxLo(myThid), myBxHi(myThid)
76			DO k=1,Nr
77			DO j=1-Oly,sNy+Oly
78			DO i=1-Olx,sNx+Olx
79			rbct0(i,j,k,bi,bj)=0. _d 0
80			rbcs0(i,j,k,bi,bj)=0. _d 0
81			rbct1(i,j,k,bi,bj)=0. _d 0
82			rbcs1(i,j,k,bi,bj)=0. _d 0
83			ENDDO
84	dimitri	1.1	ENDDO
85			ENDDO
86			ENDDO
87			ENDDO
88			#ifdef ALLOW_PTRACERS
89			DO iTracer = 1, PTRACERS_numInUse
90			DO bj = myByLo(myThid), myByHi(myThid)
91	dimitri	1.2	DO bi = myBxLo(myThid), myBxHi(myThid)
92			DO k=1,Nr
93			DO j=1-Oly,sNy+Oly
94			DO i=1-Olx,sNx+Olx
95			rbcptr0(i,j,k,bi,bj,iTracer)=0. _d 0
96			rbcptr1(i,j,k,bi,bj,iTracer)=0. _d 0
97			ENDDO
98	dimitri	1.1	ENDDO
99			ENDDO
100			ENDDO
101			ENDDO
102			ENDDO
103			#endif
104			ENDIF
105
106			C Now calculate whether it is time to update the forcing arrays
107			IF (rbcsForcingPeriod.GT.0. _d 0) THEN
108			rdt = 1. _d 0 / deltaTclock
109			nForcingPeriods = NINT(rbcsForcingCycle/rbcsForcingPeriod)
110			Imytm = NINT( (myTime-rbcsIniterdeltaTclock)rdt )
111			Ifprd = NINT(rbcsForcingPeriod*rdt)
112			Ifcyc = NINT(rbcsForcingCycle*rdt)
113			Iftm = Imytm-Ifprd/2
114			IF (rbcsForcingCycle.GT.0. _d 0)
115			& Iftm = MOD(Iftm+Ifcyc, Ifcyc)
116
117			intime0 = INT( (Iftm+Ifprd)/Ifprd )
118			IF (rbcsForcingCycle.GT.0. _d 0) THEN
119			intime1 = 1 + MOD( intime0,nForcingPeriods )
120			ELSE
121			intime1 = 1 + intime0
122			ENDIF
123			c aWght = DFLOAT( Iftm-Ifprd*(intime0 - 1) ) / DFLOAT( Ifprd )
124			aWght = FLOAT( Iftm-Ifprd*(intime0 - 1) )
125			bWght = FLOAT( Ifprd )
126			aWght = aWght / bWght
127			bWght = 1. _d 0 - aWght
128
129			ELSE
130			intime1 = 1
131			intime0 = 1
132			Iftm = 1
133			Ifprd = 0
134			aWght = .5 _d 0
135			bWght = .5 _d 0
136			ENDIF
137
138			IF (
139			& Iftm-Ifprd*(intime0-1) .EQ. 0
140			& .OR. myIter .EQ. nIter0
141			& ) THEN
142
143			C If the above condition is met then we need to read in
144			C data for the period ahead and the period behind myTime.
145			_BEGIN_MASTER(myThid)
146	dimitri	1.2	WRITE(standardMessageUnit,'(A,2I5,I10,1P1E20.12)')
147			& 'S/R RBCS_FIELDS_LOAD: Reading new data:',
148			& intime0, intime1, myIter, myTime
149	dimitri	1.1	_END_MASTER(myThid)
150
151	dimitri	1.2	IF ( useRBCtemp .AND. relaxTFile .NE. ' ' ) THEN
152	dimitri	1.1	IL=ILNBLNK( relaxTFile )
153			C assume files carry iter at end of period in name
154	dimitri	1.3	inIter = rbcsIniter + intime0*rbcsForcingPeriod/rbcsTclock
155	dimitri	1.1	WRITE(fullName,'(2a,i10.10)') relaxTFile(1:IL),'.',inIter
156			CALL READ_REC_XYZ_RS(fullName, rbct0, 1, myIter,myThid)
157			C
158	dimitri	1.3	inIter = rbcsIniter + intime1*rbcsForcingPeriod/rbcsTclock
159	dimitri	1.1	WRITE(fullName,'(2a,i10.10)') relaxTFile(1:IL),'.',inIter
160			CALL READ_REC_XYZ_RS(fullName, rbct1, 1, myIter,myThid)
161			ENDIF
162	dimitri	1.2	IF ( useRBCsalt .AND. relaxSFile .NE. ' ' ) THEN
163			IL=ILNBLNK( relaxSFile )
164	dimitri	1.3	inIter = rbcsIniter + intime0*rbcsForcingPeriod/rbcsTclock
165	dimitri	1.1	WRITE(fullName,'(2a,i10.10)') relaxSFile(1:IL),'.',inIter
166			CALL READ_REC_XYZ_RS(fullName, rbcs0, 1, myIter,myThid)
167	dimitri	1.3	inIter = rbcsIniter + intime1*rbcsForcingPeriod/rbcsTclock
168	dimitri	1.1	WRITE(fullName,'(2a,i10.10)') relaxSFile(1:IL),'.',inIter
169			CALL READ_REC_XYZ_RS(fullName, rbcs1, 1, myIter,myThid)
170			ENDIF
171
172			#ifdef ALLOW_PTRACERS
173	dimitri	1.2	IF ( usePTRACERS ) THEN
174	dimitri	1.1	DO iTracer = 1, PTRACERS_numInUse
175	dimitri	1.2	IF ( useRBCptrnum(iTracer) .AND.
176			& relaxPtracerFile(iTracer).NE. ' ' ) THEN
177			IL=ILNBLNK( relaxPtracerFile(iTracer) )
178	dimitri	1.3	inIter=rbcsIniter+intime0*rbcsForcingPeriod/rbcsTclock
179	dimitri	1.1	WRITE(fullName,'(2a,i10.10)') relaxPtracerFile(iTracer)(1:IL)
180			& ,'.',inIter
181			CALL READ_REC_XYZ_RS( fullName,
182			& rbcptr0(1-Olx,1-Oly,1,1,1,iTracer),
183			& 1, myIter, myThid )
184	dimitri	1.3	inIter=rbcsIniter+intime1*rbcsForcingPeriod/rbcsTclock
185	dimitri	1.1	WRITE(fullName,'(2a,i10.10)') relaxPtracerFile(iTracer)(1:IL)
186			& ,'.',inIter
187			CALL READ_REC_XYZ_RS( fullName,
188			& rbcptr1(1-Olx,1-Oly,1,1,1,iTracer),
189			& 1, myIter, myThid )
190			ENDIF
191			ENDDO
192			ENDIF
193			#endif
194
195	dimitri	1.2	IF ( useRBCtemp .AND. relaxTFile .NE. ' ' ) THEN
196			CALL EXCH_XYZ_RS( rbct0 , myThid )
197			CALL EXCH_XYZ_RS( rbct1 , myThid )
198			ENDIF
199			IF ( useRBCsalt .AND. relaxSFile .NE. ' ' ) THEN
200			CALL EXCH_XYZ_RS( rbcs0 , myThid )
201			CALL EXCH_XYZ_RS( rbcs1 , myThid )
202			ENDIF
203	dimitri	1.1	#ifdef ALLOW_PTRACERS
204	dimitri	1.2	IF (usePTRACERS) THEN
205	dimitri	1.1	DO iTracer = 1, PTRACERS_numInUse
206	dimitri	1.2	IF ( useRBCptrnum(iTracer) ) THEN
207			CALL EXCH_XYZ_RS( rbcptr0(1-Olx,1-Oly,1,1,1,iTracer),myThid )
208			CALL EXCH_XYZ_RS( rbcptr1(1-Olx,1-Oly,1,1,1,iTracer),myThid )
209			ENDIF
210	dimitri	1.1	ENDDO
211			ENDIF
212	dimitri	1.2	#endif /* ALLOW_PTRACERS */
213	dimitri	1.1
214			ENDIF
215
216			C-- Interpolate
217			DO bj = myByLo(myThid), myByHi(myThid)
218			DO bi = myBxLo(myThid), myBxHi(myThid)
219			DO k=1,Nr
220			DO j=1-Oly,sNy+Oly
221			DO i=1-Olx,sNx+Olx
222			RBCtemp(i,j,k,bi,bj) = bWght*rbct0(i,j,k,bi,bj)
223			& +aWght*rbct1(i,j,k,bi,bj)
224			RBCsalt(i,j,k,bi,bj) = bWght*rbcs0(i,j,k,bi,bj)
225			& +aWght*rbcs1(i,j,k,bi,bj)
226			ENDDO
227			ENDDO
228			ENDDO
229			ENDDO
230			ENDDO
231
232			#ifdef ALLOW_PTRACERS
233	dimitri	1.2	IF ( usePTRACERS ) THEN
234	dimitri	1.1	DO iTracer = 1, PTRACERS_numInUse
235			IF (useRBCptrnum(iTracer)) THEN
236			DO bj = myByLo(myThid), myByHi(myThid)
237			DO bi = myBxLo(myThid), myBxHi(myThid)
238			DO k=1,Nr
239			DO j=1-Oly,sNy+Oly
240			DO i=1-Olx,sNx+Olx
241			RBC_ptracers(i,j,k,bi,bj,iTracer) =
242			& bWght*rbcptr0(i,j,k,bi,bj,iTracer)
243			& +aWght*rbcptr1(i,j,k,bi,bj,iTracer)
244			ENDDO
245			ENDDO
246			ENDDO
247			ENDDO
248			ENDDO
249			ENDIF
250			ENDDO
251			ENDIF
252	dimitri	1.2	#endif /* ALLOW_PTRACERS */
253	dimitri	1.1
254			CALL TIMER_STOP ('RBCS_FIELDS_LOAD [I/O]', myThid)
255
256			#endif /* ALLOW_RBCS */
257
258			RETURN
259			END