/[MITgcm]/MITgcm/pkg/dic/dic_surfforcing.F
ViewVC logotype

Annotation of /MITgcm/pkg/dic/dic_surfforcing.F

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


Revision 1.7 - (hide annotations) (download)
Sun Jul 18 01:13:50 2004 UTC (19 years, 11 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint54d_post, checkpoint54c_post
Changes since 1.6: +4 -4 lines
replace surfaceTendencyPtr by surfaceForcingPtr

1 jmc 1.7 C $Header: /u/gcmpack/MITgcm/pkg/dic/dic_surfforcing.F,v 1.6 2004/07/13 18:03:31 jmc Exp $
2 jmc 1.6 C $Name: $
3    
4 edhill 1.4 #include "DIC_OPTIONS.h"
5 stephd 1.1 #include "PTRACERS_OPTIONS.h"
6     #include "GCHEM_OPTIONS.h"
7    
8 stephd 1.5 CBOP
9     C !ROUTINE: DIC_SURFFORCING
10    
11     C !INTERFACE: ==========================================================
12 stephd 1.1 SUBROUTINE DIC_SURFFORCING( PTR_CO2 , GDC,
13     I bi,bj,imin,imax,jmin,jmax,
14     I myIter,myTime,myThid)
15    
16 stephd 1.5 C !DESCRIPTION:
17     C Calculate the carbon air-sea flux terms
18     C following external_forcing_dic.F (OCMIP run) from Mick
19    
20     C !USES: ===============================================================
21 stephd 1.1 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 "DIC_ABIOTIC.h"
29     #ifdef DIC_BIOTIC
30 jmc 1.6 #include "PTRACERS_SIZE.h"
31 stephd 1.1 #include "PTRACERS.h"
32     #endif
33    
34 stephd 1.5 C !INPUT PARAMETERS: ===================================================
35     C myThid :: thread number
36     C myIter :: current timestep
37     C myTime :: current time
38     c PTR_CO2 :: DIC tracer field
39 stephd 1.1 INTEGER myIter, myThid
40     _RL myTime
41     _RL PTR_CO2(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
42 stephd 1.5 INTEGER iMin,iMax,jMin,jMax, bi, bj
43    
44     C !OUTPUT PARAMETERS: ===================================================
45 jmc 1.7 c GDC :: surface flux due to air-sea exchange
46 stephd 1.1 _RL GDC(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
47    
48     #ifdef ALLOW_PTRACERS
49 stephd 1.5
50     C !LOCAL VARIABLES: ====================================================
51 stephd 1.2 INTEGER I,J, kLev, it
52 stephd 1.1 C Number of iterations for pCO2 solvers...
53     C Solubility relation coefficients
54     _RL SchmidtNoDIC(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
55     _RL pCO2sat(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
56     _RL Kwexch(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
57     C local variables for carbon chem
58     _RL surfalk(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
59     _RL surfphos(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
60     _RL surfsi(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
61     _RL VirtualFlux(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
62 stephd 1.5 CEOP
63 stephd 1.1
64     cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
65    
66     kLev=1
67    
68     C PRE-INDUSTRIAL STEADY STATE pCO2 = 278.0 ppmv
69     DO j=1-OLy,sNy+OLy
70     DO i=1-OLx,sNx+OLx
71     AtmospCO2(i,j,bi,bj)=278.0d-6
72     ENDDO
73     ENDDO
74    
75    
76     C =================================================================
77     C determine inorganic carbon chem coefficients
78     DO j=1-OLy,sNy+OLy
79     DO i=1-OLx,sNx+OLx
80    
81     #ifdef DIC_BIOTIC
82     cQQQQ check ptracer numbers
83     surfalk(i,j) = PTRACER(i,j,klev,bi,bj,2)
84     & * maskC(i,j,kLev,bi,bj)
85     surfphos(i,j) = PTRACER(i,j,klev,bi,bj,3)
86     & * maskC(i,j,kLev,bi,bj)
87     #else
88     surfalk(i,j) = 2.366595 * salt(i,j,kLev,bi,bj)/gsm_s
89     & * maskC(i,j,kLev,bi,bj)
90     surfphos(i,j) = 5.1225e-4 * maskC(i,j,kLev,bi,bj)
91     #endif
92     C FOR NON-INTERACTIVE Si
93 stephd 1.3 surfsi(i,j) = SILICA(i,j,bi,bj) * maskC(i,j,kLev,bi,bj)
94 stephd 1.1 ENDDO
95     ENDDO
96    
97     CALL CARBON_COEFFS(
98     I theta,salt,
99     I bi,bj,iMin,iMax,jMin,jMax)
100     C====================================================================
101    
102     c pCO2 solver...
103 stephd 1.3 C$TAF LOOP = parallel
104 stephd 1.1 DO j=1-OLy,sNy+OLy
105 stephd 1.3 C$TAF LOOP = parallel
106 stephd 1.1 DO i=1-OLx,sNx+OLx
107    
108     IF(maskC(i,j,kLev,bi,bj) .NE. 0.)THEN
109     CALL CALC_PCO2_APPROX(
110     I theta(i,j,kLev,bi,bj),salt(i,j,kLev,bi,bj),
111     I PTR_CO2(i,j,kLev), surfphos(i,j),
112     I surfsi(i,j),surfalk(i,j),
113     I ak1(i,j,bi,bj),ak2(i,j,bi,bj),
114     I ak1p(i,j,bi,bj),ak2p(i,j,bi,bj),ak3p(i,j,bi,bj),
115     I aks(i,j,bi,bj),akb(i,j,bi,bj),akw(i,j,bi,bj),
116     I aksi(i,j,bi,bj),akf(i,j,bi,bj),ff(i,j,bi,bj),
117     I bt(i,j,bi,bj),st(i,j,bi,bj),ft(i,j,bi,bj),
118     U pH(i,j,bi,bj),pCO2(i,j,bi,bj) )
119     ELSE
120     pCO2(i,j,bi,bj)=0. _d 0
121     END IF
122     ENDDO
123     ENDDO
124    
125     DO j=1-OLy,sNy+OLy
126     DO i=1-OLx,sNx+OLx
127    
128     IF (maskC(i,j,kLev,bi,bj).NE.0.) THEN
129     C calculate SCHMIDT NO. for CO2
130     SchmidtNoDIC(i,j) =
131     & sca1
132     & + sca2 * theta(i,j,kLev,bi,bj)
133     & + sca3 * theta(i,j,kLev,bi,bj)*theta(i,j,kLev,bi,bj)
134     & + sca4 * theta(i,j,kLev,bi,bj)*theta(i,j,kLev,bi,bj)
135     & *theta(i,j,kLev,bi,bj)
136    
137     C Determine surface flux (FDIC)
138     C first correct pCO2at for surface atmos pressure
139     pCO2sat(i,j) =
140     & AtmosP(i,j,bi,bj)*AtmospCO2(i,j,bi,bj)
141     c find exchange coefficient
142     c account for schmidt number and and varible piston velocity
143     Kwexch(i,j) =
144     & pisvel(i,j,bi,bj)
145     & / sqrt(SchmidtNoDIC(i,j)/660.0)
146     c OR use a constant coeff
147     c Kwexch(i,j) = 5e-5
148     c ice influence
149     cQQ Kwexch(i,j) =(1.d0-Fice(i,j,bi,bj))*Kwexch(i,j)
150    
151    
152     C Calculate flux in terms of DIC units using K0, solubility
153     C Flux = Vp * ([CO2sat] - [CO2])
154     C CO2sat = K0*pCO2atmos*P/P0
155     C Converting pCO2 to [CO2] using ff, as in CALC_PCO2
156 stephd 1.2 FluxCO2(i,j,bi,bj) =
157 stephd 1.1 & maskC(i,j,kLev,bi,bj)*Kwexch(i,j)*(
158     & ak0(i,j,bi,bj)*pCO2sat(i,j) -
159     & ff(i,j,bi,bj)*pCO2(i,j,bi,bj)
160     & )
161     ELSE
162 stephd 1.2 FluxCO2(i,j,bi,bj) = 0.
163 stephd 1.1 ENDIF
164     C convert flux (mol kg-1 m s-1) to (mol m-2 s-1)
165 stephd 1.2 FluxCO2(i,j,bi,bj) = FluxCO2(i,j,bi,bj)/permil
166 stephd 1.1
167     IF (maskC(i,j,kLev,bi,bj).NE.0.) THEN
168     c calculate virtual flux
169     c EminusPforV = dS/dt*(1/Sglob)
170     C NOTE: Be very careful with signs here!
171     C Positive EminusPforV => loss of water to atmos and increase
172     C in salinity. Thus, also increase in other surface tracers
173     C (i.e. positive virtual flux into surface layer)
174     C ...so here, VirtualFLux = dC/dt!
175 jmc 1.7 VirtualFlux(i,j)=gsm_DIC*surfaceForcingS(i,j,bi,bj)/gsm_s
176 stephd 1.1 c OR
177     c let virtual flux be zero
178     c VirtualFlux(i,j)=0.d0
179     c
180     ELSE
181     VirtualFlux(i,j)=0. _d 0
182     ENDIF
183     ENDDO
184     ENDDO
185    
186     C update tendency
187     DO j=1-OLy,sNy+OLy
188     DO i=1-OLx,sNx+OLx
189     GDC(i,j)= maskC(i,j,kLev,bi,bj)*(
190 jmc 1.7 & FluxCO2(i,j,bi,bj)
191 stephd 1.1 & + VirtualFlux(i,j)
192     & )
193     ENDDO
194     ENDDO
195    
196     #endif
197     RETURN
198     END

  ViewVC Help
Powered by ViewVC 1.1.22