Changeset 3294 for trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/p4zflx.F90
- Timestamp:
- 2012-01-28T17:44:18+01:00 (12 years ago)
- File:
-
- 1 edited
-
trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/p4zflx.F90 (modified) (14 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/TOP_SRC/PISCES/p4zflx.F90
r2715 r3294 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 ! !!* nampisatm namelist (Atmospheric PRessure) * 50 LOGICAL, PUBLIC :: ln_presatm = .true. !: ref. pressure: global mean Patm (F) or a constant (F) 51 52 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:) :: patm ! atmospheric pressure at kt [N/m2] 53 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_patm ! structure of input fields (file informations, fields read) 54 55 37 56 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: oce_co2 !: ocean carbon flux 38 57 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: satmco2 !: atmospheric pco2 … … 41 60 REAL(wp) :: t_atm_co2_flx !: global mean of atmospheric pco2 42 61 REAL(wp) :: area !: ocean surface 43 REAL(wp) :: atcco2 = 278._wp !: pre-industrial atmospheric [co2] (ppm)44 REAL(wp) :: atcox = 0.20946_wp !:45 62 REAL(wp) :: xconv = 0.01_wp / 3600._wp !: coefficients for conversion 46 63 … … 60 77 !! ** Purpose : CALCULATES GAS EXCHANGE AND CHEMISTRY AT SEA SURFACE 61 78 !! 62 !! ** Method : - ??? 79 !! ** Method : 80 !! - Include total atm P correction via Esbensen & Kushnir (1981) 81 !! - Pressure correction NOT done for key_cpl_carbon_cycle 82 !! - Remove Wanninkhof chemical enhancement; 83 !! - Add option for time-interpolation of atcco2.txt 63 84 !!--------------------------------------------------------------------- 64 USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released65 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_768 85 ! 69 86 INTEGER, INTENT(in) :: kt ! 70 87 ! 71 INTEGER :: ji, jj, j rorr88 INTEGER :: ji, jj, jm, iind, iindm1 72 89 REAL(wp) :: ztc, ztc2, ztc3, zws, zkgwan 73 90 REAL(wp) :: zfld, zflu, zfld16, zflu16, zfact 74 91 REAL(wp) :: zph, zah2, zbot, zdic, zalk, zsch_o2, zalka, zsch_co2 92 REAL(wp) :: zyr_dec, zdco2dt 75 93 CHARACTER (len=25) :: charout 94 REAL(wp), POINTER, DIMENSION(:,:) :: zkgco2, zkgo2, zh2co3, zoflx 76 95 !!--------------------------------------------------------------------- 77 78 IF( wrk_in_use(2, 1,2,3,4,5,6,7) ) THEN 79 CALL ctl_stop('p4z_flx: requested workspace arrays unavailable') ; RETURN 80 ENDIF 96 ! 97 IF( nn_timing == 1 ) CALL timing_start('p4z_flx') 98 ! 99 CALL wrk_alloc( jpi, jpj, zkgco2, zkgo2, zh2co3, zoflx ) 100 ! 81 101 82 102 ! SURFACE CHEMISTRY (PCO2 AND [H+] IN … … 84 104 ! IS USED TO COMPUTE AIR-SEA FLUX OF CO2 85 105 106 IF( kt /= nit000 ) CALL p4z_patm( kt ) ! Get sea-level pressure (E&K [1981] climatology) for use in flux calcs 107 108 IF( ln_co2int ) THEN 109 ! Linear temporal interpolation of atmospheric pco2. atcco2.txt has annual values. 110 ! Caveats: First column of .txt must be in years, decimal years preferably. 111 ! For nn_offset, if your model year is iyy, nn_offset=(years(1)-iyy) 112 ! then the first atmospheric CO2 record read is at years(1) 113 zyr_dec = REAL( nyear + nn_offset, wp ) + REAL( nday_year, wp ) / REAL( nyear_len(1), wp ) 114 jm = 2 115 DO WHILE( jm <= nmaxrec .AND. years(jm-1) < zyr_dec .AND. years(jm) >= zyr_dec ) ; jm = jm + 1 ; END DO 116 iind = jm ; iindm1 = jm - 1 117 zdco2dt = ( atcco2h(iind) - atcco2h(iindm1) ) / ( years(iind) - years(iindm1) + rtrn ) 118 atcco2 = zdco2dt * ( zyr_dec - years(iindm1) ) + atcco2h(iindm1) 119 satmco2(:,:) = atcco2 120 ENDIF 121 86 122 #if defined key_cpl_carbon_cycle 87 123 satmco2(:,:) = atm_co2(:,:) 88 124 #endif 89 125 90 DO jrorr = 1, 10 91 126 DO jm = 1, 10 92 127 !CDIR NOVERRCHK 93 128 DO jj = 1, jpj … … 137 172 ! Compute the piston velocity for O2 and CO2 138 173 zkgwan = 0.3 * zws + 2.5 * ( 0.5246 + 0.016256 * ztc + 0.00049946 * ztc2 ) 174 zkgwan = zkgwan * xconv * ( 1.- fr_i(ji,jj) ) * tmask(ji,jj,1) 139 175 # 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) 176 zkgwan = zkgwan * facvol(ji,jj,1) 143 177 #endif 144 178 ! compute gas exchange for CO2 and O2 … … 151 185 DO ji = 1, jpi 152 186 ! 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) 187 zfld = satmco2(ji,jj) * patm(ji,jj) * tmask(ji,jj,1) * chemc(ji,jj,1) * zkgco2(ji,jj) ! (mol/L) * (m/s) 188 zflu = zh2co3(ji,jj) * tmask(ji,jj,1) * zkgco2(ji,jj) ! (mol/L) (m/s) ? 155 189 oce_co2(ji,jj) = ( zfld - zflu ) * rfact * e1e2t(ji,jj) * tmask(ji,jj,1) * 1000. 156 190 ! compute the trend … … 158 192 159 193 ! Compute O2 flux 160 zfld16 = atcox * chemc(ji,jj,2) * tmask(ji,jj,1) * zkgo2(ji,jj)194 zfld16 = atcox * patm(ji,jj) * chemc(ji,jj,2) * tmask(ji,jj,1) * zkgo2(ji,jj) ! (mol/L) * (m/s) 161 195 zflu16 = trn(ji,jj,1,jpoxy) * tmask(ji,jj,1) * zkgo2(ji,jj) 162 tra(ji,jj,1,jpoxy) = tra(ji,jj,1,jpoxy) + ( zfld16 - zflu16 ) / fse3t(ji,jj,1) 163 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 196 zoflx(ji,jj) = zfld16 - zflu16 197 tra(ji,jj,1,jpoxy) = tra(ji,jj,1,jpoxy) + zoflx(ji,jj) / fse3t(ji,jj,1) 180 198 END DO 181 199 END DO 182 200 183 t_oce_co2_flx = t_oce_co2_flx + glob_sum( oce_co2(:,:) ) 201 t_oce_co2_flx = t_oce_co2_flx + glob_sum( oce_co2(:,:) ) ! Cumulative Total Flux of Carbon 184 202 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 203 t_atm_co2_flx = glob_sum( satmco2(:,:) * patm(:,:) * e1e2t(:,:) ) ! Total atmospheric pCO2 204 ! 205 t_oce_co2_flx = (-1.) * t_oce_co2_flx * 12. / 1.e15 ! Conversion in PgC ; negative for out of the ocean 206 t_atm_co2_flx = t_atm_co2_flx / area ! global mean of atmospheric pCO2 189 207 ! 190 208 IF( lwp) THEN … … 205 223 ENDIF 206 224 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') 225 IF( ln_diatrc ) THEN 226 IF( lk_iomput ) THEN 227 CALL iom_put( "Cflx" , oce_co2(:,:) / e1e2t(:,:) / rfact ) 228 CALL iom_put( "Oflx" , zoflx(:,:) * 1000 * tmask(:,:,1) ) 229 CALL iom_put( "Kg" , zkgco2(:,:) * tmask(:,:,1) ) 230 CALL iom_put( "Dpco2", ( satmco2(:,:) * patm(:,:) - zh2co3(:,:) / ( chemc(:,:,1) + rtrn ) ) * tmask(:,:,1) ) 231 CALL iom_put( "Dpo2" , ( atcox * patm(:,:) - trn(:,:,1,jpoxy) / ( chemc(:,:,2) + rtrn ) ) * tmask(:,:,1) ) 232 ELSE 233 trc2d(:,:,jp_pcs0_2d ) = oce_co2(:,:) / e1e2t(:,:) / rfact 234 trc2d(:,:,jp_pcs0_2d + 1) = zoflx(:,:) * 1000 * tmask(:,:,1) 235 trc2d(:,:,jp_pcs0_2d + 2) = zkgco2(:,:) * tmask(:,:,1) 236 trc2d(:,:,jp_pcs0_2d + 3) = ( satmco2(:,:) * patm(:,:) - zh2co3(:,:) / ( chemc(:,:,1) + rtrn ) ) * tmask(:,:,1) 237 ENDIF 238 ENDIF 239 ! 240 CALL wrk_dealloc( jpi, jpj, zkgco2, zkgo2, zh2co3, zoflx ) 241 ! 242 IF( nn_timing == 1 ) CALL timing_stop('p4z_flx') 216 243 ! 217 244 END SUBROUTINE p4z_flx … … 225 252 !! 226 253 !! ** Method : Read the nampisext namelist and check the parameters 227 !! called at the first timestep (nit 000)254 !! called at the first timestep (nittrc000) 228 255 !! ** input : Namelist nampisext 229 256 !!---------------------------------------------------------------------- 230 NAMELIST/nampisext/ atcco2 231 !!---------------------------------------------------------------------- 232 ! 233 REWIND( numnat ) ! read numnat 234 READ ( numnat, nampisext ) 257 NAMELIST/nampisext/ln_co2int, atcco2, clname, nn_offset 258 INTEGER :: jm 259 !!---------------------------------------------------------------------- 260 ! 261 REWIND( numnatp ) ! read numnatp 262 READ ( numnatp, nampisext ) 235 263 ! 236 264 IF(lwp) THEN ! control print … … 238 266 WRITE(numout,*) ' Namelist parameters for air-sea exchange, nampisext' 239 267 WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~' 240 WRITE(numout,*) ' Atmospheric pCO2 atcco2 =', atcco2 268 WRITE(numout,*) ' Choice for reading in the atm pCO2 file or constant value, ln_co2int =', ln_co2int 269 WRITE(numout,*) ' ' 270 ENDIF 271 IF( .NOT.ln_co2int ) THEN 272 IF(lwp) THEN ! control print 273 WRITE(numout,*) ' Constant Atmospheric pCO2 value atcco2 =', atcco2 274 WRITE(numout,*) ' ' 275 ENDIF 276 satmco2(:,:) = atcco2 ! Initialisation of atmospheric pco2 277 ELSE 278 IF(lwp) THEN 279 WRITE(numout,*) ' Atmospheric pCO2 value from file clname =', TRIM( clname ) 280 WRITE(numout,*) ' Offset model-data start year nn_offset =', nn_offset 281 WRITE(numout,*) ' ' 282 ENDIF 283 CALL ctl_opn( numco2, TRIM( clname) , 'OLD', 'FORMATTED', 'SEQUENTIAL', -1 , numout, lwp ) 284 jm = 0 ! Count the number of record in co2 file 285 DO 286 READ(numco2,*,END=100) 287 jm = jm + 1 288 END DO 289 100 nmaxrec = jm - 1 290 ALLOCATE( years (nmaxrec) ) ; years (:) = 0._wp 291 ALLOCATE( atcco2h(nmaxrec) ) ; atcco2h(:) = 0._wp 292 293 REWIND(numco2) 294 DO jm = 1, nmaxrec ! get xCO2 data 295 READ(numco2, *) years(jm), atcco2h(jm) 296 IF(lwp) WRITE(numout, '(f6.0,f7.2)') years(jm), atcco2h(jm) 297 END DO 298 CLOSE(numco2) 241 299 ENDIF 242 300 ! … … 245 303 oce_co2(:,:) = 0._wp ! Initialization of Flux of Carbon 246 304 t_atm_co2_flx = 0._wp 247 !248 satmco2(:,:) = atcco2 ! Initialisation of atmospheric pco2249 305 t_oce_co2_flx = 0._wp 250 306 ! 307 CALL p4z_patm( nit000 ) 308 ! 251 309 END SUBROUTINE p4z_flx_init 252 310 311 SUBROUTINE p4z_patm( kt ) 312 313 !!---------------------------------------------------------------------- 314 !! *** ROUTINE p4z_atm *** 315 !! 316 !! ** Purpose : Read and interpolate the external atmospheric sea-levl pressure 317 !! ** Method : Read the files and interpolate the appropriate variables 318 !! 319 !!---------------------------------------------------------------------- 320 !! * arguments 321 INTEGER, INTENT( in ) :: kt ! ocean time step 322 ! 323 INTEGER :: ierr 324 CHARACTER(len=100) :: cn_dir ! Root directory for location of ssr files 325 TYPE(FLD_N) :: sn_patm ! informations about the fields to be read 326 !! 327 NAMELIST/nampisatm/ ln_presatm, sn_patm, cn_dir 328 329 ! ! -------------------- ! 330 IF( kt == nit000 ) THEN ! First call kt=nittrc000 ! 331 ! ! -------------------- ! 332 ! !* set file information (default values) 333 ! ... default values (NB: frequency positive => hours, negative => months) 334 ! ! file ! frequency ! variable ! time intep ! clim ! 'yearly' or ! weights ! rotation ! 335 ! ! name ! (hours) ! name ! (T/F) ! (T/F) ! 'monthly' ! filename ! pairs ! 336 sn_patm = FLD_N( 'pres' , 24 , 'patm' , .false. , .true. , 'yearly' , '' , '' ) 337 cn_dir = './' ! directory in which the Patm data are 338 339 REWIND( numnatp ) !* read in namlist nampisatm 340 READ ( numnatp, nampisatm ) 341 ! 342 ! 343 IF(lwp) THEN !* control print 344 WRITE(numout,*) 345 WRITE(numout,*) ' Namelist nampisatm : Atmospheric Pressure as external forcing' 346 WRITE(numout,*) ' constant atmopsheric pressure (F) or from a file (T) ln_presatm = ', ln_presatm 347 WRITE(numout,*) 348 ENDIF 349 ! 350 IF( ln_presatm ) THEN 351 ALLOCATE( sf_patm(1), STAT=ierr ) !* allocate and fill sf_patm (forcing structure) with sn_patm 352 IF( ierr > 0 ) CALL ctl_stop( 'STOP', 'p4z_flx: unable to allocate sf_patm structure' ) 353 ! 354 CALL fld_fill( sf_patm, (/ sn_patm /), cn_dir, 'p4z_flx', 'Atmospheric pressure ', 'nampisatm' ) 355 ALLOCATE( sf_patm(1)%fnow(jpi,jpj,1) ) 356 IF( sn_patm%ln_tint ) ALLOCATE( sf_patm(1)%fdta(jpi,jpj,1,2) ) 357 ENDIF 358 ! 359 IF( .NOT.ln_presatm ) patm(:,:) = 1.e0 ! Initialize patm if no reading from a file 360 ! 361 ENDIF 362 ! 363 IF( ln_presatm ) THEN 364 CALL fld_read( kt, 1, sf_patm ) !* input Patm provided at kt + 1/2 365 patm(:,:) = sf_patm(1)%fnow(:,:,1) ! atmospheric pressure 366 ENDIF 367 ! 368 END SUBROUTINE p4z_patm 253 369 254 370 INTEGER FUNCTION p4z_flx_alloc() … … 256 372 !! *** ROUTINE p4z_flx_alloc *** 257 373 !!---------------------------------------------------------------------- 258 ALLOCATE( oce_co2(jpi,jpj), satmco2(jpi,jpj), STAT=p4z_flx_alloc )374 ALLOCATE( oce_co2(jpi,jpj), satmco2(jpi,jpj), patm(jpi,jpj), STAT=p4z_flx_alloc ) 259 375 ! 260 376 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.
