pkg/dic/dic_fields_load.F

C $Header: /u/gcmpack/MITgcm/pkg/dic/dic_fields_load.F,v 1.14 2007/08/09 19:51:18 stephd Exp $
C $Name:  $

#include "DIC_OPTIONS.h"
#include "GCHEM_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 "DYNVARS.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#include "PTRACERS_SIZE.h"
#include "PTRACERS.h"
#include "GCHEM.h"
#include "DIC_ABIOTIC.h"
#include "DIC_BIOTIC.h"
#include "DIC_LOAD.h"
#ifdef ALLOW_THSICE
#include "THSICE_VARS.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 ( periodicExternalForcing ) THEN


C Now calculate whether it is time to update the forcing arrays
       rdt=1. _d 0 / deltaTclock
       nForcingPeriods=
     &  int(externForcingCycle/externForcingPeriod+0.5)
cswd QQ change for placement of chem forcing (ie. after timestep)
       Imytm=int(myTime*rdt+0.5)
       Ifprd=int(externForcingPeriod*rdt+0.5)
       Ifcyc=int(externForcingCycle*rdt+0.5)
       Iftm=mod( Imytm+Ifcyc-Ifprd/2 ,Ifcyc)


       intime0=int(Iftm/Ifprd)
       intime1=mod(intime0+1,nForcingPeriods)
       aWght=float( Iftm-Ifprd*intime0 )/float( Ifprd )
       bWght=1.-aWght

       intime0=intime0+1
       intime1=intime1+1


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


       IF ( WindFile .NE. ' '  ) THEN
         CALL READ_REC_XY_RS( WindFile,wspeed0,intime0,
     &        myIter,myThid )
         CALL READ_REC_XY_RS( WindFile,wspeed1,intime1,
     &        myIter,myThid )
       ENDIF
       IF ( AtmospFile .NE. ' '  ) THEN
         CALL READ_REC_XY_RS( AtmospFile,atmosp0,intime0,
     &        myIter,myThid )
         CALL READ_REC_XY_RS( AtmospFile,atmosp1,intime1,
     &        myIter,myThid )
       ENDIF
       IF ( SilicaFile .NE. ' '  ) THEN
         CALL READ_REC_XY_RS( SilicaFile,silica0,intime0,
     &        myIter,myThid )
         CALL READ_REC_XY_RS( SilicaFile,silica1,intime1,
     &        myIter,myThid )
       ENDIF
       IF ( IceFile .NE. ' '  ) THEN
         CALL READ_REC_XY_RS( IceFile,ice0,intime0,
     &       myIter,myThid )
         CALL READ_REC_XY_RS( 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 ( IronFile .NE. ' '  ) THEN
         CALL READ_REC_XY_RS( IronFile,feinput0,intime0,
     &       myIter,myThid )
         CALL READ_REC_XY_RS( IronFile,feinput1,intime1,
     &       myIter,myThid )
       ENDIF
#endif

       _EXCH_XY_R4(wspeed0, myThid )
       _EXCH_XY_R4(wspeed1, myThid )
       _EXCH_XY_R4(atmosp0, myThid )
       _EXCH_XY_R4(atmosp1, myThid )
       _EXCH_XY_R4(silica0, myThid )
       _EXCH_XY_R4(silica1, myThid )
       _EXCH_XY_R4(ice0, myThid )
       _EXCH_XY_R4(ice1, myThid )
#ifdef READ_PAR
       _EXCH_XY_R4(par0, myThid )
       _EXCH_XY_R4(par1, myThid )
#endif
#ifdef ALLOW_FE
       _EXCH_XY_R4(feinput0, myThid )
       _EXCH_XY_R4(feinput1, myThid )
#endif

C
       ENDIF

       DO bj = myByLo(myThid), myByHi(myThid)
        DO bi = myBxLo(myThid), myBxHi(myThid)
         DO j=1-Oly,sNy+Oly
          DO i=1-Olx,sNx+Olx
           IF ( WindFile .NE. ' '  ) THEN
             WIND(i,j,bi,bj)    = bWght*wspeed0(i,j,bi,bj)
     &                        +aWght*wspeed1(i,j,bi,bj)
           ENDIF
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
           IF ( AtmospFile .NE. ' '  ) THEN
             AtmosP(i,j,bi,bj)    = bWght*atmosp0(i,j,bi,bj)
     &                        +aWght*atmosp1(i,j,bi,bj)
           ENDIF
           IF ( SilicaFile .NE. ' '  ) THEN
             SILICA(i,j,bi,bj)    = bWght*silica0(i,j,bi,bj)
     &                        +aWght*silica1(i,j,bi,bj)
           ENDIF
#ifdef ALLOW_THSICE
           FIce(i,j,bi,bj)    = iceMask(i,j,bi,bj)
#else
           IF ( IceFile .NE. ' '  ) THEN
            FIce(i,j,bi,bj)    = bWght*ice0(i,j,bi,bj)
     &                          +aWght*ice1(i,j,bi,bj)
           ENDIF
           if (FIce(i,j,bi,bj).lt.1d-2) then
              FIce(i,j,bi,bj) = 0.d0
           endif
#endif
#ifdef READ_PAR
           IF ( Filename1 .NE. ' '  ) THEN
            PAR(i,j,bi,bj)    = bWght*par0(i,j,bi,bj)
     &                          +aWght*par1(i,j,bi,bj)
           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 ( IronFile .NE. ' '  ) THEN
            InputFe(i,j,bi,bj)   = bWght*feinput0(i,j,bi,bj)
     &                          +aWght*feinput1(i,j,bi,bj)
           ENDIF
#endif
          ENDDO
         ENDDO
        ENDDO
       ENDDO

C endif for periodicForcing
      ENDIF

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