Changeset 9285


Ignore:
Timestamp:
2018-01-26T13:27:36+01:00 (3 years ago)
Author:
jpalmier
Message:

JPALM —26-01-2018 — change MEDUSA atm co2 input method. now reads a file a variable length

Location:
branches/NERC/dev_r5518_GO6_CO2_cmip/NEMOGCM
Files:
1 added
3 edited

Legend:

Unmodified
Added
Removed
  • branches/NERC/dev_r5518_GO6_CO2_cmip/NEMOGCM/NEMO/TOP_SRC/MEDUSA/sms_medusa.F90

    r9258 r9285  
    283283#if defined key_roam 
    284284!!---------------------------------------------------------------------- 
    285 !! Atmospheric pCO2 data (1859 to 2100 inclusive) 
    286 !!---------------------------------------------------------------------- 
    287 !! 
    288    REAL(wp), DIMENSION(242)         ::   hist_pco2 !: pCO2 
    289  
    290 # if defined key_rcp26 
    291       !! UKMO, run AJKKH + KAAEC, RCP 2.6, pCO2 time evolution 
    292       DATA hist_pco2 / 286.0230, 286.1730, 286.3230, 286.4480, 286.5730, & 
    293       & 286.7230, 286.8480, 286.9480, 287.0480, 287.1730, & 
    294       & 287.3230, 287.4730, 287.6480, 287.8480, 288.0730, & 
    295       & 288.3480, 288.6480, 288.9730, 289.3470, 289.7470, & 
    296       & 290.1730, 290.6470, 291.1470, 291.6220, 292.0720, & 
    297       & 292.5220, 292.9220, 293.2470, 293.5220, 293.7470, & 
    298       & 293.9470, 294.1220, 294.2720, 294.4220, 294.5470, & 
    299       & 294.6470, 294.7470, 294.8470, 294.9710, 295.1710, & 
    300       & 295.4460, 295.7470, 296.0720, 296.4210, 296.7710, & 
    301       & 297.1460, 297.5710, 298.0210, 298.4460, 298.8460, & 
    302       & 299.2460, 299.6450, 300.0210, 300.3710, 300.7200, & 
    303       & 301.0450, 301.3460, 301.6710, 302.0200, 302.3450, & 
    304       & 302.6450, 302.9700, 303.3450, 303.7200, 304.0700, & 
    305       & 304.4700, 304.9200, 305.3440, 305.7700, 306.2450, & 
    306       & 306.7190, 307.1700, 307.6440, 308.1190, 308.5440, & 
    307       & 308.9440, 309.3440, 309.6940, 309.9440, 310.1190, & 
    308       & 310.2440, 310.3190, 310.3190, 310.2440, 310.1440, & 
    309       & 310.0690, 310.0440, 310.0690, 310.1440, 310.2690, & 
    310       & 310.4440, 310.6940, 311.0430, 311.4440, 311.8690, & 
    311       & 312.3680, 312.9430, 313.5430, 314.1680, 314.7900, & 
    312       & 315.4430, 316.2150, 317.0170, 317.7370, 318.3400, & 
    313       & 318.8680, 319.5900, 320.5890, 321.5470, 322.5770, & 
    314       & 323.8440, 324.9260, 325.7960, 327.0810, 328.6180, & 
    315       & 329.6830, 330.5250, 331.6880, 333.2120, 334.7870, & 
    316       & 336.4640, 338.2990, 339.6660, 340.7310, 342.1360, & 
    317       & 343.7200, 345.2200, 346.7350, 348.5820, 350.6740, & 
    318       & 352.4230, 353.7910, 354.9530, 355.8210, 356.7130, & 
    319       & 358.0630, 359.7720, 361.3970, 363.0900, 365.2560, & 
    320       & 367.2810, 368.7980, 370.4000, 372.4550, 374.6920, & 
    321       & 376.7440, 378.7440, 380.7580, 382.7080, 384.7300, & 
    322       & 386.9310, 389.2150, 391.4910, 393.7710, 396.0460, & 
    323       & 398.3240, 400.6080, 402.8950, 405.1780, 407.4550, & 
    324       & 409.7260, 411.9930, 414.2500, 416.4410, 418.5280, & 
    325       & 420.5250, 422.4390, 424.2720, 426.0200, 427.6750, & 
    326       & 429.2360, 430.7050, 432.0850, 433.3580, 434.5140, & 
    327       & 435.5740, 436.5490, 437.4420, 438.2550, 438.9810, & 
    328       & 439.6110, 440.1430, 440.5770, 440.9450, 441.2660, & 
    329       & 441.5410, 441.7840, 442.0050, 442.2040, 442.3780, & 
    330       & 442.5210, 442.6200, 442.6720, 442.6810, 442.6540, & 
    331       & 442.5830, 442.4670, 442.3270, 442.1680, 441.9960, & 
    332       & 441.8060, 441.5930, 441.3440, 441.0540, 440.7230, & 
    333       & 440.3510, 439.9300, 439.4650, 438.9730, 438.4630, & 
    334       & 437.9400, 437.4020, 436.8400, 436.2640, 435.6850, & 
    335       & 435.1030, 434.5160, 433.9170, 433.3060, 432.7010, & 
    336       & 432.1110, 431.5380, 430.9810, 430.4320, 429.8860, & 
    337       & 429.3370, 428.7810, 428.2220, 427.6490, 427.0660, & 
    338       & 426.4890, 425.9270, 425.3840, 424.8610, 424.3540, & 
    339       & 423.8540, 423.3540, 422.8530, 422.3510, 421.8410, & 
    340       & 421.3250, 420.8190 / 
    341 # else 
    342       !! UKMO, run AJKKH + KAAEF, RCP 8.5, pCO2 time evolution 
    343       DATA hist_pco2 / 286.0230, 286.1730, 286.3230, 286.4480, 286.5730, & 
    344       & 286.7230, 286.8480, 286.9480, 287.0480, 287.1730, & 
    345       & 287.3230, 287.4730, 287.6480, 287.8480, 288.0730, & 
    346       & 288.3480, 288.6480, 288.9730, 289.3470, 289.7470, & 
    347       & 290.1730, 290.6470, 291.1470, 291.6220, 292.0720, & 
    348       & 292.5220, 292.9220, 293.2470, 293.5220, 293.7470, & 
    349       & 293.9470, 294.1220, 294.2720, 294.4220, 294.5470, & 
    350       & 294.6470, 294.7470, 294.8470, 294.9710, 295.1710, & 
    351       & 295.4460, 295.7470, 296.0720, 296.4210, 296.7710, & 
    352       & 297.1460, 297.5710, 298.0210, 298.4460, 298.8460, & 
    353       & 299.2460, 299.6450, 300.0210, 300.3710, 300.7200, & 
    354       & 301.0450, 301.3460, 301.6710, 302.0200, 302.3450, & 
    355       & 302.6450, 302.9700, 303.3450, 303.7200, 304.0700, & 
    356       & 304.4700, 304.9200, 305.3440, 305.7700, 306.2450, & 
    357       & 306.7190, 307.1700, 307.6440, 308.1190, 308.5440, & 
    358       & 308.9440, 309.3440, 309.6940, 309.9440, 310.1190, & 
    359       & 310.2440, 310.3190, 310.3190, 310.2440, 310.1440, & 
    360       & 310.0690, 310.0440, 310.0690, 310.1440, 310.2690, & 
    361       & 310.4440, 310.6940, 311.0430, 311.4440, 311.8690, & 
    362       & 312.3680, 312.9430, 313.5430, 314.1680, 314.7900, & 
    363       & 315.4430, 316.2150, 317.0170, 317.7370, 318.3400, & 
    364       & 318.8680, 319.5900, 320.5890, 321.5470, 322.5770, & 
    365       & 323.8440, 324.9260, 325.7960, 327.0810, 328.6180, & 
    366       & 329.6830, 330.5250, 331.6880, 333.2120, 334.7870, & 
    367       & 336.4640, 338.2990, 339.6660, 340.7310, 342.1360, & 
    368       & 343.7200, 345.2200, 346.7350, 348.5820, 350.6740, & 
    369       & 352.4230, 353.7910, 354.9530, 355.8210, 356.7130, & 
    370       & 358.0630, 359.7720, 361.3970, 363.0900, 365.2560, & 
    371       & 367.2810, 368.7980, 370.4000, 372.4550, 374.6920, & 
    372       & 376.7440, 378.7440, 380.7580, 382.7080, 384.7300, & 
    373       & 386.9420, 389.2540, 391.5670, 393.9370, 396.3920, & 
    374       & 398.9320, 401.5550, 404.2550, 407.0220, 409.8530, & 
    375       & 412.7470, 415.7050, 418.7210, 421.7880, 424.9180, & 
    376       & 428.1200, 431.3970, 434.7470, 438.1650, 441.6410, & 
    377       & 445.1700, 448.7530, 452.3920, 456.0950, 459.8810, & 
    378       & 463.7680, 467.7660, 471.8750, 476.0960, 480.4210, & 
    379       & 484.8390, 489.3470, 493.9430, 498.6400, 503.4380, & 
    380       & 508.3410, 513.3630, 518.5160, 523.8050, 529.2290, & 
    381       & 534.7780, 540.4450, 546.2230, 552.1120, 558.1110, & 
    382       & 564.2110, 570.4130, 576.7390, 583.1990, 589.7980, & 
    383       & 596.5390, 603.4110, 610.4060, 617.4940, 624.6500, & 
    384       & 631.8800, 639.1750, 646.5360, 653.9800, 661.5230, & 
    385       & 669.1840, 676.9570, 684.8290, 692.7790, 700.7690, & 
    386       & 708.8050, 716.8870, 725.0020, 733.1770, 741.3900, & 
    387       & 749.6700, 758.0480, 766.5050, 775.0350, 783.6110, & 
    388       & 792.2200, 800.8740, 809.5680, 818.2760, 827.0090, & 
    389       & 835.8020, 844.6550, 853.5730, 862.5690, 871.6190, & 
    390       & 880.7020, 889.8240, 898.9590, 908.1270, 917.3080, & 
    391       & 926.4960, 935.7040 / 
    392 # endif 
     285!! JPALM -- change hist_pco2 init 
     286!!---------------------------------------------------------------------- 
     287!! 
     288   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:)   ::   hist_pco2 !: pCO2 
     289   INTEGER :: co2_yinit, co2_yend, co2_rec 
    393290#endif 
    394291 
  • branches/NERC/dev_r5518_GO6_CO2_cmip/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90

    r9257 r9285  
    8181                                            gdept_0, gdept_n,               & 
    8282                                            gdepw_0, gdepw_n,               & 
    83                                             nday_year, nsec_day, nyear,     & 
     83                                            nday_year, nsec_day,            & 
     84                                            nyear, nyear_len,               & 
    8485                                            rdt, tmask, mig, mjg, nimpp,    & 
    8586                                            njmpp  
     
    102103      !! JPALM (27-06-2016): add lk_oasis for CO2 and DMS coupling with atm 
    103104      USE sbc_oce,                    ONLY: lk_oasis 
    104       USE sms_medusa,                 ONLY: hist_pco2 
     105      USE sms_medusa,                 ONLY: hist_pco2, co2_yinit, co2_yend 
    105106      USE trc,                        ONLY: ln_rsttr, nittrc000, trn 
    106107      USE bio_medusa_init_mod,        ONLY: bio_medusa_init 
     
    294295      !! 
    295296      !! what's atmospheric pCO2 doing? (data start in 1859) 
    296       iyr1  = nyear - 1859 + 1 
     297      iyr1  = nyear  
    297298      iyr2  = iyr1 + 1 
    298       if (iyr1 .le. 1) then 
    299          !! before 1860 
    300          f_xco2a(:,:) = hist_pco2(1) 
    301       elseif (iyr2 .ge. 242) then 
    302          !! after 2099 
    303          f_xco2a(:,:) = hist_pco2(242) 
     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) 
    304305      else 
    305306         !! just right 
     
    307308         fq1 = hist_pco2(iyr2) 
    308309         fq2 = real(nsec_day) / (60.0 * 60.0 * 24.0) 
     310         fq3 = (real(nday_year) - 1.0 + fq2) / real(nyear_len(1)) 
    309311         !! AXY (14/06/12): tweaked to make more sense (and be correct) 
    310 #  if defined key_bs_axy_yrlen 
    311          !! bugfix: for 360d year with HadGEM2-ES forcing 
    312          fq3 = (real(nday_year) - 1.0 + fq2) / 360.0   
    313 #  else 
    314          !! original use of 365 days (not accounting for leap year or  
    315          !! 360d year) 
    316          fq3 = (real(nday_year) - 1.0 + fq2) / 365.0 
    317 #  endif 
     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 
    318320         fq4 = (fq0 * (1.0 - fq3)) + (fq1 * fq3) 
    319321         f_xco2a(:,:) = fq4 
     
    322324      !! OCMIP pre-industrial pCO2 
    323325      !! f_xco2a(:,:) = 284.725  !! CMIP5 pre-industrial pCO2 
    324       f_xco2a = 284.317          !! CMIP6 pre-industrial pCO2  
     326      f_xco2a(:,:) = 284.317          !! CMIP6 pre-industrial pCO2  
    325327#  endif 
    326       !! IF(lwp) WRITE(numout,*) ' MEDUSA nyear     =', nyear 
    327       !! IF(lwp) WRITE(numout,*) ' MEDUSA nsec_day  =', real(nsec_day) 
    328       !! IF(lwp) WRITE(numout,*) ' MEDUSA nday_year =', real(nday_year) 
     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) 
    329332      !! AXY (29/01/14): remove surplus diagnostics 
    330       !! IF(lwp) WRITE(numout,*) ' MEDUSA fq0       =', fq0 
    331       !! IF(lwp) WRITE(numout,*) ' MEDUSA fq1       =', fq1 
    332       !! IF(lwp) WRITE(numout,*) ' MEDUSA fq2       =', fq2 
    333       !! IF(lwp) WRITE(numout,*) ' MEDUSA fq3       =', fq3 
    334       IF(lwp) WRITE(numout,*) ' MEDUSA atm pCO2  =', f_xco2a(1,1) 
     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) 
    335338# endif 
    336339 
  • branches/NERC/dev_r5518_GO6_CO2_cmip/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcini_medusa.F90

    r9114 r9285  
    324324      IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~~' 
    325325      IF(lwp) CALL flush(numout) 
     326      !! 
     327      !!---------------------------------------------------------------------- 
     328      !! JPALM (23-01-2017): new way to initialize CO2-atm for cmip6  
     329      !!                     initially done in trcsms_medusa 
     330      !!---------------------------------------------------------------------- 
     331      !!      
     332      IF(lwp) WRITE(numout,*) ' trc_ini_medusa: initialisating atm CO2 record' 
     333      CALL trc_ini_medusa_co2atm 
     334 
    326335 
    327336   END SUBROUTINE trc_ini_medusa 
     
    480489   END SUBROUTINE trc_ini_medusa_river 
    481490    
     491   SUBROUTINE trc_ini_medusa_co2atm 
     492      !!---------------------------------------------------------------------- 
     493      !!                     ***  trc_ini_medusa_co2atm  ***   
     494      !! 
     495      !! ** Purpose :   initialization atmospheric co2 record 
     496      !! 
     497      !! ** Method  : - Read the xco2 file 
     498      !!---------------------------------------------------------------------- 
     499      INTEGER  ::  jn, jl, jm, io, ierr, inum 
     500      INTEGER  ::  iskip = 4   ! number of 1st descriptor lines 
     501      REAL(wp) ::  zyy 
     502      CHARACTER (len=10) ::   clname = 'xco2.atm'  !! atm CO2 record file 
     503      !!---------------------------------------------------------------------- 
     504 
     505      IF(lwp) WRITE(numout,*) 
     506      IF(lwp) WRITE(numout,*) ' trc_ini_medusa_co2atm: initialisation of atm CO2 historical record' 
     507      IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~' 
     508 
     509 
     510      IF(lwp) WRITE(numout,*) 'read of formatted file xco2.atm' 
     511 
     512      CALL ctl_opn( inum, clname, 'OLD', 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE. ) 
     513      REWIND(inum) 
     514 
     515      ! compute the number of year in the file 
     516      ! file starts in 1849 do jn represent the record number in the file. 
     517      ! the year is jn + yinit - 1 
     518      jn = 1 
     519      DO 
     520        READ(inum,'(1x)',END=100) 
     521        jn = jn + 1 
     522      END DO 
     523 100  co2_rec = jn - 1 - iskip 
     524      IF ( lwp) WRITE(numout,*) '    ', co2_rec ,' years read in the file' 
     525      !                                ! Allocate CO2 hist arrays 
     526      ierr = 0  
     527      ALLOCATE( hist_pco2(co2_rec), STAT=ierr ) 
     528      IF( ierr > 0 ) THEN 
     529         CALL ctl_stop( 'trc_ini_medusa_co2atm: unable to allocate  array' )   ; RETURN 
     530      ENDIF 
     531 
     532      REWIND(inum) 
     533 
     534      DO jm = 1, iskip        ! Skip over 1st six descriptor lines 
     535         READ(inum,'(1x)') 
     536      END DO 
     537      ! file starts in 1931 do jn represent the year in the century.jhh 
     538      ! Read file till the end 
     539      ! allocate start and end year of the file 
     540      jn = 1 
     541      DO 
     542        READ(inum,*, IOSTAT=io) zyy, hist_pco2(jn) 
     543        IF( io < 0 ) exit 
     544        IF(jn==1) co2_yinit = int(zyy) 
     545        jn = jn + 1 
     546      END DO 
     547      co2_yend = co2_yinit + co2_rec - 1 
     548 
     549      IF(lwp) THEN        ! Control print 
     550         WRITE(numout,*) 
     551         WRITE(numout,*) 'CO2 hist start year: ', co2_yinit 
     552         WRITE(numout,*) 'CO2 hist end   year: ', co2_yend 
     553         WRITE(numout,*) ' Year   xCO2 atm ' 
     554         DO jn = 1, co2_rec 
     555            jl = jn + co2_yinit - 1 
     556            WRITE(numout, '( 1I4, 6F9.2)') jl, hist_pco2(jn) 
     557         END DO 
     558      ENDIF 
     559 
     560   END SUBROUTINE trc_ini_medusa_co2atm 
     561 
     562 
    482563#else 
    483564   !!---------------------------------------------------------------------- 
Note: See TracChangeset for help on using the changeset viewer.