*IF DEF,OCEAN ORH1F305.465
C ******************************COPYRIGHT****************************** GTS2F400.7417
C (c) CROWN COPYRIGHT 1995, METEOROLOGICAL OFFICE, All Rights Reserved. GTS2F400.7418
C GTS2F400.7419
C Use, duplication or disclosure of this code is subject to the GTS2F400.7420
C restrictions as set forth in the contract. GTS2F400.7421
C GTS2F400.7422
C Meteorological Office GTS2F400.7423
C London Road GTS2F400.7424
C BRACKNELL GTS2F400.7425
C Berkshire UK GTS2F400.7426
C RG12 2SZ GTS2F400.7427
C GTS2F400.7428
C If no contract has been raised with this copy of the code, the use, GTS2F400.7429
C duplication or disclosure of it is strictly prohibited. Permission GTS2F400.7430
C to do so must first be obtained in writing from the Head of Numerical GTS2F400.7431
C Modelling at the above address. GTS2F400.7432
C ******************************COPYRIGHT****************************** GTS2F400.7433
C GTS2F400.7434
C*LL Subroutine PPCO2 PPCO2.3
CLL Can run on any FORTRAN 77 compiler with long lower case variables PPCO2.4
CLL PPCO2.5
CLL Author: N.K. TAYLOR PPCO2.6
CLL Date: 17 December 1993 PPCO2.7
CLL Version 3.3 PPCO2.8
CLL PPCO2.9
! Modification History: ORH1F305.4181
! Version Date Details ORH1F305.4182
! ------- ------- ------------------------------------------ ORH1F305.4183
! 3.5 16.01.95 Remove *IF dependency. R.Hill ORH1F305.4184
CLL Programming standards use Cox naming convention for Cox variables PPCO2.10
CLL with the addition that lower case variables are local to the PPCO2.11
CLL routine. PPCO2.12
CLL PPCO2.13
CLL This routine computes partial pressure of CO2 in the surface mixed PPCO2.14
CLL layer using the routine of Bacastow (1981) given in PPCO2.15
CLL SCOPE 16: "Carbon Cycle Modelling" PPCO2.16
CLL PPCO2.17
CLL Note that the total dissolved inorganic carbon (DIC, or TCO2) is PPCO2.18
CLL held as tracer #3 in units of micro-Moles per litre. PPCO2.19
CLL PPCO2.20
CLLEND ----------------------------------------------------------------- PPCO2.21
C* PPCO2.22
C*L-------------------------------- Arguments ------------------------ PPCO2.23
C PPCO2.24
SUBROUTINE PPCO2 (T, REK0, REK, RG1, RG2, RG3, 1OJP0F404.639
+ IMT, KM, NT, PPCO2.26
+ FM, PCO2) PPCO2.27
C PPCO2.28
IMPLICIT NONE PPCO2.29
C PPCO2.30
*CALL CNTLOCN 
ORH1F305.4185
*CALL OARRYSIZ ORH1F305.4186
*CALL OTRACPNT OJP0F404.640
*CALL UMSCALAR OJP0F404.641
! Call umscalar to get use of RHO_WATER_SI=1026.0 kg/m3 (set in OSETCON) OJP0F404.642
C Define constants for array sizes PPCO2.31
C PPCO2.32
INTEGER PPCO2.33
+ IMT ! IN Number of points in horizontal PPCO2.34
+, KM ! IN Number of layers in model PPCO2.35
+, NT ! IN Number of tracers PPCO2.36
C PPCO2.37
C Physical arguments PPCO2.38
C PPCO2.39
REAL PPCO2.40
+ T (IMT, KM, NT) ! IN Tracer values OJP0F404.643
+, REK0 (IMT) ! IN Solubility of CO2 in Sea Water PPCO2.42
+, RG1 (IMT) ! IN Products etc of the Equilibrium PPCO2.43
+, RG2 (IMT) ! IN constants, which are used for PPCO2.44
+, RG3 (IMT) ! IN the calculation of PCO2 PPCO2.45
+, REK (IMT) ! IN PPCO2.46
+, FM (IMT,KM) ! IN Land-sea mask PPCO2.47
+, PCO2 (IMT) ! OUT Partial pressure of CO2 in ppm PPCO2.48
! OJP0F404.644
! Units of dissocitaion constant products are OJP0F404.645
! REK unitless OJP0F404.646
! RG1 kg/mole RG2,RG3 moles/kg OJP0F404.647
! REK0 moles/(kg.atm) OJP0F404.648
! The solution of the carbonate system is done in units of moles/kg OJP0F404.649
! whereas the tracers TCO2 and ALK are held in micromoles per litre. OJP0F404.650
! The errors in pco2 due to converting moles/litre to moles/kg using OJP0F404.651
! a constant of 1/1.024 are only rder 1.024 , as long as BOTH ALK and OJP0F404.652
! TCO2 are subject to the same conversion error. OJP0F404.653
C* PPCO2.49
C PPCO2.50
C PPCO2.51
C Locally defined variables PPCO2.52
C PPCO2.53
INTEGER iter ! Number of iterations used PPCO2.54
+, i ! Horizontal loop index PPCO2.55
+, maxit ! Maximum number of iterations PPCO2.56
&, ISSW ! Start column index catering for cyclic conds ORH1F305.4187
&, IESW ! Start column index catering for cyclic conds ORH1F305.4188
C PPCO2.57
REAL xs,x2p,x1,x2,x3,xx,x ! Intermediates for computing CO2 p PPCO2.58
REAL xw,alkprime,qa,qb,qc ! PPCO2.59
REAL fxa ! PPCO2.60
REAL sal ! Surface salinity in ppt PPCO2.61
REAL alk_t ! Surface Alkalinity (moles/kg) OJP0F404.654
REAL tco2 ! Surface TCO2 (moles/kg) OJP0F404.655
REAL boron ! Surface Boron conc. (moles/kg) OJP0F404.656
C ! (proportional to salinity) PPCO2.64
C PPCO2.65
DATA maxit /20/ PPCO2.66
C PPCO2.67
! Set up start and end indexes controls dependent on whether cyclic ORH1F305.4189
! conditions are in use. ORH1F305.4190
IF (L_OCYCLIC) THEN ORH1F305.4191
ISSW = 2 ORH1F305.4192
IESW = IMT-1 ORH1F305.4193
ELSE ORH1F305.4194
ISSW = 1 ORH1F305.4195
IESW = IMT ORH1F305.4196
ENDIF ORH1F305.4197
CL Compute partial pressure along a row, using an iterative procedure PPCO2.68
CL at each grid point along the row. PPCO2.69
C PPCO2.70
DO I = ISSW, IESW ORH1F305.4198
C PPCO2.77
C Convert surface salinity from model units to parts per thousand PPCO2.78
C PPCO2.79
sal = T(I,1,2)*1000.0 + 35.0 OJP0F404.657
C PPCO2.81
! Convert TCO2 from micromoles per litre to moles per kg OJP0F404.658
! OJP0F404.659
tco2 = T(I,1,TCO2_TRACER) /(1.0e6 * (RHO_WATER_SI/1000.0)) OJP0F404.660
OJP0F404.661
! Use one of these formulations for surface Alkalinity OJP0F404.662
! 1) Use the modelled tracer value (needs BIOLOGY) converting OJP0F404.663
! from micromoles per litre to moles per kg OJP0F404.664
alk_t = T(I,1,ALK_TRACER) /(1.0e6 * (RHO_WATER_SI/1000.0)) OJP0F404.665
! 2) Set surface Alkalinity to a constant value appropriate to a globa OJP0F404.666
C average salinity of 35 ppt. PPCO2.90
! alk_t = 2300.0e-6 OJP0F404.667
! 3) Set surface Alkalinity to salinity*2311/34.78 (GEOSECS) OJP0F404.668
! Put limits on to cope with exceptionally high and low salinities OJP0F404.669
! alk_t = sal * (2374.0e-6/34.78) OJP0F404.670
! alk_t = amax1(alk_t,2000.0e-6) OJP0F404.671
! alk_t = amin1(alk_t,2500.0e-6) OJP0F404.672
C PPCO2.93
C Set Boron concentration according to the surface salinity (Boron PPCO2.94
C accounts for about 10% of the total Alkalinity) PPCO2.95
! Use relationship used by Peng et al, Tellus 39b 1987, originally OJP0F404.673
! from Culkin 1965: OJP0F404.674
C PPCO2.96
boron = 4.106e-4 * sal / 35.0 OJP0F404.675
OJP0F404.676
if(tco2 .le. 0.0) then OJP0F404.677
pco2(i)=0.0 PPCO2.99
goto 250 PPCO2.100
endif PPCO2.101
C PPCO2.102
xs = 1.0 PPCO2.103
x2p = 0.0 PPCO2.104
x3 = 0.0 PPCO2.105
xx = 0.0 PPCO2.106
x = xs PPCO2.107
iter = 0 PPCO2.108
80 x1 = x2p PPCO2.109
iter = iter + 1 PPCO2.110
x2 = x3 PPCO2.111
x2p = xx PPCO2.112
x3 = x PPCO2.113
IF (iter-2) 90,90,100 PPCO2.114
90 xx = x PPCO2.115
GOTO 110 PPCO2.116
100 xx = x2p+(x3-x2p)*(x1-x2p)/(x1-x2-x2p+x3) PPCO2.117
110 CONTINUE PPCO2.118
xw = boron/(1.0+RG1(I)/xx)+RG2(I)*xx-RG3(I)/xx PPCO2.119
alkprime = alk_t - xw PPCO2.120
C PPCO2.121
c x is given as the root of a quadratic equation PPCO2.122
c PPCO2.123
qa = 2.0*tco2 - alkprime OJP0F404.678
qb = -REK(I)*(alkprime-tco2) OJP0F404.679
qc = -alkprime PPCO2.126
x = (-qb+sqrt(qb*qb-4.0*qa*qc))/(2.0*qa) PPCO2.127
c PPCO2.128
c Finish iteration if sufficient accuracy has been achieved PPCO2.129
c PPCO2.130
IF (((abs(x-xx)/x).GT.1.0E-05) .and. (iter.LT.maxit)) GOTO 80 PPCO2.131
C PPCO2.132
C Finally get partial pressure of CO2 in water and mask out values PPCO2.133
! over land points. Multiply partial pressure by 1e6 to get it in OJP0F404.680
! microatmospheres. OJP0F404.681
c PPCO2.135
fxa = 1.0+REK(I)*x + x*x OJP0F404.682
if(fxa.gt.0.0) then PPCO2.137
PCO2(I) = FM(I,1) * 1.0e6 * tco2/fxa/REK0(I) OJP0F404.683
else PPCO2.139
PCO2(I) = 0.0 PPCO2.140
endif PPCO2.141
250 CONTINUE PPCO2.142
PPCO2.143
ENDDO PPCO2.144
C PPCO2.145
IF (L_OCYCLIC) THEN ORH1F305.4202
C Set cyclic boundary condition on PCO2 PPCO2.147
PCO2(1) = PCO2(IMT-1) PPCO2.148
PCO2(IMT) = PCO2(2) PPCO2.149
ENDIF ORH1F305.4203
C PPCO2.151
CL Return from PPCO2 PPCO2.152
C PPCO2.153
RETURN PPCO2.154
END PPCO2.155
*ENDIF PPCO2.156