*IF DEF,OCEAN ORH1F305.345
C ******************************COPYRIGHT****************************** GTS2F400.3007
C (c) CROWN COPYRIGHT 1995, METEOROLOGICAL OFFICE, All Rights Reserved. GTS2F400.3008
C GTS2F400.3009
C Use, duplication or disclosure of this code is subject to the GTS2F400.3010
C restrictions as set forth in the contract. GTS2F400.3011
C GTS2F400.3012
C Meteorological Office GTS2F400.3013
C London Road GTS2F400.3014
C BRACKNELL GTS2F400.3015
C Berkshire UK GTS2F400.3016
C RG12 2SZ GTS2F400.3017
C GTS2F400.3018
C If no contract has been raised with this copy of the code, the use, GTS2F400.3019
C duplication or disclosure of it is strictly prohibited. Permission GTS2F400.3020
C to do so must first be obtained in writing from the Head of Numerical GTS2F400.3021
C Modelling at the above address. GTS2F400.3022
C ******************************COPYRIGHT****************************** GTS2F400.3023
C GTS2F400.3024
C*LL Subroutine FLUX_CO2 FLUX_CO2.3
CLL Can run on any FORTRAN 77 compiler with long lower case variables FLUX_CO2.4
CLL FLUX_CO2.5
CLL Author: N.K. TAYLOR FLUX_CO2.6
CLL Date 17 December 1993 FLUX_CO2.7
CLL Version 3.3 FLUX_CO2.8
CLL FLUX_CO2.9
CLL Programming standards use Cox naming convention for Cox variables FLUX_CO2.10
CLL with the addition that lower case variables are local to the FLUX_CO2.11
CLL routine. FLUX_CO2.12
CLL FLUX_CO2.13
CLL Computes the flux of CO2 from the atmosphere into the ocean, from FLUX_CO2.14
CLL the difference in partial pressure between atmosphere and ocean FLUX_CO2.15
CLL The flux is only applied to ice-free points FLUX_CO2.17
CLL - when L_SEAICE is true. ORH1F305.346
CLL The Liss & Merlivat (1986) formula is used for wind speed dependance FLUX_CO2.20
CLL - when L_OLISS is true. ORH1F305.347
CLL The Wanninkhof (1992) formula is used for wind speed dependance FLUX_CO2.23
CLL - when L_OLISS is false. ORH1F305.348
! Modification History: ORH1F305.349
! Version Date Details ORH1F305.350
! ------- ------- ------------------------------------------ ORH1F305.351
! 3.5 16.01.95 Remove *IF dependency. R.Hill ORH1F305.352
CLL 4.5 1.07.98 Allow full 2D field of atmospheric CO2 OCN1F405.27
CLL to be used in flux calculations. C.D.Jones OCN1F405.28
CLL FLUX_CO2.25
CLLEND ----------------------------------------------------------------- FLUX_CO2.26
C* FLUX_CO2.27
C*L-------------------------------- Arguments ------------------------ FLUX_CO2.28
C FLUX_CO2.29
SUBROUTINE FLUX_CO2(TA, REK0, PCO2, CO2_FLUX, ATMPCO2_ROW, 1OCN1F405.29
+ INVADE,EVADE, FLUX_CO2.31
+ c14to12_atm, FLUX_CO2.33
+ WME, FLUX_CO2.36
+ AICE, FLUX_CO2.39
+ BICE, FLUX_CO2.42
+ IMT, KM, NT, FLUX_CO2.44
+ FM, DZ, TIMESTEP ) OJP0F404.573
C FLUX_CO2.46
IMPLICIT NONE FLUX_CO2.47
C FLUX_CO2.48
*CALL CNTLOCN 
ORH1F305.353
*CALL OARRYSIZ ORH1F305.354
*CALL OTRACPNT OJP0F404.574
*CALL UMSCALAR OJP0F404.575
! Use UMSCALAR to get RHO_WATER_SI = 1026.0 (set in OSETCON) OJP0F404.576
C Define constants for array sizes FLUX_CO2.49
C FLUX_CO2.50
INTEGER FLUX_CO2.51
+ IMT ! IN Number of points in horizontal FLUX_CO2.52
+, KM ! IN Number of layers in model FLUX_CO2.53
+, NT ! IN Number of tracers FLUX_CO2.54
C FLUX_CO2.55
C Physical arguments FLUX_CO2.56
C FLUX_CO2.57
REAL FLUX_CO2.58
& TA (IMT, KM, NT) ! INOUT Tracer values OJP0F404.577
+, CO2_FLUX (IMT) ! OUT Air-sea CO2 flux (Moles/m2/year) FLUX_CO2.60
+, INVADE (IMT) ! OUT Invasion rate of CO2 (Moles/m2/yr) FLUX_CO2.61
+, EVADE (IMT) ! OUT Invasion rate of CO2 (Moles/m2/yr) FLUX_CO2.62
&, ATMPCO2_ROW(IMT_CAR) ! OUT Atmospheric CO2 surface conc OCN1F405.30
+, PCO2_ATM ! OUT Atmospheric CO2 pressure in ppm FLUX_CO2.63
+, REK0 (IMT) ! IN Solubility of CO2 in Sea Water FLUX_CO2.64
+, DZ (KM) ! IN Layer thicknesses (in cm) FLUX_CO2.65
&, TIMESTEP ! IN Timestep in top layer OJP0F404.578
+, FM (IMT,KM) ! IN Land-sea mask FLUX_CO2.67
+, PCO2 (IMT) ! IN Partial pressure of CO2 in ppm FLUX_CO2.68
&, WME (IMT_CAR_MIX) ! IN Wind-mixing power (W/m2) ORH1F305.355
&, AICE(IMT_CAR_ICE) ! IN Ice fraction in grid box ORH1F305.356
&, BICE(IMT_CAR_PIC) ! IN ARRAY TO BLANK OUT FLUXES AT ORH1F305.357
! CLIMATOLOGICAL ice points ORH1F305.358
REAL c14to12_atm ! IN Atmospheric C14/C12 ratio ORH1F305.359
C* FLUX_CO2.82
C FLUX_CO2.83
C FLUX_CO2.84
C Locally defined variables FLUX_CO2.85
C FLUX_CO2.86
INTEGER FLUX_CO2.87
+ i ! Horizontal loop index FLUX_CO2.88
C FLUX_CO2.89
REAL FLUX_CO2.90
+ piston_vel(IMT) ! Piston velocity (m/s) FLUX_CO2.91
+, schmidt ! Schmidt Number FLUX_CO2.95
+, wind_speed ! Wind speed computed from WME (m/s) FLUX_CO2.96
+, rho_air ! Average surface density of air FLUX_CO2.97
&, fxb ! Intermediate constants OJP0F404.579
+, drag_coef ! Drag coefficient for conversion FLUX_CO2.99
C ! of wind-mixing power to wind speed FLUX_CO2.100
C FLUX_CO2.102
DATA rho_air,drag_coef / 1.2 , 1.3 e-3/ FLUX_CO2.104
C FLUX_CO2.107
CL Define conversion factor and intermediate constants FLUX_CO2.108
C FLUX_CO2.109
IF (L_OMIXLAY) THEN ORH1F305.360
fxb = sqrt ((RHO_WATER_SI**(0.333))/(drag_coef*rho_air)) OJP0F404.580
ENDIF ORH1F305.361
C FLUX_CO2.118
C FLUX_CO2.119
DO I = IFROM_CYC, ITO_CYC ORH1F305.362
IF (L_OMIXLAY) THEN ORH1F305.363
C FLUX_CO2.127
C Compute a typical "wind speed" using wind mixing power FLUX_CO2.128
C FLUX_CO2.129
! Mask windspeed calculation to zero over land. OJP0F404.581
wind_speed = fxb*((FM(I,1)*WME(I))**(0.333)) OJP0F404.582
C FLUX_CO2.131
C Compute Schmidt Number FLUX_CO2.132
C Formula from Wanninkhof (1992) FLUX_CO2.133
C FLUX_CO2.134
! Calculate schmidt number, with mask to make schmidt=1.0 over land. OJP0F404.583
schmidt = FM(I,1)* ( 2073.1 - 125.62*TA(I,1,1) OJP0F404.584
& + 3.6276*(TA(I,1,1)**2) - 0.043219*(TA(I,1,1)**3) ) OJP0F404.585
& + ( 1.0 - FM(I,1) ) OJP0F404.586
C FLUX_CO2.137
IF (L_OLISS) THEN ORH1F305.364
C Compute piston velocity (in m/s) (From Liss & Merlivat, 1986) FLUX_CO2.139
C FLUX_CO2.140
IF (wind_speed.LE.3.6) THEN FLUX_CO2.141
piston_vel(I) = 0.17 * wind_speed * (0.01/3600.0) OJP0F404.587
& /((schmidt/660.0)**0.667) FLUX_CO2.143
ELSE FLUX_CO2.144
IF (wind_speed.LE.13.0) THEN FLUX_CO2.145
piston_vel(I) = (2.85*wind_speed -9.65) * (0.01/3600.0) OJP0F404.588
& /(sqrt(schmidt/660.0)) FLUX_CO2.147
ELSE FLUX_CO2.148
piston_vel(I) = (5.9*wind_speed -49.3) * (0.01/3600.0) OJP0F404.589
& /(sqrt(schmidt/660.0)) FLUX_CO2.150
ENDIF FLUX_CO2.151
ENDIF FLUX_CO2.152
ELSE ORH1F305.365
! Compute piston velocity (in m/s) (From Wanninkhof 1992, in cm/hr) OJP0F404.590
C FLUX_CO2.156
piston_vel(I) = 0.39 * wind_speed * wind_speed OJP0F404.591
& / (sqrt(schmidt/660.0)) FLUX_CO2.158
& * 0.01 / 3600.0 OJP0F404.592
ENDIF ORH1F305.366
ELSE ORH1F305.367
piston_vel(I) = 5.556e-5 FLUX_CO2.162
ENDIF ORH1F305.368
ENDDO FLUX_CO2.164
C FLUX_CO2.165
! COMPUTE CO2 INVASION, EVASION AND NET FLUXES in moles/(m2.yr) OJP0F404.593
CL NB Flux is blanked out over land and ice points (except for leads) FLUX_CO2.167
C FLUX_CO2.168
IF (L_SEAICE) THEN ORH1F305.369
DO I = IFROM_CYC, ITO_CYC ORH1F305.370
INVADE(I) = REK0(I)*piston_vel(I)*atmpco2_row(i) OCN1F405.31
& * (1.0 - AICE(I)) * FM(I,1) ORH1F305.372
& * RHO_WATER_SI * 60.0*60.0*24.0*360.0 / 1.0e6 OJP0F404.595
EVADE(I) = REK0(I)*piston_vel(I)*PCO2(I) OJP0F404.596
& * (1.0 - AICE(I)) * FM(I,1) ORH1F305.374
& * RHO_WATER_SI * 60.0*60.0*24.0*360.0 / 1.0e6 OJP0F404.597
CO2_FLUX(I) = INVADE(I)-EVADE(I) ORH1F305.375
ENDDO ORH1F305.376
ELSE ORH1F305.377
IF (L_OPSEUDIC) THEN ORH1F305.378
DO I = IFROM_CYC, ITO_CYC ORH1F305.379
INVADE(I) = REK0(I)*piston_vel(I)*atmpco2_row(i) OCN1F405.32
& * BICE(I) * FM(I,1) ORH1F305.381
& * RHO_WATER_SI * 60.0*60.0*24.0*360.0 / 1.0e6 OJP0F404.599
EVADE(I) = REK0(I)*piston_vel(I)*PCO2(I) OJP0F404.600
& * BICE(I) * FM(I,1) ORH1F305.383
& * RHO_WATER_SI * 60.0*60.0*24.0*360.0 / 1.0e6 OJP0F404.601
CO2_FLUX(I) = INVADE(I)-EVADE(I) ORH1F305.384
ENDDO ORH1F305.385
ELSE ORH1F305.386
DO I = IFROM_CYC, ITO_CYC ORH1F305.387
INVADE(I) = REK0(I)*piston_vel(I)*atmpco2_row(i) OCN1F405.33
& * FM(I,1) OJP0F404.603
& * RHO_WATER_SI * 60.0*60.0*24.0*360.0 / 1.0e6 OJP0F404.604
EVADE(I) = REK0(I)*piston_vel(I)*PCO2(I) OJP0F404.605
& * FM(I,1) OJP0F404.606
& * RHO_WATER_SI * 60.0*60.0*24.0*360.0 / 1.0e6 OJP0F404.607
CO2_FLUX(I) = INVADE(I)-EVADE(I) ORH1F305.390
ENDDO ORH1F305.391
ENDIF ORH1F305.392
ORH1F305.393
ENDIF ORH1F305.394
C FLUX_CO2.197
C Now update TCO2 (tracer 3) according to the net CO2 flux FLUX_CO2.198
C FLUX_CO2.199
DO I = IFROM_CYC, ITO_CYC ORH1F305.395
IF (L_OCARB14) THEN ORH1F305.396
C FLUX_CO2.207
C Compute change in the C14/C12 ratio (ie tracer 14) FLUX_CO2.208
C The C14/C12 ratio for the pre-industrial atmosphere is FLUX_CO2.209
C fixed (arbitrarily) at 100 model units FLUX_CO2.210
IF(TA(I,1,TCO2_TRACER).GT.0.0) THEN OJP0F404.608
TA(I,1,C14_TRACER) = TA(I,1,C14_TRACER) OJP0F404.609
& + INVADE(I) * TIMESTEP * 100.0/DZ(1) OJP0F404.610
& * 1000.0 / (60.0*60.0*24.0*360.0) OJP0F404.611
& * (c14to12_atm-TA(I,1,C14_TRACER)) / TA(I,1,TCO2_TRACER) OJP0F404.612
ENDIF FLUX_CO2.214
ENDIF ORH1F305.397
TA(I,1,TCO2_TRACER) = TA(I,1,TCO2_TRACER) OJP0F404.613
& + CO2_FLUX(I) * TIMESTEP * 100.0/DZ(1) OJP0F404.614
& * 1000.0/(60.0*60.0*24.0*360.0) OJP0F404.615
ENDDO FLUX_CO2.219
C FLUX_CO2.220
C FLUX_CO2.221
IF (L_OCYCLIC) THEN ORH1F305.398
C Set cyclic boundary conditions FLUX_CO2.223
CO2_FLUX(1) = CO2_FLUX(IMT-1) FLUX_CO2.224
CO2_FLUX(IMT) = CO2_FLUX(2) FLUX_CO2.225
invade(1) = invade(IMT-1) FLUX_CO2.226
invade(IMT) = invade(2) FLUX_CO2.227
evade(1) = evade(IMT-1) FLUX_CO2.228
evade(IMT) = evade(2) FLUX_CO2.229
ENDIF ORH1F305.399
C FLUX_CO2.231
CL Return from FLUXCO2 FLUX_CO2.232
C FLUX_CO2.233
RETURN FLUX_CO2.234
END FLUX_CO2.235
*ENDIF FLUX_CO2.236