/[MITgcm]/MITgcm_contrib/dcarroll/highres_darwin/code/dic_budgetTheta.F
ViewVC logotype

Annotation of /MITgcm_contrib/dcarroll/highres_darwin/code/dic_budgetTheta.F

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


Revision 1.1 - (hide annotations) (download)
Sun Sep 22 21:23:46 2019 UTC (5 years, 10 months ago) by dcarroll
Branch: MAIN
CVS Tags: HEAD
Initial check in of high resolution Darwin simulation code

1 dcarroll 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_BUDGETTHETA
12    
13     C !INTERFACE: ==========================================================
14     SUBROUTINE DIC_BUDGETTHETA( PTR_DIC , PTR_ALK, PTR_PO4, PTR_SIL,
15     O deltaTheta,
16     I bi,bj,imin,imax,jmin,jmax,
17     I myIter,myTime,myThid)
18    
19     C !DESCRIPTION:
20     C Calculate the carbon air-sea flux terms
21     C following external_forcing_dic.F (OCMIP run) from Mick
22    
23     C !USES: ===============================================================
24     IMPLICIT NONE
25     #include "SIZE.h"
26     #include "DYNVARS.h"
27     #include "EEPARAMS.h"
28     #include "PARAMS.h"
29     #include "GRID.h"
30     #include "FFIELDS.h"
31     #include "DARWIN_SIZE.h"
32     #include "DARWIN_IO.h"
33     #include "DARWIN_FLUX.h"
34     #ifdef USE_EXFWIND
35     #include "EXF_FIELDS.h"
36     #endif
37    
38     C !INPUT PARAMETERS: ===================================================
39     C myThid :: thread number
40     C myIter :: current timestep
41     C myTime :: current time
42     c PTR_DIC :: DIC tracer field
43     INTEGER myIter, myThid
44     _RL myTime
45     _RL PTR_DIC(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
46     _RL PTR_ALK(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
47     _RL PTR_PO4(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
48     _RL PTR_SIL(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
49     INTEGER iMin,iMax,jMin,jMax, bi, bj
50    
51     C !LOCAL VARIABLES: ====================================================
52     INTEGER I,J, kLev, it
53     C Number of iterations for pCO2 solvers...
54     C Solubility relation coefficients
55     _RL SchmidtNoDIC(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
56     _RL pCO2sat(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
57     _RL Kwexch(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
58     _RL pisvel(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
59     C local variables for carbon chem
60     _RL surfdic(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
61     _RL surfalk(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
62     _RL surfphos(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
63     _RL surfsi(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
64     _RL surfsalt(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
65     _RL surftemp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
66     _RL budgetTemp1Pert(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
67     #ifdef ALLOW_OLD_VIRTUALFLUX
68     _RL VirtualFlux(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
69     #endif
70     C local variables for CO2_FLUX_BUDGET
71     _RL FluxCO2_loc(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
72     _RL deltaTheta(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
73     CEOP
74    
75     cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
76    
77     kLev=1
78    
79     cc if coupled to atmsopheric model, use the
80     cc Co2 value passed from the coupler
81     c#ifndef USE_ATMOSCO2
82     cC PRE-INDUSTRIAL STEADY STATE pCO2 = 278.0 ppmv
83     c DO j=1-OLy,sNy+OLy
84     c DO i=1-OLx,sNx+OLx
85     c AtmospCO2(i,j,bi,bj)=278.0 _d -6
86     c ENDDO
87     c ENDDO
88     c#endif
89     C =================================================================
90     C determine inorganic carbon chem coefficients
91     DO j=jmin,jmax
92     DO i=imin,imax
93     c put bounds on tracers so pH solver doesn't blow up
94     surfdic(i,j) =
95     & max(100. _d 0 , min(4000. _d 0, PTR_DIC(i,j)))*1e-3
96     & * maskC(i,j,kLev,bi,bj)
97     surfalk(i,j) =
98     & max(100. _d 0 , min(4000. _d 0, PTR_ALK(i,j)))*1e-3
99     & * maskC(i,j,kLev,bi,bj)
100     surfphos(i,j) =
101     & max(1. _d -10, min(10. _d 0, PTR_PO4(i,j)))*1e-3
102     & * maskC(i,j,kLev,bi,bj)
103     surfsi(i,j) =
104     & max(1. _d -8, min(500. _d 0, PTR_SIL(i,j)))*1e-3
105     & * maskC(i,j,kLev,bi,bj)
106     surfsalt(i,j) =
107     & max(4. _d 0, min(50. _d 0, salt(i,j,kLev,bi,bj)))
108     C theta from previous timestep
109     surftemp(i,j) = budgetTemp1(i,j,bi,bj)
110     C theta from current timestep
111     budgetTemp1(i,j,bi,bj) =
112     & max(-4. _d 0, min(39. _d 0, theta(i,j,kLev,bi,bj)))
113     if(budgetTStep1.EQ.0) then
114     C if first timestep
115     C this is problematic for restarts; to do correctly we will have to
116     C add to pickups or run simulation without interruptions
117     surftemp(i,j) = budgetTemp1(i,j,bi,bj)
118     endif
119     budgetTemp1Pert(i,j) = budgetTemp1(i,j,bi,bj) +
120     & budgetPert
121     ENDDO
122     ENDDO
123    
124     CALL CARBON_COEFFS(
125     I budgetTemp1Pert,surfsalt,
126     I bi,bj,iMin,iMax,jMin,jMax,myThid)
127     C====================================================================
128    
129     DO j=jmin,jmax
130     DO i=imin,imax
131     C Compute AtmosP and Kwexch_Pre which are re-used for flux of O2
132    
133     #ifdef USE_PLOAD
134     C Convert anomalous pressure pLoad (in Pa) from atmospheric model
135     C to total pressure (in Atm)
136     C Note: it is assumed the reference atmospheric pressure is 1Atm=1013mb
137     C rather than the actual ref. pressure from Atm. model so that on
138     C average AtmosP is about 1 Atm.
139     AtmosP(i,j,bi,bj)= 1. _d 0 + pLoad(i,j,bi,bj)/Pa2Atm
140     #endif
141    
142     C Pre-compute part of exchange coefficient: pisvel*(1-fice)
143     C Schmidt number is accounted for later
144     #ifdef USE_EXFWIND
145     pisvel(i,j)=0.337 _d 0 *wspeed(i,j,bi,bj)**2/3.6 _d 5
146     cBX linear piston velocity after Krakauer et al. (2006), Eq. 3
147     cBX using <k> = 20, n=0.5, and <u^n> = 2.6747 (as determined from 2010
148     cBX EXFwspee field from cube92 run)
149     cDc pisvel(i,j)=20 _d 0 *(wspeed(i,j,bi,bj)**0.5
150     cDc & /2.6747 _d 0) /3.6 _d 5
151     #else
152     pisvel(i,j)=0.337 _d 0 *wind(i,j,bi,bj)**2/3.6 _d 5
153     #endif
154     Kwexch_Pre(i,j,bi,bj) = pisvel(i,j)
155     & * (1. _d 0 - FIce(i,j,bi,bj))
156    
157     ENDDO
158     ENDDO
159    
160     c pCO2 solver...
161     C$TAF LOOP = parallel
162     DO j=jmin,jmax
163     C$TAF LOOP = parallel
164     DO i=imin,imax
165    
166     IF ( maskC(i,j,kLev,bi,bj).NE.0. _d 0 ) THEN
167     CALL CALC_PCO2_APPROX(
168     I budgetTemp1Pert(i,j),surfsalt(i,j),
169     I surfdic(i,j), surfphos(i,j),
170     I surfsi(i,j),surfalk(i,j),
171     I ak1(i,j,bi,bj),ak2(i,j,bi,bj),
172     I ak1p(i,j,bi,bj),ak2p(i,j,bi,bj),ak3p(i,j,bi,bj),
173     I aks(i,j,bi,bj),akb(i,j,bi,bj),akw(i,j,bi,bj),
174     I aksi(i,j,bi,bj),akf(i,j,bi,bj),
175     I ak0(i,j,bi,bj), fugf(i,j,bi,bj),
176     I ff(i,j,bi,bj),
177     I bt(i,j,bi,bj),st(i,j,bi,bj),ft(i,j,bi,bj),
178     U pH(i,j,bi,bj),pCO2(i,j,bi,bj),CO3(i,j,bi,bj),
179     I myThid )
180     ELSE
181     pH(i,j,bi,bj) = 0. _d 0
182     pCO2(i,j,bi,bj) = 0. _d 0
183     CO3(i,j,bi,bj) = 0. _d 0
184     ENDIF
185     ENDDO
186     ENDDO
187    
188     DO j=jmin,jmax
189     DO i=imin,imax
190    
191     IF ( maskC(i,j,kLev,bi,bj).NE.0. _d 0 ) THEN
192     C calculate SCHMIDT NO. for CO2
193     SchmidtNoDIC(i,j) =
194     & sca1
195     & + sca2 * surftemp(i,j)
196     & + sca3 * surftemp(i,j)*surftemp(i,j)
197     & + sca4 * surftemp(i,j)*surftemp(i,j)
198     & *surftemp(i,j)
199     c put positive bound on SCHMIT number (will go negative for temp>40)
200     SchmidtNoDIC(i,j) = max(1. _d -2, SchmidtNoDIC(i,j))
201    
202     C Determine surface flux (FDIC)
203     C first correct pCO2at for surface atmos pressure
204     pCO2sat(i,j) =
205     & AtmosP(i,j,bi,bj)*AtmospCO2(i,j,bi,bj)
206    
207     C then account for Schmidt number
208     Kwexch(i,j) = Kwexch_Pre(i,j,bi,bj)
209     & / sqrt(SchmidtNoDIC(i,j)/660.0 _d 0)
210    
211     #ifdef WATERVAP_BUG
212     C Calculate flux in terms of DIC units using K0, solubility
213     C Flux = Vp * ([CO2sat] - [CO2])
214     C CO2sat = K0*pCO2atmos*P/P0
215     C Converting pCO2 to [CO2] using ff, as in CALC_PCO2
216     FluxCO2_loc(i,j) =
217     & Kwexch(i,j)*(
218     & ak0(i,j,bi,bj)*pCO2sat(i,j) -
219     & ff(i,j,bi,bj)*pCO2(i,j,bi,bj)
220     & )
221     #else
222     C Corrected by Val Bennington Nov 2010 per G.A. McKinley's finding
223     C of error in application of water vapor correction
224     c Flux = kw*rho*(ff*pCO2atm-k0*FugFac*pCO2ocean)
225     FluxCO2_loc(i,j) =
226     & Kwexch(i,j)*(
227     & ff(i,j,bi,bj)*pCO2sat(i,j) -
228     & pCO2(i,j,bi,bj)*fugf(i,j,bi,bj)
229     & *ak0(i,j,bi,bj) )
230     &
231     #endif
232     ELSE
233     FluxCO2_loc(i,j) = 0. _d 0
234     ENDIF
235     C convert flux (mol kg-1 m s-1) to (mol m-2 s-1)
236     FluxCO2_loc(i,j) = FluxCO2_loc(i,j)/permil
237     c convert flux (mol m-2 s-1) to (mmol m-2 s-1)
238     FluxCO2_loc(i,j) = FluxCO2_loc(i,j)*1. _d 3
239    
240     #ifdef ALLOW_OLD_VIRTUALFLUX
241     IF (maskC(i,j,kLev,bi,bj).NE.0. _d 0) THEN
242     c calculate virtual flux
243     c EminusPforV = dS/dt*(1/Sglob)
244     C NOTE: Be very careful with signs here!
245     C Positive EminusPforV => loss of water to atmos and increase
246     C in salinity. Thus, also increase in other surface tracers
247     C (i.e. positive virtual flux into surface layer)
248     C ...so here, VirtualFLux = dC/dt!
249     VirtualFlux(i,j)=gsm_DIC*surfaceForcingS(i,j,bi,bj)/gsm_s
250     c OR
251     c let virtual flux be zero
252     c VirtualFlux(i,j)=0.d0
253     c
254     ELSE
255     VirtualFlux(i,j)=0. _d 0
256     ENDIF
257     #endif /* ALLOW_OLD_VIRTUALFLUX */
258     ENDDO
259     ENDDO
260    
261     C update tendency
262     DO j=jmin,jmax
263     DO i=imin,imax
264     if(budgetTStep1.EQ.0) then
265     C if first timestep
266     C this is problematic at restart; clean-up later
267     dFluxCO2Temp(i,j,bi,bj) = 0. _d 0
268     deltaTheta(i,j) = 0. _d 0
269     else
270     C at this point in code, fluxCO2_1 contains
271     C total flux for current time step
272     dFluxCO2Temp(i,j,bi,bj) = fluxCO2_1(i,j,bi,bj) -
273     & FluxCO2_loc(i,j)
274     C current value - value from last timestep
275     deltaTheta(i,j) = budgetTemp1(i,j,bi,bj) -
276     & surftemp(i,j)
277     endif
278     ENDDO
279     ENDDO
280    
281     RETURN
282     END
283     #endif /*ALLOW_CARBON*/
284    
285     #endif /*DARWIN*/
286     #endif /*ALLOW_PTRACERS*/
287     c ==================================================================

  ViewVC Help
Powered by ViewVC 1.1.22