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

Annotation of /MITgcm_contrib/dcarroll/highres_darwin/code/dic_budgetSalt.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_BUDGETSALT
12    
13     C !INTERFACE: ==========================================================
14     SUBROUTINE DIC_BUDGETSALT( PTR_DIC , PTR_ALK, PTR_PO4, PTR_SIL,
15     O deltaSalt,
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 budgetSalt1Pert(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 deltaSalt(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     C salt from previous timestep
107     surfsalt(i,j) = budgetSalt1(i,j,bi,bj)
108     C salt from current timestep
109     budgetSalt1(i,j,bi,bj) =
110     & max(4. _d 0, min(50. _d 0, salt(i,j,kLev,bi,bj)))
111     if(budgetTStep1.EQ.0) then
112     C if first timestep
113     C this is problematic for restarts; to do correctly we will have to
114     C add to pickups or run simulation without interruptions
115     surfsalt(i,j) = budgetSalt1(i,j,bi,bj)
116     endif
117     C perturbation is from current value
118     budgetSalt1Pert(i,j) = budgetSalt1(i,j,bi,bj) +
119     & budgetPert
120     surftemp(i,j) =
121     & max(-4. _d 0, min(39. _d 0, theta(i,j,kLev,bi,bj)))
122     ENDDO
123     ENDDO
124    
125     CALL CARBON_COEFFS(
126     I surftemp,budgetSalt1Pert,
127     I bi,bj,iMin,iMax,jMin,jMax,myThid)
128     C====================================================================
129    
130     DO j=jmin,jmax
131     DO i=imin,imax
132     C Compute AtmosP and Kwexch_Pre which are re-used for flux of O2
133    
134     #ifdef USE_PLOAD
135     C Convert anomalous pressure pLoad (in Pa) from atmospheric model
136     C to total pressure (in Atm)
137     C Note: it is assumed the reference atmospheric pressure is 1Atm=1013mb
138     C rather than the actual ref. pressure from Atm. model so that on
139     C average AtmosP is about 1 Atm.
140     AtmosP(i,j,bi,bj)= 1. _d 0 + pLoad(i,j,bi,bj)/Pa2Atm
141     #endif
142    
143     C Pre-compute part of exchange coefficient: pisvel*(1-fice)
144     C Schmidt number is accounted for later
145     #ifdef USE_EXFWIND
146     pisvel(i,j)=0.337 _d 0 *wspeed(i,j,bi,bj)**2/3.6 _d 5
147     cBX linear piston velocity after Krakauer et al. (2006), Eq. 3
148     cBX using <k> = 20, n=0.5, and <u^n> = 2.6747 (as determined from 2010
149     cBX EXFwspee field from cube92 run)
150     cDc pisvel(i,j)=20 _d 0 *(wspeed(i,j,bi,bj)**0.5
151     cDc & /2.6747 _d 0) /3.6 _d 5
152     #else
153     pisvel(i,j)=0.337 _d 0 *wind(i,j,bi,bj)**2/3.6 _d 5
154     #endif
155     Kwexch_Pre(i,j,bi,bj) = pisvel(i,j)
156     & * (1. _d 0 - FIce(i,j,bi,bj))
157    
158     ENDDO
159     ENDDO
160    
161     c pCO2 solver...
162     C$TAF LOOP = parallel
163     DO j=jmin,jmax
164     C$TAF LOOP = parallel
165     DO i=imin,imax
166    
167     IF ( maskC(i,j,kLev,bi,bj).NE.0. _d 0 ) THEN
168     CALL CALC_PCO2_APPROX(
169     I surftemp(i,j),budgetSalt1Pert(i,j),
170     I surfdic(i,j), surfphos(i,j),
171     I surfsi(i,j),surfalk(i,j),
172     I ak1(i,j,bi,bj),ak2(i,j,bi,bj),
173     I ak1p(i,j,bi,bj),ak2p(i,j,bi,bj),ak3p(i,j,bi,bj),
174     I aks(i,j,bi,bj),akb(i,j,bi,bj),akw(i,j,bi,bj),
175     I aksi(i,j,bi,bj),akf(i,j,bi,bj),
176     I ak0(i,j,bi,bj), fugf(i,j,bi,bj),
177     I ff(i,j,bi,bj),
178     I bt(i,j,bi,bj),st(i,j,bi,bj),ft(i,j,bi,bj),
179     U pH(i,j,bi,bj),pCO2(i,j,bi,bj),CO3(i,j,bi,bj),
180     I myThid )
181     ELSE
182     pH(i,j,bi,bj) = 0. _d 0
183     pCO2(i,j,bi,bj) = 0. _d 0
184     CO3(i,j,bi,bj) = 0. _d 0
185     ENDIF
186     ENDDO
187     ENDDO
188    
189    
190     DO j=jmin,jmax
191     DO i=imin,imax
192    
193     IF ( maskC(i,j,kLev,bi,bj).NE.0. _d 0 ) THEN
194     C calculate SCHMIDT NO. for CO2
195     SchmidtNoDIC(i,j) =
196     & sca1
197     & + sca2 * surftemp(i,j)
198     & + sca3 * surftemp(i,j)*surftemp(i,j)
199     & + sca4 * surftemp(i,j)*surftemp(i,j)
200     & *surftemp(i,j)
201     c put positive bound on SCHMIT number (will go negative for temp>40)
202     SchmidtNoDIC(i,j) = max(1. _d -2, SchmidtNoDIC(i,j))
203    
204     C Determine surface flux (FDIC)
205     C first correct pCO2at for surface atmos pressure
206     pCO2sat(i,j) =
207     & AtmosP(i,j,bi,bj)*AtmospCO2(i,j,bi,bj)
208    
209     C then account for Schmidt number
210     Kwexch(i,j) = Kwexch_Pre(i,j,bi,bj)
211     & / sqrt(SchmidtNoDIC(i,j)/660.0 _d 0)
212    
213     #ifdef WATERVAP_BUG
214     C Calculate flux in terms of DIC units using K0, solubility
215     C Flux = Vp * ([CO2sat] - [CO2])
216     C CO2sat = K0*pCO2atmos*P/P0
217     C Converting pCO2 to [CO2] using ff, as in CALC_PCO2
218     FluxCO2_loc(i,j) =
219     & Kwexch(i,j)*(
220     & ak0(i,j,bi,bj)*pCO2sat(i,j) -
221     & ff(i,j,bi,bj)*pCO2(i,j,bi,bj)
222     & )
223     #else
224     C Corrected by Val Bennington Nov 2010 per G.A. McKinley's finding
225     C of error in application of water vapor correction
226     c Flux = kw*rho*(ff*pCO2atm-k0*FugFac*pCO2ocean)
227     FluxCO2_loc(i,j) =
228     & Kwexch(i,j)*(
229     & ff(i,j,bi,bj)*pCO2sat(i,j) -
230     & pCO2(i,j,bi,bj)*fugf(i,j,bi,bj)
231     & *ak0(i,j,bi,bj) )
232     &
233     #endif
234     ELSE
235     FluxCO2_loc(i,j) = 0. _d 0
236     ENDIF
237     C convert flux (mol kg-1 m s-1) to (mol m-2 s-1)
238     FluxCO2_loc(i,j) = FluxCO2_loc(i,j)/permil
239     c convert flux (mol m-2 s-1) to (mmol m-2 s-1)
240     FluxCO2_loc(i,j) = FluxCO2_loc(i,j)*1. _d 3
241    
242     #ifdef ALLOW_OLD_VIRTUALFLUX
243     IF (maskC(i,j,kLev,bi,bj).NE.0. _d 0) THEN
244     c calculate virtual flux
245     c EminusPforV = dS/dt*(1/Sglob)
246     C NOTE: Be very careful with signs here!
247     C Positive EminusPforV => loss of water to atmos and increase
248     C in salinity. Thus, also increase in other surface tracers
249     C (i.e. positive virtual flux into surface layer)
250     C ...so here, VirtualFLux = dC/dt!
251     VirtualFlux(i,j)=gsm_DIC*surfaceForcingS(i,j,bi,bj)/gsm_s
252     c OR
253     c let virtual flux be zero
254     c VirtualFlux(i,j)=0.d0
255     c
256     ELSE
257     VirtualFlux(i,j)=0. _d 0
258     ENDIF
259     #endif /* ALLOW_OLD_VIRTUALFLUX */
260     ENDDO
261     ENDDO
262    
263     C update tendency
264     DO j=jmin,jmax
265     DO i=imin,imax
266     if(budgetTStep1.EQ.0) then
267     C if first timestep
268     C this is problematic at restart; clean-up later
269     dFluxCO2Salt(i,j,bi,bj) = 0. _d 0
270     deltaSalt(i,j) = 0. _d 0
271     else
272     C at this point in code, fluxCO2_1 contains
273     C total flux for current time step
274     dFluxCO2Salt(i,j,bi,bj) = fluxCO2_1(i,j,bi,bj) -
275     & FluxCO2_loc(i,j)
276     C current value - value from previous timestep
277     deltaSalt(i,j) = budgetSalt1(i,j,bi,bj) -
278     & surfsalt(i,j)
279     endif
280     ENDDO
281     ENDDO
282    
283     RETURN
284     END
285     #endif /*ALLOW_CARBON*/
286    
287     #endif /*DARWIN*/
288     #endif /*ALLOW_PTRACERS*/
289     c ==================================================================

  ViewVC Help
Powered by ViewVC 1.1.22