Changeset 7259

Timestamp:

2021-07-27T17:15:46+02:00 (3 years ago)

Author:

agnes.ducharne

Message:

As in r6385, so that interpolation using aggregate_p is now done in parallel. Runtime has been divided by 6 for a 5d run from scratch at 2deg offline!

Location:

branches/ORCHIDEE_2_2/ORCHIDEE

Files:

• src_driver/orchideedriver.f90 (modified) (1 diff)
• src_global/grid.f90 (modified) (1 diff)
• src_global/interpol_help.f90 (modified) (10 diffs)
• src_sechiba/intersurf.f90 (modified) (4 diffs)

branches/ORCHIDEE_2_2/ORCHIDEE/src_driver/orchideedriver.f90

@@ -412 +412 @@
   ! Variables *_g were allocated with the CALL init_grid
   !
-  IF ( is_root_prc) THEN
-     !
-     lalo_g(:,:) = lalo_glo(:,:)
-     lon_g(:,:) = lon_glo(:,:)
-     lat_g(:,:) = lat_glo(:,:)
-     !
-  ENDIF
+  !
+  lalo_g(:,:) = lalo_glo(:,:)
+  contfrac_g(:) = contfrac_glo(:)
+  lon_g(:,:) = lon_glo(:,:)
+  lat_g(:,:) = lat_glo(:,:)
   !
   !

branches/ORCHIDEE_2_2/ORCHIDEE/src_global/grid.f90

@@ -443 +443 @@
     !
     CALL grid_scatter()
+    !
+    CALL bcast(neighbours_g)
+    CALL bcast(resolution_g)
     !
     IF ( printlev >= 3 ) WRITE(numout,*) 'Leaving grid_stuff'

branches/ORCHIDEE_2_2/ORCHIDEE/src_global/interpol_help.f90

@@ -42 +42 @@
   SUBROUTINE aggregate_2d (nbpt, lalo, neighbours, resolution, contfrac, &
        &                iml, jml, lon_rel, lat_rel, mask, callsign, &
-       &                incmax, indinc, areaoverlap, ok)
+       &                incmax, indinc, areaoverlap, ok, opt_nbpt_start, opt_nbpt_end)
 
     USE grid, ONLY : global
@@ -61 +61 @@
     CHARACTER(LEN=*), INTENT(in) :: callsign             ! Allows to specify which variable is beeing treated
     INTEGER(i_std), INTENT(in)   :: incmax              ! Maximum point of the fine grid we can store.
+    INTEGER(i_std), OPTIONAL, INTENT(in)    :: opt_nbpt_start            ! Input Start grid cell interpolation 
+    INTEGER(i_std), OPTIONAL, INTENT(in)    :: opt_nbpt_end              ! Input End grid cell interpolation
     !
     ! Output
     !
-    INTEGER(i_std), INTENT(out)  :: indinc(nbpt,incmax,2)
-    REAL(r_std), INTENT(out)      :: areaoverlap(nbpt,incmax)
+    INTEGER(i_std), INTENT(out)  :: indinc(:,:,:)
+    REAL(r_std), INTENT(out)      :: areaoverlap(:,:)
     LOGICAL, OPTIONAL, INTENT(out)      :: ok            ! return code
     !
@@ -78 +80 @@
     REAL(r_std) :: domain_minlon,domain_maxlon,domain_minlat,domain_maxlat
     INTEGER(i_std) :: minLon(1), maxLon(1)
+    INTEGER(i_std) :: nbpt_start            ! Start grid cell interpolation 
+    INTEGER(i_std) :: nbpt_end              ! End grid cell interpolation
+    INTEGER(i_std) :: landpoint_idx 
 
     INTEGER                  :: ALLOC_ERR
@@ -108 +113 @@
     ALLOCATE (searchind(iml*jml,2), STAT=ALLOC_ERR)
     IF (ALLOC_ERR/=0) CALL ipslerr_p(3,'aggregate_2d', 'ERROR IN ALLOCATION of searchind','','')
+
+    nbpt_start = 1
+    nbpt_end = nbpt
+    IF (PRESENT(opt_nbpt_start) .AND. PRESENT(opt_nbpt_end)) THEN
+        nbpt_start = opt_nbpt_start
+        nbpt_end = opt_nbpt_end
+    ENDIF
 
     IF (PRESENT(ok)) ok = .TRUE.
@@ -256 +268 @@
     !
     fopt_max = -1
-    DO ib =1, nbpt
+    DO ib = nbpt_start, nbpt_end
+       landpoint_idx = ib - nbpt_start + 1 
        !
        !   Give a progress meter
@@ -355 +368 @@
                    ay = (MIN(lat_up, laup_rel(ip,jp))-MAX(lat_low,lalow_rel(ip,jp)))*pi/180. * R_Earth
                    !
-                   areaoverlap(ib, fopt) = ax*ay
-                   indinc(ib, fopt, 1) = ip
-                   indinc(ib, fopt, 2) = jp
+                   areaoverlap(landpoint_idx, fopt) = ax*ay
+                   indinc(landpoint_idx, fopt, 1) = ip
+                   indinc(landpoint_idx, fopt, 2) = jp
                    !
                    ! If this point was 100% within the grid then we can de-select it from our
@@ -403 +416 @@
     ENDDO
     !
-    DO ib=1,nbpt
+    DO ib=nbpt_start, nbpt_end
+       landpoint_idx = ib - nbpt_start + 1
        DO fopt=1,incmax
-          IF (( indinc(ib,fopt,1) == 0 .AND. indinc(ib,fopt,2) > 0) .OR.&
-               & ( indinc(ib,fopt,2) == 0 .AND. indinc(ib,fopt,1) > 0) ) THEN
+          IF (( indinc(landpoint_idx,fopt,1) == 0 .AND. indinc(landpoint_idx,fopt,2) > 0) .OR.&
+               & ( indinc(landpoint_idx,fopt,2) == 0 .AND. indinc(landpoint_idx,fopt,1) > 0) ) THEN
              WRITE(*,*) "aggregate_2d PROBLEM : point =",ib, fopt," Indicies = ", &
-                  & indinc(ib,fopt,1), indinc(ib,fopt,2), areaoverlap(ib,fopt)
+                  & indinc(landpoint_idx,fopt,1), indinc(landpoint_idx,fopt,2), areaoverlap(landpoint_idx,fopt)
           ENDIF
        ENDDO
@@ -815 +829 @@
     REAL(r_std), INTENT(out)       :: sub_area(nbpt,nbvmax) 
     LOGICAL, OPTIONAL, INTENT(out) :: ok            ! return code
-
-    INTEGER(i_std)   :: sub_index_g(nbp_glo,nbvmax,2)
-    REAL(r_std)       :: sub_area_g(nbp_glo,nbvmax)
+    INTEGER(i_std)                 :: nbp_start, nbp_end
+
+
+    INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:,:)  :: sub_index_g 
+    REAL(r_std), ALLOCATABLE, DIMENSION(:,:)       :: sub_area_g 
 
     IF ( grid_type == regular_lonlat ) THEN
-       IF (is_root_prc) CALL aggregate_2d(nbp_glo, lalo_g, neighbours_g, resolution_g, contfrac_g, &
-            &                                  iml, jml, lon_ful, lat_ful, mask, callsign,   &
-            &                                  nbvmax, sub_index_g, sub_area_g, ok)
+       nbp_start = nbp_mpi_para_begin(mpi_rank) + (nbp_omp_para_begin(omp_rank) - 1)
+#ifdef CPP_OMP
+       ! possible abstraction? 
+       nbp_end = nbp_mpi_para_begin(mpi_rank) + (nbp_omp_para_end(omp_rank) - 1) 
+#else
+       nbp_end = nbp_mpi_para_end(mpi_rank)
+#endif
+       CALL aggregate_2d(nbp_glo, lalo_g, neighbours_g, resolution_g, contfrac_g, &
+            &             iml, jml, lon_ful, lat_ful, mask, callsign,   &
+            &             nbvmax, sub_index, sub_area, ok, &
+                          nbp_start, nbp_end)
     ELSE IF ( grid_type == regular_xy ) THEN
+        ALLOCATE(sub_index_g(nbp_glo, nbvmax, 2))
+        ALLOCATE(sub_area_g(nbp_glo, nbvmax))
+
        IF ( proj_stack(1)%code > undef_int-1 ) THEN
           CALL ipslerr(3, "aggregate_2d_p", "Regular_xy projection was not intialized.", &
@@ -831 +858 @@
             &                               iml, jml, lon_ful, lat_ful, mask, callsign,   &
             &                               nbvmax, sub_index_g, sub_area_g, ok)
+
+       CALL BCAST(ok)
+       CALL scatter(sub_index_g,sub_index)
+       CALL scatter(sub_area_g,sub_area)
+
+       DEALLOCATE(sub_index_g)
+       DEALLOCATE(sub_area_g)
     ELSE
        CALL ipslerr(3, "aggregate_2d_p", "Interpolation is only possible for regular lat/lon grids for the moment.", &
@@ -836 +870 @@
     ENDIF
     !
-    CALL BCAST(ok)
-    CALL scatter(sub_index_g,sub_index)
-    CALL scatter(sub_area_g,sub_area)
    
   END SUBROUTINE aggregate_2d_p

branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/intersurf.f90

@@ -242 +242 @@
     CALL gather2D_mpi(lon,lon_g)
     CALL gather2D_mpi(lat,lat_g)
+
+    CALL bcast(lalo_g)
+    CALL bcast(contfrac_g)
     
     CALL ioipslctrl_restini(kjit, date0, xrdt, rest_id, rest_id_stom, itau_offset, date0_shifted)
@@ -955 +958 @@
     lalo(:,:) = latlon(:,:)
     CALL gather(lalo,lalo_g)
+    CALL bcast(lalo_g)
     !
     !-
@@ -962 +966 @@
     neighbours(:,:) = zneighbours(:,:)
     CALL gather(neighbours,neighbours_g)
+    CALL bcast(neighbours_g)
     !
     resolution(:,:) = zresolution(:,:)
     CALL gather(resolution,resolution_g)
+    CALL bcast(resolution_g)
     !
     IF (grid_type==regular_lonlat) area(:) = resolution(:,1)*resolution(:,2)
     CALL gather(area,area_g)
+    CALL bcast(area_g)
     !
     !- Store the fraction of the continents only once so that the user
@@ -974 +981 @@
     contfrac(:) = zcontfrac(:)
     CALL gather(contfrac,contfrac_g)
+    CALL bcast(contfrac_g)
     !
     !