main program and module mods

Small modifications (number of variables, dimension names, etc.) to accommodate adjustment of TraCE-21ka distribution files.

README.md

@@ -22,7 +22,8 @@ The current version is v1.1.  Relative to previous versions, this version includ
 - specification of the infofile path and name on the command line, so that once built locally, the programs can be run in a terminal window;
 - addition of the adjusted month lengths, and beginning, middle, and ending dates to the output file;
 - a choice of two mean-preserving interpolation methods, including the Epstein (1991) approach implemented in v1.0, as well as the Harzallah (1995) iterated-spline approach;
-- the inclusion of a subroutine, `enforce_mean()` that requires the pseudo-daily interpolated values to have the same monthly mean as the input monthly values.
+- the inclusion of a subroutine, `enforce_mean()` that requires the pseudo-daily interpolated values to have the same monthly mean as the input monthly values;
+- some modifications of I/O to accommodate (some) "non-standard" files (e.g. TraCE-21ka).
 
 ## Interpolation methods ##
 

f90/main_programs/cal_adjust.f90

@@ -133,7 +133,7 @@ program cal_adjust
 character(32)           :: calendar_type        ! CF calendar type (e.g. "noleap", "proleptic_gregorian", etc.)
 
 ! strings
-character(12)            :: activity             ! simple label (e.g. "PMIP3" or "PMIP4"), or some other label (e.g. "transient")
+character(12)           :: activity             ! simple label (e.g. "PMIP3" or "PMIP4"), or some other label (e.g. "transient")
 character(64)           :: variable             ! variable name
 character(8)            :: time_freq            ! CMIP/PMIP-style output time-frequency type (e.g. "Amon", "Aclim", "day", "Oclim" ...)
 character(14)           :: tol_string           ! enforce_mean() tolerance value as string
@@ -173,6 +173,7 @@ program cal_adjust
 integer(4)              :: iostatus             ! IOSTAT value
 integer(4)              :: infofile_len         ! length of the info file
 integer(4)              :: infofile_status      ! info file read status
+integer(4)              :: fill_status          ! _FillValue status
 
 ! file paths and names
 character(2048)         :: source_path, source_file, ncfile_in, adjusted_path, adjusted_file, ncfile_out
@@ -254,7 +255,7 @@ program cal_adjust
     case ("harzallah", "Harzallah")
         interp_method = "Harzallah"
         max_nctrl_in = 12 + 2 * npad            ! include padding
-        max_ntargs_in = 366 + 2 * npad * 31         ! include padding
+        max_ntargs_in = 366 + 2 * (npad + 1) * 31         ! include padding
     case ("none", "None")
         interp_method = "none"
     end select
@@ -454,8 +455,16 @@ program cal_adjust
     end select
 
     ! get _FillValue
-    call check( nf90_get_att(ncid_in, varid_in, '_FillValue', vfill) )
-    write (*,'("_FillValue:", g14.6)') vfill
+    !call check( nf90_get_att(ncid_in, varid_in, '_FillValue', vfill) )
+    !write (*,'("_FillValue:", g14.6)') vfill
+    !write (*,'(a,f7.2)') "Read time: ", secnds(read_secs)
+    vfill = 0.0
+    fill_status = nf90_get_att(ncid_in, varid_in, '_FillValue', vfill)
+    if (fill_status.eq.nf90_noerr) then
+        print '("  _FillValue = ",g14.6)', vfill
+    else
+        print '(" no _FillValue")'
+    end if
     write (*,'(a,f7.2)') "Read time: ", secnds(read_secs)
 
     ! Step 7:  Get calendar-adjusted values

f90/main_programs/demo_02_adjust_1yr.f90

@@ -74,7 +74,7 @@ program demo_02_adjust_1yr
     ydh, ym_int)
 
 ! reaggregate pseudo-daily values to new calendar
-call day_to_rmon_ts(ny, ndays, rmonbeg_06, rmonend_06, nd, ydh, yfill, ym_adj_cal)
+call day_to_mon_ts(ny, ndays, rmonbeg_06, rmonend_06, nd, ydh, yfill, ym_adj_cal)
 
 write (out_unit,'(/a)') "Paleo month-length adjusted values on 6ka calendar:"
 write (out_unit,'("   ym_adj_cal: ",12f9.3)') ym_adj_cal

f90/main_programs/demo_03_adjust_TraCE_ts.f90

@@ -69,7 +69,7 @@ program demo_03_adjust_TraCE_ts
 !dailyfile = "TraCE_c30r40_tas_land_monlen0ka_Jan-Dec_ts_daily.csv" ! uncomment to see daily data
 outfile = "TraCE_c30r40_tas_land_monlenadj_Jan-Dec.csv"
 
-monlenpath = "/Projects/Calendar/PaleoCalAdjust/data/month_lengths/"
+monlenpath = "/Projects/Calendar/PaleoCalAdjust/data/month_lengths/files/"
 monfile = "tr21_cal_noleap_rmonlen.csv"
 
 ! open output files and write headers
@@ -164,7 +164,7 @@ program demo_03_adjust_TraCE_ts
 ! calendar adjustment
 write (*,*) "Calendar adjustment..."
 
-call day_to_rmon_ts(ny, ndays, rmonbeg, rmonend, ndtot, ydh, yfill, ym_adj)
+call day_to_mon_ts(ny, ndays, rmonbeg, rmonend, ndtot, ydh, yfill, ym_adj)
 
 do n=1,ny
     write (3,'(f8.3,12(", ",g14.6))') yrbp(n),(ym_adj((n-1)*nm + m), m=1,nm)

f90/modules/CMIP_netCDF_subs.f90

@@ -12,7 +12,7 @@ module CMIP_netCDF_subs
 
     implicit none
 
-    integer, parameter          :: maxdims = 5, maxvars = 10, maxatts = 100, maxdimsize = 1440
+    integer, parameter          :: maxdims = 16, maxvars = 64, maxatts = 100, maxdimsize = 1440
 
     ! dimensions
     integer(4)                  :: ncid_in, ncid_out, ndim, nvar, natt, nglatt, unlimid, ncformat
@@ -22,7 +22,8 @@ module CMIP_netCDF_subs
     ! variables
     integer(4)                  :: varid_in, varid_out
     integer(4)                  :: xtype(maxvars), nvardims(maxvars), vardimids(maxvars,nf90_max_var_dims), nvaratts(maxvars)
-    character(nf90_max_name)    :: varname(maxvars)
+    !character(nf90_max_name)    :: varname(maxvars)
+    character(2048)    :: varname(maxvars)
     integer(4)                  :: dataid(maxvars), timedimid, rmonlenid, rmonmidid, rmonbegid, rmonendid
 
     ! input variable dimensions
@@ -116,7 +117,7 @@ subroutine copy_dims_and_glatts(ncid_in, ncid_out, addattname, addatt, nt, time,
     do i = 1,ndim
         dimid(i) = i
         call check( nf90_inquire_dimension(ncid_in, dimid(i), dimname(i), dimlen(i)) )
-        if (print_out) print '("   in:  dimid, dimlen, dimname = " ,2i7,1x,a)', i, dimlen(i), trim(dimname(i))
+        if (print_out) print '(" dimid, dimlen, dimname = " ,2i7,1x,a)', i, dimlen(i), trim(dimname(i))
 
         if (i .eq. unlimid) then
             call check( nf90_def_dim(ncid_out, dimname(i), nf90_unlimited, dimid(i)) )
@@ -133,13 +134,14 @@ subroutine copy_dims_and_glatts(ncid_in, ncid_out, addattname, addatt, nt, time,
     end do
 
     ! define dimension variables
-    if (print_out) print '("Define dimension variables:")'
+    if (print_out) print '("Finding dimension variables:")'
     varid_out = 0
     do i = 1,nvar
 
         varid_in = i
         call check( nf90_inquire_variable(ncid_in, varid_in, varname(i), xtype(i), ndims=nvardims(i), natts=nvaratts(i)) )
-        if (print_out) print '(" i, xtype, nvardims, nvaratts = ", 4i6, 1x, a)', i,xtype(i),nvardims(i),nvaratts(i),trim(varname(i))
+        if (print_out) print '(" i, xtype, nvardims, nvaratts, varname = ", 4i6, 1x, a)', &
+            i,xtype(i),nvardims(i),nvaratts(i),trim(varname(i))
         call check( nf90_inquire_variable(ncid_in, varid_in, dimids=vardimids(i,:nvardims(i))) )
         if (print_out) print '("    vardimids = ", 6i6)', vardimids(i,:nvardims(i))
 
@@ -210,12 +212,15 @@ subroutine copy_dims_and_glatts(ncid_in, ncid_out, addattname, addatt, nt, time,
     call check( nf90_enddef(ncid_out) )
 
     ! copy dimension variable values, replacing time and time_bnds
+
     if (print_out) print '("Copy dimension variables")'
     varid_out = 0
     do i = 1,nvar
 
         varid_in = i
+        if (print_out) print '(" i, varname(i): ", i4, 1x, a)', i, trim(varname(i))
         select case (varname(i))
+        ! add other dimension variables as needed
         case ('lon', 'lon_bnds', 'lat', 'lat_bnds', 'time', 'time_bnds', 'climatology_bnds', 'height', &
                 'plev', 'j', 'i', 'vertices', 'lon_vertices', 'lat_vertices', 'lev', 'lev_bnds')
 
@@ -389,6 +394,8 @@ subroutine define_outvar(ncid_in, ncid_out, varinname, varid_out, varoutname, ad
 
     ! find input variable number
     do i=1,nvar
+        if (print_out) print '(" i, nvardim(i), varname(i): ", 2i4, 1x, a)',  &
+            i, nvardims(i), trim(varname(i))
         if (trim(varname(i)) .eq. trim(varinname)) exit
     end do
     varid_in = i

f90/modules/pseudo_daily_interp_subs.f90

@@ -259,6 +259,106 @@ subroutine day_to_rmon_ts(ny, ndays, imonbeg, imonend, rmonbeg, rmonend, ndtot,
 
 end subroutine day_to_rmon_ts
 
+subroutine day_to_mon_ts(ny,ndays,rmonbeg,rmonend,ndtot,xd,xfill,xm_adj)
+! aggregation of pseudo- or actual daily data to months using a paleo calendar
+! used in demonstration of calendar adjustment of TraCE data
+
+    implicit none
+
+    integer(4), parameter       :: nm=12, nd=366
+    integer(4), intent(in)      :: ny, ndtot        ! number of years, total number of days
+    integer(4), intent(in)      :: ndays(ny)        ! number of days in each year
+    real(8), intent(in)         :: rmonbeg(ny,nm), rmonend(ny,nm)   ! beginning and ending days of each month
+    real(8), intent(in)         :: xd(ndtot)        ! daily values
+    real(8), intent(in)         :: xfill            ! _FillValue
+    real(8), intent(out)        :: xm_adj(ny*nm)    ! (aggregated) average monthly values
+
+    ! variables used to calculate monthly means
+    integer(4)              :: ibegday, iendday             ! beginning day and ending day of each year
+    integer(4)              :: ibeg(nm), iend(nm)           ! beginning and ending (integer) day of each month
+    integer(4)              :: ndays_in_month(nm)           ! integer number of days in month
+    real(8)                 :: xdx(-29:nd+30)               ! daily data for current year, padded by data from adjacent years
+    real(8)                 :: wgt(-29:nd+30), wsum         ! weights (for interpolating over fractional days)
+    integer(4)              :: nfill                        ! number of days with fill values
+
+    integer(4)              :: n, m, i, nn
+
+    logical                 :: debug_write_cal_effects = .false.
+
+    debug_write_cal_effects = .false.
+    if (debug_write_cal_effects) write (10,'(a)') "day_to_mon"
+    if (debug_write_cal_effects) write (10,*) ny, ndtot
+
+    ! loop over years, collecting daily data for each year, and getting monthly means
+    iendday = 0; nn = 0
+    xm_adj=0.0d0
+    do n=1,ny
+        if (debug_write_cal_effects) write (10,'("n, ndays:", 2i5)') n,ndays(n)
+        ibegday = iendday + 1
+        iendday = ibegday + ndays(n) - 1
+        if (debug_write_cal_effects) write (10,*) n,ibegday,iendday
+
+        if (ny .eq. 1) then       ! single-year Aclim data  
+            ! wrap the input daily data
+            xdx(-29:0)=xd(ndays(n)-30+1:ndays(n))
+            xdx(1:ndays(n))=xd(1:ndays(n))
+            xdx(ndays(n)+1:ndays(n)+30)=xd(1:30)
+        else 
+            ! copy current year into xdx
+            xdx(1:ndays(n)) = xd(ibegday:iendday)
+            ! pad beginning and end of xdx
+            if (n .eq. 1) then
+                xdx(-29:0) = xd(ndays(n)-30+1:ndays(n))
+                xdx(ndays(n)+1:ndays(n)+30) = xd(iendday+1:iendday+30)
+            elseif (n .eq. ny) then
+                xdx(-29:0) = xd(ibegday-30:ibegday-1)
+                xdx(ndays(n)+1:ndays(n)+30) = xd(ibegday+1:ibegday+30)
+            else
+                xdx(-29:0) = xd(ibegday-30:ibegday-1)
+                xdx(ndays(n)+1:ndays(n)+30) = xd(iendday+1:iendday+30)
+            end if
+        end if
+
+        ! integer beginning and end of each month, and number of days in each month
+        ! ndays_in_month should be equal to the integer month length + 1
+        ibeg=ceiling(rmonbeg(n,:)); iend=ceiling(rmonend(n,:)); ndays_in_month=(iend-ibeg+1)
+        if (debug_write_cal_effects) write (10,'("rmonbeg:  ",12f14.8)') rmonbeg(n,:)
+        if (debug_write_cal_effects) write (10,'("rmonend:  ",12f14.8)') rmonend(n,:)
+        if (debug_write_cal_effects) write (10,'("ibeg:  ",12i4)') ibeg
+        if (debug_write_cal_effects) write (10,'("iend:  ",12i4)') iend
+        if (debug_write_cal_effects) write (10,'("ndays: ",13i4)') ndays_in_month, sum(ndays_in_month)
+
+        ! monthly means
+        do m=1,nm
+            nn = nn + 1
+            nfill = 0; wgt=1.0d0; wsum=0.0d0
+            wgt(ibeg(m))=abs(rmonbeg(n,m)-dble(ibeg(m)))
+            wgt(iend(m))=abs(rmonend(n,m)-dble(iend(m)-1))
+            do i=ibeg(m),iend(m)
+                if (xdx(i) .ne. xfill) then
+                    xm_adj(nn)=xm_adj(nn)+xdx(i)*wgt(i)
+                    wsum=wsum+wgt(i)
+                    if (debug_write_cal_effects) &
+                    write (10,'("m, i, xd(i), xm_adj(nn), wgt(i), wsum: ",2i4,2f12.6,2f12.6)') &
+                        m,i,xdx(i),xm_adj(nn),wgt(i),wsum
+                else
+                    nfill = nfill + 1
+                end if
+            end do
+            if (wsum .ne. 0.0d0 .and. nfill .eq. 0) then
+                xm_adj(nn)=xm_adj(nn)/wsum
+            else
+                xm_adj(nn)=xfill
+            end if
+            if (debug_write_cal_effects) write(10,'("m, xm_adj(nn): ",i4,2f12.6)') &
+                m,xm_adj(nn),sngl(xm_adj(nn))
+        end do
+    end do 
+    if (debug_write_cal_effects) write (10,'(12g14.6)') xm_adj
+
+end subroutine day_to_mon_ts
+
+
 subroutine dayinterp(nm,nd,monlen,ym,yd)
 ! Interpolate pseudo-daily values of monthly data.  Not mean-preserving.