/[MITgcm]/MITgcm_contrib/ecco_darwin/v4_3deg/code/dic_budgetApCO2.F
ViewVC logotype

Annotation of /MITgcm_contrib/ecco_darwin/v4_3deg/code/dic_budgetApCO2.F

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


Revision 1.2 - (hide annotations) (download)
Mon Jan 27 06:31:50 2020 UTC (6 years, 7 months ago) by dimitri
Branch: MAIN
CVS Tags: HEAD
Changes since 1.1: +0 -0 lines
FILE REMOVED
where possible, using code and input files from va_llc270

1 dimitri 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_BUDGETAPCO2
12    
13     C !INTERFACE: ==========================================================
14     SUBROUTINE DIC_BUDGETAPCO2( PTR_DIC , PTR_ALK, PTR_PO4, PTR_SIL,
15     O deltaApCO2,
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     #ifdef ALLOW_OLD_VIRTUALFLUX
67     _RL VirtualFlux(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
68     #endif
69     C local variables for CO2_FLUX_BUDGET
70     _RL deltaApCO2(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
71     CEOP
72    
73     cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
74    
75     kLev=1
76    
77     cc if coupled to atmsopheric model, use the
78     cc Co2 value passed from the coupler
79     c#ifndef USE_ATMOSCO2
80     cC PRE-INDUSTRIAL STEADY STATE pCO2 = 278.0 ppmv
81     c DO j=1-OLy,sNy+OLy
82     c DO i=1-OLx,sNx+OLx
83     c AtmospCO2(i,j,bi,bj)=278.0 _d -6
84     c ENDDO
85     c ENDDO
86     c#endif
87     C =================================================================
88     C determine inorganic carbon chem coefficients
89     DO j=jmin,jmax
90     DO i=imin,imax
91     c put bounds on tracers so pH solver doesn't blow up
92     surfdic(i,j) =
93     & max(100. _d 0 , min(4000. _d 0, PTR_DIC(i,j)))*1e-3
94     & * maskC(i,j,kLev,bi,bj)
95     surfalk(i,j) =
96     & max(100. _d 0 , min(4000. _d 0, PTR_ALK(i,j)))*1e-3
97     & * maskC(i,j,kLev,bi,bj)
98     surfphos(i,j) =
99     & max(1. _d -10, min(10. _d 0, PTR_PO4(i,j)))*1e-3
100     & * maskC(i,j,kLev,bi,bj)
101     surfsi(i,j) =
102     & max(1. _d -8, min(500. _d 0, PTR_SIL(i,j)))*1e-3
103     & * maskC(i,j,kLev,bi,bj)
104     surfsalt(i,j) =
105     & max(4. _d 0, min(50. _d 0, salt(i,j,kLev,bi,bj)))
106     surftemp(i,j) =
107     & max(-4. _d 0, min(39. _d 0, theta(i,j,kLev,bi,bj)))
108     ENDDO
109     ENDDO
110    
111     CALL CARBON_COEFFS(
112     I surftemp,surfsalt,
113     I bi,bj,iMin,iMax,jMin,jMax,myThid)
114     C====================================================================
115    
116     DO j=jmin,jmax
117     DO i=imin,imax
118     C Compute AtmosP and Kwexch_Pre which are re-used for flux of O2
119    
120     #ifdef USE_PLOAD
121     C Convert anomalous pressure pLoad (in Pa) from atmospheric model
122     C to total pressure (in Atm)
123     C Note: it is assumed the reference atmospheric pressure is 1Atm=1013mb
124     C rather than the actual ref. pressure from Atm. model so that on
125     C average AtmosP is about 1 Atm.
126     AtmosP(i,j,bi,bj)= 1. _d 0 + pLoad(i,j,bi,bj)/Pa2Atm
127     #endif
128    
129     C Pre-compute part of exchange coefficient: pisvel*(1-fice)
130     C Schmidt number is accounted for later
131     #ifdef USE_EXFWIND
132     pisvel(i,j)=0.337 _d 0 *wspeed(i,j,bi,bj)**2/3.6 _d 5
133     cBX linear piston velocity after Krakauer et al. (2006), Eq. 3
134     cBX using <k> = 20, n=0.5, and <u^n> = 2.6747 (as determined from 2010
135     cBX EXFwspee field from cube92 run)
136     cDc pisvel(i,j)=20 _d 0 *(wspeed(i,j,bi,bj)**0.5
137     cDc & /2.6747 _d 0) /3.6 _d 5
138     #else
139     pisvel(i,j)=0.337 _d 0 *wind(i,j,bi,bj)**2/3.6 _d 5
140     #endif
141     Kwexch_Pre(i,j,bi,bj) = pisvel(i,j)
142     & * (1. _d 0 - FIce(i,j,bi,bj))
143    
144     ENDDO
145     ENDDO
146    
147     c pCO2 solver...
148     C$TAF LOOP = parallel
149     DO j=jmin,jmax
150     C$TAF LOOP = parallel
151     DO i=imin,imax
152    
153     IF ( maskC(i,j,kLev,bi,bj).NE.0. _d 0 ) THEN
154     CALL CALC_PCO2_APPROX(
155     I surftemp(i,j),surfsalt(i,j),
156     I surfdic(i,j), surfphos(i,j),
157     I surfsi(i,j),surfalk(i,j),
158     I ak1(i,j,bi,bj),ak2(i,j,bi,bj),
159     I ak1p(i,j,bi,bj),ak2p(i,j,bi,bj),ak3p(i,j,bi,bj),
160     I aks(i,j,bi,bj),akb(i,j,bi,bj),akw(i,j,bi,bj),
161     I aksi(i,j,bi,bj),akf(i,j,bi,bj),
162     I ak0(i,j,bi,bj), fugf(i,j,bi,bj),
163     I ff(i,j,bi,bj),
164     I bt(i,j,bi,bj),st(i,j,bi,bj),ft(i,j,bi,bj),
165     U pH(i,j,bi,bj),pCO2(i,j,bi,bj),CO3(i,j,bi,bj),
166     I myThid )
167     ELSE
168     pH(i,j,bi,bj) = 0. _d 0
169     pCO2(i,j,bi,bj) = 0. _d 0
170     CO3(i,j,bi,bj) = 0. _d 0
171     ENDIF
172     ENDDO
173     ENDDO
174    
175     DO j=jmin,jmax
176     DO i=imin,imax
177     IF ( maskC(i,j,kLev,bi,bj).NE.0. _d 0 ) THEN
178     pCO2sat(i,j) = AtmosP(i,j,bi,bj)*AtmospCO2(i,j,bi,bj)
179     ENDIF
180     ENDDO
181     ENDDO
182    
183     C update tendency
184     DO j=jmin,jmax
185     DO i=imin,imax
186     if(budgetTStep1.EQ.0) then
187     C if first timestep
188     C this is problematic at restart; clean-up later
189     deltaApCO2(i,j) = 0. _d 0
190     previousApCO2(i,j,bi,bj) = pCO2sat(i,j)
191     else
192     C current value - value from last timestep
193     deltaApCO2(i,j) = pCO2sat(i,j) - previousApCO2(i,j,bi,bj)
194     previousApCO2(i,j,bi,bj) = pCO2sat(i,j)
195     endif
196     ENDDO
197     ENDDO
198    
199     RETURN
200     END
201     #endif /*ALLOW_CARBON*/
202    
203     #endif /*DARWIN*/
204     #endif /*ALLOW_PTRACERS*/
205     c ==================================================================

  ViewVC Help
Powered by ViewVC 1.1.22