Fix serial-MPI non-reproducibility for gswp3 CASA-CNP configuration

src/offline/cable_mpicommon.F90

@@ -99,7 +99,7 @@ MODULE cable_mpicommon
   ! MPI: number of final casa result matrices and vectors to receive
   ! by the master for casa_poolout and casa_fluxout
   INTEGER, PARAMETER :: ncasa_mat = 37    ! add three more wood product variables
-  INTEGER, PARAMETER :: ncasa_vec = 58   ! vh changed on 5-feb-2016 for adding sapwood area and frac_sapwood
+  INTEGER, PARAMETER :: ncasa_vec = 66
   ! 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

src/offline/cable_mpimaster.F90

@@ -248,7 +248,9 @@ SUBROUTINE mpidrv_master (comm, trunk_sumbal, dels, koffset, kend, PLUME, CRU)
          ctime = 0,  &  ! day count for casacnp
          YYYY,       &  !
          LOY,        &  ! Length of Year
-         maxdiff(2)     ! location of maximum in convergence test
+         maxdiff(2), &  ! location of maximum in convergence test
+         count_sum_casa ! number of time steps over which casa pools &
+                                !and fluxes are aggregated (for output)
 
     CHARACTER :: dum*9, str1*9, str2*9, str3*9  ! dummy char for fileName generation
 
@@ -602,8 +604,8 @@ SUBROUTINE mpidrv_master (comm, trunk_sumbal, dels, koffset, kend, PLUME, CRU)
             CALL master_restart_types (comm, canopy, air, bgc)
           END IF
 
-          !  CALL zero_sum_casa(sum_casapool, sum_casaflux)
-          !  count_sum_casa = 0
+          CALL zero_sum_casa(sum_casapool, sum_casaflux)
+          count_sum_casa = 0
 
           ! CALL master_sumcasa_types(comm, sum_casapool, sum_casaflux)
           IF( icycle>0 .AND. spincasa) THEN
@@ -766,6 +768,13 @@ SUBROUTINE mpidrv_master (comm, trunk_sumbal, dels, koffset, kend, PLUME, CRU)
               ! receive casa update from worker
               CALL master_receive (ocomm, oktau, casa_ts)
 
+              IF(MOD((oktau-kstart+1),ktauday)==0) THEN
+                ! update time-aggregates of casa pools and fluxes
+                CALL update_sum_casa(sum_casapool, sum_casaflux, casapool, casaflux, &
+                    & .TRUE. , .FALSE., 1)
+                count_sum_casa = count_sum_casa + 1
+              END IF
+
               CALL MPI_Waitall (wnp, recv_req, recv_stats, ierr)
               ! receive casa dump requirements from worker
               IF ( ((.NOT.spinup).OR.(spinup.AND.spinConv)) .AND.   &
@@ -835,11 +844,12 @@ SUBROUTINE mpidrv_master (comm, trunk_sumbal, dels, koffset, kend, PLUME, CRU)
               IF ( IS_CASA_TIME("write", yyyy, oktau, kstart, &
                                 koffset, kend, ktauday, logn) ) THEN
                 ctime = ctime +1
-
-
-                CALL WRITE_CASA_OUTPUT_NC (veg, casamet, casapool, casabal, casaflux, &
-                     CASAONLY, ctime, &
-                     ( ktau.EQ.kend .AND. YYYY .EQ.cable_user%YearEnd ) )
+                CALL update_sum_casa(sum_casapool, sum_casaflux, casapool, casaflux, &
+                    .FALSE. , .TRUE. , count_sum_casa)
+                CALL WRITE_CASA_OUTPUT_NC (veg, casamet, sum_casapool, casabal, sum_casaflux, &
+                    CASAONLY, ctime, ( oktau == kend .AND. YYYY == cable_user%YearEnd ) )
+                count_sum_casa = 0
+                CALL zero_sum_casa(sum_casapool, sum_casaflux)
               ENDIF
             ENDIF
 
@@ -952,6 +962,13 @@ SUBROUTINE mpidrv_master (comm, trunk_sumbal, dels, koffset, kend, PLUME, CRU)
 
             CALL master_receive (ocomm, oktau, casa_ts)
 
+            IF(MOD((oktau-kstart+1),ktauday)==0) THEN
+              ! update time-aggregates of casa pools and fluxes
+              CALL update_sum_casa(sum_casapool, sum_casaflux, casapool, casaflux, &
+                  & .TRUE. , .FALSE., 1)
+              count_sum_casa = count_sum_casa + 1
+            END IF
+
             IF ( ((.NOT.spinup).OR.(spinup.AND.spinConv)) .AND. &
                  ( IS_CASA_TIME("dwrit", yyyy, oktau, kstart, &
                  koffset, kend, ktauday, logn) ) ) THEN
@@ -1024,9 +1041,12 @@ SUBROUTINE mpidrv_master (comm, trunk_sumbal, dels, koffset, kend, PLUME, CRU)
         IF((.NOT.spinup).OR.(spinup.AND.spinConv)) THEN
           IF(icycle >0) THEN
             ctime = ctime +1
-            CALL WRITE_CASA_OUTPUT_NC (veg, casamet, casapool, casabal, casaflux, &
-                 CASAONLY, ctime, ( ktau.EQ.kend .AND. YYYY .EQ.               &
-                 cable_user%YearEnd ) )
+            CALL update_sum_casa(sum_casapool, sum_casaflux, casapool, casaflux, &
+                .FALSE. , .TRUE. , count_sum_casa)
+            CALL WRITE_CASA_OUTPUT_NC (veg, casamet, sum_casapool, casabal, sum_casaflux, &
+                CASAONLY, ctime, ( oktau == kend .AND. YYYY == cable_user%YearEnd ) )
+            count_sum_casa = 0
+            CALL zero_sum_casa(sum_casapool, sum_casaflux)
             IF ( cable_user%CALL_POP ) THEN
 
               ! CALL master_receive_pop(POP, ocomm)
@@ -6656,6 +6676,38 @@ SUBROUTINE master_casa_types (comm, casapool, casaflux, &
        CALL MPI_Get_address (casaflux%Cplant_turnover_resource_limitation(off), displs(bidx), ierr)
        blocks(bidx) = r2len
 
+       bidx = bidx + 1
+       CALL MPI_Get_address (casaflux%Pupland(off), displs(bidx), ierr)
+       blocks(bidx) = r2len
+
+       bidx = bidx + 1
+       CALL MPI_Get_address (casaflux%Plittermin(off), displs(bidx), ierr)
+       blocks(bidx) = r2len
+
+       bidx = bidx + 1
+       CALL MPI_Get_address (casaflux%Psmin(off), displs(bidx), ierr)
+       blocks(bidx) = r2len
+
+       bidx = bidx + 1
+       CALL MPI_Get_address (casaflux%Psimm(off), displs(bidx), ierr)
+       blocks(bidx) = r2len
+
+       bidx = bidx + 1
+       CALL MPI_Get_address (casaflux%kplab(off), displs(bidx), ierr)
+       blocks(bidx) = r2len
+
+       bidx = bidx + 1
+       CALL MPI_Get_address (casaflux%kpsorb(off), displs(bidx), ierr)
+       blocks(bidx) = r2len
+
+       bidx = bidx + 1
+       CALL MPI_Get_address (casaflux%kpocc(off), displs(bidx), ierr)
+       blocks(bidx) = r2len
+
+       bidx = bidx + 1
+       CALL MPI_Get_address (casaflux%FluxCtoco2(off), displs(bidx), ierr)
+       blocks(bidx) = r2len
+
        types(last2d+1:bidx) = MPI_BYTE
 
        ! MPI: sanity check

src/offline/cable_mpiworker.F90

@@ -6055,6 +6055,38 @@ SUBROUTINE worker_casa_type (comm, casapool,casaflux, &
     CALL MPI_Get_address (casaflux%Cplant_turnover_resource_limitation(off), displs(bidx), ierr)
     blocks(bidx) = r2len
 
+    bidx = bidx + 1
+    CALL MPI_Get_address (casaflux%Pupland(off), displs(bidx), ierr)
+    blocks(bidx) = r2len
+
+    bidx = bidx + 1
+    CALL MPI_Get_address (casaflux%Plittermin(off), displs(bidx), ierr)
+    blocks(bidx) = r2len
+
+    bidx = bidx + 1
+    CALL MPI_Get_address (casaflux%Psmin(off), displs(bidx), ierr)
+    blocks(bidx) = r2len
+
+    bidx = bidx + 1
+    CALL MPI_Get_address (casaflux%Psimm(off), displs(bidx), ierr)
+    blocks(bidx) = r2len
+
+    bidx = bidx + 1
+    CALL MPI_Get_address (casaflux%kplab(off), displs(bidx), ierr)
+    blocks(bidx) = r2len
+
+    bidx = bidx + 1
+    CALL MPI_Get_address (casaflux%kpsorb(off), displs(bidx), ierr)
+    blocks(bidx) = r2len
+
+    bidx = bidx + 1
+    CALL MPI_Get_address (casaflux%kpocc(off), displs(bidx), ierr)
+    blocks(bidx) = r2len
+
+    bidx = bidx + 1
+    CALL MPI_Get_address (casaflux%FluxCtoco2(off), displs(bidx), ierr)
+    blocks(bidx) = r2len
+
     ! MPI: sanity check
     IF (bidx /= ntyp) THEN
        WRITE (*,*) 'worker: invalid number of casa fields, fix it!'

src/science/casa-cnp/casa_phenology.F90

@@ -46,21 +46,24 @@ MODULE phenvariable
 CONTAINS
 
   SUBROUTINE alloc_phenvariable(phen,arraysize)
+    !* Allocate phen derived type instance.
+    ! Allocated arrays are initialised to zero.
 
-    IMPLICIT NONE
     TYPE(phen_variable), INTENT(INOUT) :: phen
-    INTEGER,             INTENT(IN) :: arraysize
+    INTEGER,             INTENT(IN   ) :: arraysize
 
-    ALLOCATE(phen%Tkshed(mvtype))
-    ALLOCATE(phen%phase(arraysize),         &
-         phen%doyphase(arraysize,mphase))
-    ALLOCATE(phen%phen(arraysize), &
-         phen%aphen(arraysize), &
-         phen%phasespin(arraysize,mdyear), &
-         phen%doyphasespin_1(arraysize,mdyear), &
-         phen%doyphasespin_2(arraysize,mdyear), &
-         phen%doyphasespin_3(arraysize,mdyear), &
-         phen%doyphasespin_4(arraysize,mdyear))
+    ALLOCATE(phen%Tkshed(mvtype), source=0.0_r_2)
+    ALLOCATE(phen%phen(arraysize), phen%aphen(arraysize), source=0.0)
+    ALLOCATE(                                &
+      phen%phase(arraysize),                 &
+      phen%doyphase(arraysize,mphase),       &
+      phen%phasespin(arraysize,mdyear),      &
+      phen%doyphasespin_1(arraysize,mdyear), &
+      phen%doyphasespin_2(arraysize,mdyear), &
+      phen%doyphasespin_3(arraysize,mdyear), &
+      phen%doyphasespin_4(arraysize,mdyear), &
+      source=0                               &
+    )
   END SUBROUTINE alloc_phenvariable
 
 END MODULE phenvariable