pkg/dic/dic_fields_load.F

C $Header: /u/gcmpack/MITgcm/pkg/dic/dic_fields_load.F,v 1.27 2008/04/24 21:35:53 gforget Exp $
C $Name:  $

#include "DIC_OPTIONS.h"

CBOP
C !ROUTINE: DIC_FIELDS_LOAD

C !INTERFACE: ==========================================================
      SUBROUTINE DIC_FIELDS_LOAD (
     I           myIter,myTime,myThid)

C !DESCRIPTION:
C  Read in fields needed for CO2,O2 fluxterms, silica for pH calculation

C !USES: ===============================================================
      IMPLICIT NONE
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#include "DIC_VARS.h"
#include "DIC_LOAD.h"
#ifdef ALLOW_THSICE
#include "THSICE_VARS.h"
#endif
#ifdef ALLOW_SEAICE
#include "SEAICE.h"
#endif

C !INPUT PARAMETERS: ===================================================
C  myThid               :: thread number
C  myIter               :: current timestep
C  myTime               :: current time
      INTEGER myIter
      _RL myTime
      INTEGER myThid

#ifdef ALLOW_PTRACERS

c !LOCAL VARIABLES: ===================================================
       INTEGER bi,bj,i,j,intime0,intime1
      _RL aWght,bWght,rdt
      INTEGER nForcingPeriods,Imytm,Ifprd,Ifcyc,Iftm
      CHARACTER*(MAX_LEN_MBUF) msgBuf
CEOP

c
      IF (  DIC_forcingCycle.gt.0. _d 0 ) THEN


C Now calculate whether it is time to update the forcing arrays
       rdt = 1. _d 0 / deltaTclock
       nForcingPeriods = NINT(DIC_forcingCycle/DIC_forcingPeriod)
cswd QQ change for placement of chem forcing (ie. after timestep)
       Imytm = NINT(myTime*rdt)
       Ifprd = NINT(DIC_forcingPeriod*rdt)
       Ifcyc = NINT(DIC_forcingCycle*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

       IF (
     &   Iftm-Ifprd*(intime0-1).EQ. 0
     &   .OR. myIter .EQ. nIter0
     &    ) THEN
C-     this is time to load new fields

        _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(standardMessageUnit,'(A,2I5,I10,1P1E20.12)')
     &  'S/R DIC_FIELDS_LOAD: Reading new dic data:',
     &  intime0, intime1, myIter, myTime
        _END_MASTER(myThid)


        IF ( DIC_windFile .NE. ' '  ) THEN
         CALL READ_REC_XY_RS( DIC_windFile,dicwind0,intime0,
     &        myIter,myThid )
         CALL READ_REC_XY_RS( DIC_windFile,dicwind1,intime1,
     &        myIter,myThid )
        ENDIF
        IF ( DIC_atmospFile .NE. ' '  ) THEN
         CALL READ_REC_XY_RS( DIC_atmospFile,atmosp0,intime0,
     &        myIter,myThid )
         CALL READ_REC_XY_RS( DIC_atmospFile,atmosp1,intime1,
     &        myIter,myThid )
        ENDIF
        IF ( DIC_silicaFile .NE. ' '  ) THEN
         CALL READ_REC_XY_RS( DIC_silicaFile,silica0,intime0,
     &        myIter,myThid )
         CALL READ_REC_XY_RS( DIC_silicaFile,silica1,intime1,
     &        myIter,myThid )
        ENDIF
        IF ( DIC_iceFile .NE. ' '  ) THEN
         CALL READ_REC_XY_RS( DIC_iceFile,ice0,intime0,
     &       myIter,myThid )
         CALL READ_REC_XY_RS( DIC_iceFile,ice1,intime1,
     &       myIter,myThid )
        ENDIF
#ifdef READ_PAR
        IF ( Filename1 .NE. ' '  ) THEN
         CALL READ_REC_XY_RS( Filename1,par0,intime0,
     &       myIter,myThid )
         CALL READ_REC_XY_RS( Filename1,par1,intime1,
     &       myIter,myThid )
        ENDIF
#endif
#ifdef ALLOW_FE
        IF ( DIC_ironFile .NE. ' '  ) THEN
         CALL READ_REC_XY_RS( DIC_ironFile,feinput0,intime0,
     &       myIter,myThid )
         CALL READ_REC_XY_RS( DIC_ironFile,feinput1,intime1,
     &       myIter,myThid )
        ENDIF
#endif
        _BARRIER

C-     end if-bloc (time to load new fields)
       ENDIF

       DO bj = myByLo(myThid), myByHi(myThid)
        DO bi = myBxLo(myThid), myBxHi(myThid)
         IF ( DIC_windFile .NE. ' '  ) THEN
           DO j=1-Oly,sNy+Oly
            DO i=1-Olx,sNx+Olx
             WIND(i,j,bi,bj) = bWght*dicwind0(i,j,bi,bj)
     &                       + aWght*dicwind1(i,j,bi,bj)
            ENDDO
           ENDDO
C calculate piston velocity
C QQ: note - we should have wind speed variance in here
C QQ         also need to check units, and conversion factors
c          pisvel(i,j,bi,bj)  =0.337*wind(i,j,bi,bj)**2/3.6d5    !QQQQ
         ENDIF
#ifndef USE_PLOAD
         IF ( DIC_atmospFile .NE. ' '  ) THEN
           DO j=1-Oly,sNy+Oly
            DO i=1-Olx,sNx+Olx
             AtmosP(i,j,bi,bj) = bWght*atmosp0(i,j,bi,bj)
     &                         + aWght*atmosp1(i,j,bi,bj)
            ENDDO
           ENDDO
         ENDIF
#endif
         IF ( DIC_silicaFile .NE. ' '  ) THEN
           DO j=1-Oly,sNy+Oly
            DO i=1-Olx,sNx+Olx
             SILICA(i,j,bi,bj) = bWght*silica0(i,j,bi,bj)
     &                         + aWght*silica1(i,j,bi,bj)
            ENDDO
           ENDDO
         ENDIF
         IF ( DIC_iceFile .NE. ' '  ) THEN
           DO j=1-Oly,sNy+Oly
            DO i=1-Olx,sNx+Olx
             FIce(i,j,bi,bj) = bWght*ice0(i,j,bi,bj)
     &                       + aWght*ice1(i,j,bi,bj)
            ENDDO
           ENDDO
         ENDIF

#ifdef READ_PAR
         IF ( Filename1 .NE. ' '  ) THEN
           DO j=1-Oly,sNy+Oly
            DO i=1-Olx,sNx+Olx
             PAR(i,j,bi,bj) = bWght*par0(i,j,bi,bj)
     &                      + aWght*par1(i,j,bi,bj)
            ENDDO
           ENDDO
         ELSE
            WRITE(msgBuf,'(2A)')
     &       ' DIC_FIELDS_LOAD: You need to provide ',
     &       ' a file if you want to use READ_PAR'
            CALL PRINT_ERROR( msgBuf, myThid )
            STOP 'ABNORMAL END: S/R DIC_FIELDS_LOAD'
         ENDIF
#endif
#ifdef ALLOW_FE
         IF ( DIC_ironFile .NE. ' '  ) THEN
           DO j=1-Oly,sNy+Oly
            DO i=1-Olx,sNx+Olx
             InputFe(i,j,bi,bj) = bWght*feinput0(i,j,bi,bj)
     &                          + aWght*feinput1(i,j,bi,bj)
            ENDDO
           ENDDO
         ENDIF
#endif
        ENDDO
       ENDDO

C endif for DIC_forcingCycle
      ENDIF

      DO bj = myByLo(myThid), myByHi(myThid)
       DO bi = myBxLo(myThid), myBxHi(myThid)
#ifdef ALLOW_SEAICE
         IF ( useSEAICE ) THEN
           DO j=1-Oly,sNy+Oly
            DO i=1-Olx,sNx+Olx
             FIce(i,j,bi,bj) = AREA(i,j,bi,bj)
            ENDDO
           ENDDO
         ENDIF
#endif
#ifdef ALLOW_THSICE
         IF ( useThSIce ) THEN
           DO j=1-Oly,sNy+Oly
            DO i=1-Olx,sNx+Olx
             FIce(i,j,bi,bj) = iceMask(i,j,bi,bj)
            ENDDO
           ENDDO
         ENDIF
#endif
       ENDDO
      ENDDO

#endif /* ALLOW_PTRACERS */
      RETURN
      END
1	dimitri	1.28	C $Header: /u/gcmpack/MITgcm/pkg/dic/dic_fields_load.F,v 1.27 2008/04/24 21:35:53 gforget Exp $
2	jmc	1.8	C $Name: $
3
4	edhill	1.4	#include "DIC_OPTIONS.h"
5	stephd	1.1
6	stephd	1.7	CBOP
7			C !ROUTINE: DIC_FIELDS_LOAD
8
9			C !INTERFACE: ==========================================================
10	stephd	1.1	SUBROUTINE DIC_FIELDS_LOAD (
11			I myIter,myTime,myThid)
12
13	stephd	1.7	C !DESCRIPTION:
14			C Read in fields needed for CO2,O2 fluxterms, silica for pH calculation
15
16			C !USES: ===============================================================
17	stephd	1.1	IMPLICIT NONE
18			#include "SIZE.h"
19			#include "EEPARAMS.h"
20			#include "PARAMS.h"
21			#include "GRID.h"
22	dfer	1.22	#include "DIC_VARS.h"
23	stephd	1.3	#include "DIC_LOAD.h"
24	stephd	1.14	#ifdef ALLOW_THSICE
25			#include "THSICE_VARS.h"
26			#endif
27	stephd	1.19	#ifdef ALLOW_SEAICE
28			#include "SEAICE.h"
29			#endif
30	stephd	1.1
31	stephd	1.7	C !INPUT PARAMETERS: ===================================================
32			C myThid :: thread number
33			C myIter :: current timestep
34			C myTime :: current time
35	stephd	1.1	INTEGER myIter
36			_RL myTime
37			INTEGER myThid
38
39			#ifdef ALLOW_PTRACERS
40	stephd	1.7
41			c !LOCAL VARIABLES: ===================================================
42	stephd	1.1	INTEGER bi,bj,i,j,intime0,intime1
43			_RL aWght,bWght,rdt
44			INTEGER nForcingPeriods,Imytm,Ifprd,Ifcyc,Iftm
45	stephd	1.13	CHARACTER*(MAX_LEN_MBUF) msgBuf
46	stephd	1.7	CEOP
47
48	stephd	1.1	c
49	jmc	1.25	IF ( DIC_forcingCycle.gt.0. _d 0 ) THEN
50	stephd	1.1
51
52			C Now calculate whether it is time to update the forcing arrays
53	jmc	1.26	rdt = 1. _d 0 / deltaTclock
54			nForcingPeriods = NINT(DIC_forcingCycle/DIC_forcingPeriod)
55	stephd	1.1	cswd QQ change for placement of chem forcing (ie. after timestep)
56	jmc	1.26	Imytm = NINT(myTime*rdt)
57			Ifprd = NINT(DIC_forcingPeriod*rdt)
58			Ifcyc = NINT(DIC_forcingCycle*rdt)
59			Iftm = MOD( Imytm+Ifcyc-Ifprd/2, Ifcyc)
60
61			intime0 = 1 + INT(Iftm/Ifprd)
62			intime1 = 1 + MOD(intime0,nForcingPeriods)
63			c aWght = DFLOAT( Iftm-Ifprd*(intime0 - 1) ) / DFLOAT( Ifprd )
64			aWght = FLOAT( Iftm-Ifprd*(intime0 - 1) )
65			bWght = FLOAT( Ifprd )
66			aWght = aWght / bWght
67			bWght = 1. _d 0 - aWght
68	stephd	1.1
69			IF (
70			& Iftm-Ifprd*(intime0-1).EQ. 0
71	stephd	1.3	& .OR. myIter .EQ. nIter0
72	stephd	1.1	& ) THEN
73	jmc	1.26	C- this is time to load new fields
74	stephd	1.1
75	jmc	1.11	_BARRIER
76	stephd	1.1
77			C If the above condition is met then we need to read in
78			C data for the period ahead and the period behind myTime.
79	jmc	1.11	_BEGIN_MASTER(myThid)
80			WRITE(standardMessageUnit,'(A,2I5,I10,1P1E20.12)')
81			& 'S/R DIC_FIELDS_LOAD: Reading new dic data:',
82			& intime0, intime1, myIter, myTime
83			_END_MASTER(myThid)
84
85	stephd	1.1
86	jmc	1.26	IF ( DIC_windFile .NE. ' ' ) THEN
87	gforget	1.27	CALL READ_REC_XY_RS( DIC_windFile,dicwind0,intime0,
88	stephd	1.1	& myIter,myThid )
89	gforget	1.27	CALL READ_REC_XY_RS( DIC_windFile,dicwind1,intime1,
90	stephd	1.1	& myIter,myThid )
91	jmc	1.26	ENDIF
92			IF ( DIC_atmospFile .NE. ' ' ) THEN
93	jmc	1.25	CALL READ_REC_XY_RS( DIC_atmospFile,atmosp0,intime0,
94	stephd	1.1	& myIter,myThid )
95	jmc	1.25	CALL READ_REC_XY_RS( DIC_atmospFile,atmosp1,intime1,
96	stephd	1.1	& myIter,myThid )
97	jmc	1.26	ENDIF
98			IF ( DIC_silicaFile .NE. ' ' ) THEN
99	jmc	1.25	CALL READ_REC_XY_RS( DIC_silicaFile,silica0,intime0,
100	stephd	1.2	& myIter,myThid )
101	jmc	1.25	CALL READ_REC_XY_RS( DIC_silicaFile,silica1,intime1,
102	stephd	1.2	& myIter,myThid )
103	jmc	1.26	ENDIF
104			IF ( DIC_iceFile .NE. ' ' ) THEN
105	jmc	1.25	CALL READ_REC_XY_RS( DIC_iceFile,ice0,intime0,
106	stephd	1.1	& myIter,myThid )
107	jmc	1.25	CALL READ_REC_XY_RS( DIC_iceFile,ice1,intime1,
108	stephd	1.1	& myIter,myThid )
109	jmc	1.26	ENDIF
110	stephd	1.9	#ifdef READ_PAR
111	jmc	1.26	IF ( Filename1 .NE. ' ' ) THEN
112	stephd	1.9	CALL READ_REC_XY_RS( Filename1,par0,intime0,
113			& myIter,myThid )
114			CALL READ_REC_XY_RS( Filename1,par1,intime1,
115			& myIter,myThid )
116	jmc	1.26	ENDIF
117	stephd	1.9	#endif
118	stephd	1.3	#ifdef ALLOW_FE
119	jmc	1.26	IF ( DIC_ironFile .NE. ' ' ) THEN
120	jmc	1.25	CALL READ_REC_XY_RS( DIC_ironFile,feinput0,intime0,
121	stephd	1.1	& myIter,myThid )
122	jmc	1.25	CALL READ_REC_XY_RS( DIC_ironFile,feinput1,intime1,
123	stephd	1.1	& myIter,myThid )
124	jmc	1.26	ENDIF
125	stephd	1.1	#endif
126	jmc	1.26	_BARRIER
127	stephd	1.1
128	jmc	1.26	C- end if-bloc (time to load new fields)
129	stephd	1.1	ENDIF
130
131			DO bj = myByLo(myThid), myByHi(myThid)
132			DO bi = myBxLo(myThid), myBxHi(myThid)
133	jmc	1.25	IF ( DIC_windFile .NE. ' ' ) THEN
134	dfer	1.16	DO j=1-Oly,sNy+Oly
135			DO i=1-Olx,sNx+Olx
136	gforget	1.27	WIND(i,j,bi,bj) = bWght*dicwind0(i,j,bi,bj)
137			& + aWght*dicwind1(i,j,bi,bj)
138	dfer	1.16	ENDDO
139			ENDDO
140	jmc	1.25	C calculate piston velocity
141			C QQ: note - we should have wind speed variance in here
142			C QQ also need to check units, and conversion factors
143	stephd	1.10	c pisvel(i,j,bi,bj) =0.337wind(i,j,bi,bj)*2/3.6d5 !QQQQ
144	dfer	1.16	ENDIF
145			#ifndef USE_PLOAD
146	jmc	1.25	IF ( DIC_atmospFile .NE. ' ' ) THEN
147	dfer	1.16	DO j=1-Oly,sNy+Oly
148			DO i=1-Olx,sNx+Olx
149			AtmosP(i,j,bi,bj) = bWght*atmosp0(i,j,bi,bj)
150			& + aWght*atmosp1(i,j,bi,bj)
151			ENDDO
152			ENDDO
153			ENDIF
154			#endif
155	jmc	1.25	IF ( DIC_silicaFile .NE. ' ' ) THEN
156	dfer	1.16	DO j=1-Oly,sNy+Oly
157			DO i=1-Olx,sNx+Olx
158			SILICA(i,j,bi,bj) = bWght*silica0(i,j,bi,bj)
159			& + aWght*silica1(i,j,bi,bj)
160			ENDDO
161			ENDDO
162			ENDIF
163	jmc	1.25	IF ( DIC_iceFile .NE. ' ' ) THEN
164	dfer	1.16	DO j=1-Oly,sNy+Oly
165			DO i=1-Olx,sNx+Olx
166			FIce(i,j,bi,bj) = bWght*ice0(i,j,bi,bj)
167			& + aWght*ice1(i,j,bi,bj)
168			ENDDO
169			ENDDO
170			ENDIF
171	stephd	1.20
172	stephd	1.9	#ifdef READ_PAR
173	dfer	1.16	IF ( Filename1 .NE. ' ' ) THEN
174			DO j=1-Oly,sNy+Oly
175			DO i=1-Olx,sNx+Olx
176			PAR(i,j,bi,bj) = bWght*par0(i,j,bi,bj)
177			& + aWght*par1(i,j,bi,bj)
178			ENDDO
179			ENDDO
180			ELSE
181	stephd	1.13	WRITE(msgBuf,'(2A)')
182			& ' DIC_FIELDS_LOAD: You need to provide ',
183			& ' a file if you want to use READ_PAR'
184			CALL PRINT_ERROR( msgBuf, myThid )
185			STOP 'ABNORMAL END: S/R DIC_FIELDS_LOAD'
186	dfer	1.16	ENDIF
187	stephd	1.9	#endif
188	stephd	1.3	#ifdef ALLOW_FE
189	jmc	1.25	IF ( DIC_ironFile .NE. ' ' ) THEN
190	dfer	1.16	DO j=1-Oly,sNy+Oly
191			DO i=1-Olx,sNx+Olx
192			InputFe(i,j,bi,bj) = bWght*feinput0(i,j,bi,bj)
193			& + aWght*feinput1(i,j,bi,bj)
194			ENDDO
195			ENDDO
196			ENDIF
197	stephd	1.1	#endif
198			ENDDO
199			ENDDO
200
201	jmc	1.25	C endif for DIC_forcingCycle
202	jmc	1.11	ENDIF
203	stephd	1.1
204	jmc	1.26	DO bj = myByLo(myThid), myByHi(myThid)
205			DO bi = myBxLo(myThid), myBxHi(myThid)
206			#ifdef ALLOW_SEAICE
207			IF ( useSEAICE ) THEN
208	stephd	1.20	DO j=1-Oly,sNy+Oly
209			DO i=1-Olx,sNx+Olx
210	dimitri	1.28	FIce(i,j,bi,bj) = AREA(i,j,bi,bj)
211	stephd	1.20	ENDDO
212			ENDDO
213	jmc	1.26	ENDIF
214			#endif
215			#ifdef ALLOW_THSICE
216			IF ( useThSIce ) THEN
217	stephd	1.20	DO j=1-Oly,sNy+Oly
218			DO i=1-Olx,sNx+Olx
219	jmc	1.26	FIce(i,j,bi,bj) = iceMask(i,j,bi,bj)
220	stephd	1.20	ENDDO
221			ENDDO
222	jmc	1.26	ENDIF
223	stephd	1.20	#endif
224	jmc	1.26	ENDDO
225			ENDDO
226	stephd	1.20
227	dfer	1.16	#endif /* ALLOW_PTRACERS */
228	stephd	1.1	RETURN
229			END