Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
MODULE p4zflx
!!======================================================================
!! *** MODULE p4zflx ***
!! TOP : PISCES CALCULATES GAS EXCHANGE AND CHEMISTRY AT SEA SURFACE
!!======================================================================
!! History : - ! 1988-07 (E. MAIER-REIMER) Original code
!! - ! 1998 (O. Aumont) additions
!! - ! 1999 (C. Le Quere) modifications
!! 1.0 ! 2004 (O. Aumont) modifications
!! 2.0 ! 2007-12 (C. Ethe, G. Madec) F90
!! ! 2011-02 (J. Simeon, J. Orr) Include total atm P correction
!!----------------------------------------------------------------------
!! p4z_flx : CALCULATES GAS EXCHANGE AND CHEMISTRY AT SEA SURFACE
!! p4z_flx_init : Read the namelist
!! p4z_patm : Read sfc atm pressure [atm] for each grid cell
!!----------------------------------------------------------------------
USE oce_trc ! shared variables between ocean and passive tracers
USE trc ! passive tracers common variables
USE sms_pisces ! PISCES Source Minus Sink variables
USE p4zche ! Chemical model
USE prtctl ! print control for debugging
USE iom ! I/O manager
USE fldread ! read input fields
IMPLICIT NONE
PRIVATE
PUBLIC p4z_flx
PUBLIC p4z_flx_init
PUBLIC p4z_flx_alloc
! !!** Namelist nampisext **
REAL(wp) :: atcco2 !: pre-industrial atmospheric [co2] (ppm)
LOGICAL :: ln_co2int !: flag to read in a file and interpolate atmospheric pco2 or not
CHARACTER(len=34) :: clname !: filename of pco2 values
INTEGER :: nn_offset !: Offset model-data start year (default = 0)
!! Variables related to reading atmospheric CO2 time history
INTEGER :: nmaxrec, numco2 !
REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: atcco2h, years !
! !!* nampisatm namelist (Atmospheric PRessure) *
LOGICAL, PUBLIC :: ln_presatm !: ref. pressure: global mean Patm (F) or a constant (F)
LOGICAL, PUBLIC :: ln_presatmco2 !: accounting for spatial atm CO2 in the compuation of carbon flux (T) or not (F)
REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:) :: patm ! atmospheric pressure at kt [N/m2]
TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_patm ! structure of input fields (file informations, fields read)
TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_atmco2 ! structure of input fields (file informations, fields read)
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: satmco2 !: atmospheric pco2
REAL(wp) :: xconv = 0.01_wp / 3600._wp !: coefficients for conversion
!! * Substitutions
# include "do_loop_substitute.h90"
# include "domzgr_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/TOP 4.0 , NEMO Consortium (2018)
!! $Id: p4zflx.F90 15459 2021-10-29 08:19:18Z cetlod $
!! Software governed by the CeCILL license (see ./LICENSE)
!!----------------------------------------------------------------------
CONTAINS
SUBROUTINE p4z_flx ( kt, knt, Kbb, Kmm, Krhs )
!!---------------------------------------------------------------------
!! *** ROUTINE p4z_flx ***
!!
!! ** Purpose : CALCULATES GAS EXCHANGE AND CHEMISTRY AT SEA SURFACE
!!
!! ** Method :
!! - Include total atm P correction via Esbensen & Kushnir (1981)
!! - Remove Wanninkhof chemical enhancement;
!! - Add option for time-interpolation of atcco2.txt
!!---------------------------------------------------------------------
INTEGER, INTENT(in) :: kt, knt !
INTEGER, INTENT(in) :: Kbb, Kmm, Krhs ! time level indices
!
INTEGER :: ji, jj, jm, iind, iindm1
REAL(wp) :: ztc, ztc2, ztc3, ztc4, zws, zkgwan
REAL(wp) :: zfld, zflu, zfld16, zflu16, zfact
REAL(wp) :: zvapsw, zsal, zfco2, zxc2, xCO2approx, ztkel, zfugcoeff
REAL(wp) :: zph, zdic, zsch_o2, zsch_co2
REAL(wp) :: zyr_dec, zdco2dt
CHARACTER (len=25) :: charout
REAL(wp), DIMENSION(jpi,jpj) :: zkgco2, zkgo2, zh2co3, zoflx, zpco2atm, zpco2oce
!!---------------------------------------------------------------------
!
IF( ln_timing ) CALL timing_start('p4z_flx')
!
! SURFACE CHEMISTRY (PCO2 AND [H+] IN
! SURFACE LAYER); THE RESULT OF THIS CALCULATION
! IS USED TO COMPUTE AIR-SEA FLUX OF CO2
IF( kt /= nit000 .AND. .NOT.l_co2cpl .AND. knt == 1 ) CALL p4z_patm( kt ) ! Get sea-level pressure (E&K [1981] climatology) for use in flux calcs
IF( ln_co2int .AND. .NOT.ln_presatmco2 .AND. .NOT.l_co2cpl ) THEN
! Linear temporal interpolation of atmospheric pco2. atcco2.txt has annual values.
! Caveats: First column of .txt must be in years, decimal years preferably.
! For nn_offset, if your model year is iyy, nn_offset=(years(1)-iyy)
! then the first atmospheric CO2 record read is at years(1)
zyr_dec = REAL( nyear + nn_offset, wp ) + REAL( nday_year, wp ) / REAL( nyear_len(1), wp )
jm = 1
DO WHILE( jm <= nmaxrec .AND. years(jm) < zyr_dec ) ; jm = jm + 1 ; END DO
iind = jm ; iindm1 = jm - 1
zdco2dt = ( atcco2h(iind) - atcco2h(iindm1) ) / ( years(iind) - years(iindm1) + rtrn )
atcco2 = zdco2dt * ( zyr_dec - years(iindm1) ) + atcco2h(iindm1)
satmco2(:,:) = atcco2
ENDIF
IF( l_co2cpl ) satmco2(:,:) = atm_co2(:,:)
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
! DUMMY VARIABLES FOR DIC, H+, AND BORATE
zfact = rhop(ji,jj,1) / 1000. + rtrn
zdic = tr(ji,jj,1,jpdic,Kbb)
zph = MAX( hi(ji,jj,1), 1.e-10 ) / zfact
! CALCULATE [H2CO3]
zh2co3(ji,jj) = zdic/(1. + ak13(ji,jj,1)/zph + ak13(ji,jj,1)*ak23(ji,jj,1)/zph**2)
END_2D
! --------------
! COMPUTE FLUXES
! --------------
! FIRST COMPUTE GAS EXCHANGE COEFFICIENTS
! -------------------------------------------
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
ztc = MIN( 35., ts(ji,jj,1,jp_tem,Kmm) )
ztc2 = ztc * ztc
ztc3 = ztc * ztc2
ztc4 = ztc2 * ztc2
! Compute the schmidt Number both O2 and CO2
zsch_co2 = 2116.8 - 136.25 * ztc + 4.7353 * ztc2 - 0.092307 * ztc3 + 0.0007555 * ztc4
zsch_o2 = 1920.4 - 135.6 * ztc + 5.2122 * ztc2 - 0.109390 * ztc3 + 0.0009377 * ztc4
! wind speed
zws = wndm(ji,jj) * wndm(ji,jj)
! Compute the piston velocity for O2 and CO2
zkgwan = 0.251 * zws
zkgwan = zkgwan * xconv * ( 1.- fr_i(ji,jj) ) * tmask(ji,jj,1)
! compute gas exchange for CO2 and O2
zkgco2(ji,jj) = zkgwan * SQRT( 660./ zsch_co2 )
zkgo2 (ji,jj) = zkgwan * SQRT( 660./ zsch_o2 )
END_2D
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
ztkel = tempis(ji,jj,1) + 273.15
zsal = salinprac(ji,jj,1) + ( 1.- tmask(ji,jj,1) ) * 35.
zvapsw = EXP(24.4543 - 67.4509*(100.0/ztkel) - 4.8489*LOG(ztkel/100) - 0.000544*zsal)
zpco2atm(ji,jj) = satmco2(ji,jj) * ( patm(ji,jj) - zvapsw )
zxc2 = ( 1.0 - zpco2atm(ji,jj) * 1E-6 )**2
zfugcoeff = EXP( patm(ji,jj) * (chemc(ji,jj,2) + 2.0 * zxc2 * chemc(ji,jj,3) ) &
& / ( 82.05736 * ztkel ))
zfco2 = zpco2atm(ji,jj) * zfugcoeff
! Compute CO2 flux for the sea and air
zfld = zfco2 * chemc(ji,jj,1) * zkgco2(ji,jj) ! (mol/L) * (m/s)
zflu = zh2co3(ji,jj) * zkgco2(ji,jj) ! (mol/L) (m/s) ?
zpco2oce(ji,jj) = zh2co3(ji,jj) / ( chemc(ji,jj,1) * zfugcoeff + rtrn )
oce_co2(ji,jj) = ( zfld - zflu ) * tmask(ji,jj,1)
! compute the trend
tr(ji,jj,1,jpdic,Krhs) = tr(ji,jj,1,jpdic,Krhs) + oce_co2(ji,jj) * rfact2 / e3t(ji,jj,1,Kmm)
! Compute O2 flux
zfld16 = patm(ji,jj) * chemo2(ji,jj,1) * zkgo2(ji,jj) ! (mol/L) * (m/s)
zflu16 = tr(ji,jj,1,jpoxy,Kbb) * zkgo2(ji,jj)
zoflx(ji,jj) = ( zfld16 - zflu16 ) * tmask(ji,jj,1)
tr(ji,jj,1,jpoxy,Krhs) = tr(ji,jj,1,jpoxy,Krhs) + zoflx(ji,jj) * rfact2 / e3t(ji,jj,1,Kmm)
END_2D
IF( iom_use("tcflx") .OR. iom_use("tcflxcum") .OR. kt == nitrst &
& .OR. (ln_check_mass .AND. kt == nitend) ) &
t_oce_co2_flx = glob_sum( 'p4zflx', oce_co2(:,:) * e1e2t(:,:) * 1000. ) ! Total Flux of Carbon
t_oce_co2_flx_cum = t_oce_co2_flx_cum + t_oce_co2_flx ! Cumulative Total Flux of Carbon
! t_atm_co2_flx = glob_sum( 'p4zflx', satmco2(:,:) * e1e2t(:,:) ) ! Total atmospheric pCO2
t_atm_co2_flx = atcco2 ! Total atmospheric pCO2
IF(sn_cfctl%l_prttrc) THEN ! print mean trends (used for debugging)
WRITE(charout, FMT="('flx ')")
CALL prt_ctl_info( charout, cdcomp = 'top' )
CALL prt_ctl(tab4d_1=tr(:,:,:,:,Krhs), mask1=tmask, clinfo=ctrcnm)
ENDIF
IF( lk_iomput .AND. knt == nrdttrc ) THEN
CALL iom_put( "AtmCo2" , satmco2(:,:) * tmask(:,:,1) ) ! Atmospheric CO2 concentration
CALL iom_put( "Cflx" , oce_co2(:,:) * 1000. )
CALL iom_put( "Oflx" , zoflx(:,:) * 1000. )
CALL iom_put( "Kg" , zkgco2(:,:) * tmask(:,:,1) )
CALL iom_put( "Dpco2" , ( zpco2atm(:,:) - zpco2oce(:,:) ) * tmask(:,:,1) )
CALL iom_put( "pCO2sea" , zpco2oce(:,:) * tmask(:,:,1) )
CALL iom_put( "Dpo2" , ( atcox * patm(:,:) - atcox * tr(:,:,1,jpoxy,Kbb) / ( chemo2(:,:,1) + rtrn ) ) * tmask(:,:,1) )
CALL iom_put( "tcflx" , t_oce_co2_flx ) ! molC/s
CALL iom_put( "tcflxcum", t_oce_co2_flx_cum ) ! molC
ENDIF
!
IF( ln_timing ) CALL timing_stop('p4z_flx')
!
END SUBROUTINE p4z_flx
SUBROUTINE p4z_flx_init
!!----------------------------------------------------------------------
!! *** ROUTINE p4z_flx_init ***
!!
!! ** Purpose : Initialization of atmospheric conditions
!!
!! ** Method : Read the nampisext namelist and check the parameters
!! called at the first timestep (nittrc000)
!!
!! ** input : Namelist nampisext
!!----------------------------------------------------------------------
INTEGER :: jm, ios ! Local integer
!!
NAMELIST/nampisext/ln_co2int, atcco2, clname, nn_offset
!!----------------------------------------------------------------------
IF(lwp) THEN
WRITE(numout,*)
WRITE(numout,*) ' p4z_flx_init : atmospheric conditions for air-sea flux calculation'
WRITE(numout,*) ' ~~~~~~~~~~~~'
ENDIF
!
READ ( numnatp_ref, nampisext, IOSTAT = ios, ERR = 901)
901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nampisext in reference namelist' )
READ ( numnatp_cfg, nampisext, IOSTAT = ios, ERR = 902 )
902 IF( ios > 0 ) CALL ctl_nam ( ios , 'nampisext in configuration namelist' )
IF(lwm) WRITE ( numonp, nampisext )
!
IF(lwp) THEN ! control print
WRITE(numout,*) ' Namelist : nampisext --- parameters for air-sea exchange'
WRITE(numout,*) ' reading in the atm pCO2 file or constant value ln_co2int =', ln_co2int
ENDIF
!
CALL p4z_patm( nit000 )
!
IF( .NOT.ln_co2int .AND. .NOT.ln_presatmco2 ) THEN
IF(lwp) THEN ! control print
WRITE(numout,*) ' Constant Atmospheric pCO2 value atcco2 =', atcco2
ENDIF
satmco2(:,:) = atcco2 ! Initialisation of atmospheric pco2
ELSEIF( ln_co2int .AND. .NOT.ln_presatmco2 ) THEN
IF(lwp) THEN
WRITE(numout,*) ' Constant Atmospheric pCO2 value atcco2 =', atcco2
WRITE(numout,*) ' Atmospheric pCO2 value from file clname =', TRIM( clname )
WRITE(numout,*) ' Offset model-data start year nn_offset =', nn_offset
ENDIF
CALL ctl_opn( numco2, TRIM( clname) , 'OLD', 'FORMATTED', 'SEQUENTIAL', -1 , numout, lwp )
jm = 0 ! Count the number of record in co2 file
DO
READ(numco2,*,END=100)
jm = jm + 1
END DO
100 nmaxrec = jm - 1
ALLOCATE( years (nmaxrec) ) ; years (:) = 0._wp
ALLOCATE( atcco2h(nmaxrec) ) ; atcco2h(:) = 0._wp
!
REWIND(numco2)
DO jm = 1, nmaxrec ! get xCO2 data
READ(numco2, *) years(jm), atcco2h(jm)
IF(lwp) WRITE(numout, '(f6.0,f7.2)') years(jm), atcco2h(jm)
END DO
CLOSE(numco2)
ELSEIF( .NOT.ln_co2int .AND. ln_presatmco2 ) THEN
IF(lwp) WRITE(numout,*) ' Spatialized Atmospheric pCO2 from an external file'
ELSE
IF(lwp) WRITE(numout,*) ' Spatialized Atmospheric pCO2 from an external file'
ENDIF
!
oce_co2(:,:) = 0._wp ! Initialization of Flux of Carbon
t_oce_co2_flx = 0._wp
t_atm_co2_flx = 0._wp
!
END SUBROUTINE p4z_flx_init
SUBROUTINE p4z_patm( kt )
!!----------------------------------------------------------------------
!! *** ROUTINE p4z_atm ***
!!
!! ** Purpose : Read and interpolate the external atmospheric sea-level pressure
!! ** Method : Read the files and interpolate the appropriate variables
!!
!!----------------------------------------------------------------------
INTEGER, INTENT(in) :: kt ! ocean time step
!
INTEGER :: ierr, ios ! Local integer
CHARACTER(len=100) :: cn_dir ! Root directory for location of ssr files
TYPE(FLD_N) :: sn_patm ! informations about the fields to be read
TYPE(FLD_N) :: sn_atmco2 ! informations about the fields to be read
!!
NAMELIST/nampisatm/ ln_presatm, ln_presatmco2, sn_patm, sn_atmco2, cn_dir
!!----------------------------------------------------------------------
!
IF( kt == nit000 ) THEN !== First call kt=nittrc000 ==!
!
IF(lwp) THEN
WRITE(numout,*)
WRITE(numout,*) ' p4z_patm : sea-level atmospheric pressure'
WRITE(numout,*) ' ~~~~~~~~'
ENDIF
!
READ ( numnatp_ref, nampisatm, IOSTAT = ios, ERR = 901)
901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nampisatm in reference namelist' )
READ ( numnatp_cfg, nampisatm, IOSTAT = ios, ERR = 902 )
902 IF( ios > 0 ) CALL ctl_nam ( ios , 'nampisatm in configuration namelist' )
IF(lwm) WRITE ( numonp, nampisatm )
!
!
IF(lwp) THEN !* control print
WRITE(numout,*) ' Namelist : nampisatm --- Atmospheric Pressure as external forcing'
WRITE(numout,*) ' constant atmopsheric pressure (F) or from a file (T) ln_presatm = ', ln_presatm
WRITE(numout,*) ' spatial atmopsheric CO2 for flux calcs ln_presatmco2 = ', ln_presatmco2
ENDIF
!
IF( ln_presatm ) THEN
ALLOCATE( sf_patm(1), STAT=ierr ) !* allocate and fill sf_patm (forcing structure) with sn_patm
IF( ierr > 0 ) CALL ctl_stop( 'STOP', 'p4z_flx: unable to allocate sf_patm structure' )
!
CALL fld_fill( sf_patm, (/ sn_patm /), cn_dir, 'p4z_flx', 'Atmospheric pressure ', 'nampisatm' )
ALLOCATE( sf_patm(1)%fnow(jpi,jpj,1) )
IF( sn_patm%ln_tint ) ALLOCATE( sf_patm(1)%fdta(jpi,jpj,1,2) )
ENDIF
!
IF( ln_presatmco2 ) THEN
ALLOCATE( sf_atmco2(1), STAT=ierr ) !* allocate and fill sf_atmco2 (forcing structure) with sn_atmco2
IF( ierr > 0 ) CALL ctl_stop( 'STOP', 'p4z_flx: unable to allocate sf_atmco2 structure' )
!
CALL fld_fill( sf_atmco2, (/ sn_atmco2 /), cn_dir, 'p4z_flx', 'Atmospheric co2 partial pressure ', 'nampisatm' )
ALLOCATE( sf_atmco2(1)%fnow(jpi,jpj,1) )
IF( sn_atmco2%ln_tint ) ALLOCATE( sf_atmco2(1)%fdta(jpi,jpj,1,2) )
ENDIF
!
IF( .NOT.ln_presatm ) patm(:,:) = 1._wp ! Initialize patm if no reading from a file
!
ENDIF
!
IF( ln_presatm ) THEN
CALL fld_read( kt, 1, sf_patm ) !* input Patm provided at kt + 1/2
patm(:,:) = sf_patm(1)%fnow(:,:,1)/101325.0 ! atmospheric pressure
ENDIF
!
IF( ln_presatmco2 ) THEN
CALL fld_read( kt, 1, sf_atmco2 ) !* input atmco2 provided at kt + 1/2
satmco2(:,:) = sf_atmco2(1)%fnow(:,:,1) ! atmospheric pressure
ELSE
satmco2(:,:) = atcco2 ! Initialize atmco2 if no reading from a file
ENDIF
!
END SUBROUTINE p4z_patm
INTEGER FUNCTION p4z_flx_alloc()
!!----------------------------------------------------------------------
!! *** ROUTINE p4z_flx_alloc ***
!!----------------------------------------------------------------------
ALLOCATE( satmco2(jpi,jpj), patm(jpi,jpj), STAT=p4z_flx_alloc )
!
IF( p4z_flx_alloc /= 0 ) CALL ctl_stop( 'STOP', 'p4z_flx_alloc : failed to allocate arrays' )
!
END FUNCTION p4z_flx_alloc
!!======================================================================
END MODULE p4zflx