pkg/rbcs/rbcs_fields_load.F

C $Header: /u/gcmpack/MITgcm/pkg/rbcs/rbcs_fields_load.F,v 1.7 2007/11/05 19:10:09 jmc Exp $
C $Name:  $

#include "CPP_OPTIONS.h"
#include "PACKAGES_CONFIG.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     !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
      INTEGER iTracer
      INTEGER  IFNBLNK, ILNBLNK
      EXTERNAL IFNBLNK, ILNBLNK

#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 (rbcsForcingCycle.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  = MOD( Imytm+Ifcyc-Ifprd/2, Ifcyc)

       intime0 = 1 + INT( Iftm/Ifprd )
       intime1 = 1 + MOD( intime0,nForcingPeriods )
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

       _BARRIER

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(*,*)
     &  'S/R RBCS_FIELDS_LOAD: Reading new data',myTime,myIter
     &            , nIter0, intime0,intime1
       _END_MASTER(myThid)

       IF ( relaxTFile      .NE. ' '  ) THEN
        CALL READ_REC_XYZ_RS(relaxTFile, rbct0, intime0, myIter,myThid)
        CALL READ_REC_XYZ_RS(relaxTFile, rbct1, intime1, myIter,myThid)
       ENDIF
       IF ( relaxSFile      .NE. ' '  ) THEN
        CALL READ_REC_XYZ_RS(relaxSFile, rbcs0, intime0, myIter,myThid)
        CALL READ_REC_XYZ_RS(relaxSFile, rbcs1, intime1, myIter,myThid)
       ENDIF

#ifdef ALLOW_PTRACERS
       IF (useRBCptracers) THEN
        DO iTracer = 1, PTRACERS_numInUse
         IF (useRBCptrnum(iTracer)) THEN
           CALL READ_REC_XYZ_RS( relaxPtracerFile(iTracer),
     &                           rbcptr0(1-Olx,1-Oly,1,1,1,iTracer),
     &                           intime0, myIter, myThid )
           CALL READ_REC_XYZ_RS( relaxPtracerFile(iTracer),
     &                           rbcptr1(1-Olx,1-Oly,1,1,1,iTracer),
     &                           intime1, myIter, myThid )
         ENDIF
        ENDDO
       ENDIF
#endif

       _EXCH_XYZ_R4(rbct0 , myThid )
       _EXCH_XYZ_R4(rbct1 , myThid )
       _EXCH_XYZ_R4(rbcs0 , myThid )
       _EXCH_XYZ_R4(rbcs1 , myThid )
#ifdef ALLOW_PTRACERS
       IF (useRBCptracers) THEN
        DO iTracer = 1, PTRACERS_numInUse
         _EXCH_XYZ_R4(rbcptr0(1-Olx,1-Oly,1,1,1,iTracer),myThid)
         _EXCH_XYZ_R4(rbcptr1(1-Olx,1-Oly,1,1,1,iTracer),myThid)
        ENDDO
       ENDIF
#endif

      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 (useRBCptracers) 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

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

#endif /* ALLOW_RBCS */

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