- Timestamp:
- 2011-10-22T15:46:41+02:00 (13 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2011/dev_LOCEAN_2011/NEMOGCM/NEMO/TOP_SRC/PISCES/p4zflx.F90
r2715 r2977 9 9 !! 1.0 ! 2004 (O. Aumont) modifications 10 10 !! 2.0 ! 2007-12 (C. Ethe, G. Madec) F90 11 !! ! 2011-02 (J. Simeon, J. Orr) Include total atm P correction 11 12 !!---------------------------------------------------------------------- 12 13 #if defined key_pisces … … 16 17 !! p4z_flx : CALCULATES GAS EXCHANGE AND CHEMISTRY AT SEA SURFACE 17 18 !! p4z_flx_init : Read the namelist 18 !!---------------------------------------------------------------------- 19 USE trc 20 USE oce_trc ! 21 USE trc 22 USE sms_pisces 23 USE prtctl_trc 24 USE p4zche 25 USE iom 19 !! p4z_patm : Read sfc atm pressure [atm] for each grid cell 20 !!---------------------------------------------------------------------- 21 USE oce_trc ! shared variables between ocean and passive tracers 22 USE trc ! passive tracers common variables 23 USE sms_pisces ! PISCES Source Minus Sink variables 24 USE p4zche ! Chemical model 25 USE prtctl_trc ! print control for debugging 26 USE iom ! I/O manager 27 USE fldread ! read input fields 26 28 #if defined key_cpl_carbon_cycle 27 USE sbc_oce , ONLY : atm_co229 USE sbc_oce, ONLY : atm_co2 ! atmospheric pCO2 28 30 #endif 29 31 … … 35 37 PUBLIC p4z_flx_alloc 36 38 39 ! !!** Namelist nampisext ** 40 REAL(wp) :: atcco2 = 278._wp !: pre-industrial atmospheric [co2] (ppm) 41 LOGICAL :: ln_co2int = .FALSE. !: flag to read in a file and interpolate atmospheric pco2 or not 42 CHARACTER(len=34) :: clname = 'atcco2.txt' !: filename of pco2 values 43 INTEGER :: nn_offset = 0 !: Offset model-data start year (default = 0) 44 45 !! Variables related to reading atmospheric CO2 time history 46 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: atcco2h, years 47 INTEGER :: nmaxrec, numco2 48 49 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:) :: patm ! atmospheric pressure at kt [N/m2] 50 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_patm ! structure of input fields (file informations, fields read) 51 52 37 53 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: oce_co2 !: ocean carbon flux 38 54 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: satmco2 !: atmospheric pco2 … … 41 57 REAL(wp) :: t_atm_co2_flx !: global mean of atmospheric pco2 42 58 REAL(wp) :: area !: ocean surface 43 REAL(wp) :: atcco2 = 278._wp !: pre-industrial atmospheric [co2] (ppm)44 REAL(wp) :: atcox = 0.20946_wp !:45 59 REAL(wp) :: xconv = 0.01_wp / 3600._wp !: coefficients for conversion 46 60 … … 60 74 !! ** Purpose : CALCULATES GAS EXCHANGE AND CHEMISTRY AT SEA SURFACE 61 75 !! 62 !! ** Method : - ??? 76 !! ** Method : 77 !! - Include total atm P correction via Esbensen & Kushnir (1981) 78 !! - Pressure correction NOT done for key_cpl_carbon_cycle 79 !! - Remove Wanninkhof chemical enhancement; 80 !! - Add option for time-interpolation of atcco2.txt 63 81 !!--------------------------------------------------------------------- 64 82 USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released 65 USE wrk_nemo, ONLY: zkgco2 => wrk_2d_1 , zkgo2 => wrk_2d_2 , zh2co3 => wrk_2d_366 USE wrk_nemo, ONLY: zoflx => wrk_2d_ 4 , zkg => wrk_2d_567 USE wrk_nemo, ONLY: zdpco2 => wrk_2d_ 6 , zdpo2 => wrk_2d_783 USE wrk_nemo, ONLY: zkgco2 => wrk_2d_11 , zkgo2 => wrk_2d_12 , zh2co3 => wrk_2d_13 84 USE wrk_nemo, ONLY: zoflx => wrk_2d_14 , zkg => wrk_2d_15 85 USE wrk_nemo, ONLY: zdpco2 => wrk_2d_16 , zdpo2 => wrk_2d_17 68 86 ! 69 87 INTEGER, INTENT(in) :: kt ! 70 88 ! 71 INTEGER :: ji, jj, j rorr89 INTEGER :: ji, jj, jm, iind, iindm1 72 90 REAL(wp) :: ztc, ztc2, ztc3, zws, zkgwan 73 91 REAL(wp) :: zfld, zflu, zfld16, zflu16, zfact 74 92 REAL(wp) :: zph, zah2, zbot, zdic, zalk, zsch_o2, zalka, zsch_co2 93 REAL(wp) :: zyr_dec, zdco2dt 75 94 CHARACTER (len=25) :: charout 76 95 !!--------------------------------------------------------------------- 77 96 78 IF( wrk_in_use(2, 1 ,2,3,4,5,6,7) ) THEN97 IF( wrk_in_use(2, 11,12,13,14,15,16,17) ) THEN 79 98 CALL ctl_stop('p4z_flx: requested workspace arrays unavailable') ; RETURN 80 99 ENDIF … … 84 103 ! IS USED TO COMPUTE AIR-SEA FLUX OF CO2 85 104 105 CALL p4z_patm( kt ) ! Get sea-level pressure (E&K [1981] climatology) for use in flux calcs 106 107 IF( ln_co2int ) THEN 108 ! Linear temporal interpolation of atmospheric pco2. atcco2.txt has annual values. 109 ! Caveats: First column of .txt must be in years, decimal years preferably. 110 ! For nn_offset, if your model year is iyy, nn_offset=(years(1)-iyy) 111 ! then the first atmospheric CO2 record read is at years(1) 112 zyr_dec = REAL( nyear + nn_offset, wp ) + REAL( nday_year, wp ) / REAL( nyear_len(1), wp ) 113 jm = 2 114 DO WHILE( jm <= nmaxrec .AND. years(jm-1) < zyr_dec .AND. years(jm) >= zyr_dec ) ; jm = jm + 1 ; END DO 115 iind = jm ; iindm1 = jm - 1 116 zdco2dt = ( atcco2h(iind) - atcco2h(iindm1) ) / ( years(iind) - years(iindm1) + rtrn ) 117 atcco2 = zdco2dt * ( zyr_dec - years(iindm1) ) + atcco2h(iindm1) 118 satmco2(:,:) = atcco2 119 ENDIF 120 86 121 #if defined key_cpl_carbon_cycle 87 122 satmco2(:,:) = atm_co2(:,:) 88 123 #endif 89 124 90 DO jrorr = 1, 10 91 125 DO jm = 1, 10 92 126 !CDIR NOVERRCHK 93 127 DO jj = 1, jpj … … 137 171 ! Compute the piston velocity for O2 and CO2 138 172 zkgwan = 0.3 * zws + 2.5 * ( 0.5246 + 0.016256 * ztc + 0.00049946 * ztc2 ) 173 zkgwan = zkgwan * xconv * ( 1.- fr_i(ji,jj) ) * tmask(ji,jj,1) 139 174 # if defined key_degrad 140 zkgwan = zkgwan * xconv * ( 1.- fr_i(ji,jj) ) * tmask(ji,jj,1) * facvol(ji,jj,1) 141 #else 142 zkgwan = zkgwan * xconv * ( 1.- fr_i(ji,jj) ) * tmask(ji,jj,1) 175 zkgwan = zkgwan * facvol(ji,jj,1) 143 176 #endif 144 177 ! compute gas exchange for CO2 and O2 … … 151 184 DO ji = 1, jpi 152 185 ! Compute CO2 flux for the sea and air 153 zfld = satmco2(ji,jj) * tmask(ji,jj,1) * chemc(ji,jj,1) * zkgco2(ji,jj)154 zflu = zh2co3(ji,jj) * tmask(ji,jj,1) * zkgco2(ji,jj) 186 zfld = satmco2(ji,jj) * patm(ji,jj) * tmask(ji,jj,1) * chemc(ji,jj,1) * zkgco2(ji,jj) ! (mol/L) * (m/s) 187 zflu = zh2co3(ji,jj) * tmask(ji,jj,1) * zkgco2(ji,jj) ! (mol/L) (m/s) ? 155 188 oce_co2(ji,jj) = ( zfld - zflu ) * rfact * e1e2t(ji,jj) * tmask(ji,jj,1) * 1000. 156 189 ! compute the trend … … 158 191 159 192 ! Compute O2 flux 160 zfld16 = atcox * chemc(ji,jj,2) * tmask(ji,jj,1) * zkgo2(ji,jj)193 zfld16 = atcox * patm(ji,jj) * chemc(ji,jj,2) * tmask(ji,jj,1) * zkgo2(ji,jj) ! (mol/L) * (m/s) 161 194 zflu16 = trn(ji,jj,1,jpoxy) * tmask(ji,jj,1) * zkgo2(ji,jj) 162 195 tra(ji,jj,1,jpoxy) = tra(ji,jj,1,jpoxy) + ( zfld16 - zflu16 ) / fse3t(ji,jj,1) 163 196 164 #if defined key_diatrc 165 ! Save diagnostics 166 # if ! defined key_iomput 167 zfact = 1. / e1e2t(ji,jj) / rfact 168 trc2d(ji,jj,jp_pcs0_2d ) = oce_co2(ji,jj) * zfact 169 trc2d(ji,jj,jp_pcs0_2d + 1) = ( zfld16 - zflu16 ) * 1000. * tmask(ji,jj,1) 170 trc2d(ji,jj,jp_pcs0_2d + 2) = zkgco2(ji,jj) * tmask(ji,jj,1) 171 trc2d(ji,jj,jp_pcs0_2d + 3) = ( satmco2(ji,jj) - zh2co3(ji,jj) / ( chemc(ji,jj,1) + rtrn ) ) & 172 & * tmask(ji,jj,1) 173 # else 174 zoflx(ji,jj) = ( zfld16 - zflu16 ) * 1000. * tmask(ji,jj,1) 175 zkg (ji,jj) = zkgco2(ji,jj) * tmask(ji,jj,1) 176 zdpco2(ji,jj) = ( satmco2(ji,jj) - zh2co3(ji,jj) / ( chemc(ji,jj,1) + rtrn ) ) * tmask(ji,jj,1) 177 zdpo2 (ji,jj) = ( atcox - trn(ji,jj,1,jpoxy) / ( chemc(ji,jj,2) + rtrn ) ) * tmask(ji,jj,1) 178 # endif 179 #endif 197 IF( ln_diatrc ) THEN ! Save diagnostics 198 IF( lk_iomput ) THEN 199 zoflx(ji,jj) = ( zfld16 - zflu16 ) * 1000. * tmask(ji,jj,1) 200 zkg (ji,jj) = zkgco2(ji,jj) * tmask(ji,jj,1) 201 zdpco2(ji,jj) = ( satmco2(ji,jj) * patm(ji,jj) - zh2co3(ji,jj) / ( chemc(ji,jj,1) + rtrn ) ) * tmask(ji,jj,1) 202 zdpo2 (ji,jj) = ( atcox * patm(ji,jj) - trn(ji,jj,1,jpoxy) / ( chemc(ji,jj,2) + rtrn ) ) * tmask(ji,jj,1) 203 ELSE 204 zfact = 1. / e1e2t(ji,jj) / rfact 205 trc2d(ji,jj,jp_pcs0_2d ) = oce_co2(ji,jj) * zfact 206 trc2d(ji,jj,jp_pcs0_2d + 1) = ( zfld16 - zflu16 ) * 1000. * tmask(ji,jj,1) 207 trc2d(ji,jj,jp_pcs0_2d + 2) = zkgco2(ji,jj) * tmask(ji,jj,1) 208 trc2d(ji,jj,jp_pcs0_2d + 3) = ( satmco2(ji,jj) * patm(ji,jj) - zh2co3(ji,jj) / ( chemc(ji,jj,1) + rtrn ) ) & 209 & * tmask(ji,jj,1) 210 ENDIF 211 ENDIF 180 212 END DO 181 213 END DO 182 214 183 t_oce_co2_flx = t_oce_co2_flx + glob_sum( oce_co2(:,:) ) 215 t_oce_co2_flx = t_oce_co2_flx + glob_sum( oce_co2(:,:) ) ! Cumulative Total Flux of Carbon 184 216 IF( kt == nitend ) THEN 185 t_atm_co2_flx = glob_sum( satmco2(:,:) * e1e2t(:,:) ) ! Total atmospheric pCO2186 ! 187 t_oce_co2_flx = (-1.) * t_oce_co2_flx * 12. / 1.e15 188 t_atm_co2_flx = t_atm_co2_flx / area 217 t_atm_co2_flx = glob_sum( satmco2(:,:) * patm(:,:) * e1e2t(:,:) ) ! Total atmospheric pCO2 218 ! 219 t_oce_co2_flx = (-1.) * t_oce_co2_flx * 12. / 1.e15 ! Conversion in PgC ; negative for out of the ocean 220 t_atm_co2_flx = t_atm_co2_flx / area ! global mean of atmospheric pCO2 189 221 ! 190 222 IF( lwp) THEN … … 205 237 ENDIF 206 238 207 # if defined key_diatrc && defined key_iomput 208 CALL iom_put( "Cflx" , oce_co2(:,:) / e1e2t(:,:) / rfact ) 209 CALL iom_put( "Oflx" , zoflx ) 210 CALL iom_put( "Kg" , zkg ) 211 CALL iom_put( "Dpco2", zdpco2 ) 212 CALL iom_put( "Dpo2" , zdpo2 ) 213 #endif 214 ! 215 IF( wrk_not_released(2, 1,2,3,4,5,6,7) ) CALL ctl_stop('p4z_flx: failed to release workspace arrays') 239 IF( ln_diatrc ) THEN 240 CALL iom_put( "Cflx" , oce_co2(:,:) / e1e2t(:,:) / rfact ) 241 CALL iom_put( "Oflx" , zoflx ) 242 CALL iom_put( "Kg" , zkg ) 243 CALL iom_put( "Dpco2", zdpco2 ) 244 CALL iom_put( "Dpo2" , zdpo2 ) 245 ENDIF 246 ! 247 IF( wrk_not_released(2, 11,12,13,14,15,16,17) ) & 248 & CALL ctl_stop('p4z_flx: failed to release workspace arrays') 216 249 ! 217 250 END SUBROUTINE p4z_flx … … 228 261 !! ** input : Namelist nampisext 229 262 !!---------------------------------------------------------------------- 230 NAMELIST/nampisext/ atcco2 231 !!---------------------------------------------------------------------- 232 ! 233 REWIND( numnat ) ! read numnat 234 READ ( numnat, nampisext ) 263 NAMELIST/nampisext/ln_co2int, atcco2, clname, nn_offset 264 INTEGER :: jm 265 !!---------------------------------------------------------------------- 266 ! 267 REWIND( numnatp ) ! read numnatp 268 READ ( numnatp, nampisext ) 235 269 ! 236 270 IF(lwp) THEN ! control print … … 238 272 WRITE(numout,*) ' Namelist parameters for air-sea exchange, nampisext' 239 273 WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~' 240 WRITE(numout,*) ' Atmospheric pCO2 atcco2 =', atcco2 274 WRITE(numout,*) ' Choice for reading in the atm pCO2 file or constant value, ln_co2int =', ln_co2int 275 WRITE(numout,*) ' ' 276 ENDIF 277 IF( .NOT.ln_co2int ) THEN 278 IF(lwp) THEN ! control print 279 WRITE(numout,*) ' Constant Atmospheric pCO2 value atcco2 =', atcco2 280 WRITE(numout,*) ' ' 281 ENDIF 282 satmco2(:,:) = atcco2 ! Initialisation of atmospheric pco2 283 ELSE 284 IF(lwp) THEN 285 WRITE(numout,*) ' Atmospheric pCO2 value from file clname =', TRIM( clname ) 286 WRITE(numout,*) ' Offset model-data start year nn_offset =', nn_offset 287 WRITE(numout,*) ' ' 288 ENDIF 289 CALL ctl_opn( numco2, TRIM( clname) , 'OLD', 'FORMATTED', 'SEQUENTIAL', -1 , numout, lwp ) 290 jm = 0 ! Count the number of record in co2 file 291 DO 292 READ(numco2,*,END=100) 293 jm = jm + 1 294 END DO 295 100 nmaxrec = jm - 1 296 ALLOCATE( years (nmaxrec) ) ; years (:) = 0._wp 297 ALLOCATE( atcco2h(nmaxrec) ) ; atcco2h(:) = 0._wp 298 299 REWIND(numco2) 300 DO jm = 1, nmaxrec ! get xCO2 data 301 READ(numco2, *) years(jm), atcco2h(jm) 302 IF(lwp) WRITE(numout, '(f6.0,f7.2)') years(jm), atcco2h(jm) 303 END DO 304 CLOSE(numco2) 241 305 ENDIF 242 306 ! … … 245 309 oce_co2(:,:) = 0._wp ! Initialization of Flux of Carbon 246 310 t_atm_co2_flx = 0._wp 247 !248 satmco2(:,:) = atcco2 ! Initialisation of atmospheric pco2249 311 t_oce_co2_flx = 0._wp 250 312 ! 251 313 END SUBROUTINE p4z_flx_init 252 314 315 SUBROUTINE p4z_patm( kt ) 316 317 !!---------------------------------------------------------------------- 318 !! *** ROUTINE p4z_atm *** 319 !! 320 !! ** Purpose : Read and interpolate the external atmospheric sea-levl pressure 321 !! ** Method : Read the files and interpolate the appropriate variables 322 !! 323 !!---------------------------------------------------------------------- 324 !! * arguments 325 INTEGER, INTENT( in ) :: kt ! ocean time step 326 ! 327 INTEGER :: ierr 328 CHARACTER(len=100) :: cn_dir ! Root directory for location of ssr files 329 TYPE(FLD_N) :: sn_patm ! informations about the fields to be read 330 !! 331 NAMELIST/nampisatm/ sn_patm, cn_dir 332 333 ! ! -------------------- ! 334 IF( kt == nit000 ) THEN ! First call kt=nit000 ! 335 ! ! -------------------- ! 336 ! !* set file information (default values) 337 ! ... default values (NB: frequency positive => hours, negative => months) 338 ! ! file ! frequency ! variable ! time intep ! clim ! 'yearly' or ! weights ! rotation ! 339 ! ! name ! (hours) ! name ! (T/F) ! (T/F) ! 'monthly' ! filename ! pairs ! 340 sn_patm = FLD_N( 'pres' , 24 , 'patm' , .false. , .true. , 'yearly' , '' , '' ) 341 cn_dir = './' ! directory in which the Patm data are 342 343 REWIND( numnatp ) !* read in namlist nampisatm 344 READ ( numnatp, nampisatm ) 345 ! 346 ALLOCATE( sf_patm(1), STAT=ierr ) !* allocate and fill sf_patm (forcing structure) with sn_patm 347 IF( ierr > 0 ) CALL ctl_stop( 'STOP', 'p4z_flx: unable to allocate sf_patm structure' ) 348 ! 349 CALL fld_fill( sf_patm, (/ sn_patm /), cn_dir, 'p4z_flx', 'Atmospheric pressure ', 'nampisatm' ) 350 ALLOCATE( sf_patm(1)%fnow(jpi,jpj,1) ) 351 IF( sn_patm%ln_tint ) ALLOCATE( sf_patm(1)%fdta(jpi,jpj,1,2) ) 352 ! 353 ENDIF 354 ! 355 CALL fld_read( kt, 1, sf_patm ) !* input Patm provided at kt + 1/2 356 patm(:,:) = sf_patm(1)%fnow(:,:,1) ! atmospheric pressure 357 358 END SUBROUTINE p4z_patm 253 359 254 360 INTEGER FUNCTION p4z_flx_alloc() … … 256 362 !! *** ROUTINE p4z_flx_alloc *** 257 363 !!---------------------------------------------------------------------- 258 ALLOCATE( oce_co2(jpi,jpj), satmco2(jpi,jpj), STAT=p4z_flx_alloc )364 ALLOCATE( oce_co2(jpi,jpj), satmco2(jpi,jpj), patm(jpi,jpj), STAT=p4z_flx_alloc ) 259 365 ! 260 366 IF( p4z_flx_alloc /= 0 ) CALL ctl_warn('p4z_flx_alloc : failed to allocate arrays')
Note: See TracChangeset
for help on using the changeset viewer.
