/[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.4 - (hide annotations) (download)
Wed May 11 18:43:22 2011 UTC (14 years, 7 months ago) by stephd
Branch: MAIN
CVS Tags: ctrb_darwin2_ckpt64k_20130723, ctrb_darwin2_ckpt63l_20120405, ctrb_darwin2_ckpt64h_20130528, ctrb_darwin2_ckpt64m_20130820, ctrb_darwin2_ckpt64f_20130405, ctrb_darwin2_ckpt63f_20111201, ctrb_darwin2_ckpt64a_20121116, ctrb_darwin2_ckpt64n_20130826, ctrb_darwin2_ckpt62y_20110526, ctrb_darwin2_ckpt64i_20130622, ctrb_darwin2_ckpt62x_20110513, ctrb_darwin2_ckpt63o_20120629, ctrb_darwin2_ckpt64e_20130305, ctrb_darwin2_ckpt63c_20111011, ctrb_darwin2_ckpt63i_20120124, ctrb_darwin2_ckpt63m_20120506, ctrb_darwin2_ckpt63s_20120908, ctrb_darwin2_ckpt63e_20111107, ctrb_darwin2_ckpt63b_20110830, ctrb_darwin2_ckpt63j_20120217, ctrb_darwin2_ckpt63r_20120817, ctrb_darwin2_ckpt64g_20130503, ctrb_darwin2_ckpt64l_20130806, ctrb_darwin2_ckpt63g_20111220, ctrb_darwin2_ckpt64c_20130120, ctrb_darwin2_ckpt63a_20110804, ctrb_darwin2_ckpt64j_20130704, ctrb_darwin2_ckpt63h_20111230, ctrb_darwin2_ckpt63p_20120707, ctrb_darwin2_ckpt63d_20111107, ctrb_darwin2_ckpt63q_20120731, ctrb_darwin2_ckpt63_20110728, ctrb_darwin2_ckpt64b_20121224, ctrb_darwin2_ckpt64d_20130219, ctrb_darwin2_ckpt64_20121012, ctrb_darwin2_ckpt63n_20120604, ctrb_darwin2_ckpt63k_20120317, ctrb_darwin2_ckpt62z_20110622
Changes since 1.3: +0 -1 lines
o fix bug: extra line in subroutine call. BAD!

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     pH(i,j,bi,bj) = 8. _d 0
177     ENDDO
178     ENDDO
179     ENDDO
180     ENDDO
181    
182     DO bj=myByLo(myThid),myByHi(myThid)
183     DO bi=myBxLo(myThid),myBxHi(myThid)
184     DO j=1-OLy,sNy+OLy
185     DO i=1-OLx,sNx+OLx
186     ak0(i,j,bi,bj)=0. _d 0
187     ak1(i,j,bi,bj)=0. _d 0
188     ak2(i,j,bi,bj)=0. _d 0
189     akw(i,j,bi,bj)=0. _d 0
190     akb(i,j,bi,bj)=0. _d 0
191     akf(i,j,bi,bj)=0. _d 0
192     ak1p(i,j,bi,bj)=0. _d 0
193     ak2p(i,j,bi,bj)=0. _d 0
194     ak3p(i,j,bi,bj)=0. _d 0
195     aksi(i,j,bi,bj)=0. _d 0
196 stephd 1.2 fugf(i,j,bi,bj)=0. _d 0
197 jahn 1.1 ff(i,j,bi,bj)=0. _d 0
198     ft(i,j,bi,bj)=0. _d 0
199     st(i,j,bi,bj)=0. _d 0
200     bt(i,j,bi,bj)=0. _d 0
201     ENDDO
202     ENDDO
203     ENDDO
204     ENDDO
205    
206     pH_isLoaded = .FALSE.
207     IF ( nIter0.GT.PTRACERS_Iter0 ) THEN
208     C Read pH from a pickup file if needed
209     CALL DIC_READ_PICKUP(
210     O pH_isLoaded,
211     I nIter0, myThid )
212     ENDIF
213    
214     DO bj=myByLo(myThid),myByHi(myThid)
215     DO bi=myBxLo(myThid),myBxHi(myThid)
216    
217     C determine inorganic carbon chem coefficients
218     DO j=jMin,jMax
219     DO i=iMin,iMax
220    
221     c use surface layer values and convert to mol/m3
222 stephd 1.3 c and put bounds on tracers so pH solver doesn't blow up
223     surfdic(i,j) =
224     & max(10. _d 0 ,
225     & min(4000. _d 0, Ptracer(i,j,1,bi,bj,iDIC)))*1e-3
226 jahn 1.1 & * maskC(i,j,kLev,bi,bj)
227 stephd 1.3 surfalk(i,j) =
228     & max(10. _d 0 ,
229     & min(4000. _d 0, Ptracer(i,j,1,bi,bj,iALK)))*1e-3
230 jahn 1.1 & * maskC(i,j,kLev,bi,bj)
231 stephd 1.3 surfphos(i,j) =
232     & max(1. _d -10,
233     & min(10. _d 0, Ptracer(i,j,1,bi,bj,iPO4)))*1e-3
234 jahn 1.1 & * maskC(i,j,kLev,bi,bj)
235 stephd 1.3 surfsi(i,j) =
236     & max(1. _d -8,
237     & min(500. _d 0, Ptracer(i,j,1,bi,bj,iSi)))*1e-3
238 jahn 1.1 & * maskC(i,j,kLev,bi,bj)
239 stephd 1.3 surfsalt(i,j) =
240     & max(4. _d 0, min(50. _d 0, salt(i,j,kLev,bi,bj)))
241     surftemp(i,j) =
242     & max(-4. _d 0, min(39. _d 0, theta(i,j,kLev,bi,bj)))
243 jahn 1.1 c
244     WIND(i,j,bi,bj) = 5. _d 0*maskC(i,j,1,bi,bj)
245     AtmosP(i,j,bi,bj) = 1. _d 0*maskC(i,j,1,bi,bj)
246     ENDDO
247     ENDDO
248    
249     CALL CARBON_COEFFS(
250 stephd 1.3 I surftemp,surfsalt,
251 jahn 1.1 I bi,bj,iMin,iMax,jMin,jMax,myThid)
252    
253     C====================================================================
254    
255     IF ( .NOT.pH_isLoaded ) THEN
256     C set guess of pH for first step here
257    
258     print*,'QQ: pCO2 approximation method'
259     c first approximation
260     C$TAF LOOP = parallel
261     DO j=jMin,jMax
262     C$TAF LOOP = parallel
263     DO i=iMin,iMax
264     IF ( maskC(i,j,kLev,bi,bj) .NE. 0. _d 0) THEN
265     C$TAF init dic_surf = static, 10
266     DO it=1,10
267     C$TAF STORE pH(i,j,bi,bj), PTR_CO2(i,j,kLev) = dic_surf
268     C$TAF STORE surfalk(i,j), surfphos(i,j), surfsi(i,j) = dic_surf
269     CALL CALC_PCO2_APPROX(
270 stephd 1.3 I surftemp(i,j),surfsalt(i,j),
271 jahn 1.1 I surfdic(i,j), surfphos(i,j),
272     I surfsi(i,j),surfalk(i,j),
273     I ak1(i,j,bi,bj),ak2(i,j,bi,bj),
274     I ak1p(i,j,bi,bj),ak2p(i,j,bi,bj),ak3p(i,j,bi,bj),
275     I aks(i,j,bi,bj),akb(i,j,bi,bj),akw(i,j,bi,bj),
276 stephd 1.2 I aksi(i,j,bi,bj),akf(i,j,bi,bj),
277     I ak0(i,j,bi,bj), fugf(i,j,bi,bj),
278     I ff(i,j,bi,bj),
279 jahn 1.1 I bt(i,j,bi,bj),st(i,j,bi,bj),ft(i,j,bi,bj),
280     U pH(i,j,bi,bj),pCO2(i,j,bi,bj),
281     I myThid )
282     ENDDO
283     ENDIF
284     ENDDO
285     ENDDO
286     iprt = MIN(20,sNx)
287     jprt = MIN(20,sNy)
288     print*,'QQ first guess pH', pH(iprt,jprt,bi,bj),
289     & theta(iprt,jprt,1,bi,bj), salt(iprt,jprt,1,bi,bj),
290     & surfdic(iprt,jprt), surfphos(iprt,jprt),
291     & surfsi(iprt,jprt),surfalk(iprt,jprt)
292    
293     ENDIF
294    
295     C end bi,bj loops
296     ENDDO
297     ENDDO
298    
299     RETURN
300     END
301     #endif /*ALLOW_CARBON*/
302    
303     #endif /*DARWIN*/
304     #endif /*ALLOW_PTRACERS*/
305     c ==================================================================

  ViewVC Help
Powered by ViewVC 1.1.22