Add serial driver gswp3 support and fix serial-MPI non-reproducibility (#544)

SeanBryan51 · web-flow · commit 284e1ef129eb · 2025-03-07T09:39:02.000+11:00
Currently the offline serial driver does not support the namelist option `cable_user%MetType = 'gswp3'` - the option is only supported in the MPI drivers. This PR adds support for the namelist option to the serial driver so we can use these configurations for testing the MPI implementation with parallel-IO (which will be based off the serial driver). This PR also fixes several bugs found in the MPI drivers such that we have bitwise reproducibility between serial and MPI for the [gswp3 configuration](https://github.com/SeanBryan51/crujra_accessN96_1h) for outputs and restarts (run for 1 year on a reduced domain). Fixes #388 ## Type of change - [x] Bug fix ## Checklist - [x] The new content is accessible and located in the appropriate section - [x] I have checked that links are valid and point to the intended content - [x] I have checked my code/text and corrected any misspellings ## Testing - [x] Are the changes non bitwise-compatible with the main branch because of a bug fix or a feature being newly implemented or improved? If yes, add the link to the modelevaluation.org analysis versus the main branch or equivalent results below this line. - Changes are bitwise compatible with main for site runs but not for spatial MPI runs due to the bug fixes in the MPI drivers: ``` 2025-03-05 17:37:00,335 - INFO - benchcab.benchcab.py:380 - Running comparison tasks... 2025-03-05 17:37:00,358 - INFO - benchcab.benchcab.py:381 - tasks: 168 (models: 2, sites: 42, science configurations: 4) 2025-03-05 17:40:02,866 - INFO - benchcab.benchcab.py:391 - 0 failed, 168 passed ```  ---- 📚 Documentation preview 📚: https://cable--544.org.readthedocs.build/en/544/
diff --git a/src/offline/cable_mpicommon.F90 b/src/offline/cable_mpicommon.F90
@@ -29,7 +29,7 @@ MODULE cable_mpicommon
 
   ! base number of input fields: must correspond to CALLS to 
   ! MPI_address (field ) in *_mpimaster/ *_mpiworker
-  INTEGER, PARAMETER :: nparam = 340
+  INTEGER, PARAMETER :: nparam = 341
    
   ! MPI: extra params sent only if nsoilparmnew is true
   INTEGER, PARAMETER :: nsoilnew = 1
@@ -103,7 +103,7 @@ MODULE cable_mpicommon
   ! MPI: number of fields included in restart_t type for data
   ! that is returned only for creating a restart file at the end of the run
   ! MPI: gol124: canopy%rwater removed when Bernard ported to CABLE_r491
-  INTEGER, PARAMETER :: nrestart = 15
+  INTEGER, PARAMETER :: nrestart = 17
   INTEGER, PARAMETER :: nsumcasaflux = 62
   INTEGER, PARAMETER :: nsumcasapool = 40
   INTEGER, PARAMETER :: nclimate = 30
diff --git a/src/offline/cable_mpimaster.F90 b/src/offline/cable_mpimaster.F90
@@ -500,6 +500,77 @@ SUBROUTINE mpidrv_master (comm, trunk_sumbal, dels, koffset, kend, PLUME, CRU)
           CALL alloc_cbm_var (imet, mp)
           CALL alloc_cbm_var (iveg, mp)
 
+          imet%ca(:) = met%ca(:)
+          imet%year(:) = met%year(:)
+          imet%moy(:) = met%moy(:)
+          imet%doy(:) = met%doy(:)
+          imet%hod(:) = met%hod(:)
+          imet%fsd(:,:) = met%fsd(:,:)
+          imet%ofsd(:) = met%ofsd(:)
+          imet%fld(:) = met%fld(:)
+          imet%precip(:) = met%precip(:)
+          imet%precip_sn(:) = met%precip_sn(:)
+          imet%tk(:) = met%tk(:)
+          imet%tvair(:) = met%tvair(:)
+          imet%tvrad(:) = met%tvrad(:)
+          imet%pmb(:) = met%pmb(:)
+          imet%ua(:) = met%ua(:)
+          imet%qv(:) = met%qv(:)
+          imet%qvair(:) = met%qvair(:)
+          imet%da(:) = met%da(:)
+          imet%dva(:) = met%dva(:)
+          imet%coszen(:) = met%coszen(:)
+          imet%Ndep(:) = met%Ndep(:)
+          imet%Pdep(:) = met%Pdep(:)
+
+          iveg%canst1(:) =  veg%canst1(:)
+          iveg%dleaf(:) = veg%dleaf(:)
+          iveg%ejmax(:) = veg%ejmax(:)
+          iveg%iveg(:) = veg%iveg(:)
+          iveg%iLU(:) = veg%iLU(:)
+          iveg%meth(:) = veg%meth(:)
+          iveg%frac4(:) = veg%frac4(:)
+          iveg%hc(:) = veg%hc(:)
+          iveg%vlai(:) = veg%vlai(:)
+          iveg%xalbnir(:) = veg%xalbnir(:)
+          iveg%rp20(:) = veg%rp20(:)
+          iveg%rpcoef(:) = veg%rpcoef(:)
+          iveg%rs20(:) = veg%rs20(:)
+          iveg%shelrb(:) = veg%shelrb(:)
+          iveg%vegcf(:) = veg%vegcf(:)
+          iveg%tminvj(:) = veg%tminvj(:)
+          iveg%toptvj(:) = veg%toptvj(:)
+          iveg%tmaxvj(:) = veg%tmaxvj(:)
+          iveg%vbeta(:) = veg%vbeta(:)
+          iveg%vcmax(:) = veg%vcmax(:)
+          iveg%xfang(:) = veg%xfang(:)
+          iveg%extkn(:) = veg%extkn(:)
+          iveg%wai(:) = veg%wai(:)
+          iveg%deciduous(:) = veg%deciduous(:)
+          iveg%froot(:,:) = veg%froot(:,:)
+          iveg%refl(:,:) = veg%refl(:,:)
+          iveg%taul(:,:) = veg%taul(:,:)
+          iveg%vlaimax(:) = veg%vlaimax(:)
+          iveg%a1gs(:) = veg%a1gs(:)
+          iveg%d0gs(:) = veg%d0gs(:)
+          iveg%alpha(:) = veg%alpha(:)
+          iveg%convex(:) = veg%convex(:)
+          iveg%cfrd(:) = veg%cfrd(:)
+          iveg%gswmin(:) = veg%gswmin(:)
+          iveg%conkc0(:) = veg%conkc0(:)
+          iveg%conko0(:) = veg%conko0(:)
+          iveg%ekc(:) = veg%ekc(:)
+          iveg%eko(:) = veg%eko(:)
+          iveg%g0(:) = veg%g0(:)
+          iveg%g1(:) = veg%g1(:)
+          iveg%rootbeta(:) = veg%rootbeta(:)
+          iveg%gamma(:) = veg%gamma(:)
+          iveg%F10(:) = veg%F10(:)
+          iveg%ZR(:) = veg%ZR(:)
+          iveg%clitt(:) = veg%clitt(:)
+          iveg%disturbance_interval(:,:) = veg%disturbance_interval(:,:)
+          iveg%disturbance_intensity(:,:) = veg%disturbance_intensity(:,:)
+
           ! MPI: create inp_t types to scatter input data to the workers
           ! at the start of every timestep
           !CALL master_intypes (comm,met,veg)
@@ -528,7 +599,7 @@ SUBROUTINE mpidrv_master (comm, trunk_sumbal, dels, koffset, kend, PLUME, CRU)
           ! MPI: create type to send restart data back to the master
           ! only if restart file is to be created
           IF(output%restart) THEN
-            CALL master_restart_types (comm, canopy, air)
+            CALL master_restart_types (comm, canopy, air, bgc)
           END IF
 
           !  CALL zero_sum_casa(sum_casapool, sum_casaflux)
@@ -680,6 +751,8 @@ SUBROUTINE mpidrv_master (comm, trunk_sumbal, dels, koffset, kend, PLUME, CRU)
 !$          CALL read_casa_dump( ncfile, casamet, casaflux, casa_it, kend, .FALSE. )
 !$        ENDIF
 
+          ! Zero out lai where there is no vegetation acc. to veg. index
+          WHERE ( iveg%iveg(:) .GE. 14 ) iveg%vlai = 0.
 
 !$        ! At first time step of year, set tile area according to updated LU areas
 !$        IF (ktau == 1 .and. CABLE_USER%POPLUC) THEN
@@ -750,9 +823,6 @@ SUBROUTINE mpidrv_master (comm, trunk_sumbal, dels, koffset, kend, PLUME, CRU)
           met%ofsd = met%fsd(:,1) + met%fsd(:,2)
           canopy%oldcansto=canopy%cansto
 
-          ! Zero out lai where there is no vegetation acc. to veg. index
-          WHERE ( iveg%iveg(:) .GE. 14 ) iveg%vlai = 0.
-
           ! Write time step's output to file if either: we're not spinning up
           ! or we're spinning up and the spinup has converged:
           ! MPI: TODO: pull mass and energy balance calculation from write_output
@@ -1815,6 +1885,12 @@ SUBROUTINE master_cable_params (comm,met,air,ssnow,veg,bgc,soil,canopy,&
             &                             types(bidx), ierr)
        blen(bidx) = 1
 
+       bidx = bidx + 1
+       CALL MPI_Get_address (ssnow%wbliq(off,1), displs(bidx), ierr)
+       CALL MPI_Type_create_hvector (ms, r2len, r2stride, MPI_BYTE, &
+            &                             types(bidx), ierr)
+       blen(bidx) = 1
+
        ! additional  for sli
        bidx = bidx + 1
        CALL MPI_Get_address (ssnow%S(off,1), displs(bidx), ierr)
@@ -7095,7 +7171,7 @@ END SUBROUTINE master_climate_types
   !CLNEND SUBROUTINE master_casa_restart_types
 
   ! MPI: creates datatype handles to receive restart data from workers
-  SUBROUTINE master_restart_types (comm, canopy, air)
+  SUBROUTINE master_restart_types (comm, canopy, air, bgc)
 
     USE mpi
 
@@ -7109,6 +7185,7 @@ SUBROUTINE master_restart_types (comm, canopy, air)
 
     TYPE(canopy_type), INTENT(IN) :: canopy
     TYPE (air_type),INTENT(IN)        :: air
+    TYPE (bgc_pool_type), INTENT(IN) :: bgc
     !  TYPE (casa_pool),           INTENT(INOUT) :: casapool
     !  TYPE (casa_flux),           INTENT(INOUT) :: casaflux
     !  TYPE (casa_met),            INTENT(INOUT) :: casamet
@@ -7170,6 +7247,18 @@ SUBROUTINE master_restart_types (comm, canopy, air)
             &                             types(bidx), ierr)
        blocks(bidx) = 1
 
+       bidx = bidx + 1
+       CALL MPI_Get_address (bgc%cplant(off,1), displs(bidx), ierr)
+       CALL MPI_Type_create_hvector (ncp, r1len, r1stride, MPI_BYTE, &
+            &                             types(bidx), ierr)
+       blocks(bidx) = 1
+
+       bidx = bidx + 1
+       CALL MPI_Get_address (bgc%csoil(off,1), displs(bidx), ierr)
+       CALL MPI_Type_create_hvector (ncs, r1len, r1stride, MPI_BYTE, &
+            &                             types(bidx), ierr)
+       blocks(bidx) = 1
+
        last2d = bidx
 
        ! ------------- 1D vectors -------------
diff --git a/src/offline/cable_mpiworker.F90 b/src/offline/cable_mpiworker.F90
@@ -365,7 +365,7 @@ SUBROUTINE mpidrv_worker (comm)
           ! only if restart file is to be created
           IF(output%restart) THEN
 
-            CALL worker_restart_type (comm, canopy, air)
+            CALL worker_restart_type (comm, canopy, air, bgc)
 
           END IF
 
@@ -1147,6 +1147,10 @@ SUBROUTINE worker_cable_params (comm,met,air,ssnow,veg,bgc,soil,canopy,&
     CALL MPI_Get_address (ssnow%wblf, displs(bidx), ierr)
     blen(bidx) = ms * r2len
 
+    bidx = bidx + 1
+    CALL MPI_Get_address (ssnow%wbliq, displs(bidx), ierr)
+    blen(bidx) = ms * r2len
+
     ! additional  for sli
     bidx = bidx + 1
     CALL MPI_Get_address (ssnow%S, displs(bidx), ierr)
@@ -6350,7 +6354,7 @@ END SUBROUTINE worker_climate_types
   ! MPI: creates restart_t type to send to the master the fields
   ! that are only required for the restart file but not included in the
   ! results sent at the end of each time step
-  SUBROUTINE worker_restart_type (comm, canopy, air)
+  SUBROUTINE worker_restart_type (comm, canopy, air, bgc)
 
     USE mpi
 
@@ -6362,6 +6366,7 @@ SUBROUTINE worker_restart_type (comm, canopy, air)
 
     TYPE(canopy_type), INTENT(IN) :: canopy
     TYPE (air_type),INTENT(IN)     :: air
+    TYPE (bgc_pool_type), INTENT(IN) :: bgc
 
     ! MPI: temp arrays for marshalling all types into a struct
     INTEGER, ALLOCATABLE, DIMENSION(:) :: blocks
@@ -6404,6 +6409,14 @@ SUBROUTINE worker_restart_type (comm, canopy, air)
     ! MPI: gol124: changed to r1 when Bernard ported to CABLE_r491
     blocks(bidx) = r1len * ms
 
+    bidx = bidx + 1
+    CALL MPI_Get_address (bgc%cplant(off,1), displs(bidx), ierr)
+    blocks(bidx) = r1len * ncp
+
+    bidx = bidx + 1
+    CALL MPI_Get_address (bgc%csoil(off,1), displs(bidx), ierr)
+    blocks(bidx) = r1len * ncs
+
     bidx = bidx + 1
     CALL MPI_Get_address (canopy%cduv(off), displs(bidx), ierr)
     blocks(bidx) = r1len
diff --git a/src/offline/cable_serial.F90 b/src/offline/cable_serial.F90
@@ -338,6 +338,10 @@ SUBROUTINE serialdrv(trunk_sumbal, NRRRR, dels, koffset, kend, GSWP_MID, PLUME,
           LOY = 365
           kend = NINT(24.0*3600.0/dels) * LOY
 
+        CASE ('gswp3')
+          ncciy = CurYear
+          CALL open_met_file( dels, koffset, kend, spinup, CTFRZ )
+
         CASE ('site')
           ! site experiment eg AmazonFace (spinup or  transient run type)
           kend = NINT(24.0*3600.0/dels) * LOY
@@ -532,7 +536,7 @@ SUBROUTINE serialdrv(trunk_sumbal, NRRRR, dels, koffset, kend, GSWP_MID, PLUME,
             met%Pdep = site%Pdep  *1000./10000./365. ! kg ha-1 y-1 > g m-2 d-1
             met%fsd = MAX(met%fsd,0.0)
 
-          CASE ('')
+          CASE DEFAULT
             CALL get_met_data( spinup, spinConv, met, soil,            &
               rad, veg, kend, dels, CTFRZ, ktau+koffset,                 &
               kstart+koffset )
@@ -702,7 +706,7 @@ SUBROUTINE serialdrv(trunk_sumbal, NRRRR, dels, koffset, kend, GSWP_MID, PLUME,
           IF ( (.NOT. CASAONLY) .AND. spinConv ) THEN
             !mpidiff
             SELECT CASE (TRIM(cable_user%MetType))
-            CASE ('plum', 'cru', 'bios', 'gswp', 'site')
+            CASE ('plum', 'cru', 'bios', 'gswp', 'gswp3', 'site')
               CALL write_output( dels, ktau_tot, met, canopy, casaflux, casapool, casamet, &
                    ssnow, rad, bal, air, soil, veg, CSBOLTZ, CEMLEAF, CEMSOIL )
             CASE DEFAULT
@@ -911,6 +915,7 @@ SUBROUTINE serialdrv(trunk_sumbal, NRRRR, dels, koffset, kend, GSWP_MID, PLUME,
           IF ( YYYY .EQ. CABLE_USER%YearEnd .AND. &
                NRRRR .GT. 1 ) DEALLOCATE ( GSWP_MID )
         ENDIF
+        IF (TRIM(cable_user%MetType) == "gswp3") CALL close_met_file
 
         IF ((icycle.GT.0).AND.(.NOT.casaonly)) THEN
           ! re-initalise annual flux sums
@@ -1025,6 +1030,7 @@ SUBROUTINE serialdrv(trunk_sumbal, NRRRR, dels, koffset, kend, GSWP_MID, PLUME,
 
 
   IF ( TRIM(cable_user%MetType) .NE. "gswp" .AND. &
+       TRIM(cable_user%MetType) .NE. "gswp3" .AND. &
        TRIM(cable_user%MetType) .NE. "plum" .AND. &
        TRIM(cable_user%MetType) .NE. "cru" ) CALL close_met_file
 
