Merge pull request #43 from COSIMA/feat-155-JRA55-v1p4

Feat 155 jra55 v1p4

.gitignore

@@ -2,13 +2,15 @@
 *.swp
 __pycache__
 *.pyc
+.swo
 .cache
 build*
 building-errors.log
 debug.root.01
 debug.01.000000
 debug.02.000000
 debug.03.000000
+debug.root.*
 nout.000000
 a2i.nc
 accessom2_restart.nml

.gitmodules

@@ -0,0 +1,9 @@
+[submodule "json-fortran"]
+	path = json-fortran
+	url = https://github.com/jacobwilliams/json-fortran.git
+[submodule "datetime-fortran"]
+	path = datetime-fortran
+	url = https://github.com/wavebitscientific/datetime-fortran.git
+[submodule "oasis3-mct"]
+	path = oasis3-mct
+	url = https://github.com/COSIMA/oasis3-mct.git

CMakeLists.txt

@@ -61,8 +61,8 @@ include(ExternalProject)
 # We use a forked version because of this behaviour:
 # https://github.com/wavebitscientific/datetime-fortran/issues/51
 ExternalProject_Add(datetime-fortran
-    GIT_REPOSITORY   https://github.com/nicjhan/datetime-fortran.git
-    GIT_TAG          3e7755c585cdc62ed262e91d57e7d54f82046087
+    DOWNLOAD_COMMAND ""
+    SOURCE_DIR ${CMAKE_SOURCE_DIR}/datetime-fortran
     INSTALL_COMMAND cmake -E echo "Skipping datetime-fortran install step."
 )
 ExternalProject_Get_property(datetime-fortran BINARY_DIR)
@@ -73,8 +73,8 @@ set_property(TARGET datetime PROPERTY IMPORTED_LOCATION ${BINARY_DIR}/lib/libdat
 
 # json-fortran external project
 ExternalProject_Add(json-fortran
-    GIT_REPOSITORY    https://github.com/jacobwilliams/json-fortran.git
-    GIT_TAG          37de3ade555e84779ebf8a26c23e86a3ae268390
+    DOWNLOAD_COMMAND ""
+    SOURCE_DIR ${CMAKE_SOURCE_DIR}/json-fortran
     INSTALL_COMMAND cmake -E echo "Skipping json-fortran install step."
 )
 ExternalProject_Get_property(json-fortran BINARY_DIR)
@@ -86,8 +86,8 @@ set_property(TARGET jsonfortran PROPERTY IMPORTED_LOCATION ${BINARY_DIR}/lib/lib
 # oasis external project
 # json-fortran external project
 ExternalProject_Add(oasis3-mct
-    GIT_REPOSITORY    https://github.com/COSIMA/oasis3-mct.git
-    GIT_TAG           87a873aa7b545a6ee46a42f811f5c7a5dfc81a01
+    DOWNLOAD_COMMAND ""
+    SOURCE_DIR ${CMAKE_SOURCE_DIR}/oasis3-mct
     CONFIGURE_COMMAND ""
     BUILD_IN_SOURCE 1
     BUILD_COMMAND make ${PLATFORM}

README.md

@@ -60,13 +60,20 @@ cd libaccessom2
 
 # Run tests on Gadi (NCI)
 
-First do build as above. Then:
+First do build as above. Then to get some computer resources:
+
+```{bash}
+qsub -I -P x77 -q normal -lncpus=4 -lmem=16Gb -lwalltime=3:00:00 -lstorage=gdata/ua8+gdata/qv56+gdata/hh5+gdata/ik11
+
+/g/data1b/qv56/
+```
+
 
 ```{bash}
 export LIBACCESSOM2_DIR=$(pwd)
 module load openmpi
 cd tests/
-./copy_test_data_from_gadi.sh
+./copy_test_data.sh
 cd JRA55_IAF
 rm -rf log ; mkdir log ; rm -f accessom2_restart_datetime.nml ; cp ../test_data/i2o.nc ./ ; cp ../test_data/o2i.nc ./
 mpirun -np 1 $LIBACCESSOM2_DIR/build/bin/yatm.exe : -np 1 $LIBACCESSOM2_DIR/build/bin/ice_stub.exe : -np 1 $LIBACCESSOM2_DIR/build/bin/ocean_stub.exe

atm/src/atm.F90

@@ -2,7 +2,7 @@ program atm
 
     use mod_oasis, only : OASIS_IN, OASIS_OUT
     use forcing_mod, only : forcing_type => forcing
-    use field_mod, only : field_type => field
+    use field_mod, only : field_type => field, FIELD_DOMAIN_LAND
     use coupler_mod, only : coupler_type => coupler
     use error_handler, only : assert
     use ice_grid_proxy_mod, only : ice_grid_type => ice_grid_proxy
@@ -22,13 +22,16 @@ program atm
     type(ice_grid_type) :: ice_grid
     type(runoff_type) :: runoff
     type(field_type), dimension(:), allocatable :: fields
-    type(field_type) :: runoff_field
+    integer, dimension(:), allocatable :: to_runoff_map
+    ! Liquid (river) and solid (iceberg) runoff
+    type(field_type), dimension(:), allocatable :: runoff_fields
     character(len=MAX_FILE_NAME_LEN) :: forcing_file, accessom2_config_dir
     character(len=9) :: calendar
     integer, dimension(2) :: ice_shape
-    integer :: i, err, tmp_unit
+    integer :: i, ri, err, tmp_unit
     logical :: file_exists
     integer :: num_atm_to_ice_fields, dt, cur_runtime_in_seconds
+    integer :: num_land_fields
 
     type(simple_timer_type) :: field_read_timer, ice_wait_timer
     type(simple_timer_type) :: init_runoff_timer, remap_runoff_timer
@@ -51,11 +54,18 @@ program atm
     call accessom2%init('matmxx', config_dir=trim(accessom2_config_dir))
     call accessom2%print_version_info()
 
-    ! Initialise forcing object and fields, involves reading details of each
-    ! field from disk.
+    ! Initialise forcing object, this reads config and
+    ! tells us how man atm-to-ice fields there are.
     call forcing%init(forcing_file, accessom2%logger, num_atm_to_ice_fields)
+
+    ! Initialise forcing fields, involves reading details of each from disk,
+    ! and allocating necessary memory.
     allocate(fields(num_atm_to_ice_fields))
-    call forcing%init_fields(fields, accessom2%get_cur_forcing_date(), dt, calendar)
+    call forcing%init_fields(fields, accessom2%get_cur_forcing_date(), &
+                             dt, calendar, num_land_fields)
+    ! Create intermediate fields for runoff,
+    ! these are a copy/variation of the forcing fields
+    allocate(runoff_fields(num_land_fields))
 
     ! Initialise the coupler.
     call coupler%init_begin('matmxx', accessom2%logger, &
@@ -68,9 +78,11 @@ program atm
     ! Synchronise accessom2 'state' (i.e. configuration) between all PEs of all models.
     call accessom2%sync_config(coupler)
 
-    ! Get information about the ice grid needed for runoff remapping.
+    ! Initialise ice grid proxy and get information about it,
+    ! this is needed for local remapping.
     call ice_grid%init(coupler%ice_root)
     call ice_grid%recv()
+    ice_shape = ice_grid%get_shape()
 
     ! Initialise timers
     call field_read_timer%init('field_read', accessom2%logger, &
@@ -88,21 +100,34 @@ program atm
     call init_runoff_timer%start()
     call runoff%init(ice_grid)
     call init_runoff_timer%stop()
-    ice_shape = ice_grid%get_shape()
 
-    ! Initialise OASIS3-MCT fields. Runoff done seperately for now.
+    ! Create a little map to go from atm_to_ice field indices to
+    ! runoff field indices, simplifies the code below
+    allocate(to_runoff_map(num_atm_to_ice_fields))
+    ri = 1
+    do i=1, num_atm_to_ice_fields
+        if (fields(i)%domain == FIELD_DOMAIN_LAND) then
+            to_runoff_map(i) = ri
+            ri = ri + 1
+        else
+            to_runoff_map(i) = 0
+        endif
+    enddo
+
+    ! Initialise coupling fields, runoff fields need special treatment.
     do i=1, num_atm_to_ice_fields
-        if (index(fields(i)%name, 'runof') /= 0) then
-            call assert(.not. allocated(runoff_field%data_array), &
-                        'Runoff already associated')
-            runoff_field%name = fields(i)%name
-            runoff_field%timestamp = fields(i)%timestamp
-            allocate(runoff_field%data_array(ice_shape(1), ice_shape(2)))
-            call coupler%init_field(runoff_field, OASIS_OUT)
+        if (to_runoff_map(i) /= 0) then
+            ri = to_runoff_map(i)
+            runoff_fields(ri)%name = fields(i)%name
+            runoff_fields(ri)%domain = fields(i)%domain
+            runoff_fields(ri)%timestamp = fields(i)%timestamp
+            allocate(runoff_fields(ri)%data_array(ice_shape(1), ice_shape(2)))
+            call coupler%init_field(runoff_fields(ri), OASIS_OUT)
         else
             call coupler%init_field(fields(i), OASIS_OUT)
         endif
     enddo
+
     ! Finish coupler initialisation. Tell oasis how long the run is and the
     ! coupling timesteps.
     call coupler%init_end(accessom2%get_total_runtime_in_seconds(), &
@@ -114,22 +139,24 @@ program atm
 
         ! Send each forcing field
         do i=1, num_atm_to_ice_fields
+            ri = to_runoff_map(i)
+
             if (mod(cur_runtime_in_seconds, fields(i)%dt) == 0) then
                 call field_read_timer%start()
                 call forcing%update_field(fields(i), &
                                           accessom2%get_cur_forcing_date())
                 call field_read_timer%stop()
-                if (index(fields(i)%name, 'runof') /= 0) then
+                if (ri /= 0) then
                     call remap_runoff_timer%start()
                     call runoff%remap(fields(i)%data_array, &
-                                      runoff_field%data_array, ice_grid%mask)
+                                      runoff_fields(ri)%data_array, ice_grid%mask)
                     call remap_runoff_timer%stop()
                 endif
             endif
 
             call coupler_put_timer%start()
-            if (index(fields(i)%name, 'runof') /= 0) then
-                call coupler%put(runoff_field, cur_runtime_in_seconds, err)
+            if (ri /= 0) then
+                call coupler%put(runoff_fields(ri), cur_runtime_in_seconds, err)
             else
                 call coupler%put(fields(i), cur_runtime_in_seconds, err)
             endif

atm/src/forcing.F90

@@ -11,7 +11,6 @@ module forcing_mod
 use, intrinsic :: iso_fortran_env, only : stderr=>error_unit
 
 implicit none
-
 private
 
 type, public :: forcing
@@ -57,26 +56,28 @@ subroutine forcing_init(self, config, logger, nfields)
 endsubroutine forcing_init
 
 !> Parse forcing file into a dictionary.
-subroutine forcing_init_fields(self, fields, forcing_date, min_dt, calendar)
+subroutine forcing_init_fields(self, fields, forcing_date, &
+                               min_dt, calendar, num_land_fields)
 
     class(forcing), intent(inout) :: self
     type(field_type), dimension(:), intent(inout) :: fields
     type(datetime), intent(in) :: forcing_date
-    integer, intent(out) :: min_dt
+    integer, intent(out) :: min_dt, num_land_fields
     character(len=9), intent(out) :: calendar
 
     type(json_value), pointer :: fp
     integer :: i
-    character(kind=CK, len=:), allocatable :: cname, fieldname
+    character(kind=CK, len=:), allocatable :: cname, fieldname, domain
     character(kind=CK, len=:), allocatable :: filename_template
     character(kind=CK, len=:), allocatable :: scaling_filename
     character(len=1024) :: filename
     character(len=9) :: calendar_str
-    logical :: found, scaling_found
+    logical :: found, scaling_found, domain_found
 
     min_dt = huge(min_dt)
     calendar_str = ''
 
+    num_land_fields = 0
     do i=1, size(fields)
         call self%core%get_child(self%inputs, i, fp, found)
         call assert(found, "Input not found in forcing config.")
@@ -92,15 +93,27 @@ subroutine forcing_init_fields(self, fields, forcing_date, min_dt, calendar)
         call self%core%get(fp, "cname", cname, found)
         call assert(found, "Entry 'cname' not found in forcing config.")
 
+        call self%core%get(fp, "domain", domain, domain_found)
+        if (domain_found) then
+            call assert(domain == "land" .or. domain == "atmosphere", &
+                        "forcing_init_fields: invalid domain value.")
+        else
+            domain = "atmosphere"
+        endif
+
+        if (domain == "land") then
+            num_land_fields = num_land_fields + 1
+        endif
+
         ! Get the shape of forcing fields
         filename = filename_for_year(filename_template, forcing_date%getYear())
         ! Initialise a new field object.
         if (scaling_found) then
             call fields(i)%init(cname, fieldname, filename_template, &
-                                filename, self%logger, scaling_filename)
+                                filename, domain, self%logger, scaling_filename)
         else
             call fields(i)%init(cname, fieldname, filename_template, &
-                                filename, self%logger)
+                                filename, domain, self%logger)
         endif
 
         if (fields(i)%dt < min_dt) then

ice_stub/src/ice.F90

@@ -89,7 +89,6 @@ program ice
     call ice_grid%init(trim(ice_grid_file), trim(ice_mask_file), resolution)
     call ice_grid%send()
 
-
     ! Set up coupling fields
     do i=1, size(from_atm_fields)
         from_atm_fields(i)%name = trim(from_atm_field_names(i))

libcouple/src/field.F90

@@ -7,11 +7,17 @@ module field_mod
 use logger_mod, only : logger_type => logger, LOG_DEBUG
 
 implicit none
-
 private
 
+! Fields can be either atmospheric or ocean. YATM will treat the differently,
+! see atm.F90 for details.
+integer, parameter, public :: FIELD_DOMAIN_NONE = 0
+integer, parameter, public :: FIELD_DOMAIN_ATMOSPHERE = 10
+integer, parameter, public :: FIELD_DOMAIN_LAND = 20
+
 type, public :: field
     character(len=64) :: name
+    integer :: domain
     character(len=1024) :: filename_template
     character(len=1024) :: scaling_filename
     type(datetime) :: timestamp
@@ -34,10 +40,11 @@ module field_mod
 contains
 
 subroutine field_init(self, name, ncname, filename_template, &
-                      filename, logger, scaling_filename)
+                      filename, domain, logger, scaling_filename)
     class(field), intent(inout) :: self
     character(len=*), intent(in) :: name, ncname
     character(len=*), intent(in) :: filename_template, filename
+    character(len=*), intent(in) :: domain
     type(logger_type), intent(in) :: logger
     character(len=*), intent(in), optional :: scaling_filename
 
@@ -59,6 +66,16 @@ subroutine field_init(self, name, ncname, filename_template, &
         self%scaling_data_array(:, :) = HUGE(1.0)
     endif
 
+    ! Set the field domain
+    self%domain = FIELD_DOMAIN_NONE
+    if (domain == 'atmosphere') then
+        self%domain = FIELD_DOMAIN_ATMOSPHERE
+    elseif (domain == 'land') then
+        self%domain = FIELD_DOMAIN_LAND
+    endif
+    call assert(self%domain /= FIELD_DOMAIN_NONE, &
+                'field_init: invalid field domain')
+
 end subroutine
 
 subroutine field_update_data(self, filename, forcing_date)

libcouple/src/ncvar.F90

@@ -67,7 +67,7 @@ subroutine ncvar_refresh(self, filename)
     call ncheck(nf90_open(trim(self%filename), NF90_NOWRITE, self%ncid), &
                 'ncvar opening '//trim(self%filename))
     call ncheck(nf90_inq_varid(self%ncid, trim(self%name), self%varid), &
-               'ncvar inquire: '//trim(self%name))
+               'ncvar inquire: '//trim(self%name)//' in '//trim(self%filename))
 
     ! Initialise dimensions
     call get_var_dims(self%ncid, self%varid, ndims, self%nx, self%ny, num_times)

ocean_stub/src/ocean.F90

@@ -78,23 +78,26 @@ program ocean
     call coupler%init_end(accessom2%get_total_runtime_in_seconds(), &
                           accessom2%get_coupling_field_timesteps())
 
-    do while (.not. accessom2%run_finished())
-        cur_runtime_in_seconds = accessom2%get_cur_runtime_in_seconds()
-
+    cur_runtime_in_seconds = accessom2%get_cur_runtime_in_seconds()
+    do
         ! Get fields from ice
         do i=1, size(in_fields)
             call coupler%get(in_fields(i), cur_runtime_in_seconds, err)
         enddo
 
         ! Do work, i.e. use the in_fields and populate the out_fields
 
+        call accessom2%progress_date(dt)
+        if (accessom2%run_finished()) then
+            exit
+        endif
+        cur_runtime_in_seconds = accessom2%get_cur_runtime_in_seconds()
+
         ! Send fields to ice
         do i=1, size(out_fields)
             out_fields(i)%data_array(:, :) = 0.0
             call coupler%put(out_fields(i), cur_runtime_in_seconds, err)
         enddo
-
-        call accessom2%progress_date(dt)
     enddo
 
     ! Write out restart.