Changeset 9566


Ignore:
Timestamp:
2018-05-10T15:50:33+02:00 (2 years ago)
Author:
dancopsey
Message:

Use complex numbers in area calculations.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/UKMO/dev_r5518_new_runoff_coupling/NEMOGCM/NEMO/OPA_SRC/SBC/cpl_rnf_1d.F90

    r9515 r9566  
    2222   USE dom_oce         ! Domain sizes (for grid box area e1e2t) 
    2323   USE sbc_oce         ! Surface boundary condition: ocean fields 
     24   USE lib_fortran,    ONLY: DDPDD 
    2425    
    2526   IMPLICIT NONE 
     
    3233      INTEGER, POINTER, DIMENSION(:,:)    ::   river_number       !: River outflow number 
    3334      REAL(wp), POINTER, DIMENSION(:)     ::   river_area         ! 1D array listing areas of each river outflow (m2) 
     35      COMPLEX(wp), POINTER, DIMENSION(:)  ::   river_area_c       ! Comlex version of river_area for use in bit reproducible sums (m2) 
    3436   END TYPE RIVERS_DATA 
    3537    
     
    5658      INTEGER                                   ::   ios                 ! Local integer output status for namelist read 
    5759      INTEGER                                   ::   inum 
    58       INTEGER                                   ::   ii, jj              !: Loop indices 
     60      INTEGER                                   ::   ii, jj, rr          !: Loop indices 
    5961      INTEGER                                   ::   max_river 
    6062      REAL(wp), POINTER, DIMENSION(:,:)         ::   river_number        ! 2D array containing the river outflow numbers 
     
    111113      ! Get the area of each river outflow 
    112114      ALLOCATE( rivers%river_area( nn_cpl_river ) ) 
    113       rivers%river_area(:) = 0.0 
     115      ALLOCATE( rivers%river_area_c( nn_cpl_river ) ) 
     116      rivers%river_area_c(:) = CMPLX( 0.e0, 0.e0, wp ) 
    114117      DO ii = nldi, nlei       
    115118        DO jj = nldj, nlej 
    116           IF ( rivers%river_number(ii,jj) > 0 .AND. rivers%river_number(ii,jj) <= nn_cpl_river ) THEN 
    117             rivers%river_area(rivers%river_number(ii,jj)) = rivers%river_area(rivers%river_number(ii,jj)) + e1e2t(ii,jj) 
     119          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) 
     120            IF ( rivers%river_number(ii,jj) > 0 .AND. rivers%river_number(ii,jj) <= nn_cpl_river ) THEN 
     121              ! Add the area of each grid box (e1e2t) into river_area_c using DDPDD which should maintain bit reproducibility (needs to be checked) 
     122              CALL DDPDD( CMPLX( e1e2t(ii,jj), 0.e0, wp ), rivers%river_area_c(rivers%river_number(ii,jj) ) ) 
     123            END IF 
     124          ELSE 
     125            IF ( jj == nldj .AND. tmask(ii,jj) > 0.5 ) THEN 
     126              WRITE(numout,*) 'At duplicated point at ii = ',ii 
     127            END IF 
    118128          END IF 
    119129        END DO 
     
    121131       
    122132      ! Use mpp_sum to add together river areas on other processors 
    123       CALL mpp_sum( rivers%river_area, nn_cpl_river ) 
     133      CALL mpp_sum( rivers%river_area_c, nn_cpl_river ) 
     134       
     135      ! Convert from complex number to real 
     136      ! DO rr = 1, nn_cpl_river 
     137      !    rivers%river_area(rr) = rivers%river_area_c(rr) 
     138      ! END DO 
     139      rivers%river_area(:) = REAL(rivers%river_area_c(:),wp) 
     140       
    124141      IF ( ln_print_river_info ) THEN 
    125142        WRITE(numout,*) 'Area of rivers 1 to 10 are ',rivers%river_area(1:10) 
Note: See TracChangeset for help on using the changeset viewer.