Ignore:
Timestamp:
2018-02-01T12:58:33+01:00 (3 years ago)
Author:
jpalmier
Message:

JPALM — correct the CO2 interpolation for ocean-only; atm CO2 read from file in ocean-only whatever record length; report corrctly what is read in ocean.output

Location:
branches/NERC/dev_r5518_GO6_CO2_cmip/NEMOGCM
Files:
4 edited

Legend:

Unmodified
Added
Removed
  • branches/NERC/dev_r5518_GO6_CO2_cmip/NEMOGCM/CONFIG/SHARED/xco2.atm

    r9285 r9298  
    33% Year  xCO2 
    44%       (ppm) 
    5 1849.5  284.3200000 
     51849.5  284.3170000 
    661850.5  284.3169999 
    771851.5  284.4509998 
  • branches/NERC/dev_r5518_GO6_CO2_cmip/NEMOGCM/NEMO/TOP_SRC/MEDUSA/sms_medusa.F90

    r9285 r9298  
    287287!! 
    288288   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:)   ::   hist_pco2 !: pCO2 
    289    INTEGER :: co2_yinit, co2_yend, co2_rec 
     289   INTEGER  :: co2_rec 
     290   REAL(wp) :: co2_yinit, co2_yend 
    290291#endif 
    291292 
  • branches/NERC/dev_r5518_GO6_CO2_cmip/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90

    r9285 r9298  
    182182      !!  
    183183      !! 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 
    185186      !! 
    186187      !! T and S check temporary variable 
     
    295296      !! 
    296297      !! 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 
    306320         !! 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)) 
    321325         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 
    323335#  if defined key_axy_pi_co2 
    324336      !! OCMIP pre-industrial pCO2 
     
    326338      f_xco2a(:,:) = 284.317          !! CMIP6 pre-industrial pCO2  
    327339#  endif 
    328        IF(lwp) WRITE(numout,*) ' MEDUSA nyear     =', nyear 
    329        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) 
    332340      !! 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) 
    338343# endif 
    339344 
  • branches/NERC/dev_r5518_GO6_CO2_cmip/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcini_medusa.F90

    r9285 r9298  
    497497      !! ** Method  : - Read the xco2 file 
    498498      !!---------------------------------------------------------------------- 
    499       INTEGER  ::  jn, jl, jm, io, ierr, inum 
     499      INTEGER  ::  jn, jm, io, ierr, inum 
    500500      INTEGER  ::  iskip = 4   ! number of 1st descriptor lines 
    501       REAL(wp) ::  zyy 
     501      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:)   ::   zyy !: xCO2 record years 
    502502      CHARACTER (len=10) ::   clname = 'xco2.atm'  !! atm CO2 record file 
    503503      !!---------------------------------------------------------------------- 
     
    525525      !                                ! Allocate CO2 hist arrays 
    526526      ierr = 0  
    527       ALLOCATE( hist_pco2(co2_rec), STAT=ierr ) 
     527      ALLOCATE( hist_pco2(co2_rec),zyy(co2_rec), STAT=ierr ) 
    528528      IF( ierr > 0 ) THEN 
    529529         CALL ctl_stop( 'trc_ini_medusa_co2atm: unable to allocate  array' )   ; RETURN 
     
    540540      jn = 1 
    541541      DO 
    542         READ(inum,*, IOSTAT=io) zyy, hist_pco2(jn) 
     542        READ(inum,*, IOSTAT=io) zyy(jn), hist_pco2(jn) 
    543543        IF( io < 0 ) exit 
    544         IF(jn==1) co2_yinit = int(zyy) 
     544        IF(jn==1) co2_yinit = zyy(jn) 
    545545        jn = jn + 1 
    546546      END DO 
    547       co2_yend = co2_yinit + co2_rec - 1 
     547      co2_yend = co2_yinit + real(co2_rec) - 1. 
    548548 
    549549      IF(lwp) THEN        ! Control print 
     
    553553         WRITE(numout,*) ' Year   xCO2 atm ' 
    554554         DO jn = 1, co2_rec 
    555             jl = jn + co2_yinit - 1 
    556             WRITE(numout, '( 1I4, 6F9.2)') jl, hist_pco2(jn) 
     555            WRITE(numout, '( 5F7.1, 6F9.2)') zyy(jn), hist_pco2(jn) 
    557556         END DO 
    558557      ENDIF 
Note: See TracChangeset for help on using the changeset viewer.