Changeset 9285 for branches/NERC/dev_r5518_GO6_CO2_cmip/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90

Timestamp:

2018-01-26T13:27:36+01:00 (6 years ago)

Author:

jpalmier

Message:

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

File:

• branches/NERC/dev_r5518_GO6_CO2_cmip/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90 (modified) (5 diffs)

branches/NERC/dev_r5518_GO6_CO2_cmip/NEMOGCM/NEMO/TOP_SRC/MEDUSA/trcbio_medusa.F90

@@ -81 +81 @@
                                             gdept_0, gdept_n,               &
                                             gdepw_0, gdepw_n,               &
-                                            nday_year, nsec_day, nyear,     &
+                                            nday_year, nsec_day,            &
+                                            nyear, nyear_len,               &
                                             rdt, tmask, mig, mjg, nimpp,    &
                                             njmpp 
@@ -102 +103 @@
       !! JPALM (27-06-2016): add lk_oasis for CO2 and DMS coupling with atm
       USE sbc_oce,                    ONLY: lk_oasis
-      USE sms_medusa,                 ONLY: hist_pco2
+      USE sms_medusa,                 ONLY: hist_pco2, co2_yinit, co2_yend
       USE trc,                        ONLY: ln_rsttr, nittrc000, trn
       USE bio_medusa_init_mod,        ONLY: bio_medusa_init
@@ -294 +295 @@
       !!
       !! what's atmospheric pCO2 doing? (data start in 1859)
-      iyr1  = nyear - 1859 + 1
+      iyr1  = nyear 
       iyr2  = iyr1 + 1
-      if (iyr1 .le. 1) then
-         !! before 1860
-         f_xco2a(:,:) = hist_pco2(1)
-      elseif (iyr2 .ge. 242) then
-         !! after 2099
-         f_xco2a(:,:) = hist_pco2(242)
+      if (iyr1 .le. co2_yinit) then
+         !! before 1850
+         f_xco2a(:,:) = 284.317  
+      elseif (iyr2 .ge. co2_yend) then
+         !! after 2010
+         f_xco2a(:,:) = hist_pco2(co2_yend - co2_yinit +1)
       else
          !! just right
@@ -307 +308 @@
          fq1 = hist_pco2(iyr2)
          fq2 = real(nsec_day) / (60.0 * 60.0 * 24.0)
+         fq3 = (real(nday_year) - 1.0 + fq2) / real(nyear_len(1))
          !! AXY (14/06/12): tweaked to make more sense (and be correct)
-#  if defined key_bs_axy_yrlen
-         !! bugfix: for 360d year with HadGEM2-ES forcing
-         fq3 = (real(nday_year) - 1.0 + fq2) / 360.0  
-#  else
-         !! original use of 365 days (not accounting for leap year or 
-         !! 360d year)
-         fq3 = (real(nday_year) - 1.0 + fq2) / 365.0
-#  endif
+!!#  if defined key_bs_axy_yrlen
+!!         !! bugfix: for 360d year with HadGEM2-ES forcing
+!!         fq3 = (real(nday_year) - 1.0 + fq2) / 360.0  
+!!#  else
+!!         !! original use of 365 days (not accounting for leap year or 
+!!         !! 360d year)
+!!         fq3 = (real(nday_year) - 1.0 + fq2) / 365.0
+!!#  endif
          fq4 = (fq0 * (1.0 - fq3)) + (fq1 * fq3)
          f_xco2a(:,:) = fq4
@@ -322 +324 @@
       !! OCMIP pre-industrial pCO2
       !! f_xco2a(:,:) = 284.725  !! CMIP5 pre-industrial pCO2
-      f_xco2a = 284.317          !! CMIP6 pre-industrial pCO2 
+      f_xco2a(:,:) = 284.317          !! CMIP6 pre-industrial pCO2 
 #  endif
-      !! IF(lwp) WRITE(numout,*) ' MEDUSA nyear     =', nyear
-      !! IF(lwp) WRITE(numout,*) ' MEDUSA nsec_day  =', real(nsec_day)
-      !! IF(lwp) WRITE(numout,*) ' MEDUSA nday_year =', real(nday_year)
+       IF(lwp) WRITE(numout,*) ' MEDUSA nyear     =', nyear
+       IF(lwp) WRITE(numout,*) ' MEDUSA nyear     =', real(nyear_len(1))
+       IF(lwp) WRITE(numout,*) ' MEDUSA nsec_day  =', real(nsec_day)
+       IF(lwp) WRITE(numout,*) ' MEDUSA nday_year =', real(nday_year)
       !! AXY (29/01/14): remove surplus diagnostics
-      !! IF(lwp) WRITE(numout,*) ' MEDUSA fq0       =', fq0
-      !! IF(lwp) WRITE(numout,*) ' MEDUSA fq1       =', fq1
-      !! IF(lwp) WRITE(numout,*) ' MEDUSA fq2       =', fq2
-      !! IF(lwp) WRITE(numout,*) ' MEDUSA fq3       =', fq3
-      IF(lwp) WRITE(numout,*) ' MEDUSA atm pCO2  =', f_xco2a(1,1)
+       IF(lwp) WRITE(numout,*) ' MEDUSA fq0       =', fq0
+       IF(lwp) WRITE(numout,*) ' MEDUSA fq1       =', fq1
+       IF(lwp) WRITE(numout,*) ' MEDUSA fq2       =', fq2
+       IF(lwp) WRITE(numout,*) ' MEDUSA fq3       =', fq3
+       IF(lwp) WRITE(numout,*) ' MEDUSA atm pCO2  =', f_xco2a(1,1)
 # endif