Changeset 9298 for branches/NERC/dev_r5518_GO6_CO2_cmip/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90
- Timestamp:
- 2018-02-01T12:58:33+01:00 (6 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/NERC/dev_r5518_GO6_CO2_cmip/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90
r9285 r9298 182 182 !! 183 183 !! temporary variables 184 REAL(wp) :: fq0,fq1,fq2,fq3,fq4 184 REAL(wp) :: fq3,fq4 185 REAL(wp) :: this_y, this_d, this_s, fyear 185 186 !! 186 187 !! T and S check temporary variable … … 295 296 !! 296 297 !! what's atmospheric pCO2 doing? (data start in 1859) 297 iyr1 = nyear 298 iyr2 = iyr1 + 1 299 if (iyr1 .le. co2_yinit) then 300 !! before 1850 301 f_xco2a(:,:) = 284.317 302 elseif (iyr2 .ge. co2_yend) then 303 !! after 2010 304 f_xco2a(:,:) = hist_pco2(co2_yend - co2_yinit +1) 305 else 298 !! AXY - JPALM new interpolation scheme usinf nyear_len 299 this_y = real(nyear) 300 this_d = real(nday_year) 301 this_s = real(nsec_day) 302 !! 303 fyear = this_y + ((this_d - 1) + (this_s / (60. * 60. * 24.))) / real(nyear_len(1)) 304 !! 305 IF ( ( kt == nittrc000 ) .AND. lwp ) THEN 306 WRITE(numout,*) ' MEDUSA year =', this_y 307 WRITE(numout,*) ' Year length =', real(nyear_len(1)) 308 WRITE(numout,*) ' MEDUSA nday_year =', this_d 309 WRITE(numout,*) ' MEDUSA nsec_day =', this_s 310 ENDIF 311 !! 312 !! different case test 313 IF (fyear .LE. co2_yinit) THEN 314 !! before first record -- pre-industrial value 315 f_xco2a(:,:) = hist_pco2(1) 316 ELSEIF (fyear .GE. co2_yend) THEN 317 !! after last record - continue to use the last value 318 f_xco2a(:,:) = hist_pco2(int(co2_yend - co2_yinit + 1.) ) 319 ELSE 306 320 !! just right 307 fq0 = hist_pco2(iyr1) 308 fq1 = hist_pco2(iyr2) 309 fq2 = real(nsec_day) / (60.0 * 60.0 * 24.0) 310 fq3 = (real(nday_year) - 1.0 + fq2) / real(nyear_len(1)) 311 !! AXY (14/06/12): tweaked to make more sense (and be correct) 312 !!# if defined key_bs_axy_yrlen 313 !! !! bugfix: for 360d year with HadGEM2-ES forcing 314 !! fq3 = (real(nday_year) - 1.0 + fq2) / 360.0 315 !!# else 316 !! !! original use of 365 days (not accounting for leap year or 317 !! !! 360d year) 318 !! fq3 = (real(nday_year) - 1.0 + fq2) / 365.0 319 !!# endif 320 fq4 = (fq0 * (1.0 - fq3)) + (fq1 * fq3) 321 iyr1 = int(fyear - co2_yinit) + 1 322 iyr2 = iyr1 + 1 323 fq3 = fyear - real(iyr1) - co2_yinit + 1. 324 fq4 = ((1 - fq3) * hist_pco2(iyr1)) + (fq3 * hist_pco2(iyr2)) 321 325 f_xco2a(:,:) = fq4 322 endif 326 !! 327 IF ( ( kt == nittrc000 ) .AND. lwp ) THEN 328 WRITE(numout,*) ' MEDUSA year1 =', iyr1 329 WRITE(numout,*) ' MEDUSA year2 =', iyr2 330 WRITE(numout,*) ' xCO2 year1 =', hist_pco2(iyr1) 331 WRITE(numout,*) ' xCO2 year2 =', hist_pco2(iyr2) 332 WRITE(numout,*) ' Year2 weight =', fq3 333 ENDIF 334 ENDIF 323 335 # if defined key_axy_pi_co2 324 336 !! OCMIP pre-industrial pCO2 … … 326 338 f_xco2a(:,:) = 284.317 !! CMIP6 pre-industrial pCO2 327 339 # endif 328 IF(lwp) WRITE(numout,*) ' MEDUSA nyear =', nyear329 IF(lwp) WRITE(numout,*) ' MEDUSA nyear =', real(nyear_len(1))330 IF(lwp) WRITE(numout,*) ' MEDUSA nsec_day =', real(nsec_day)331 IF(lwp) WRITE(numout,*) ' MEDUSA nday_year =', real(nday_year)332 340 !! AXY (29/01/14): remove surplus diagnostics 333 IF(lwp) WRITE(numout,*) ' MEDUSA fq0 =', fq0 334 IF(lwp) WRITE(numout,*) ' MEDUSA fq1 =', fq1 335 IF(lwp) WRITE(numout,*) ' MEDUSA fq2 =', fq2 336 IF(lwp) WRITE(numout,*) ' MEDUSA fq3 =', fq3 337 IF(lwp) WRITE(numout,*) ' MEDUSA atm pCO2 =', f_xco2a(1,1) 341 IF ( ( kt == nittrc000 ) .AND. lwp ) & 342 WRITE(numout,*) ' final atm xCO2 =', f_xco2a(1,1) 338 343 # endif 339 344
Note: See TracChangeset
for help on using the changeset viewer.