/[MITgcm]/MITgcm_contrib/darwin2/pkg/darwin/dic_surfforcing_init.F
ViewVC logotype

Annotation of /MITgcm_contrib/darwin2/pkg/darwin/dic_surfforcing_init.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.6 - (hide annotations) (download)
Wed Oct 9 17:14:37 2013 UTC (12 years, 2 months ago) by stephd
Branch: MAIN
CVS Tags: ctrb_darwin2_ckpt64o_20131024, ctrb_darwin2_ckpt64p_20131118, ctrb_darwin2_ckpt64q_20131118, ctrb_darwin2_ckpt64p_20131024
Changes since 1.5: +29 -5 lines
o changes so that pH and pCO2 can be calculated for full water column
  by defining pH_3D in DARWIN_OPTIONS.h. Includes 3D diags and pickup

1 jahn 1.1 #include "CPP_OPTIONS.h"
2     #include "PTRACERS_OPTIONS.h"
3     #include "DARWIN_OPTIONS.h"
4    
5     #ifdef ALLOW_PTRACERS
6     #ifdef ALLOW_DARWIN
7    
8     #ifdef ALLOW_CARBON
9    
10     CBOP
11     C !ROUTINE: DIC_SURFFORCING_INIT
12    
13     C !INTERFACE: ==========================================================
14     SUBROUTINE DIC_SURFFORCING_INIT(
15     I myThid)
16    
17     C !DESCRIPTION:
18     C Calculate first guess of pH
19    
20     C !USES: ===============================================================
21     IMPLICIT NONE
22     #include "SIZE.h"
23     #include "DYNVARS.h"
24     #include "EEPARAMS.h"
25     #include "PARAMS.h"
26     #include "GRID.h"
27     #include "FFIELDS.h"
28     #include "PTRACERS_SIZE.h"
29     #include "PTRACERS_PARAMS.h"
30     #include "PTRACERS_FIELDS.h"
31     #include "DARWIN_SIZE.h"
32     #include "DARWIN_IO.h"
33     #include "DARWIN_FLUX.h"
34     #include "DIC_ATMOS.h"
35    
36     C !INPUT PARAMETERS: ===================================================
37     C myThid :: thread number
38     INTEGER myThid
39    
40    
41     C !LOCAL VARIABLES: ====================================================
42     INTEGER i,j, k, kLev, it
43     INTEGER intime0,intime1
44     _RL otime
45     _RL aWght,bWght,rdt
46     INTEGER nForcingPeriods,Imytm,Ifprd,Ifcyc,Iftm
47     C Number of iterations for pCO2 solvers...
48     C Solubility relation coefficients
49     C local variables for carbon chem
50     INTEGER iMin,iMax,jMin,jMax, bi, bj
51     _RL surfdic(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
52     _RL surfalk(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
53     _RL surfphos(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
54     _RL surfsi(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
55     _RL surfsalt(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
56 stephd 1.3 _RL surftemp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
57 jahn 1.1 INTEGER iprt,jprt
58     CHARACTER*(MAX_LEN_MBUF) msgBuf
59     INTEGER iUnit
60     CHARACTER*(MAX_LEN_FNAM) iName
61     integer ilo,ihi
62     integer ilnblnk,ifnblnk
63     external ilnblnk,ifnblnk
64     LOGICAL pH_isLoaded
65     CEOP
66    
67     cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
68    
69     kLev=1
70    
71     CALL DIC_ATMOS(0, startTime, nIter0, myThid )
72    
73     _BEGIN_MASTER(myThid)
74    
75     C set up coefficients for DIC chemistry
76     C define Schmidt no. coefficients for CO2
77     sca1 = 2073.1 _d 0
78     sca2 = -125.62 _d 0
79     sca3 = 3.6276 _d 0
80     sca4 = -0.043219 _d 0
81     C define Schmidt no. coefficients for O2
82     C based on Keeling et al [GBC, 12, 141, (1998)]
83     sox1 = 1638.0 _d 0
84     sox2 = -81.83 _d 0
85     sox3 = 1.483 _d 0
86     sox4 = -0.008004 _d 0
87    
88     C coefficients for determining saturation O2
89     oA0= 2.00907 _d 0
90     oA1= 3.22014 _d 0
91     oA2= 4.05010 _d 0
92     oA3= 4.94457 _d 0
93     oA4= -2.56847 _d -1
94     oA5= 3.88767 _d 0
95     oB0= -6.24523 _d -3
96     oB1= -7.37614 _d -3
97     oB2= -1.03410 _d -2
98     oB3= -8.17083 _d -3
99     oC0= -4.88682 _d -7
100     C Set other constant/flag
101    
102     #ifndef USE_ATMOSCO2
103    
104 stephd 1.3 #ifndef USE_EXFCO2
105 jahn 1.1 if (dic_int1.eq.2) then
106     call mdsfindunit( iUnit, mythid )
107     open(UNIT=iUnit,FILE='co2atmos.dat',STATUS='old')
108     do k=1,dic_int2
109     read(iUnit,*) co2atmos(k)
110     print*,'co2atmos',co2atmos(k)
111     enddo
112     close(iUnit)
113     endif
114    
115     if (dic_int1.eq.3) then
116     write(iName,'(A,I10.10)') 'dic_atmos.',nIter0
117     ilo = ifnblnk(iName)
118     ihi = ilnblnk(iName)
119     call mdsfindunit( iUnit, mythid )
120     open(UNIT=iUnit,FILE=iname(ilo:ihi),STATUS='old')
121     read(iUnit,*) total_atmos_carbon_ini,
122     & atpco2_ini
123     close(iUnit)
124     endif
125 stephd 1.3 #endif
126 jahn 1.1
127     #endif
128     _END_MASTER(myThid)
129    
130     ccccccccccccccccccccccccccccccccccccccccc
131     IF ( periodicExternalForcing ) THEN
132    
133    
134     rdt = 1. _d 0 / deltaTclock
135     nForcingPeriods = NINT(externForcingCycle/externForcingPeriod)
136     cswd QQ change for placement of chem forcing (ie. after timestep)
137     Imytm = NINT(startTime*rdt)
138     Ifprd = NINT(externForcingPeriod*rdt)
139     Ifcyc = NINT(externForcingCycle*rdt)
140     Iftm = MOD( Imytm+Ifcyc-Ifprd/2, Ifcyc)
141    
142     intime0 = 1 + INT(Iftm/Ifprd)
143     intime1 = 1 + MOD(intime0,nForcingPeriods)
144     c aWght = DFLOAT( Iftm-Ifprd*(intime0 - 1) ) / DFLOAT( Ifprd )
145     aWght = FLOAT( Iftm-Ifprd*(intime0 - 1) )
146     bWght = FLOAT( Ifprd )
147     aWght = aWght / bWght
148     bWght = 1. _d 0 - aWght
149    
150     _BARRIER
151     _BEGIN_MASTER(myThid)
152    
153     _END_MASTER(myThid)
154    
155     #ifdef ALLOW_OFFLINE
156     IF ( useOffLine ) THEN
157     otime=nIter0*deltaTclock
158     CALL OFFLINE_FIELDS_LOAD( otime, nIter0, myThid )
159     ENDIF
160     #endif
161    
162     c end periodicExternalForcing
163     ENDIF
164    
165     C =================================================================
166    
167     jMin=1
168     jMax=sNy
169     iMin=1
170     iMax=sNx
171    
172     DO bj=myByLo(myThid),myByHi(myThid)
173     DO bi=myBxLo(myThid),myBxHi(myThid)
174     DO j=1-OLy,sNy+OLy
175     DO i=1-Olx,sNx+OLx
176 stephd 1.6 #ifdef pH_3D
177     DO k=1,Nr
178     pH(i,j,k,bi,bj) = 8. _d 0
179     ENDDO
180     #else
181 jahn 1.1 pH(i,j,bi,bj) = 8. _d 0
182 stephd 1.6 #endif
183 jahn 1.1 ENDDO
184     ENDDO
185     ENDDO
186     ENDDO
187    
188     DO bj=myByLo(myThid),myByHi(myThid)
189     DO bi=myBxLo(myThid),myBxHi(myThid)
190     DO j=1-OLy,sNy+OLy
191     DO i=1-OLx,sNx+OLx
192     ak0(i,j,bi,bj)=0. _d 0
193     ak1(i,j,bi,bj)=0. _d 0
194     ak2(i,j,bi,bj)=0. _d 0
195     akw(i,j,bi,bj)=0. _d 0
196     akb(i,j,bi,bj)=0. _d 0
197     akf(i,j,bi,bj)=0. _d 0
198     ak1p(i,j,bi,bj)=0. _d 0
199     ak2p(i,j,bi,bj)=0. _d 0
200     ak3p(i,j,bi,bj)=0. _d 0
201     aksi(i,j,bi,bj)=0. _d 0
202 stephd 1.2 fugf(i,j,bi,bj)=0. _d 0
203 jahn 1.1 ff(i,j,bi,bj)=0. _d 0
204     ft(i,j,bi,bj)=0. _d 0
205     st(i,j,bi,bj)=0. _d 0
206     bt(i,j,bi,bj)=0. _d 0
207     ENDDO
208     ENDDO
209     ENDDO
210     ENDDO
211    
212     pH_isLoaded = .FALSE.
213     IF ( nIter0.GT.PTRACERS_Iter0 ) THEN
214     C Read pH from a pickup file if needed
215     CALL DIC_READ_PICKUP(
216     O pH_isLoaded,
217     I nIter0, myThid )
218     ENDIF
219    
220     DO bj=myByLo(myThid),myByHi(myThid)
221     DO bi=myBxLo(myThid),myBxHi(myThid)
222    
223 stephd 1.6 #ifdef pH_3D
224     DO kLev=1,Nr
225     #endif
226 jahn 1.1 C determine inorganic carbon chem coefficients
227     DO j=jMin,jMax
228     DO i=iMin,iMax
229    
230     c use surface layer values and convert to mol/m3
231 stephd 1.3 c and put bounds on tracers so pH solver doesn't blow up
232     surfdic(i,j) =
233     & max(10. _d 0 ,
234 stephd 1.6 & min(4000. _d 0, Ptracer(i,j,kLev,bi,bj,iDIC)))*1e-3
235 jahn 1.1 & * maskC(i,j,kLev,bi,bj)
236 stephd 1.3 surfalk(i,j) =
237     & max(10. _d 0 ,
238 stephd 1.6 & min(4000. _d 0, Ptracer(i,j,kLev,bi,bj,iALK)))*1e-3
239 jahn 1.1 & * maskC(i,j,kLev,bi,bj)
240 stephd 1.3 surfphos(i,j) =
241     & max(1. _d -10,
242 stephd 1.6 & min(10. _d 0, Ptracer(i,j,kLev,bi,bj,iPO4)))*1e-3
243 jahn 1.1 & * maskC(i,j,kLev,bi,bj)
244 stephd 1.3 surfsi(i,j) =
245     & max(1. _d -8,
246 stephd 1.6 & min(500. _d 0, Ptracer(i,j,kLev,bi,bj,iSi)))*1e-3
247 jahn 1.1 & * maskC(i,j,kLev,bi,bj)
248 stephd 1.3 surfsalt(i,j) =
249     & max(4. _d 0, min(50. _d 0, salt(i,j,kLev,bi,bj)))
250     surftemp(i,j) =
251     & max(-4. _d 0, min(39. _d 0, theta(i,j,kLev,bi,bj)))
252 jahn 1.1 c
253     WIND(i,j,bi,bj) = 5. _d 0*maskC(i,j,1,bi,bj)
254     AtmosP(i,j,bi,bj) = 1. _d 0*maskC(i,j,1,bi,bj)
255     ENDDO
256     ENDDO
257    
258     CALL CARBON_COEFFS(
259 stephd 1.3 I surftemp,surfsalt,
260 stephd 1.5 I bi,bj,iMin,iMax,jMin,jMax,
261     I kLev,myThid)
262 jahn 1.1
263     C====================================================================
264    
265     IF ( .NOT.pH_isLoaded ) THEN
266     C set guess of pH for first step here
267    
268     print*,'QQ: pCO2 approximation method'
269     c first approximation
270     C$TAF LOOP = parallel
271     DO j=jMin,jMax
272     C$TAF LOOP = parallel
273     DO i=iMin,iMax
274     IF ( maskC(i,j,kLev,bi,bj) .NE. 0. _d 0) THEN
275     C$TAF init dic_surf = static, 10
276     DO it=1,10
277     C$TAF STORE pH(i,j,bi,bj), PTR_CO2(i,j,kLev) = dic_surf
278     C$TAF STORE surfalk(i,j), surfphos(i,j), surfsi(i,j) = dic_surf
279     CALL CALC_PCO2_APPROX(
280 stephd 1.3 I surftemp(i,j),surfsalt(i,j),
281 jahn 1.1 I surfdic(i,j), surfphos(i,j),
282     I surfsi(i,j),surfalk(i,j),
283     I ak1(i,j,bi,bj),ak2(i,j,bi,bj),
284     I ak1p(i,j,bi,bj),ak2p(i,j,bi,bj),ak3p(i,j,bi,bj),
285     I aks(i,j,bi,bj),akb(i,j,bi,bj),akw(i,j,bi,bj),
286 stephd 1.2 I aksi(i,j,bi,bj),akf(i,j,bi,bj),
287     I ak0(i,j,bi,bj), fugf(i,j,bi,bj),
288     I ff(i,j,bi,bj),
289 jahn 1.1 I bt(i,j,bi,bj),st(i,j,bi,bj),ft(i,j,bi,bj),
290 stephd 1.6 #ifdef pH_3D
291     U pH(i,j,kLev,bi,bj),pCO2(i,j,kLev,bi,bj),
292     #else
293 jahn 1.1 U pH(i,j,bi,bj),pCO2(i,j,bi,bj),
294 stephd 1.6 #endif
295 jahn 1.1 I myThid )
296     ENDDO
297     ENDIF
298     ENDDO
299     ENDDO
300     iprt = MIN(20,sNx)
301     jprt = MIN(20,sNy)
302 stephd 1.6 #ifdef pH_3D
303     print*,'QQ first guess pH',kLev, pH(iprt,jprt,kLev,bi,bj),
304     & theta(iprt,jprt,kLev,bi,bj), salt(iprt,jprt,kLev,bi,bj),
305     & surfdic(iprt,jprt), surfphos(iprt,jprt),
306     & surfsi(iprt,jprt),surfalk(iprt,jprt)
307     #else
308 jahn 1.1 print*,'QQ first guess pH', pH(iprt,jprt,bi,bj),
309     & theta(iprt,jprt,1,bi,bj), salt(iprt,jprt,1,bi,bj),
310     & surfdic(iprt,jprt), surfphos(iprt,jprt),
311     & surfsi(iprt,jprt),surfalk(iprt,jprt)
312 stephd 1.6 #endif
313     ENDIF
314    
315     #ifdef pH_3D
316     ENDDO
317     #endif
318 jahn 1.1
319    
320     C end bi,bj loops
321     ENDDO
322     ENDDO
323    
324     RETURN
325     END
326     #endif /*ALLOW_CARBON*/
327    
328     #endif /*DARWIN*/
329     #endif /*ALLOW_PTRACERS*/
330     c ==================================================================

  ViewVC Help
Powered by ViewVC 1.1.22