Allow multiple forcing file inputs per coupling field. Coupling field…

… can be calculated on the fly. COSIMA/access-om2#242

atm/src/atm.F90

@@ -79,6 +79,7 @@ program atm
     ! tells us how man atm-to-ice fields there are.
     call forcing_config%init(forcing_file, accessom2%logger, &
                              num_atm_to_ice_fields)
+    print*, 'atm, num_atm_to_ice_fields: ', num_atm_to_ice_fields
 
     ! Initialise forcing fields, involves reading details of each from
     ! config file and from netcdf files on disk, and allocating
@@ -163,6 +164,7 @@ program atm
         do i=1, num_atm_to_ice_fields
             ri = to_runoff_map(i)
 
+            print*, 'forcing_fields(i)%dt: ', forcing_fields(i)%dt
             if (mod(cur_runtime_in_seconds, forcing_fields(i)%dt) == 0) then
                 call field_read_timer%start()
                 call forcing_fields(i)%update(accessom2%get_cur_forcing_date(), &

libforcing/forcing.json

@@ -0,0 +1,21 @@
+
+{
+  "description": "JRA55-do V1.3 RYF 1990-91 forcing",
+  "inputs": [
+    {
+      "filename": "/g/data/ua8/JRA55-do/RYF/v1-3/RYF.rsds.1990_1991.nc",
+      "fieldname": "rsds",
+      "cname": "swfld_ai",
+      "pertubations": [
+
+{
+    "type": "offset",
+    "dimension": "constant",
+    "value": "84",
+    "calendar": "forcing"
+}
+
+      ]
+    }
+  ]
+}

libforcing/src/forcing_config.F90

@@ -66,6 +66,7 @@ subroutine forcing_config_init(self, config, loggerin, num_fields)
 
     self%num_inputs = self%core%count(inputs)
     num_fields = self%num_inputs
+    print*, 'forcing_config_init num_fields: ', num_fields
     self%logger => loggerin
 
 endsubroutine forcing_config_init
@@ -83,6 +84,7 @@ subroutine forcing_config_parse(self, fields, start_date, &
     type(json_value), pointer :: field_jv_ptr
     integer :: i, dt
     character(len=9) :: calendar_str
+    character(kind=CK, len=:), allocatable :: product_name
 
     type(json_value), pointer :: root, inputs
     logical :: found, is_land_field
@@ -94,14 +96,17 @@ subroutine forcing_config_parse(self, fields, start_date, &
     call self%core%get_child(root, "inputs", inputs, found)
     call assert(found, "No inputs found in forcing config.")
 
+    call self%core%get(root, "forcing_product_name", product_name, found)
+    call assert(found, "No forcing_product_name found in forcing config.")
+
     min_dt = HUGE(1)
     calendar = ''
     num_land_fields = 0
     do i=1, self%num_inputs
         call self%core%get_child(inputs, i, field_jv_ptr, found)
         call assert(found, "No inputs found in forcing config.")
         call self%parse_field(field_jv_ptr, fields(i), start_date, &
-                              dt, calendar_str, is_land_field)
+                              product_name, dt, calendar_str, is_land_field)
         if (dt < min_dt) then
             min_dt = dt
         endif
@@ -120,13 +125,14 @@ subroutine forcing_config_parse(self, fields, start_date, &
 
 
 subroutine forcing_config_parse_field(self, field_jv_ptr, field_ptr, &
-                                      start_date, dt, forcing_calendar, &
+                                      start_date, product_name, dt, forcing_calendar, &
                                       is_land_field)
 
     class(forcing_config), intent(inout) :: self
     type(json_value), pointer :: field_jv_ptr
     type(forcing_field) :: field_ptr
     type(datetime), intent(in) :: start_date
+    character(len=*), intent(in) :: product_name
     integer, intent(out) :: dt
     character(len=9), intent(out) :: forcing_calendar
     logical, intent(out) :: is_land_field
@@ -138,21 +144,48 @@ subroutine forcing_config_parse_field(self, field_jv_ptr, field_ptr, &
     character(kind=CK, len=:), allocatable :: dimension_type
     character(kind=CK, len=:), allocatable :: perturbation_calendar
 
+    character(len=256), dimension(:), allocatable :: fieldname_list, filename_list
+
+    integer :: num_fieldnames, num_filenames
     integer :: perturbation_constant_value
     logical :: found, domain_found
     integer :: num_perturbations, num_fields
     integer :: i, j
 
+    type(json_value), pointer :: fieldname_jv_list, filename_jv_list
+    type(json_value), pointer :: fieldname_jv_ptr, filename_jv_ptr
     type(json_value), pointer :: perturbation_jv_ptr
     type(json_value), pointer :: dimension_jv_ptr, value_jv_ptr
     type(json_value), pointer :: perturbation_list, dimension_list
     type(json_value), pointer :: value_list
 
-    call self%core%get(field_jv_ptr, "fieldname", fieldname, found)
-    call assert(found, "Entry 'fieldname' not found in forcing config.")
-
-    call self%core%get(field_jv_ptr, "filename", filename, found)
-    call assert(found, "Entry 'filename' not found in forcing config.")
+    ! Allow there to be multiple 
+    call self%core%get_child(field_jv_ptr, "fieldnames", fieldname_jv_list, found)
+    call assert(found, "Entry 'fieldnames' not found in forcing config.")
+    num_fieldnames = self%core%count(fieldname_jv_list)
+
+    call self%core%get_child(field_jv_ptr, "filenames", filename_jv_list, found)
+    call assert(found, "Entry 'filenames' not found in forcing config.")
+    num_filenames = self%core%count(filename_jv_list)
+
+    call assert(num_fieldnames == num_filenames, &
+                'Number of fieldnames does not match number of filenames')
+    allocate(fieldname_list(num_fieldnames))
+    allocate(filename_list(num_filenames))
+
+    do i=1, num_filenames
+        call self%core%get_child(fieldname_jv_list, i, &
+                                 fieldname_jv_ptr, found)
+        call assert(found, "Expected to find fieldname entry.")
+        call self%core%get(fieldname_jv_ptr, value=fieldname)
+        fieldname_list(i) = trim(fieldname)
+
+        call self%core%get_child(filename_jv_list, i, &
+                                 filename_jv_ptr, found)
+        call assert(found, "Expected to find filename entry.")
+        call self%core%get(filename_jv_ptr, value=filename)
+        filename_list(i) = trim(filename)
+    enddo
 
     call self%core%get(field_jv_ptr, "cname", cname, found)
     call assert(found, "Entry 'cname' not found in forcing config.")
@@ -170,8 +203,8 @@ subroutine forcing_config_parse_field(self, field_jv_ptr, field_ptr, &
         is_land_field = .true.
     endif
 
-    call field_ptr%init(fieldname, filename, cname, domain_str, start_date, &
-                        self%logger, dt, forcing_calendar)
+    call field_ptr%init(fieldname_list, filename_list, cname, domain_str, start_date, &
+                        product_name, self%logger, dt, forcing_calendar)
 
     call self%core%get_child(field_jv_ptr, "perturbations", perturbation_list, found)
     if (.not. found) then

libforcing/src/forcing_field.F90

@@ -19,15 +19,16 @@ module forcing_field_mod
 integer, parameter, public :: FORCING_FIELD_DOMAIN_LAND = 20
 
 type, public :: forcing_field
-    character(len=64) :: name
+    character(len=64), dimension(:), allocatable :: names
     character(len=64) :: coupling_name
-    character(len=1024) :: filename_template
+    character(len=1024), dimension(:), allocatable :: filename_templates
     integer :: domain
 
     integer :: dt
     character(len=9) :: calendar
+    character(len=64) :: product_name
 
-    type(ncvar_type) :: ncvar
+    type(ncvar_type), dimension(:), allocatable :: ncvars
     real, dimension(:, :), allocatable :: data_array
     type(forcing_perturbation_type), dimension(:), allocatable :: perturbations
     type(forcing_perturbation_type), dimension(:), allocatable :: &
@@ -37,30 +38,40 @@ module forcing_field_mod
 contains
     procedure, pass(self), public :: init => forcing_field_init
     procedure, pass(self), public :: update => forcing_field_update
+    procedure, pass(self), private :: calculate => forcing_field_calculate
     procedure, pass(self), private :: apply_perturbations => &
                 forcing_field_apply_perturbations
     procedure, pass(self), public :: get_shape
 endtype forcing_field
 
 contains
 
-subroutine forcing_field_init(self, name, filename_template, cname, domain, &
-                              start_date, loggerin, dt, calendar)
+subroutine forcing_field_init(self, name_list, filename_template_list, cname, domain, &
+                              start_date, product_name, loggerin, dt, calendar)
     class(forcing_field), intent(inout) :: self
-    character(len=*), intent(in) :: name
+    character(len=*), dimension(:), intent(in) :: name_list
+    character(len=*), dimension(:), intent(in) :: filename_template_list
     character(len=*), intent(in) :: cname
-    character(len=*), intent(in) :: filename_template
     character(len=*), intent(in) :: domain
     type(datetime), intent(in) :: start_date
+    character(len=*), intent(in) :: product_name
     type(logger_type), target, intent(in) :: loggerin
     integer, intent(out) :: dt
     character(len=9), intent(out) :: calendar
 
     character(len=1024) :: filename
+    integer :: num_file_inputs, i
 
-    self%name = name
+    num_file_inputs = size(name_list)
+    print*, 'num_file_inputs: ', num_file_inputs
+
+    allocate(self%names(num_file_inputs))
+    allocate(self%filename_templates(num_file_inputs))
+    allocate(self%ncvars(num_file_inputs))
+
+    self%names(:) = name_list(:)
+    self%filename_templates(:) = filename_template_list(:)
     self%coupling_name = cname
-    self%filename_template = filename_template
     if (domain == 'atmosphere') then
         self%domain = FORCING_FIELD_DOMAIN_ATMOSPHERE
     else
@@ -69,37 +80,64 @@ subroutine forcing_field_init(self, name, filename_template, cname, domain, &
         self%domain = FORCING_FIELD_DOMAIN_LAND
     endif
 
+    self%product_name = trim(product_name)
     self%logger => loggerin
 
-    filename = filename_for_year(self%filename_template, start_date%getYear())
-    call self%ncvar%init(self%name, filename)
-    allocate(self%data_array(self%ncvar%nx, self%ncvar%ny))
+    do i=1, num_file_inputs
+        filename = filename_for_year(self%filename_templates(i), start_date%getYear())
+        call self%ncvars(i)%init(self%names(i), filename)
+    enddo
+
+    ! Check that the metadata is the same on all input filenames
+    if (num_file_inputs > 1) then
+        do i=2,num_file_inputs
+            call assert(self%ncvars(i)%dt == self%ncvars(1)%dt, &
+             'File metadata (dt) for multi-field forcing does not match.')
+            call assert(self%ncvars(i)%calendar == self%ncvars(1)%calendar, &
+             'File metadata (calendar) for multi-field forcing does not match.')
+            call assert(self%ncvars(i)%nx == self%ncvars(1)%nx, &
+             'File metadata (nx) for multi-field forcing does not match.')
+            call assert(self%ncvars(i)%ny == self%ncvars(1)%ny, &
+             'File metadata (ny) for multi-field forcing does not match.')
+        enddo
+    endif
+
+    allocate(self%data_array(self%ncvars(1)%nx, self%ncvars(1)%ny))
     self%data_array(:, :) = HUGE(1.0)
-    self%dt = self%ncvar%dt
+
+    self%dt = self%ncvars(1)%dt
+    self%calendar = self%ncvars(1)%calendar
     dt = self%dt
-    self%calendar = self%ncvar%calendar
     calendar = self%calendar
 
+    print*, 'dt, calendar, nx, ny:', dt, calendar, self%ncvars(1)%nx, self%ncvars(1)%ny
+
 endsubroutine forcing_field_init
 
+
 subroutine forcing_field_update(self, forcing_date, experiment_date)
     class(forcing_field), intent(inout) :: self
     type(datetime), intent(in) :: forcing_date, experiment_date
 
     character(len=1024) :: filename
     character(len=10) :: int_str
-    integer :: indx
+    integer :: indx, test_indx
+    integer :: num_file_inputs, i
 
-    filename = filename_for_year(self%filename_template, forcing_date%getYear())
-    call assert(trim(filename) /= '', "File not found: "//filename)
-    if (trim(filename) /= trim(self%ncvar%filename)) then
-        call self%ncvar%refresh(filename)
-    endif
+    num_file_inputs = size(self%ncvars)
 
-    indx = self%ncvar%get_index_for_datetime(forcing_date)
+    do i=1, num_file_inputs
+        filename = filename_for_year(self%filename_templates(i), forcing_date%getYear())
+        call assert(trim(filename) /= '', "File not found: "//filename)
+        if (trim(filename) /= trim(self%ncvars(i)%filename)) then
+            call self%ncvars(i)%refresh(filename)
+        endif
+    enddo
+
+    indx = self%ncvars(1)%get_index_for_datetime(forcing_date)
     if (indx == -1) then
         ! Search from the beginning before failing
-        indx = self%ncvar%get_index_for_datetime(forcing_date, .true.)
+        indx = self%ncvars(1)%get_index_for_datetime(forcing_date, .true.)
         call self%logger%write(LOG_DEBUG, &
                  'forcing_field_update: long forcing index search')
     endif
@@ -112,13 +150,51 @@ subroutine forcing_field_update(self, forcing_date, experiment_date)
     call self%logger%write(LOG_DEBUG, '{ "forcing_field_update-index" : '// &
                                        trim(int_str)//' }')
 
-    call self%ncvar%read_data(indx, self%data_array)
+    ! If there are other ncvars then we need to check that they have
+    ! give us the same index. This should be fast in the no-error case.
+    if (num_file_inputs > 1) then
+        do i=2,num_file_inputs
+            test_indx = self%ncvars(i)%get_index_for_datetime(forcing_date, &
+                                                             guess=indx)
+            call assert(test_indx == indx, &
+                    'Time indices on multi-file forcing do not match.')
+        enddo
+    endif
 
+    call self%calculate(indx, self%data_array)
     call self%apply_perturbations(forcing_date, experiment_date)
 
 end subroutine forcing_field_update
 
 
+subroutine forcing_field_calculate(self, file_index, result_array)
+
+    class(forcing_field), intent(inout) :: self
+    integer, intent(in) :: file_index
+    real, dimension(:, :), intent(inout) :: result_array
+
+    if (trim(self%product_name) == 'ERA5') then
+        if (trim(self%coupling_name) == 'rain_ai') then
+            ! Rain is calculated as total precipitation - snow
+            ! FIXME: do calculation with field name checks
+            call self%ncvars(1)%read_data(file_index, result_array)
+        elseif (trim(self%coupling_name) == 'qair_ai') then
+            ! Humidity is calculated using surface temperature and dew point
+            ! FIXME: do calculation with field name checks
+            call self%ncvars(1)%read_data(file_index, result_array)
+        else
+            call self%ncvars(1)%read_data(file_index, result_array)
+        endif
+    elseif (trim(self%product_name) == 'JRA55-do') then
+        call self%ncvars(1)%read_data(file_index, result_array)
+    else
+        call assert(.false., &
+            'Unsupported forcing_product_name. Valid names are: "ERA5", "JRA55-do"')
+    endif
+
+endsubroutine forcing_field_calculate
+
+
 ! Iterate through perturbations and apply to base field in self%data
 subroutine forcing_field_apply_perturbations(self, forcing_date, experiment_date)
     class(forcing_field), intent(inout) :: self
@@ -135,9 +211,9 @@ subroutine forcing_field_apply_perturbations(self, forcing_date, experiment_date
         return
     endif
 
-    allocate(pertub_array(self%ncvar%nx, self%ncvar%ny))
-    allocate(tmp(self%ncvar%nx, self%ncvar%ny))
-    allocate(another_tmp(self%ncvar%nx, self%ncvar%ny))
+    allocate(pertub_array(self%ncvars(1)%nx, self%ncvars(1)%ny))
+    allocate(tmp(self%ncvars(1)%nx, self%ncvars(1)%ny))
+    allocate(another_tmp(self%ncvars(1)%nx, self%ncvars(1)%ny))
     pertub_array(:, :) = 1.0
 
     ! First iterate over all of the scaling fields
@@ -240,6 +316,9 @@ function get_shape(self)
     class(forcing_field), intent(in) :: self
     integer, dimension(2) :: get_shape
 
+    print*, 'cname: ', self%coupling_name
+    print*, 'shape: ', shape(self%data_array)
+
     get_shape = shape(self%data_array)
 endfunction
 

libforcing/src/ncvar.F90

@@ -189,16 +189,21 @@ subroutine get_start_date_and_calendar(self, time_varid, start_date, calendar)
 endsubroutine get_start_date_and_calendar
 
 !> Return the time index of a particular date.
-function get_index_for_datetime(self, target_date, from_beginning)
+function get_index_for_datetime(self, target_date, from_beginning, guess)
     class(ncvar), intent(inout) :: self
     type(datetime), intent(in) :: target_date
     logical, optional, intent(in) :: from_beginning
+    integer, optional, intent(in) :: guess
 
     integer :: i, get_index_for_datetime
     integer :: days, seconds
 
     type(timedelta) :: td, td_before, td_after
 
+    if (present(guess)) then
+        self%idx_guess = guess
+    endif
+
     if (present(from_beginning)) then
         if (from_beginning) then
             self%idx_guess = 1