Changeset 9566

Timestamp:

2018-05-10T15:50:33+02:00 (6 years ago)

Author:

dancopsey

Message:

Use complex numbers in area calculations.

File:

• branches/UKMO/dev_r5518_new_runoff_coupling/NEMOGCM/NEMO/OPA_SRC/SBC/cpl_rnf_1d.F90 (modified) (5 diffs)

branches/UKMO/dev_r5518_new_runoff_coupling/NEMOGCM/NEMO/OPA_SRC/SBC/cpl_rnf_1d.F90

@@ -22 +22 @@
    USE dom_oce         ! Domain sizes (for grid box area e1e2t)
    USE sbc_oce         ! Surface boundary condition: ocean fields
+   USE lib_fortran,    ONLY: DDPDD
    
    IMPLICIT NONE
@@ -32 +33 @@
       INTEGER, POINTER, DIMENSION(:,:)    ::   river_number       !: River outflow number
       REAL(wp), POINTER, DIMENSION(:)     ::   river_area         ! 1D array listing areas of each river outflow (m2)
+      COMPLEX(wp), POINTER, DIMENSION(:)  ::   river_area_c       ! Comlex version of river_area for use in bit reproducible sums (m2)
    END TYPE RIVERS_DATA
    
@@ -56 +58 @@
       INTEGER                                   ::   ios                 ! Local integer output status for namelist read
       INTEGER                                   ::   inum
-      INTEGER                                   ::   ii, jj              !: Loop indices
+      INTEGER                                   ::   ii, jj, rr          !: Loop indices
       INTEGER                                   ::   max_river
       REAL(wp), POINTER, DIMENSION(:,:)         ::   river_number        ! 2D array containing the river outflow numbers
@@ -111 +113 @@
       ! Get the area of each river outflow
       ALLOCATE( rivers%river_area( nn_cpl_river ) )
-      rivers%river_area(:) = 0.0
+      ALLOCATE( rivers%river_area_c( nn_cpl_river ) )
+      rivers%river_area_c(:) = CMPLX( 0.e0, 0.e0, wp )
       DO ii = nldi, nlei      
         DO jj = nldj, nlej
-          IF ( rivers%river_number(ii,jj) > 0 .AND. rivers%river_number(ii,jj) <= nn_cpl_river ) THEN
-            rivers%river_area(rivers%river_number(ii,jj)) = rivers%river_area(rivers%river_number(ii,jj)) + e1e2t(ii,jj)
+          IF ( tmask_i(ii,jj) > 0.5 ) THEN  ! This makes sure we are not at a duplicated point (at north fold or east-west cyclic point)
+            IF ( rivers%river_number(ii,jj) > 0 .AND. rivers%river_number(ii,jj) <= nn_cpl_river ) THEN
+              ! Add the area of each grid box (e1e2t) into river_area_c using DDPDD which should maintain bit reproducibility (needs to be checked)
+              CALL DDPDD( CMPLX( e1e2t(ii,jj), 0.e0, wp ), rivers%river_area_c(rivers%river_number(ii,jj) ) )
+            END IF
+          ELSE
+            IF ( jj == nldj .AND. tmask(ii,jj) > 0.5 ) THEN
+              WRITE(numout,*) 'At duplicated point at ii = ',ii
+            END IF
           END IF
         END DO
@@ -121 +131 @@
       
       ! Use mpp_sum to add together river areas on other processors
-      CALL mpp_sum( rivers%river_area, nn_cpl_river )
+      CALL mpp_sum( rivers%river_area_c, nn_cpl_river )
+      
+      ! Convert from complex number to real
+      ! DO rr = 1, nn_cpl_river
+      !    rivers%river_area(rr) = rivers%river_area_c(rr)
+      ! END DO
+      rivers%river_area(:) = REAL(rivers%river_area_c(:),wp)
+      
       IF ( ln_print_river_info ) THEN
         WRITE(numout,*) 'Area of rivers 1 to 10 are ',rivers%river_area(1:10)