Changeset 953

Timestamp:

07/15/19 12:29:31 (5 years ago)

Author:

adurocher

Message:

trunk : GPU implementation with OpenACC ( merge from glcp.idris.fr )

Location:

codes/icosagcm/trunk

Files:

• arch/arch-HPE_INTEL.env (added)
• arch/arch-HPE_INTEL.fcm (added)
• arch/arch-HPE_INTEL.path (added)
• arch/arch-HPE_PGI.env (added)
• arch/arch-HPE_PGI.fcm (added)
• arch/arch-HPE_PGI.path (added)
• arch/arch-HPE_PGI_ACC.env (added)
• arch/arch-HPE_PGI_ACC.fcm (added)
• arch/arch-HPE_PGI_ACC.path (added)
• bld.cfg (modified) (1 diff)
• src/base/abort.F90 (modified) (1 diff)
• src/base/field.f90 (modified) (13 diffs)
• src/diagnostics/check_conserve.f90 (modified) (4 diffs)
• src/diagnostics/diagflux.F90 (modified) (1 diff)
• src/diagnostics/observable.f90 (modified) (3 diffs)
• src/dissip/dissip_gcm.F90 (moved) (moved from codes/icosagcm/trunk/src/dissip/dissip_gcm.f90) (50 diffs)
• src/dissip/sponge.f90 (modified) (2 diffs)
• src/dynamics/caldyn_gcm.F90 (modified) (4 diffs)
• src/dynamics/caldyn_hevi.F90 (moved) (moved from codes/icosagcm/trunk/src/dynamics/caldyn_hevi.f90) (7 diffs)
• src/dynamics/caldyn_kernels_base.F90 (modified) (18 diffs)
• src/dynamics/caldyn_kernels_hevi.F90 (modified) (22 diffs)
• src/kernels/advect_horiz.k90 (modified) (2 diffs)
• src/output/xios_mod.F90 (modified) (1 diff)
• src/parallel/mpipara.F90 (modified) (5 diffs)
• src/parallel/transfert_mpi.f90 (modified) (35 diffs)
• src/physics/physics.f90 (modified) (3 diffs)
• src/sphere/geometry.f90 (modified) (2 diffs)
• src/time/euler_scheme.F90 (moved) (moved from codes/icosagcm/trunk/src/time/euler_scheme.f90) (11 diffs)
• src/time/hevi_scheme.F90 (moved) (moved from codes/icosagcm/trunk/src/time/hevi_scheme.f90) (8 diffs)
• src/time/timeloop_gcm.F90 (moved) (moved from codes/icosagcm/trunk/src/time/timeloop_gcm.f90) (16 diffs)
• src/transport/advect.F90 (modified) (14 diffs)
• src/transport/advect_tracer.F90 (moved) (moved from codes/icosagcm/trunk/src/transport/advect_tracer.f90) (16 diffs)
• src/vertical/disvert.f90 (modified) (6 diffs)

codes/icosagcm/trunk/bld.cfg

@@ -52 +52 @@
 bld::excl_dep        use::netcdf
 bld::excl_dep        use::omp_lib
+bld::excl_dep        use::openacc
 bld::excl_dep        use::mpi
 bld::excl_dep        inc::mpif.h

codes/icosagcm/trunk/src/base/abort.F90

@@ -17 +17 @@
     !$omp end single
   end subroutine
+
+  !!!Abort execution when openacc is on
+  subroutine abort_acc( message )
+    use mpi_mod
+    implicit none
+    character(len=*), optional, intent(in) :: message  
+#ifdef _OPENACC
+    call dynamico_abort( "Not tested with OpenACC ! " // message )
+#endif
+  end subroutine
 end module

codes/icosagcm/trunk/src/base/field.f90

@@ -13 +13 @@
   TYPE t_field
     CHARACTER(30)      :: name
-    REAL(rstd),POINTER :: rval2d(:)
-    REAL(rstd),POINTER :: rval3d(:,:)
-    REAL(rstd),POINTER :: rval4d(:,:,:)
+    REAL(rstd),POINTER :: rval2d(:) => null()
+    REAL(rstd),POINTER :: rval3d(:,:) => null()
+    REAL(rstd),POINTER :: rval4d(:,:,:) => null()
 
     INTEGER,POINTER :: ival2d(:)
@@ -30 +30 @@
     INTEGER :: dim3
     INTEGER :: dim4
+    
+    LOGICAL :: ondevice !< flag if field is allocated on device as well
   END TYPE t_field   
 
@@ -48 +50 @@
 CONTAINS
 
-  SUBROUTINE allocate_field(field,field_type,data_type,dim1,dim2,name)
+  SUBROUTINE allocate_field(field,field_type,data_type,dim1,dim2,name,ondevice)
   USE domain_mod
   USE omp_para
@@ -56 +58 @@
     INTEGER,OPTIONAL :: dim1,dim2
     CHARACTER(*), OPTIONAL :: name
+    LOGICAL, INTENT(IN), OPTIONAL :: ondevice 
 !$OMP BARRIER
 !$OMP MASTER
-    ALLOCATE(field(ndomain))
+    ALLOCATE(field(ndomain))    
 !$OMP END MASTER
 !$OMP BARRIER
-    CALL allocate_field_(field,field_type,data_type,dim1,dim2,name)
+
+    CALL allocate_field_(field,field_type,data_type,dim1,dim2,name,ondevice)
+    
   END SUBROUTINE allocate_field
 
-  SUBROUTINE allocate_fields(nfield,field,field_type,data_type,dim1,dim2,name)
+  SUBROUTINE allocate_fields(nfield,field,field_type,data_type,dim1,dim2,name, ondevice)
   USE domain_mod
   USE omp_para
@@ -73 +78 @@
     INTEGER,OPTIONAL :: dim1,dim2
     CHARACTER(*), OPTIONAL :: name
+    LOGICAL, INTENT(IN), OPTIONAL :: ondevice
     INTEGER :: i
 !$OMP BARRIER
@@ -80 +86 @@
 !$OMP BARRIER
     DO i=1,nfield
-       CALL allocate_field_(field(:,i),field_type,data_type,dim1,dim2,name)
+       CALL allocate_field_(field(:,i),field_type,data_type,dim1,dim2,name,ondevice)
     END DO
   END SUBROUTINE allocate_fields
 
-  SUBROUTINE allocate_field_(field,field_type,data_type,dim1,dim2,name)
+  SUBROUTINE allocate_field_(field,field_type,data_type,dim1,dim2,name,ondevice)
   USE domain_mod
   USE omp_para
@@ -93 +99 @@
     INTEGER,OPTIONAL :: dim1,dim2
     CHARACTER(*), OPTIONAL :: name
+    LOGICAL, INTENT(IN), OPTIONAL :: ondevice
     INTEGER :: ind
     INTEGER :: ii_size,jj_size
@@ -119 +126 @@
       field(ind)%data_type=data_type
       field(ind)%field_type=field_type
-    
+
       IF (field_type==field_T) THEN 
         jj_size=domain(ind)%jjm
@@ -131 +138 @@
         
       IF (field(ind)%ndim==4) THEN
-        IF (data_type==type_integer) ALLOCATE(field(ind)%ival4d(ii_size*jj_size,dim1,dim2))
-        IF (data_type==type_real)    ALLOCATE(field(ind)%rval4d(ii_size*jj_size,dim1,dim2))
-        IF (data_type==type_logical) ALLOCATE(field(ind)%lval4d(ii_size*jj_size,dim1,dim2))
+         IF (data_type==type_integer) ALLOCATE(field(ind)%ival4d(ii_size*jj_size,dim1,dim2))
+         IF (data_type==type_real)    ALLOCATE(field(ind)%rval4d(ii_size*jj_size,dim1,dim2))
+         IF (data_type==type_logical) ALLOCATE(field(ind)%lval4d(ii_size*jj_size,dim1,dim2))
+
       ELSE IF (field(ind)%ndim==3) THEN
-        IF (data_type==type_integer) ALLOCATE(field(ind)%ival3d(ii_size*jj_size,dim1))
-        IF (data_type==type_real)    ALLOCATE(field(ind)%rval3d(ii_size*jj_size,dim1))
-        IF (data_type==type_logical) ALLOCATE(field(ind)%lval3d(ii_size*jj_size,dim1))
+         IF (data_type==type_integer) ALLOCATE(field(ind)%ival3d(ii_size*jj_size,dim1))
+         IF (data_type==type_real)    ALLOCATE(field(ind)%rval3d(ii_size*jj_size,dim1))
+         IF (data_type==type_logical) ALLOCATE(field(ind)%lval3d(ii_size*jj_size,dim1))
+
       ELSE IF (field(ind)%ndim==2) THEN
-        IF (data_type==type_integer) ALLOCATE(field(ind)%ival2d(ii_size*jj_size))
-        IF (data_type==type_real)    ALLOCATE(field(ind)%rval2d(ii_size*jj_size))
-        IF (data_type==type_logical) ALLOCATE(field(ind)%lval2d(ii_size*jj_size))
-      ENDIF
-      
+         IF (data_type==type_integer) ALLOCATE(field(ind)%ival2d(ii_size*jj_size))
+         IF (data_type==type_real)    ALLOCATE(field(ind)%rval2d(ii_size*jj_size))
+         IF (data_type==type_logical) ALLOCATE(field(ind)%lval2d(ii_size*jj_size))
+
+      ENDIF
+
+      field(ind)%ondevice = .FALSE.
+      IF (PRESENT(ondevice)) THEN
+         IF (ondevice) CALL create_device_field(field(ind))
+      END IF
+   
    ENDDO
 !$OMP BARRIER
@@ -160 +175 @@
     INTEGER :: ii_size,jj_size
 
-    ALLOCATE(field(ndomain_glo))    
+    ALLOCATE(field(ndomain_glo)) 
 
     DO ind=1,ndomain_glo
@@ -184 +199 @@
       field(ind)%field_type=field_type
     
+      field(ind)%ondevice = .FALSE.
+
       IF (field_type==field_T) THEN 
         jj_size=domain_glo(ind)%jjm
@@ -251 +268 @@
     INTEGER :: ind
     DO ind=1,ndomain
-      IF (.NOT. assigned_domain(ind)  .OR. .NOT. is_omp_level_master) CYCLE
-
-      data_type=field(ind)%data_type
-        
-      IF (field(ind)%ndim==4) THEN
-        IF (data_type==type_integer) DEALLOCATE(field(ind)%ival4d)
-        IF (data_type==type_real)    DEALLOCATE(field(ind)%rval4d)
-        IF (data_type==type_logical) DEALLOCATE(field(ind)%lval4d)
-      ELSE IF (field(ind)%ndim==3) THEN
-        IF (data_type==type_integer) DEALLOCATE(field(ind)%ival3d)
-        IF (data_type==type_real)    DEALLOCATE(field(ind)%rval3d)
-        IF (data_type==type_logical) DEALLOCATE(field(ind)%lval3d)
-      ELSE IF (field(ind)%ndim==2) THEN
-        IF (data_type==type_integer) DEALLOCATE(field(ind)%ival2d)
-        IF (data_type==type_real)    DEALLOCATE(field(ind)%rval2d)
-        IF (data_type==type_logical) DEALLOCATE(field(ind)%lval2d)
-      ENDIF
-      
-   ENDDO
+       IF (.NOT. assigned_domain(ind)  .OR. .NOT. is_omp_level_master) CYCLE
+
+       data_type=field(ind)%data_type
+
+       IF (field(ind)%ndim==4) THEN
+          IF (data_type==type_integer) THEN
+             DEALLOCATE(field(ind)%ival4d)
+             IF (field(ind)%ondevice) THEN
+                !$acc exit data delete(field(ind)%ival4d)
+                CONTINUE
+             END IF
+          END IF
+
+          IF (data_type==type_real) THEN
+             DEALLOCATE(field(ind)%rval4d)
+             IF (field(ind)%ondevice) THEN
+                !$acc exit data delete(field(ind)%rval4d)
+                CONTINUE
+             END IF
+          END IF
+
+          IF (data_type==type_logical) THEN
+             DEALLOCATE(field(ind)%lval4d)
+             IF (field(ind)%ondevice) THEN
+                !$acc exit data delete(field(ind)%lval4d)
+                CONTINUE
+             END IF
+          END IF
+
+       ELSE IF (field(ind)%ndim==3) THEN
+          IF (data_type==type_integer) THEN
+             DEALLOCATE(field(ind)%ival3d)
+             IF (field(ind)%ondevice) THEN
+                !$acc exit data delete(field(ind)%ival3d)
+                CONTINUE
+             END IF
+          END IF
+
+          IF (data_type==type_real) THEN
+             DEALLOCATE(field(ind)%rval3d)
+             IF (field(ind)%ondevice) THEN
+                !$acc exit data delete(field(ind)%rval3d)
+                CONTINUE
+             END IF
+          END IF
+
+          IF (data_type==type_logical) THEN
+             DEALLOCATE(field(ind)%lval3d)
+             IF (field(ind)%ondevice) THEN
+                !$acc exit data delete(field(ind)%lval3d)
+                CONTINUE
+             END IF
+          END IF
+
+       ELSE IF (field(ind)%ndim==2) THEN
+          IF (data_type==type_integer) THEN
+             DEALLOCATE(field(ind)%ival2d)
+             IF (field(ind)%ondevice) THEN
+                !$acc exit data delete(field(ind)%ival2d)
+                CONTINUE
+             END IF
+          END IF
+
+          IF (data_type==type_real) THEN
+             DEALLOCATE(field(ind)%rval2d)
+             IF (field(ind)%ondevice) THEN
+                !$acc exit data delete(field(ind)%rval2d)
+                CONTINUE
+             END IF
+          END IF
+
+          IF (data_type==type_logical) THEN
+             DEALLOCATE(field(ind)%lval2d)
+             IF (field(ind)%ondevice) THEN
+                !$acc exit data delete(field(ind)%lval2d)
+                CONTINUE
+             END IF
+          END IF
+
+       ENDIF
+
+    ENDDO
   END SUBROUTINE deallocate_field_
 
@@ -460 +540 @@
   END SUBROUTINE  getval_l4d    
 
+
+  SUBROUTINE update_device_field(field)
+  USE domain_mod
+  USE omp_para
+  IMPLICIT NONE
+    TYPE(t_field) :: field(:)
+    INTEGER :: ind
+
+    DO ind=1,ndomain
+      IF (.NOT. assigned_domain(ind)  .OR. .NOT. is_omp_level_master) CYCLE
+
+      IF (.NOT. field(ind)%ondevice) CALL create_device_field(field(ind))
+
+      IF (field(ind)%ndim==4) THEN
+         IF (field(ind)%data_type==type_integer) THEN
+            !$acc update device(field(ind)%ival4d(:,:,:))
+            CONTINUE
+         END IF
+
+         IF (field(ind)%data_type==type_real) THEN
+            !$acc update device(field(ind)%rval4d(:,:,:))
+            CONTINUE
+         END IF
+
+         IF (field(ind)%data_type==type_logical) THEN
+            !$acc update device(field(ind)%lval4d(:,:,:))
+            CONTINUE
+         END IF
+
+      ELSE IF (field(ind)%ndim==3) THEN
+         IF (field(ind)%data_type==type_integer) THEN
+            !$acc update device(field(ind)%ival3d(:,:))
+            CONTINUE
+         END IF
+
+         IF (field(ind)%data_type==type_real) THEN
+            !$acc update device(field(ind)%rval3d(:,:))
+            CONTINUE
+         END IF
+
+         IF (field(ind)%data_type==type_logical) THEN
+            !$acc update device(field(ind)%lval3d(:,:))
+            CONTINUE
+         END IF
+
+      ELSE IF (field(ind)%ndim==2) THEN
+         IF (field(ind)%data_type==type_integer) THEN
+            !$acc update device(field(ind)%ival2d(:))
+            CONTINUE
+         END IF
+
+         IF (field(ind)%data_type==type_real) THEN
+            !$acc update device(field(ind)%rval2d(:))
+            CONTINUE
+         END IF
+
+         IF (field(ind)%data_type==type_logical) THEN
+            !$acc update device(field(ind)%lval2d(:))
+            CONTINUE
+         END IF
+      ENDIF
+   ENDDO
+   !$OMP BARRIER
+ END SUBROUTINE update_device_field
+ 
+  SUBROUTINE update_host_field(field)
+  USE domain_mod
+  USE omp_para
+  IMPLICIT NONE
+    TYPE(t_field) :: field(:)
+    INTEGER :: ind
+
+    DO ind=1,ndomain
+      IF (.NOT. assigned_domain(ind)  .OR. .NOT. is_omp_level_master) CYCLE
+
+      IF (field(ind)%ondevice) THEN
+        
+         IF (field(ind)%ndim==4) THEN
+            IF (field(ind)%data_type==type_integer) THEN
+               !$acc update host(field(ind)%ival4d(:,:,:)) wait
+               CONTINUE
+            END IF
+
+            IF (field(ind)%data_type==type_real) THEN
+               !$acc update host(field(ind)%rval4d(:,:,:)) wait
+               CONTINUE
+            END IF
+
+            IF (field(ind)%data_type==type_logical) THEN
+               !$acc update host(field(ind)%lval4d(:,:,:)) wait
+               CONTINUE
+            END IF
+         
+         ELSE IF (field(ind)%ndim==3) THEN
+            IF (field(ind)%data_type==type_integer) THEN
+               !$acc update host(field(ind)%ival3d(:,:)) wait
+               CONTINUE
+            END IF
+
+            IF (field(ind)%data_type==type_real) THEN
+               !$acc update host(field(ind)%rval3d(:,:)) wait
+               CONTINUE
+            END IF
+
+            IF (field(ind)%data_type==type_logical) THEN
+               !$acc update host(field(ind)%lval3d(:,:)) wait
+               CONTINUE
+            END IF
+
+         ELSE IF (field(ind)%ndim==2) THEN
+            IF (field(ind)%data_type==type_integer) THEN
+               !$acc update host(field(ind)%ival2d(:)) wait
+               CONTINUE
+            END IF
+
+            IF (field(ind)%data_type==type_real) THEN
+               !$acc update host(field(ind)%rval2d(:)) wait
+               CONTINUE
+            END IF
+
+            IF (field(ind)%data_type==type_logical) THEN
+               !$acc update host(field(ind)%lval2d(:)) wait
+               CONTINUE
+            END IF
+         ENDIF
+      END IF
+   ENDDO
+   !$OMP BARRIER
+ END SUBROUTINE update_host_field
+
+ SUBROUTINE create_device_field(field)
+    TYPE(t_field) :: field
+
+    IF (field%ondevice) THEN
+       PRINT *, "Field is already on device !"
+       STOP 1
+    END IF
+    IF (field%ndim==4) THEN
+       IF (field%data_type==type_integer) THEN
+          !$acc enter data create(field%ival4d(:,:,:))
+       END IF
+
+       IF (field%data_type==type_real) THEN
+          !$acc enter data create(field%rval4d(:,:,:))
+       END IF
+
+       IF (field%data_type==type_logical) THEN
+          !$acc enter data create(field%lval4d(:,:,:))
+       END IF
+
+    ELSE IF (field%ndim==3) THEN
+       IF (field%data_type==type_integer) THEN
+          !$acc enter data create(field%ival3d(:,:))
+       END IF
+
+       IF (field%data_type==type_real) THEN
+          !$acc enter data create(field%rval3d(:,:))
+       END IF
+
+       IF (field%data_type==type_logical) THEN
+          !$acc enter data create(field%lval3d(:,:))
+       END IF
+
+    ELSE IF (field%ndim==2) THEN
+       IF (field%data_type==type_integer) THEN
+          !$acc enter data create(field%ival2d(:))
+       END IF
+
+       IF (field%data_type==type_real) THEN
+          !$acc enter data create(field%rval2d(:))
+       END IF
+
+       IF (field%data_type==type_logical) THEN
+          !$acc enter data create(field%lval2d(:))
+       END IF
+    ENDIF
+    field%ondevice = .TRUE.
+  END SUBROUTINE create_device_field
+ 
 END MODULE field_mod   

codes/icosagcm/trunk/src/diagnostics/check_conserve.f90

@@ -1 +1 @@
 MODULE check_conserve_mod
   USE icosa 
+  USE abort_mod
   IMPLICIT NONE 
 
@@ -28 +29 @@
     USE getin_mod
     USE omp_para, ONLY : is_master
+    USE abort_mod
     CHARACTER(LEN=255) :: check_type_str
     CALL allocate_field(f_pk,field_t,type_real,llm)
@@ -47 +49 @@
        STOP
     END SELECT
+
+    IF (check_type /= check_basic) THEN
+       CALL abort_acc("check_conservation /= 'basic'")
+    END IF
   END SUBROUTINE init_check_conserve
 
@@ -179 +185 @@
     
     IF(check_type == check_detailed) THEN
-
+       CALL abort_acc("!check_detailed")
        CALL transfert_request(f_ue,req_e1_vect)
        CALL pression(f_ps,f_p)

codes/icosagcm/trunk/src/diagnostics/diagflux.F90

@@ -21 +21 @@
   SUBROUTINE init_diagflux
     USE getin_mod
+    USE abort_mod
     INTEGER :: ll
     diagflux_on = .FALSE.
     CALL getin("diagflux", diagflux_on)
+    IF (diagflux_on) THEN
+       CALL abort_acc("diagflux /= .FALSE.")
+    END IF
+
     ll = MERGE(llm,1,diagflux_on)
     CALL allocate_field(f_masst,         field_t,type_real,ll,       name="masst")

codes/icosagcm/trunk/src/diagnostics/observable.f90

@@ -30 +30 @@
     CALL allocate_field(f_buf_s,   field_t,type_real, name="buf_s")
 
-    CALL allocate_field(f_theta, field_t,type_real,llm,nqdyn,  name='theta') ! potential temperature
+    CALL allocate_field(f_theta, field_t,type_real,llm,nqdyn,  name='theta', ondevice=.TRUE.) ! potential temperature
     CALL allocate_field(f_pmid,  field_t,type_real,llm,  name='pmid')       ! mid layer pressure
   END SUBROUTINE init_observable
@@ -56 +56 @@
 
     CALL transfert_request(f_ps,req_i1)
+    CALL update_host_field(f_ps)
     
     IF(init) THEN
@@ -109 +110 @@
     
     CALL progonostic_vel_to_horiz(f_geopot, f_ps, f_mass, f_u, f_W, f_buf_uh, f_buf_i)
-    CALL transfert_request(f_buf_uh,req_e1_vect) 
+    CALL transfert_request(f_buf_uh,req_e1_vect)
     CALL un2ulonlat(f_buf_uh, f_buf_ulon, f_buf_ulat)
     IF(init) THEN

codes/icosagcm/trunk/src/dissip/dissip_gcm.F90

@@ -1 +1 @@
 MODULE dissip_gcm_mod
   USE icosa
-
+  USE abort_mod
   PRIVATE
 
@@ -48 +48 @@
   USE icosa
   IMPLICIT NONE  
-    CALL allocate_field(f_due_diss1,field_u,type_real,llm)
-    CALL allocate_field(f_due_diss2,field_u,type_real,llm)
+    CALL allocate_field(f_due_diss1,field_u,type_real,llm,ondevice=.TRUE.)
+    CALL allocate_field(f_due_diss2,field_u,type_real,llm,ondevice=.TRUE.)
     CALL allocate_field(f_dtheta_diss,field_t,type_real,llm)
-    CALL allocate_field(f_dtheta_rhodz_diss,field_t,type_real,llm)
+    CALL allocate_field(f_dtheta_rhodz_diss,field_t,type_real,llm,ondevice=.TRUE.)
     ALLOCATE(tau_graddiv(llm))
     ALLOCATE(tau_gradrot(llm))    
     ALLOCATE(tau_divgrad(llm))
+    !$acc enter data create(tau_graddiv(:),tau_gradrot(:),tau_divgrad(:)) async
   END SUBROUTINE allocate_dissip
  
@@ -66 +67 @@
   USE transfert_omp_mod
   USE omp_para
+  USE abort_mod
   IMPLICIT NONE
   
@@ -98 +100 @@
       IF (is_master) PRINT *, 'No Rayleigh friction'
    CASE('dcmip2_schaer_noshear')
+      CALL abort_acc("rayleigh_friction_type /= 'none'")
       rayleigh_friction_type=1
       rayleigh_shear=0
       IF (is_master) PRINT *, 'Rayleigh friction : Schaer-like mountain without shear DCMIP2.1'
    CASE('dcmip2_schaer_shear')
+      CALL abort_acc("rayleigh_friction_type /= 'none'")
       rayleigh_shear=1
       rayleigh_friction_type=2
       IF (is_master) PRINT *, 'Rayleigh friction : Schaer-like mountain with shear DCMIP2.2'
    CASE('giant_liu_schneider') 
+      CALL abort_acc("rayleigh_friction_type /= 'none'")
       rayleigh_friction_type=99
       IF (is_master) PRINT *, 'Rayleigh friction : giant planets Liu Schneider 2010'
@@ -130 +135 @@
 
     CALL allocate_dissip
-    CALL allocate_field(f_u,field_u,type_real)
-    CALL allocate_field(f_du,field_u,type_real)
-    CALL allocate_field(f_theta,field_t,type_real)
-    CALL allocate_field(f_dtheta,field_t,type_real)
+    CALL allocate_field(f_u,field_u,type_real,ondevice=.TRUE.)
+    CALL allocate_field(f_du,field_u,type_real,ondevice=.TRUE.)
+    CALL allocate_field(f_theta,field_t,type_real,ondevice=.TRUE.)
+    CALL allocate_field(f_dtheta,field_t,type_real,ondevice=.TRUE.)
     
     CALL init_message(f_due_diss1,req_e1_vect,req_due)
@@ -173 +178 @@
       CALL RANDOM_SEED(put=(/(i,i=1,M)/))
 
+      ! This cannot be ported on GPU due to compiler limitations
       DO j=jj_begin,jj_end
         DO i=ii_begin,ii_end
@@ -192 +198 @@
       dumax=0
       DO iter=1,nitergdiv
+        CALL update_device_field(f_u)
         CALL transfert_request(f_u,req_e1_vect)
+         
         DO ind=1,ndomain
           IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE
@@ -200 +208 @@
           du=f_du(ind)
           CALL compute_gradiv_inplace(u,1,1)
+          ! This should be ported on GPU but we are running into compiler issues...
+          !$acc update host(u(:)) wait
           du=u
         ENDDO
       ENDDO
 
+      CALL update_device_field(f_du)
       CALL transfert_request(f_du,req_e1_vect)
-      
+      CALL update_host_field(f_du)
+
       DO ind=1,ndomain
         IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE
         CALL swap_dimensions(ind)
         CALL swap_geometry(ind)
-        u=f_u(ind)
         du=f_du(ind)
          
+        ! Not ported on GPU because all the other kernels cannot be ported
         DO j=jj_begin,jj_end
           DO i=ii_begin,ii_end
@@ -240 +252 @@
         u=f_u(ind)
         du=f_du(ind)
+        ! This should be ported on GPU but we are running into compiler issues...
         u=du/dumax
       ENDDO
@@ -265 +278 @@
       CALL RANDOM_SEED(put=(/(i,i=1,M)/))
  
+      ! This cannot be ported on GPU due to compiler limitations
        DO j=jj_begin,jj_end
         DO i=ii_begin,ii_end
@@ -280 +294 @@
 !$OMP BARRIER
 
-
     DO it=1,20
  
       dumax=0
       DO iter=1,nitergrot
+        CALL update_device_field(f_u)
         CALL transfert_request(f_u,req_e1_vect)
         DO ind=1,ndomain
@@ -293 +307 @@
           du=f_du(ind)
           CALL compute_gradrot_inplace(u,1,1)
+          ! This should be ported on GPU but we are running into compiler issues...
+          !$acc update host(u(:)) wait
           du=u
         ENDDO
       ENDDO
-
+      
+      CALL update_device_field(f_du)
       CALL transfert_request(f_du,req_e1_vect)
+      CALL update_host_field(f_du)
       
       DO ind=1,ndomain
@@ -303 +321 @@
         CALL swap_dimensions(ind)
         CALL swap_geometry(ind)
-        u=f_u(ind)
         du=f_du(ind)
         
+        ! Not ported on GPU because all the other kernels cannot be ported
         DO j=jj_begin,jj_end
           DO i=ii_begin,ii_end
@@ -334 +352 @@
         u=f_u(ind)
         du=f_du(ind)
+        ! This should be ported on GPU but we are running into compiler issues...
         u=du/dumax
       ENDDO
@@ -358 +377 @@
       CALL RANDOM_SEED(put=(/(i,i=1,M)/))
 
+      ! This cannot be ported on GPU due to compiler limitations
       DO j=jj_begin,jj_end
         DO i=ii_begin,ii_end
@@ -373 +393 @@
       dthetamax=0
       DO iter=1,niterdivgrad
+        CALL update_device_field(f_theta)
         CALL transfert_request(f_theta,req_i1)
         DO ind=1,ndomain
@@ -381 +402 @@
           dtheta=f_dtheta(ind)
           CALL compute_divgrad_inplace(theta,1,1)
+          ! This should be ported on GPU but we are running into compiler issues...
+          !$acc update host(theta(:)) wait
           dtheta=theta
         ENDDO
       ENDDO
 
+      CALL update_device_field(f_dtheta)
       CALL transfert_request(f_dtheta,req_i1)
+      CALL update_host_field(f_dtheta)
       
       DO ind=1,ndomain
@@ -391 +416 @@
         CALL swap_dimensions(ind)
         CALL swap_geometry(ind)
-        theta=f_theta(ind)
         dtheta=f_dtheta(ind)
         
+        ! Not ported on GPU because all the other kernels cannot be ported
         DO j=jj_begin,jj_end
           DO i=ii_begin,ii_end
@@ -421 +446 @@
         theta=f_theta(ind)
         dtheta=f_dtheta(ind)
+        ! This should be ported on GPU but we are running into compiler issues...
         theta=dtheta/dthetamax
       ENDDO
@@ -489 +515 @@
     ENDIF
 
+    !$acc update device(tau_graddiv(:),tau_gradrot(:),tau_divgrad(:)) async
+
   END SUBROUTINE init_dissip 
  
@@ -501 +529 @@
   USE time_mod
   USE omp_para
+  USE abort_mod
   IMPLICIT NONE
     TYPE(t_field),POINTER :: f_ps(:), f_mass(:), f_phis(:), f_geopot(:)
@@ -534 +563 @@
       dtheta_rhodz_diss=f_dtheta_rhodz_diss(ind)
 
+      !$acc parallel loop collapse(2) present(due(:,:), dtheta_rhodz(:,:,:), due_diss1(:,:), due_diss2(:,:), dtheta_rhodz_diss(:,:), tau_graddiv(:), tau_gradrot(:), tau_divgrad(:)) async
       DO l=ll_begin,ll_end
 !DIR$ SIMD
@@ -545 +575 @@
         ENDDO
       ENDDO
+      !$acc end parallel loop
 
 !      dtheta_rhodz=0
@@ -550 +581 @@
 
       IF(rayleigh_friction_type>0) THEN
+        CALL abort_acc("dissip/rayleigh_friction_type>0")
        IF(rayleigh_friction_type<99) THEN
          phi=f_geopot(ind)
@@ -562 +594 @@
        ELSE
          ue=f_ue(ind)
+            !$acc parallel loop present(ue(:,:), due(:,:), lat_e(:))
             DO ij=ij_begin,ij_end
               nn = ij+u_right
@@ -577 +610 @@
               ENDIF
             ENDDO
+            !$acc end parallel loop
        ENDIF
       END IF
    END DO
-
    CALL trace_end("dissip")
 
-!   CALL write_dissip_tendencies
+   !CALL write_dissip_tendencies
 !$OMP BARRIER
 
@@ -653 +686 @@
       ue=f_ue(ind)
       due=f_due(ind) 
+      !$acc parallel loop present(ue(:,:), due(:,:)) async
       DO  l = ll_begin, ll_end
+        !$acc loop
 !DIR$ SIMD
         DO ij=ij_begin,ij_end
@@ -664 +699 @@
 
     DO it=1,nitergdiv
-        
       CALL send_message(f_due,req_due)
       CALL wait_message(req_due)
-        
+
       DO ind=1,ndomain
         IF (.NOT. assigned_domain(ind)) CYCLE
         CALL swap_dimensions(ind)
         CALL swap_geometry(ind)
-        due=f_due(ind) 
+        due=f_due(ind)
         CALL compute_gradiv_inplace(due,ll_begin,ll_end)
       ENDDO
@@ -702 +736 @@
       ue=f_ue(ind)
       due=f_due(ind) 
+      !$acc parallel loop present(ue(:,:), due(:,:)) async
       DO  l = ll_begin, ll_end
+         !$acc loop
 !DIR$ SIMD
         DO ij=ij_begin,ij_end
@@ -713 +749 @@
 
     DO it=1,nitergrot
-        
       CALL send_message(f_due,req_due)
       CALL wait_message(req_due)
@@ -740 +775 @@
     REAL(rstd),POINTER    :: theta(:,:)
     REAL(rstd),POINTER    :: dtheta(:,:)
-    INTEGER :: ind
+    INTEGER :: ind, l, ij
     INTEGER :: it
 
@@ -751 +786 @@
       theta=f_theta(ind)
       dtheta=f_dtheta(ind) 
-      dtheta=theta
+      ! Replace Fortran 90 construct dtheta=theta because it confuses PGI acc kernels...
+      !$acc parallel loop collapse(2) present(theta(:,:), dtheta(:,:)) async
+      DO l=ll_begin,ll_end
+!DIR$ SIMD
+        DO ij=1,iim*jjm
+          dtheta(ij,l)=theta(ij,l)
+        ENDDO
+      ENDDO
+      !$acc end parallel loop
     ENDDO
 
     DO it=1,niterdivgrad
-        
       CALL transfert_request(f_dtheta,req_i1)
         
@@ -797 +839 @@
       theta_rhodz=f_theta_rhodz(ind)
       dtheta_rhodz=f_dtheta_rhodz(ind) 
+      !$acc parallel loop present(mass(:,:), theta_rhodz(:,:,:), dtheta_rhodz(:,:)) async
       DO  l = ll_begin, ll_end
+        !$acc loop
 !DIR$ SIMD
         DO ij=ij_begin,ij_end
@@ -806 +850 @@
 
     DO it=1,niterdivgrad
-        
       CALL send_message(f_dtheta_rhodz,req_dtheta)
       CALL wait_message(req_dtheta)
@@ -826 +869 @@
       mass=f_mass(ind)
     
+      !$acc parallel loop collapse(2) present(mass(:,:), dtheta_rhodz(:,:)) async
       DO  l = ll_begin, ll_end
 !DIR$ SIMD
@@ -832 +876 @@
         ENDDO
       ENDDO
+      !$acc end parallel loop
     ENDDO
 
@@ -837 +882 @@
     CALL trace_end("divgrad")
 
-  END SUBROUTINE divgrad_theta_rhodz  
-  
+  END SUBROUTINE divgrad_theta_rhodz
+
   SUBROUTINE compute_gradiv(ue,gradivu_e,llb,lle)
-  USE icosa
-  IMPLICIT NONE
     INTEGER,INTENT(IN)     :: llb
     INTEGER,INTENT(IN)     :: lle
+    REAL(rstd),INTENT(OUT) :: gradivu_e(iim*3*jjm,llm)
     REAL(rstd),INTENT(IN)  :: ue(iim*3*jjm,llm)
-    REAL(rstd),INTENT(OUT) :: gradivu_e(iim*3*jjm,llm)
-    
+
     gradivu_e = ue
-    CALL compute_gradiv_inplace(gradivu_e, llb, lle)
+    CALL compute_gradiv_inplace(gradivu_e,llb,lle)
+
   END SUBROUTINE compute_gradiv
   
   SUBROUTINE compute_gradiv_inplace(ue_gradivu_e,llb,lle)
-  USE icosa
-  IMPLICIT NONE
+    USE geometry, ONLY : Ai, ne, le, de
     INTEGER,INTENT(IN)     :: llb
     INTEGER,INTENT(IN)     :: lle
@@ -861 +904 @@
     INTEGER :: ij,l
 
+    ! ue and gradivu_e second dimension is not always llm so use the bounds explicitly
+    !$acc data present( ue_gradivu_e(:,llb:lle), Ai(:), ne(:,:), le(:), de(:)) create(divu_i(:,:)) async
+
+    !$acc parallel loop collapse(2) async
     DO l=llb,lle
       !DIR$ SIMD
@@ -872 +919 @@
       ENDDO
     ENDDO
+    !$acc end parallel loop
     
+    !$acc parallel loop collapse(2) async
     DO l=llb,lle
       !DIR$ SIMD
@@ -881 +930 @@
       ENDDO
     ENDDO
-
+    !$acc end parallel loop
+
+    !$acc parallel loop collapse(2) async
     DO l=llb,lle
       !DIR$ SIMD
@@ -889 +940 @@
           ue_gradivu_e(ij+u_ldown,l)=-ue_gradivu_e(ij+u_ldown,l)*cgraddiv
       ENDDO
-    ENDDO   
+    ENDDO
+    !$acc end parallel loop
+    !$acc end data   
   END SUBROUTINE compute_gradiv_inplace
-  
+
   SUBROUTINE compute_divgrad(theta,divgrad_i,llb,lle)
-  USE icosa
-  IMPLICIT NONE
     INTEGER,INTENT(IN)     :: llb
     INTEGER,INTENT(IN)     :: lle
     REAL(rstd),INTENT(IN) :: theta(iim*jjm,1:lle)
     REAL(rstd),INTENT(OUT) :: divgrad_i(iim*jjm,1:lle)
-    
+
     divgrad_i = theta
-    CALL compute_divgrad_inplace(divgrad_i, llb, lle)
-  END SUBROUTINE
+    CALL compute_divgrad_inplace(divgrad_i,llb,lle)
+  END SUBROUTINE compute_divgrad
   
   SUBROUTINE compute_divgrad_inplace(theta_divgrad_i,llb,lle)
-  USE icosa
-  IMPLICIT NONE
+    USE geometry, ONLY : Ai, ne, le, de
     INTEGER,INTENT(IN)     :: llb
     INTEGER,INTENT(IN)     :: lle
@@ -912 +962 @@
     REAL(rstd)  :: grad_e(3*iim*jjm,llb:lle)
 
-    INTEGER :: ij,l    
+    INTEGER :: ij,l
+
+    ! theta and divgrad_i second dimension is not always llm so use the bounds explicitly
+    !$acc data present(theta_divgrad_i(:,llb:lle), Ai(:), de(:), ne(:,:), le(:)) create(grad_e(:,:)) async
+       
+    !$acc parallel loop collapse(2) async
     DO l=llb,lle
       !DIR$ SIMD
@@ -920 +975 @@
           grad_e(ij+u_ldown,l)=-1/de(ij+u_ldown)*(ne(ij,ldown)*theta_divgrad_i(ij,l)+ne(ij+t_ldown,rup)*theta_divgrad_i(ij+t_ldown,l) )
       ENDDO
-    ENDDO    
+    ENDDO
+    !$acc end parallel loop
     
+    
+    !$acc parallel loop collapse(2) async
     DO l=llb,lle
       !DIR$ SIMD
@@ -933 +991 @@
       ENDDO
     ENDDO
+    !$acc end parallel loop
    
+    !$acc parallel loop collapse(2) async
     DO l=llb,lle
       DO ij=ij_begin,ij_end
@@ -939 +999 @@
       ENDDO
     ENDDO
+    !$acc end parallel loop
+    !$acc end data
   END SUBROUTINE compute_divgrad_inplace
 
@@ -946 +1008 @@
     REAL(rstd),INTENT(IN) :: ue(iim*3*jjm,lle)
     REAL(rstd),INTENT(OUT) :: gradrot_e(iim*3*jjm,lle)
-    
+
     gradrot_e = ue
     CALL compute_gradrot_inplace(gradrot_e,llb,lle)
-  END SUBROUTINE
-    
+  END SUBROUTINE compute_gradrot
+
   SUBROUTINE compute_gradrot_inplace(ue_gradrot_e,llb,lle)
-  USE icosa
-  IMPLICIT NONE
+    USE geometry, ONLY : Av, ne, le, de
     INTEGER,INTENT(IN)     :: llb
     INTEGER,INTENT(IN)     :: lle
@@ -961 +1022 @@
     INTEGER :: ij,l
 
+    ! ue and gradrot_e second dimension is not always llm so use the bounds explicitly
+    ! gradrot_e should be copyout but using copy instead allows to compare the output
+    ! more easily as the code sometimes uses unintialed values
+    !$acc data present(ue_gradrot_e(:,llb:lle), Av(:), ne(:,:), de(:), le(:)) create(rot_v(:,:)) async
+     
+    !$acc parallel loop collapse(2) async
     DO l=llb,lle
 !DIR$ SIMD
@@ -972 +1039 @@
       ENDDO
     ENDDO                              
-  
+    !$acc end parallel loop
+  
+    !$acc parallel loop collapse(2) async
     DO l=llb,lle
 !DIR$ SIMD
@@ -981 +1050 @@
       ENDDO
     ENDDO
-
+    !$acc end parallel loop
+
+    !$acc parallel loop collapse(2) async
     DO l=llb,lle
 !DIR$ SIMD
@@ -989 +1060 @@
           ue_gradrot_e(ij+u_ldown,l)=-ue_gradrot_e(ij+u_ldown,l)*cgradrot        
       ENDDO
-    ENDDO  
+    ENDDO
+    !$acc end parallel loop
+    !$acc end data   
   END SUBROUTINE compute_gradrot_inplace
 
 
 END MODULE dissip_gcm_mod
-       

codes/icosagcm/trunk/src/dissip/sponge.f90

@@ -24 +24 @@
   USE omp_para
   USE mpipara, ONLY: is_mpi_master
+  USE abort_mod
   IMPLICIT NONE
     INTEGER :: l
@@ -43 +44 @@
       RETURN
     ENDIF
+
+    IF (iflag_sponge > 0) THEN
+       CALL abort_acc("iflag_sponge > 0")
+    END IF
 
 !$OMP MASTER

codes/icosagcm/trunk/src/dynamics/caldyn_gcm.F90

@@ -15 +15 @@
   
   SUBROUTINE init_caldyn
+    USE abort_mod
     CHARACTER(len=255) :: def
     INTEGER            :: ind
@@ -21 +22 @@
     hydrostatic=.TRUE.
     CALL getin("hydrostatic",hydrostatic)
-  
+
+    IF (.NOT. hydrostatic) THEN
+       CALL abort_acc("hydrostatic /= .TRUE.")
+    END IF
+
     dysl=.NOT.hydrostatic ! dysl defaults to .FALSE. if hydrostatic, .TRUE. if NH                                              
     CALL getin("dysl",dysl)
@@ -34 +39 @@
     dysl_caldyn_coriolis=dysl
     CALL getin("dysl_caldyn_coriolis",dysl_caldyn_coriolis)
+
+    ! dysl is not supported with openACC
+    IF (dysl_geopot .OR. dysl_caldyn_fast .OR. dysl_slow_hydro .OR. dysl_pvort_only .OR. dysl_caldyn_coriolis) THEN
+       CALL abort_acc("dysl not supported")
+    END IF
 
     def='energy'
@@ -122 +132 @@
   SUBROUTINE allocate_caldyn
     CALL allocate_field(f_out_u,field_u,type_real,llm) 
-    CALL allocate_field(f_qu,field_u,type_real,llm) 
-    CALL allocate_field(f_qv,field_z,type_real,llm) 
-    CALL allocate_field(f_pk,    field_t,type_real,llm,  name='pk')
-    CALL allocate_field(f_wwuu,  field_u,type_real,llm+1,name='wwuu')    
+    CALL allocate_field(f_qu,field_u,type_real,llm, ondevice=.TRUE.) 
+    CALL allocate_field(f_qv,field_z,type_real,llm, ondevice=.TRUE.) 
+    CALL allocate_field(f_pk,    field_t,type_real,llm,  name='pk', ondevice=.TRUE.)
+    CALL allocate_field(f_wwuu,  field_u,type_real,llm+1,name='wwuu', ondevice=.TRUE.)    
     CALL allocate_field(f_planetvel,  field_u,type_real, name='planetvel') ! planetary velocity at r=a
     IF(.NOT.hydrostatic) THEN

codes/icosagcm/trunk/src/dynamics/caldyn_hevi.F90

@@ -26 +26 @@
     USE output_field_mod
     USE checksum_mod
+    USE abort_mod
     IMPLICIT NONE
     LOGICAL,INTENT(IN)    :: write_out
@@ -86 +87 @@
        CALL wait_message(req_ps) ! COM00
     ELSE
+       CALL abort_acc("HEVI_scheme/!eta_mass")
        CALL send_message(f_mass,req_mass) ! COM00
        CALL wait_message(req_mass) ! COM00
@@ -93 +95 @@
 
     IF(.NOT.hydrostatic) THEN
+       CALL abort_acc("HEVI_scheme/!hydrostatic")
        CALL send_message(f_geopot,req_geopot) ! COM03
        CALL wait_message(req_geopot) ! COM03
@@ -112 +115 @@
        du=f_du_fast(ind)
        IF(hydrostatic) THEN
-          du(:,:)=0.
+          !$acc kernels present(du) async
+          du(:,:)=0.0d0
+          !$acc end kernels
           CALL compute_geopot(mass,theta, ps,pk,geopot)
        ELSE
+          CALL abort_acc("HEVI_scheme/!hydrostatic")
           phis = f_phis(ind)
           W = f_W(ind)
@@ -161 +167 @@
           CALL compute_caldyn_slow_hydro(u,mass,hflux,du, .TRUE.)
        ELSE
+          CALL abort_acc("HEVI_scheme/!hydrostatic")
           W = f_W(ind)
           dW = f_dW_slow(ind)
@@ -170 +177 @@
           CALL compute_caldyn_slow_NH(u,mass,geopot,W, F_el,gradPhi2,w_il, hflux,du,dPhi,dW)
        END IF
-       CALL compute_caldyn_Coriolis(hflux,theta,qu, convm,dtheta_rhodz,du)
+       CALL compute_caldyn_Coriolis(hflux,theta,qu,convm,dtheta_rhodz,du)
+       
        IF(caldyn_eta==eta_mass) THEN
           wflux=f_wflux(ind)
@@ -177 +185 @@
           CALL compute_caldyn_vert(u,theta,mass,convm, wflux,wwuu, dps, dtheta_rhodz, du)
           IF(.NOT.hydrostatic) THEN
+             CALL abort_acc("HEVI_scheme/!hydrostatic")
              W_etadot=f_Wetadot(ind)
              CALL compute_caldyn_vert_NH(mass,geopot,W,wflux, W_etadot, du,dPhi,dW)

codes/icosagcm/trunk/src/dynamics/caldyn_kernels_base.F90

@@ -5 +5 @@
   USE omp_para
   USE trace
+  USE abort_mod
   IMPLICIT NONE
   PRIVATE
@@ -25 +26 @@
 
   !**************************** Geopotential *****************************
-
   SUBROUTINE compute_geopot(rhodz,theta, ps,pk,geopot) 
     REAL(rstd),INTENT(IN)    :: rhodz(iim*jjm,llm)
@@ -46 +46 @@
 
     IF(dysl_geopot) THEN
+      CALL abort_acc("HEVI_scheme/!dysl_geopot")
 #include "../kernels/compute_geopot.k90"
     ELSE
@@ -52 +53 @@
     ! Works also when caldyn_eta=eta_mass          
 
+    !$acc data present(rhodz(:,:), ps(:), pk(:,:), geopot(:,:), theta(:,:,:)) async
+
     IF(boussinesq) THEN ! compute geopotential and pk=Lagrange multiplier
+       CALL abort_acc("HEVI_scheme/boussinesq")
        ! specific volume 1 = dphi/g/rhodz
        !         IF (is_omp_level_master) THEN ! no openMP on vertical due to dependency
        DO l = 1,llm
+          !$acc parallel loop
           !DIR$ SIMD
           DO ij=ij_omp_begin_ext,ij_omp_end_ext         
              geopot(ij,l+1) = geopot(ij,l) + g*rhodz(ij,l)
           ENDDO
+          !$acc end parallel loop
        ENDDO
        ! use hydrostatic balance with theta*rhodz to find pk (Lagrange multiplier=pressure) 
        ! uppermost layer
+       !$acc parallel loop
        !DIR$ SIMD
        DO ij=ij_begin_ext,ij_end_ext         
           pk(ij,llm) = ptop + (.5*g)*theta(ij,llm,1)*rhodz(ij,llm)
        END DO
+       !$acc end parallel loop
        ! other layers
        DO l = llm-1, 1, -1
           !          !$OMP DO SCHEDULE(STATIC) 
+          !$acc parallel loop
           !DIR$ SIMD
           DO ij=ij_begin_ext,ij_end_ext         
              pk(ij,l) = pk(ij,l+1) + (.5*g)*(theta(ij,l,1)*rhodz(ij,l)+theta(ij,l+1,1)*rhodz(ij,l+1))
           END DO
+          !$acc end parallel loop
        END DO
        ! now pk contains the Lagrange multiplier (pressure)
@@ -81 +91 @@
        SELECT CASE(caldyn_thermo)
           CASE(thermo_theta, thermo_entropy)
+             !$acc parallel loop async
              !DIR$ SIMD
              DO ij=ij_omp_begin_ext,ij_omp_end_ext
                 pk(ij,llm) = ptop + (.5*g)*rhodz(ij,llm)
              END DO
+
              ! other layers
+             ! We use kernels here instead of "loop seq" + "loop gang vector"
+             ! to be sure the code really abides the standard. In practice,
+             ! it seems like the compiler interchanges the loops.
+             !$acc kernels async
              DO l = llm-1, 1, -1
                 !DIR$ SIMD
@@ -92 +108 @@
                 END DO
              END DO
+             !$acc end kernels
              ! surface pressure (for diagnostics)
              IF(caldyn_eta==eta_lag) THEN
+                !$acc parallel loop async
                 DO ij=ij_omp_begin_ext,ij_omp_end_ext         
                    ps(ij) = pk(ij,1) + (.5*g)*rhodz(ij,1)
@@ -99 +117 @@
              END IF
           CASE(thermo_moist) ! theta(ij,l,2) = qv = mv/md
+             CALL abort_acc("HEVI_scheme/thermo_moist")
+             !$acc parallel loop
              !DIR$ SIMD
              DO ij=ij_omp_begin_ext,ij_omp_end_ext
                 pk(ij,llm) = ptop + (.5*g)*rhodz(ij,llm)*(1.+theta(ij,l,2))
              END DO
+             !$acc end parallel loop
+
              ! other layers
              DO l = llm-1, 1, -1
+                !$acc parallel loop
                 !DIR$ SIMD
                 DO ij=ij_omp_begin_ext,ij_omp_end_ext         
@@ -111 +134 @@
                         rhodz(ij,l+1)*(1.+theta(ij,l+1,2)) )
                 END DO
+                !$acc end parallel loop
              END DO
              ! surface pressure (for diagnostics)
              IF(caldyn_eta==eta_lag) THEN
+                !$acc parallel loop
                 DO ij=ij_omp_begin_ext,ij_omp_end_ext         
                    ps(ij) = pk(ij,1) + (.5*g)*rhodz(ij,1)*(1.+theta(ij,l,2))
                 END DO
+                !$acc end parallel loop
              END IF
           END SELECT
 
-       DO l = 1,llm
-          SELECT CASE(caldyn_thermo)
+        SELECT CASE(caldyn_thermo)
           CASE(thermo_theta)
-             !DIR$ SIMD
-             DO ij=ij_omp_begin_ext,ij_omp_end_ext
-                p_ik = pk(ij,l)
-                exner_ik = cpp * (p_ik/preff) ** kappa
-                pk(ij,l) = exner_ik
-                ! specific volume v = kappa*theta*pi/p = dphi/g/rhodz
-                geopot(ij,l+1) = geopot(ij,l) + (g*kappa)*rhodz(ij,l)*theta(ij,l,1)*exner_ik/p_ik
-             ENDDO
+            ! We use kernels here instead of "loop seq" + "loop gang vector"
+            ! to be sure the code really abides the standard. In practice,
+            ! it seems like the compiler interchanges the loops.
+            !$acc kernels async
+            DO l = 1,llm
+               !DIR$ SIMD
+               DO ij=ij_omp_begin_ext,ij_omp_end_ext
+                  p_ik = pk(ij,l)
+                  exner_ik = cpp * (p_ik/preff) ** kappa
+                  pk(ij,l) = exner_ik
+                  ! specific volume v = kappa*theta*pi/p = dphi/g/rhodz
+                  geopot(ij,l+1) = geopot(ij,l) + (g*kappa)*rhodz(ij,l)*theta(ij,l,1)*exner_ik/p_ik
+               ENDDO
+            ENDDO
+            !$acc end kernels
+
           CASE(thermo_entropy) ! theta is in fact entropy = cpp*log(theta/Treff) = cpp*log(T/Treff) - Rd*log(p/preff)
-             !DIR$ SIMD
-             DO ij=ij_omp_begin_ext,ij_omp_end_ext
-                p_ik = pk(ij,l)
-                temp_ik = Treff*exp((theta(ij,l,1) + Rd*log(p_ik/preff))/cpp)
-                pk(ij,l) = temp_ik
-                ! specific volume v = Rd*T/p = dphi/g/rhodz
-                geopot(ij,l+1) = geopot(ij,l) + (g*Rd)*rhodz(ij,l)*temp_ik/p_ik
+             CALL abort_acc("HEVI_scheme/thermo_entropy")
+             DO l = 1,llm
+               !$acc parallel loop
+               !DIR$ SIMD
+               DO ij=ij_omp_begin_ext,ij_omp_end_ext
+                  p_ik = pk(ij,l)
+                  temp_ik = Treff*exp((theta(ij,l,1) + Rd*log(p_ik/preff))/cpp)
+                  pk(ij,l) = temp_ik
+                  ! specific volume v = Rd*T/p = dphi/g/rhodz
+                  geopot(ij,l+1) = geopot(ij,l) + (g*Rd)*rhodz(ij,l)*temp_ik/p_ik
+               ENDDO
+               !$acc end parallel loop
              ENDDO
           CASE(thermo_moist) ! theta is moist pseudo-entropy per dry air mass
-             DO ij=ij_omp_begin_ext,ij_omp_end_ext
-                p_ik = pk(ij,l)
-                qv = theta(ij,l,2) ! water vaper mixing ratio = mv/md
-                Rmix = Rd+qv*Rv
-                chi = ( theta(ij,l,1) + Rmix*log(p_ik/preff) ) / (cpp + qv*cppv) ! log(T/Treff)
-                temp_ik = Treff*exp(chi)
-                pk(ij,l) = temp_ik
-                ! specific volume v = R*T/p = dphi/g/rhodz
-                ! R = (Rd + qv.Rv)/(1+qv)
-                geopot(ij,l+1) = geopot(ij,l) + g*Rmix*rhodz(ij,l)*temp_ik/(p_ik*(1+qv))
+             CALL abort_acc("HEVI_scheme/thermo_moist")
+             DO l = 1,llm
+               !$acc parallel loop
+               DO ij=ij_omp_begin_ext,ij_omp_end_ext
+                  p_ik = pk(ij,l)
+                  qv = theta(ij,l,2) ! water vaper mixing ratio = mv/md
+                  Rmix = Rd+qv*Rv
+                  chi = ( theta(ij,l,1) + Rmix*log(p_ik/preff) ) / (cpp + qv*cppv) ! log(T/Treff)
+                  temp_ik = Treff*exp(chi)
+                  pk(ij,l) = temp_ik
+                  ! specific volume v = R*T/p = dphi/g/rhodz
+                  ! R = (Rd + qv.Rv)/(1+qv)
+                  geopot(ij,l+1) = geopot(ij,l) + g*Rmix*rhodz(ij,l)*temp_ik/(p_ik*(1+qv))
+               ENDDO
+               !$acc end parallel loop
              ENDDO
           CASE DEFAULT
              STOP
-          END SELECT
-       ENDDO
+        END SELECT
     END IF
-
+    !$acc end data
     END IF ! dysl
 
@@ -167 +209 @@
   END SUBROUTINE compute_geopot
 
-  SUBROUTINE compute_caldyn_vert(u,theta,rhodz,convm, wflux,wwuu, dps,dtheta_rhodz,du)
+  SUBROUTINE compute_caldyn_vert(u, theta, rhodz, convm, wflux, wwuu, dps, dtheta_rhodz, du)
+    USE disvert_mod, ONLY : bp
     REAL(rstd),INTENT(IN)  :: u(iim*3*jjm,llm)
     REAL(rstd),INTENT(IN)  :: theta(iim*jjm,llm,nqdyn)
@@ -193 +236 @@
     CALL trace_start("compute_caldyn_vert")
 
+    !$acc data async &
+    !$acc present(rhodz(:,:), u(:,:), wwuu(:,:), wflux(:,:), dps(:), convm(:,:), du(:,:), dtheta_rhodz(:,:,:), theta(:,:,:), bp(:))
+
+
     !$OMP BARRIER   
 !!! cumulate mass flux convergence from top to bottom
     !  IF (is_omp_level_master) THEN
+    ! We use kernels here instead of "loop seq" + "loop gang vector"
+    ! to be sure the code really abides the standard. In practice,
+    ! it seems like the compiler interchanges the loops.
+    !$acc kernels async
     DO  l = llm-1, 1, -1
        !    IF (caldyn_conserv==energy) CALL test_message(req_qu) 
@@ -205 +256 @@
        ENDDO
     ENDDO
+    !$acc end kernels
     !  ENDIF
 
@@ -213 +265 @@
     ! compute dps
     IF (is_omp_first_level) THEN
+       !$acc parallel loop  async
        !DIR$ SIMD
        DO ij=ij_begin,ij_end         
@@ -221 +274 @@
 
 !!! Compute vertical mass flux (l=1,llm+1 done by caldyn_BC)
+    !$acc parallel loop collapse(2) async
     DO l=ll_beginp1,ll_end
        !    IF (caldyn_conserv==energy) CALL test_message(req_qu) 
@@ -233 +287 @@
     !--> flush wflux  
     !$OMP BARRIER
-
+    !$acc parallel loop collapse(3) async
     DO iq=1,nqdyn
        DO l=ll_begin,ll_endm1
@@ -242 +296 @@
           END DO
        END DO
+    END DO
+    
+    !$acc parallel loop collapse(3) async
+    DO iq=1,nqdyn
        DO l=ll_beginp1,ll_end
           !DIR$ SIMD
@@ -252 +310 @@
 
     ! Compute vertical transport
+    !$acc parallel loop collapse(2) async
     DO l=ll_beginp1,ll_end
        !DIR$ SIMD
@@ -265 +324 @@
 
     ! Add vertical transport to du
+    !$acc parallel loop collapse(2) async
     DO l=ll_begin,ll_end
        !DIR$ SIMD
@@ -283 +343 @@
     !   ENDDO
 
+    !$acc end data
     CALL trace_end("compute_caldyn_vert")
 

codes/icosagcm/trunk/src/dynamics/caldyn_kernels_hevi.F90

@@ -6 +6 @@
   USE transfert_mod
   USE caldyn_kernels_base_mod
+  USE abort_mod
   IMPLICIT NONE
   PRIVATE
@@ -20 +21 @@
 
   SUBROUTINE compute_theta(ps,theta_rhodz, rhodz,theta)
+    USE disvert_mod, ONLY : mass_dbk, mass_dak
     REAL(rstd),INTENT(IN)    :: ps(iim*jjm)
     REAL(rstd),INTENT(IN)    :: theta_rhodz(iim*jjm,llm,nqdyn)
@@ -26 +28 @@
     INTEGER :: ij,l,iq
     REAL(rstd) :: m
+
     CALL trace_start("compute_theta")
 
+    !$acc data async &
+    !$acc present(ps(:), theta_rhodz(:,:,:), rhodz(:,:), theta(:,:,:), mass_dak(:), mass_dbk(:))
+
     IF(caldyn_eta==eta_mass) THEN ! Compute mass
+       !$acc parallel loop collapse(2) async
        DO l = ll_begin,ll_end
           !DIR$ SIMD
@@ -38 +45 @@
     END IF
 
-    DO l = ll_begin,ll_end
-       DO iq=1,nqdyn
+    !$acc parallel loop collapse(3) async
+    DO iq=1,nqdyn
+       DO l = ll_begin,ll_end
           !DIR$ SIMD
           DO ij=ij_begin_ext,ij_end_ext
@@ -47 +55 @@
     END DO
 
+    !$acc end data
     CALL trace_end("compute_theta")
   END SUBROUTINE compute_theta
 
   SUBROUTINE compute_pvort_only(u,rhodz,qu,qv)
+    USE geometry, ONLY : Av, Riv2, fv
     REAL(rstd),INTENT(IN)  :: u(iim*3*jjm,llm)
     REAL(rstd),INTENT(INOUT) :: rhodz(iim*jjm,llm)
@@ -61 +71 @@
     CALL trace_start("compute_pvort_only")  
 !!! Compute shallow-water potential vorticity
+
     IF(dysl_pvort_only) THEN
+      CALL abort_acc("HEVI_scheme/dysl_pvort_only")
 #include "../kernels/pvort_only.k90"
     ELSE
 
+    !$acc data async &
+    !$acc present(rhodz(:,:), u(:,:), qu(:,:), qv(:,:), Av(:), Riv2(:,:), fv(:))
+    
     radius_m2=radius**(-2)
+    !$acc parallel loop collapse(2) async
     DO l = ll_begin,ll_end
        !DIR$ SIMD
@@ -85 +101 @@
           qv(ij+z_down,l) =( etav+fv(ij+z_down) )/hv
        ENDDO
-
+    END DO
+
+    !$acc parallel loop collapse(2) async
+    DO l = ll_begin,ll_end
        !DIR$ SIMD
        DO ij=ij_begin,ij_end
@@ -94 +113 @@
 
     ENDDO
+
+    !$acc end data
    
     END IF ! dysl
+
     CALL trace_end("compute_pvort_only")
 
@@ -409 +431 @@
 
     IF(dysl_caldyn_fast) THEN
+       CALL abort_acc("HEVI_scheme/dysl_caldyn_fast")
 #include "../kernels/caldyn_fast.k90"
     ELSE
 
-    ! Compute Bernoulli term
-    IF(boussinesq) THEN
-       DO l=ll_begin,ll_end
-          !DIR$ SIMD
-          DO ij=ij_begin,ij_end         
-             berni(ij,l) = pk(ij,l)
-             ! from now on pk contains the vertically-averaged geopotential
-             pk(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1))
+
+       ! Default case : ported to openacc
+       IF(.not. boussinesq .and. caldyn_thermo /= thermo_moist) THEN
+          !$acc data async &
+          !$acc present(rhodz(:,:), u(:,:),pk(:,:), geopot(:,:), theta(:,:,:), du(:,:))
+
+#define BERNI(ij) .5*(geopot(ij,l)+geopot(ij,l+1))
+          !$acc parallel loop collapse(2) async
+          DO l=ll_begin,ll_end
+             !DIR$ SIMD
+             DO ij=ij_begin,ij_end
+                due_right =  0.5*(theta(ij,l,1)+theta(ij+t_right,l,1))  &
+                     *(pk(ij+t_right,l)-pk(ij,l))        &
+                     +  BERNI(ij+t_right)-BERNI(ij)
+                due_lup = 0.5*(theta(ij,l,1)+theta(ij+t_lup,l,1))    &
+                     *(pk(ij+t_lup,l)-pk(ij,l))          &
+                     +  BERNI(ij+t_lup)-BERNI(ij)
+                due_ldown = 0.5*(theta(ij,l,1)+theta(ij+t_ldown,l,1)) &
+                     *(pk(ij+t_ldown,l)-pk(ij,l))      &
+                     +   BERNI(ij+t_ldown)-BERNI(ij)
+                du(ij+u_right,l) = du(ij+u_right,l) - ne_right*due_right
+                du(ij+u_lup,l)   = du(ij+u_lup,l)   - ne_lup*due_lup
+                du(ij+u_ldown,l) = du(ij+u_ldown,l) - ne_ldown*due_ldown
+                u(ij+u_right,l) = u(ij+u_right,l) + tau*du(ij+u_right,l)
+                u(ij+u_lup,l)   = u(ij+u_lup,l)   + tau*du(ij+u_lup,l)
+                u(ij+u_ldown,l) = u(ij+u_ldown,l) + tau*du(ij+u_ldown,l)
+             END DO
           END DO
-       END DO
-    ELSE ! compressible
-
-       DO l=ll_begin,ll_end
-          SELECT CASE(caldyn_thermo)
-          CASE(thermo_theta) ! vdp = theta.dpi => B = Phi
-             !DIR$ SIMD
-             DO ij=ij_begin,ij_end         
-                berni(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1))
-             END DO
-          CASE(thermo_entropy) ! vdp = dG + sdT => B = Phi + G, G=h-Ts=T*(cpp-s) 
-             !DIR$ SIMD
-             DO ij=ij_begin,ij_end         
-                berni(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1)) &
-                     + pk(ij,l)*(cpp-theta(ij,l,1)) ! pk=temperature, theta=entropy
-             END DO
-          CASE(thermo_moist) 
-             !DIR$ SIMD
-             DO ij=ij_begin,ij_end         
-                ! du/dt = grad(Bd)+rv.grad(Bv)+s.grad(T)
-                ! Bd = Phi + gibbs_d
-                ! Bv = Phi + gibbs_v
-                ! pk=temperature, theta=entropy
-                qv = theta(ij,l,2)
-                temp = pk(ij,l)
-                chi = log(temp/Treff)
-                nu = (chi*(cpp+qv*cppv)-theta(ij,l,1))/(Rd+qv*Rv) ! log(p/preff)
-                berni(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1)) &
-                     + temp*(cpp*(1.-chi)+Rd*nu)
-                berniv(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1)) &
-                     + temp*(cppv*(1.-chi)+Rv*nu)
-            END DO
-          END SELECT
-       END DO
-
-    END IF ! Boussinesq/compressible
+
+#undef BERNI
+
+          !$acc end data
+       ELSE ! Not default case : not ported to openacc
+          CALL abort_acc("HEVI_scheme/compute_caldyn_fast")
+          !!$acc parallel loop gang private(berni(:,:), berniv(:,:)) async
+          DO l=ll_begin,ll_end
+             ! Compute Bernoulliterm
+             IF(boussinesq) THEN
+                !!$acc loop vector
+                !DIR$ SIMD
+                DO ij=ij_begin,ij_end
+                   berni(ij,l) = pk(ij,l)
+                   ! from now on pk contains the vertically-averaged geopotential
+                   pk(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1))
+                END DO
+             ELSE ! compressible
+                SELECT CASE(caldyn_thermo)
+                CASE(thermo_theta) ! vdp = theta.dpi => B = Phi
+                   CALL dynamico_abort("Case already dealt with")
+                CASE(thermo_entropy) ! vdp = dG + sdT => B = Phi + G, G=h-Ts=T*(cpp-s)
+                   !$acc loop vector
+                   !DIR$ SIMD
+                   DO ij=ij_begin,ij_end
+                      berni(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1)) &
+                           + pk(ij,l)*(cpp-theta(ij,l,1)) ! pk=temperature, theta=entropy
+                   END DO
+                CASE(thermo_moist)
+!!$acc loop vector
+                   !DIR$ SIMD
+                   DO ij=ij_begin,ij_end
+                      ! du/dt = grad(Bd)+rv.grad(Bv)+s.grad(T)
+                      ! Bd = Phi + gibbs_d
+                      ! Bv = Phi + gibbs_v
+                      ! pk=temperature, theta=entropy
+                      qv = theta(ij,l,2)
+                      temp = pk(ij,l)
+                      chi = log(temp/Treff)
+                      nu = (chi*(cpp+qv*cppv)-theta(ij,l,1))/(Rd+qv*Rv) ! log(p/preff)
+                      berni(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1)) &
+                           + temp*(cpp*(1.-chi)+Rd*nu)
+                      berniv(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1)) &
+                           + temp*(cppv*(1.-chi)+Rv*nu)
+                   END DO
+                END SELECT
+             END IF ! Boussinesq/compressible
 
 !!! u:=u+tau*du, du = -grad(B)-theta.grad(pi)
-    DO l=ll_begin,ll_end
-       IF(caldyn_thermo == thermo_moist) THEN
-          !DIR$ SIMD
-          DO ij=ij_begin,ij_end         
-             due_right =  berni(ij+t_right,l)-berni(ij,l)       &
-                  + 0.5*(theta(ij,l,1)+theta(ij+t_right,l,1))   &
-                       *(pk(ij+t_right,l)-pk(ij,l))             &
-                  + 0.5*(theta(ij,l,2)+theta(ij+t_right,l,2))   &
-                       *(berniv(ij+t_right,l)-berniv(ij,l))
-             
-             due_lup = berni(ij+t_lup,l)-berni(ij,l)            &
-                  + 0.5*(theta(ij,l,1)+theta(ij+t_lup,l,1))     &
-                       *(pk(ij+t_lup,l)-pk(ij,l))               &
-                  + 0.5*(theta(ij,l,2)+theta(ij+t_lup,l,2))     &
-                       *(berniv(ij+t_lup,l)-berniv(ij,l))
-             
-             due_ldown = berni(ij+t_ldown,l)-berni(ij,l)        &
-                  + 0.5*(theta(ij,l,1)+theta(ij+t_ldown,l,1))   &
-                       *(pk(ij+t_ldown,l)-pk(ij,l))             &
-                  + 0.5*(theta(ij,l,2)+theta(ij+t_ldown,l,2))   &
-                       *(berniv(ij+t_ldown,l)-berniv(ij,l))
-             
-             du(ij+u_right,l) = du(ij+u_right,l) - ne_right*due_right
-             du(ij+u_lup,l)   = du(ij+u_lup,l)   - ne_lup*due_lup
-             du(ij+u_ldown,l) = du(ij+u_ldown,l) - ne_ldown*due_ldown
-             u(ij+u_right,l) = u(ij+u_right,l) + tau*du(ij+u_right,l)
-             u(ij+u_lup,l)   = u(ij+u_lup,l)   + tau*du(ij+u_lup,l)
-             u(ij+u_ldown,l) = u(ij+u_ldown,l) + tau*du(ij+u_ldown,l)
-          END DO
-       ELSE
-          !DIR$ SIMD
-          DO ij=ij_begin,ij_end         
-             due_right =  0.5*(theta(ij,l,1)+theta(ij+t_right,l,1))  &
-                  *(pk(ij+t_right,l)-pk(ij,l))        &
-                  +  berni(ij+t_right,l)-berni(ij,l)
-             due_lup = 0.5*(theta(ij,l,1)+theta(ij+t_lup,l,1))    &
-                  *(pk(ij+t_lup,l)-pk(ij,l))          &
-                  +  berni(ij+t_lup,l)-berni(ij,l)
-             due_ldown = 0.5*(theta(ij,l,1)+theta(ij+t_ldown,l,1)) &
-                  *(pk(ij+t_ldown,l)-pk(ij,l))      &
-                  +   berni(ij+t_ldown,l)-berni(ij,l)
-             du(ij+u_right,l) = du(ij+u_right,l) - ne_right*due_right
-             du(ij+u_lup,l)   = du(ij+u_lup,l)   - ne_lup*due_lup
-             du(ij+u_ldown,l) = du(ij+u_ldown,l) - ne_ldown*due_ldown
-             u(ij+u_right,l) = u(ij+u_right,l) + tau*du(ij+u_right,l)
-             u(ij+u_lup,l)   = u(ij+u_lup,l)   + tau*du(ij+u_lup,l)
-             u(ij+u_ldown,l) = u(ij+u_ldown,l) + tau*du(ij+u_ldown,l)
+             IF(caldyn_thermo == thermo_moist) THEN
+!!$acc loop vector
+                !DIR$ SIMD
+                DO ij=ij_begin,ij_end
+                   due_right =  berni(ij+t_right,l)-berni(ij,l)       &
+                        + 0.5*(theta(ij,l,1)+theta(ij+t_right,l,1))   &
+                        *(pk(ij+t_right,l)-pk(ij,l))             &
+                        + 0.5*(theta(ij,l,2)+theta(ij+t_right,l,2))   &
+                        *(berniv(ij+t_right,l)-berniv(ij,l))
+
+                   due_lup = berni(ij+t_lup,l)-berni(ij,l)            &
+                        + 0.5*(theta(ij,l,1)+theta(ij+t_lup,l,1))     &
+                        *(pk(ij+t_lup,l)-pk(ij,l))               &
+                        + 0.5*(theta(ij,l,2)+theta(ij+t_lup,l,2))     &
+                        *(berniv(ij+t_lup,l)-berniv(ij,l))
+
+                   due_ldown = berni(ij+t_ldown,l)-berni(ij,l)        &
+                        + 0.5*(theta(ij,l,1)+theta(ij+t_ldown,l,1))   &
+                        *(pk(ij+t_ldown,l)-pk(ij,l))             &
+                        + 0.5*(theta(ij,l,2)+theta(ij+t_ldown,l,2))   &
+                        *(berniv(ij+t_ldown,l)-berniv(ij,l))
+
+                   du(ij+u_right,l) = du(ij+u_right,l) - ne_right*due_right
+                   du(ij+u_lup,l)   = du(ij+u_lup,l)   - ne_lup*due_lup
+                   du(ij+u_ldown,l) = du(ij+u_ldown,l) - ne_ldown*due_ldown
+                   u(ij+u_right,l) = u(ij+u_right,l) + tau*du(ij+u_right,l)
+                   u(ij+u_lup,l)   = u(ij+u_lup,l)   + tau*du(ij+u_lup,l)
+                   u(ij+u_ldown,l) = u(ij+u_ldown,l) + tau*du(ij+u_ldown,l)
+                END DO
+             ELSE
+                CALL dynamico_abort("Case already dealt with")
+             END IF
           END DO
        END IF
-    END DO
-
-    END IF ! dysl
+    END IF
     CALL trace_end("compute_caldyn_fast")
 
   END SUBROUTINE compute_caldyn_fast
-
+    
   SUBROUTINE compute_caldyn_Coriolis(hflux,theta,qu, convm,dtheta_rhodz,du)
+    USE geometry, ONLY : Ai, wee
     REAL(rstd),INTENT(IN)  :: hflux(3*iim*jjm,llm)  ! hflux in kg/s
     REAL(rstd),INTENT(IN)  :: theta(iim*jjm,llm,nqdyn) ! active scalars
@@ -534 +571 @@
 
     ELSE
-#define FTHETA(ij) Ftheta(ij,1)
-
-    DO l=ll_begin, ll_end
-       ! compute theta flux
-       DO iq=1,nqdyn
-       !DIR$ SIMD
+#define FTHETA(ij) Ftheta(ij,l)
+
+    !$acc data async &
+    !$acc create(Ftheta(3*iim*jjm,llm)) &
+    !$acc present(qu(:,:), hflux(:,:), convm(:,:), du(:,:), dtheta_rhodz(:,:,:), theta(:,:,:), Ai(:), wee(:,:,:))
+
+    ! compute theta flux
+    DO iq=1,nqdyn
+       !$acc parallel loop collapse(2) present(Ftheta, theta) async
+       DO l=ll_begin, ll_end
+          !DIR$ SIMD
           DO ij=ij_begin_ext,ij_end_ext      
              FTHETA(ij+u_right) = 0.5*(theta(ij,l,iq)+theta(ij+t_right,l,iq)) &
@@ -548 +590 @@
                   * hflux(ij+u_ldown,l)
           END DO
+       END DO
+
+       !$acc parallel loop collapse(2) present(Ftheta) async
+       DO l=ll_begin, ll_end
           ! horizontal divergence of fluxes
-       !DIR$ SIMD
+          !DIR$ SIMD
           DO ij=ij_begin,ij_end         
              ! dtheta_rhodz =  -div(flux.theta)
@@ -561 +607 @@
           END DO
        END DO
-
+    END DO
+
+    !$acc parallel loop collapse(2) present(Ai) async
+    DO l=ll_begin, ll_end
        !DIR$ SIMD
        DO ij=ij_begin,ij_end         
@@ -578 +627 @@
 
     CASE(energy) ! energy-conserving TRiSK
-
+       !$acc parallel loop collapse(2) present(qu, wee, du) async
        DO l=ll_begin,ll_end
           !DIR$ SIMD
@@ -622 +671 @@
 
     CASE(enstrophy) ! enstrophy-conserving TRiSK
-
+       
+       !$acc parallel loop collapse(2) present(wee, du) async
        DO l=ll_begin,ll_end
           !DIR$ SIMD
@@ -665 +715 @@
           END DO
        END DO
+       !$acc end parallel loop
     CASE DEFAULT
        STOP
     END SELECT
+
+    !$acc end data
+
 #undef FTHETA
 
@@ -676 +730 @@
   END SUBROUTINE compute_caldyn_Coriolis
 
-  SUBROUTINE compute_caldyn_slow_hydro(u,rhodz,hflux,du, zero)
+  SUBROUTINE compute_caldyn_slow_hydro(u,rhodz,hflux,du,zero)
+    USE geometry, ONLY : Ai, le_de
     LOGICAL, INTENT(IN) :: zero
     REAL(rstd),INTENT(IN)  :: u(3*iim*jjm,llm)    ! prognostic "velocity"
@@ -684 +739 @@
     
     REAL(rstd) :: berni(iim*jjm,llm)  ! Bernoulli function
-    REAL(rstd) :: berni1(iim*jjm)  ! Bernoulli function
     REAL(rstd) :: uu_right, uu_lup, uu_ldown, ke, uu
     INTEGER :: ij,l
@@ -698 +752 @@
      ELSE
 
-#define BERNI(ij) berni1(ij)
+#define BERNI(ij) berni(ij,l)
 
     DO l = ll_begin, ll_end
-       !  Compute mass fluxes
        IF (caldyn_conserv==energy) CALL test_message(req_qu) 
+    END DO
+    
+    !$acc data async &
+    !$acc create(berni(:,ll_begin:ll_end)) &
+    !$acc present(rhodz(:,ll_begin:ll_end), u(:,ll_begin:ll_end), hflux(:,:), du(:,ll_begin:ll_end), Ai(:), le_de(:))
+
+    !  Compute mass fluxes
+    !$acc parallel loop collapse(2) async
+    DO l = ll_begin, ll_end
        !DIR$ SIMD
        DO ij=ij_begin_ext,ij_end_ext
@@ -715 +777 @@
           hflux(ij+u_ldown,l)=uu_ldown
        ENDDO
-       ! Compute Bernoulli=kinetic energy
+    END DO
+
+    ! Compute Bernoulli=kinetic energy
+    !$acc parallel loop collapse(2) async
+    DO l = ll_begin, ll_end
        !DIR$ SIMD
        DO ij=ij_begin,ij_end         
@@ -726 +792 @@
                le_de(ij+u_rdown)*u(ij+u_rdown,l)**2 )  
        ENDDO
-       ! Compute du=-grad(Bernoulli)
-       IF(zero) THEN
+    END DO
+
+    ! Compute du=-grad(Bernoulli)
+    IF(zero) THEN
+       !$acc parallel loop collapse(2) async
+       DO l = ll_begin, ll_end
           !DIR$ SIMD
           DO ij=ij_begin,ij_end
@@ -734 +804 @@
              du(ij+u_ldown,l) = ne_ldown*(BERNI(ij)-BERNI(ij+t_ldown))
           END DO
-       ELSE
+       END DO
+    ELSE
+       !$acc parallel loop collapse(2) async
+       DO l = ll_begin, ll_end
           !DIR$ SIMD
           DO ij=ij_begin,ij_end
@@ -744 +817 @@
                   ne_ldown*(BERNI(ij)-BERNI(ij+t_ldown))
           END DO
-       END IF
-    END DO
+       END DO
+    END IF
+
+    !$acc end data
 
 #undef BERNI
 
     END IF ! dysl
+
     CALL trace_end("compute_caldyn_slow_hydro")    
   END SUBROUTINE compute_caldyn_slow_hydro

codes/icosagcm/trunk/src/kernels/advect_horiz.k90

@@ -1 +1 @@
    !--------------------------------------------------------------------------
    !---------------------------- advect_horiz ----------------------------------
+
+   !$acc data create(qflux(:,:)) present(qi(:,:), cc(:,:,:), gradq3d(:,:,:), mass(:,:), hfluxt(:,:), Ai(:), xyz_i(:,:)) async
+
    ! evaluate tracer field at point cc using piecewise linear reconstruction
    ! q(cc)= q0 + gradq.(cc-xyz_i) with xi centroid of hexagon
    ! sign*hfluxt>0 iff outgoing
+   !$acc parallel loop collapse(2) async
    DO l = ll_begin, ll_end
       !DIR$ SIMD
       DO ij=ij_begin_ext, ij_end_ext
          IF(ne_right*hfluxt(ij+u_right,l)>0.) THEN
-            qe = qi(ij,l)
-            qe = qe + (cc(ij+u_right,l,1)-xyz_i(ij,1))*gradq3d(ij,l,1)
-            qe = qe + (cc(ij+u_right,l,2)-xyz_i(ij,2))*gradq3d(ij,l,2)
-            qe = qe + (cc(ij+u_right,l,3)-xyz_i(ij,3))*gradq3d(ij,l,3)
+            ij_tmp = ij
          ELSE
-            qe = qi(ij+t_right,l)
-            qe = qe + (cc(ij+u_right,l,1)-xyz_i(ij+t_right,1))*gradq3d(ij+t_right,l,1)
-            qe = qe + (cc(ij+u_right,l,2)-xyz_i(ij+t_right,2))*gradq3d(ij+t_right,l,2)
-            qe = qe + (cc(ij+u_right,l,3)-xyz_i(ij+t_right,3))*gradq3d(ij+t_right,l,3)
+            ij_tmp = ij+t_right
          END IF
+
+         qe = qi(ij_tmp,l)
+         qe = qe + (cc(ij+u_right,l,1)-xyz_i(ij_tmp,1))*gradq3d(ij_tmp,l,1)
+         qe = qe + (cc(ij+u_right,l,2)-xyz_i(ij_tmp,2))*gradq3d(ij_tmp,l,2)
+         qe = qe + (cc(ij+u_right,l,3)-xyz_i(ij_tmp,3))*gradq3d(ij_tmp,l,3)
+
          qflux(ij+u_right,l) = hfluxt(ij+u_right,l)*qe
-         IF(diagflux_on) qfluxt(ij+u_right,l) = qfluxt(ij+u_right,l)+qflux(ij+u_right,l)
+
          IF(ne_lup*hfluxt(ij+u_lup,l)>0.) THEN
-            qe = qi(ij,l)
-            qe = qe + (cc(ij+u_lup,l,1)-xyz_i(ij,1))*gradq3d(ij,l,1)
-            qe = qe + (cc(ij+u_lup,l,2)-xyz_i(ij,2))*gradq3d(ij,l,2)
-            qe = qe + (cc(ij+u_lup,l,3)-xyz_i(ij,3))*gradq3d(ij,l,3)
+            ij_tmp = ij
          ELSE
-            qe = qi(ij+t_lup,l)
-            qe = qe + (cc(ij+u_lup,l,1)-xyz_i(ij+t_lup,1))*gradq3d(ij+t_lup,l,1)
-            qe = qe + (cc(ij+u_lup,l,2)-xyz_i(ij+t_lup,2))*gradq3d(ij+t_lup,l,2)
-            qe = qe + (cc(ij+u_lup,l,3)-xyz_i(ij+t_lup,3))*gradq3d(ij+t_lup,l,3)
+            ij_tmp = ij+t_lup
          END IF
+
+         qe = qi(ij_tmp,l)
+         qe = qe + (cc(ij+u_lup,l,1)-xyz_i(ij_tmp,1))*gradq3d(ij_tmp,l,1)
+         qe = qe + (cc(ij+u_lup,l,2)-xyz_i(ij_tmp,2))*gradq3d(ij_tmp,l,2)
+         qe = qe + (cc(ij+u_lup,l,3)-xyz_i(ij_tmp,3))*gradq3d(ij_tmp,l,3)
+
          qflux(ij+u_lup,l) = hfluxt(ij+u_lup,l)*qe
-         IF(diagflux_on) qfluxt(ij+u_lup,l) = qfluxt(ij+u_lup,l)+qflux(ij+u_lup,l)
+
          IF(ne_ldown*hfluxt(ij+u_ldown,l)>0.) THEN
-            qe = qi(ij,l)
-            qe = qe + (cc(ij+u_ldown,l,1)-xyz_i(ij,1))*gradq3d(ij,l,1)
-            qe = qe + (cc(ij+u_ldown,l,2)-xyz_i(ij,2))*gradq3d(ij,l,2)
-            qe = qe + (cc(ij+u_ldown,l,3)-xyz_i(ij,3))*gradq3d(ij,l,3)
+            ij_tmp = ij
          ELSE
-            qe = qi(ij+t_ldown,l)
-            qe = qe + (cc(ij+u_ldown,l,1)-xyz_i(ij+t_ldown,1))*gradq3d(ij+t_ldown,l,1)
-            qe = qe + (cc(ij+u_ldown,l,2)-xyz_i(ij+t_ldown,2))*gradq3d(ij+t_ldown,l,2)
-            qe = qe + (cc(ij+u_ldown,l,3)-xyz_i(ij+t_ldown,3))*gradq3d(ij+t_ldown,l,3)
+            ij_tmp = ij+t_ldown
          END IF
+
+         qe = qi(ij_tmp,l)
+         qe = qe + (cc(ij+u_ldown,l,1)-xyz_i(ij_tmp,1))*gradq3d(ij_tmp,l,1)
+         qe = qe + (cc(ij+u_ldown,l,2)-xyz_i(ij_tmp,2))*gradq3d(ij_tmp,l,2)
+         qe = qe + (cc(ij+u_ldown,l,3)-xyz_i(ij_tmp,3))*gradq3d(ij_tmp,l,3)
+
          qflux(ij+u_ldown,l) = hfluxt(ij+u_ldown,l)*qe
-         IF(diagflux_on) qfluxt(ij+u_ldown,l) = qfluxt(ij+u_ldown,l)+qflux(ij+u_ldown,l)
       END DO
    END DO
+
+   IF(diagflux_on) THEN
+     !$acc parallel loop collapse(2) copy(qfluxt(:,:)) async
+     DO l = ll_begin, ll_end
+        !DIR$ SIMD
+        DO ij=ij_begin_ext, ij_end_ext
+           qfluxt(ij+u_right,l) = qfluxt(ij+u_right,l)+qflux(ij+u_right,l)
+           qfluxt(ij+u_lup,l) = qfluxt(ij+u_lup,l)+qflux(ij+u_lup,l)
+           qfluxt(ij+u_ldown,l) = qfluxt(ij+u_ldown,l)+qflux(ij+u_ldown,l)
+        END DO
+     END DO
+   END IF
+
    ! update q and, if update_mass, update mass
+   !$acc parallel loop collapse(2) async
    DO l = ll_begin, ll_end
       !DIR$ SIMD
@@ -71 +88 @@
       END DO
    END DO
+   
+   !$acc end data
+
    !---------------------------- advect_horiz ----------------------------------
    !--------------------------------------------------------------------------

codes/icosagcm/trunk/src/output/xios_mod.F90

@@ -1120 +1120 @@
    CHARACTER(LEN=*),INTENT(IN) :: name
    REAL(rstd), INTENT(OUT) :: field
+   field=0
  END SUBROUTINE
  

codes/icosagcm/trunk/src/parallel/mpipara.F90

@@ -14 +14 @@
   INTEGER,SAVE :: id_mpi ! id for profiling
 
+  INTEGER, SAVE :: device_id
+
   INTERFACE allocate_mpi_buffer
     MODULE PROCEDURE allocate_mpi_buffer_r2, allocate_mpi_buffer_r3,allocate_mpi_buffer_r4
@@ -46 +48 @@
   USE xios
 #endif
+  USE abort_mod
   IMPLICIT NONE
     CHARACTER(LEN=256) :: required_mode_str
@@ -75 +78 @@
       END SELECT
       
+      IF (required_mode==MPI_THREAD_SERIALIZED .OR. required_mode==MPI_THREAD_MULTIPLE) THEN
+        CALL abort_acc("mpi_threading_mode /= 'single' .AND. mpi_threading_mode /= 'funneled'")
+      ENDIF
 
       IF (required_mode==MPI_THREAD_SINGLE)     PRINT*,'MPI_INIT_THREAD : MPI_SINGLE_THREAD required'
@@ -120 +126 @@
     ENDIF
     
+
+#ifdef _OPENACC
+    device_id = setDevice(mpi_size, mpi_rank)
+    PRINT *, 'GPU device ', device_id
+#else
+     device_id = -1
+#endif
+
   END SUBROUTINE  init_mpipara
 
@@ -233 +247 @@
   END SUBROUTINE free_mpi_buffer_r4
    
+#ifdef _OPENACC
+  FUNCTION setDevice(nprocs, myrank)
+    use iso_c_binding
+    use openacc
+    USE mpi_mod
+    implicit none
+    
+    interface 
+       function gethostid() bind(C)
+       use iso_c_binding
+       integer(C_INT) :: gethostid
+     end function gethostid
+  end interface
+  
+  integer, intent(in) :: nprocs, myrank
+  integer :: hostids(nprocs), localprocs(nprocs)
+  integer :: hostid, ierr, numdev, mydev, i, numlocal
+  integer :: setDevice
+
+  ! get the hostids so we can determine what other processes are on this node
+  hostid = gethostid()
+  call mpi_allgather(hostid,1,MPI_INTEGER,hostids,1,MPI_INTEGER, MPI_COMM_WORLD, ierr)
+
+  ! determine which processors are on this node
+  numlocal = 0
+  localprocs(:) = 0
+  do i = 1, nprocs
+     if (hostid == hostids(i)) then
+        localprocs(i) = numlocal
+        numlocal = numlocal + 1
+     end if
+  end do
+
+  ! get the number of device on this node
+  numdev = acc_get_num_devices(ACC_DEVICE_NVIDIA)
+
+  if (numdev < 1) then
+     print *, "Error: there are no devices available on this host. ABORTING", myrank
+     stop
+  end if
+
+  ! print a warning if the number of devices is less than the number of processes on this node. Having multiple processes share a devices is not recommended
+  if (numdev < numlocal) then
+     if (localprocs(myrank+1) == 1) then
+        ! print warning message only once per node
+        print *, "WARNING: the number of process is greater than the number of GPUs.", myrank
+     end if
+     mydev = mod(localprocs(myrank+1), numdev)
+  else
+     mydev = localprocs(myrank+1)
+  end if
+
+  call acc_set_device_num(mydev,ACC_DEVICE_NVIDIA)
+  call acc_init(ACC_DEVICE_NVIDIA)
+  setDevice = acc_get_device_num(ACC_DEVICE_NVIDIA)
+END FUNCTION setDevice
+
+#endif
+
+     
+     
 END MODULE mpipara

codes/icosagcm/trunk/src/parallel/transfert_mpi.f90

@@ -5 +5 @@
   
   TYPE array
-    INTEGER,POINTER :: value(:)
-    INTEGER,POINTER :: sign(:)
+    INTEGER,POINTER :: value(:)=>null()
+    INTEGER,POINTER :: sign(:)=>null()
     INTEGER         :: domain
     INTEGER         :: rank
@@ -13 +13 @@
     INTEGER         :: offset
     INTEGER         :: ireq
-    INTEGER,POINTER :: buffer(:)
-    REAL,POINTER    :: buffer_r(:)
-    INTEGER,POINTER :: src_value(:)
+    INTEGER,POINTER :: buffer(:)=>null()
+    REAL,POINTER    :: buffer_r(:)=>null()
+    INTEGER,POINTER :: src_value(:)=>null()
   END TYPE array
   
@@ -77 +77 @@
     LOGICAL :: open      ! for debug
     INTEGER :: number
+    LOGICAL :: ondevice=.false. !<Ready to transfer ondevice field
   END TYPE t_message
 
@@ -942 +943 @@
 !$OMP BARRIER
 !$OMP MASTER
+    !init off-device
+    message%ondevice=.false.
+
     IF(PRESENT(name)) THEN
        message%name = TRIM(name)
@@ -1058 +1062 @@
     TYPE(t_message) :: message
 
-    INTEGER :: ireq
+    INTEGER :: ireq, ibuff
 
 !$OMP BARRIER
@@ -1084 +1088 @@
         ENDDO
 
-      ENDIF      
+      ENDIF
     ENDIF
-    
+
+    !deallocate device data if ondevice
+    if(message%ondevice) then
+      do ireq=1, ndomain
+        do ibuff=1,message%request(ireq)%nsend
+          !$acc exit data delete(message%request(ireq)%send(ibuff)%buffer(:))
+          !$acc exit data delete(message%request(ireq)%send(ibuff)%buffer_r(:))
+          !$acc exit data delete(message%request(ireq)%send(ibuff)%sign(:))
+          !$acc exit data delete(message%request(ireq)%send(ibuff)%src_value(:))
+          !$acc exit data delete(message%request(ireq)%send(ibuff)%value(:))
+        end do
+        do ibuff=1,message%request(ireq)%nrecv
+          !$acc exit data delete(message%request(ireq)%recv(ibuff)%buffer(:))
+          !$acc exit data delete(message%request(ireq)%recv(ibuff)%buffer_r(:))
+          !$acc exit data delete(message%request(ireq)%recv(ibuff)%sign(:))
+          !$acc exit data delete(message%request(ireq)%recv(ibuff)%src_value(:))
+          !$acc exit data delete(message%request(ireq)%recv(ibuff)%value(:))
+        end do
+      end do
+      DO ireq=1,message%nreq
+        !$acc exit data delete(message%buffers(ireq)%r)
+      ENDDO
+      message%ondevice=.false.
+    end if
+
     DEALLOCATE(message%mpi_req)
     DEALLOCATE(message%buffers)
@@ -1120 +1148 @@
 
 
+  !!!Update buffers on device for 'message'
+  !!! does create_device_message when not already ondevice
+  SUBROUTINE update_device_message_mpi(message)
+    USE domain_mod
+    IMPLICIT NONE
+    TYPE(t_message), intent(inout) :: message
+    INTEGER :: ireq, ibuff
+
+    !if(.not. message%ondevice) call create_device_message_mpi(message)
+
+    do ireq=1, ndomain
+      do ibuff=1,message%request(ireq)%nsend
+        !device buffers updated even if pointers not attached :
+        !non allocated buffers in 'message' must be set to NULL()
+        !$acc enter data copyin(message%request(ireq)%send(ibuff)%buffer(:)) async
+        !$acc enter data copyin(message%request(ireq)%send(ibuff)%buffer_r(:)) async
+        !$acc enter data copyin(message%request(ireq)%send(ibuff)%sign(:)) async
+        !$acc enter data copyin(message%request(ireq)%send(ibuff)%src_value(:)) async
+        !$acc enter data copyin(message%request(ireq)%send(ibuff)%value(:)) async
+      end do
+      do ibuff=1,message%request(ireq)%nrecv
+        !$acc enter data copyin(message%request(ireq)%recv(ibuff)%buffer(:)) async
+        !$acc enter data copyin(message%request(ireq)%recv(ibuff)%buffer_r(:)) async
+        !$acc enter data copyin(message%request(ireq)%recv(ibuff)%sign(:)) async
+        !$acc enter data copyin(message%request(ireq)%recv(ibuff)%src_value(:)) async
+        !$acc enter data copyin(message%request(ireq)%recv(ibuff)%value(:)) async
+      end do
+    end do
+    DO ireq=1,message%nreq
+      !$acc enter data copyin(message%buffers(ireq)%r) async
+    ENDDO
+    message%ondevice=.true.
+  END SUBROUTINE
+
+  !TODO : add openacc with multiple process
   SUBROUTINE send_message_mpi(field,message)
   USE abort_mod
@@ -1129 +1192 @@
   USE omp_para
   USE trace
+  USE abort_mod
   IMPLICIT NONE
     TYPE(t_field),POINTER :: field(:)
@@ -1145 +1209 @@
     INTEGER :: dim3,dim4,d3,d4
     INTEGER,POINTER :: src_value(:)
-    INTEGER :: offset,msize,rank
+    INTEGER :: offset,msize,moffset,rank
     INTEGER :: lbegin, lend
     INTEGER :: max_req
@@ -1152 +1216 @@
 
     CALL enter_profile(id_mpi)
+
+    !Prepare 'message' for on-device copies if field is on device
+    if( field(1)%ondevice .AND. .NOT. message%ondevice ) call update_device_message_mpi(message)
 
 !$OMP BARRIER
@@ -1195 +1262 @@
               buffer_r=>message%buffers(ireq)%r
               offset=send%offset
-                                            
-              DO n=1,send%size
+              msize=send%size
+              !$acc parallel loop default(present) async if (field(ind)%ondevice)
+              DO n=1,msize
                 buffer_r(offset+n)=rval2d(value(n))
               ENDDO
 
               IF (mpi_threading_mode==MPI_THREAD_SERIALIZED) THEN
+                CALL abort_acc("mpi_threading_mode==MPI_THREAD_SERIALIZED")
                 !$OMP CRITICAL            
                 CALL MPI_ISEND(buffer_r,send%size,MPI_REAL8,send%rank,               &
@@ -1206 +1275 @@
                 !$OMP END CRITICAL
               ELSE IF (mpi_threading_mode==MPI_THREAD_MULTIPLE) THEN
+                CALL abort_acc("mpi_threading_mode==MPI_THREAD_MULTIPLE")
                 CALL MPI_ISEND(buffer_r,send%size,MPI_REAL8,send%rank,               &
                   send%tag+1000*message%number,comm_icosa, message%mpi_req(ireq),ierr)
@@ -1228 +1298 @@
               src_rval2d=>field(recv%domain)%rval2d
               sgn=>recv%sign
-
-              DO n=1,recv%size
+              msize=recv%size
+              !$acc parallel loop default(present) async if (field(ind)%ondevice)
+              DO n=1,msize
                 rval2d(value(n))=src_rval2d(src_value(n))*sgn(n)
               ENDDO
@@ -1239 +1310 @@
               buffer_r=>message%buffers(ireq)%r
               IF (mpi_threading_mode==MPI_THREAD_SERIALIZED) THEN
+                CALL abort_acc("mpi_threading_mode==MPI_THREAD_SERIALIZED")
                !$OMP CRITICAL            
                 CALL MPI_IRECV(buffer_r,recv%size,MPI_REAL8,recv%rank,               &
@@ -1244 +1316 @@
                !$OMP END CRITICAL
               ELSE IF (mpi_threading_mode==MPI_THREAD_MULTIPLE) THEN
+                 CALL abort_acc("mpi_threading_mode==MPI_THREAD_MULTIPLE")
                  CALL MPI_IRECV(buffer_r,recv%size,MPI_REAL8,recv%rank,              &
                    recv%tag+1000*message%number,comm_icosa, message%mpi_req(ireq),ierr)
@@ -1278 +1351 @@
               buffer_r=>message%buffers(ireq)%r
 
-              offset=send%offset*dim3 + (lbegin-1)*send%size
-          
+              msize=send%size
+              moffset=send%offset
               CALL trace_start("copy_to_buffer")
 
+              !$acc parallel loop default(present) async if (field(ind)%ondevice)
               DO d3=lbegin,lend
-                DO n=1,send%size
+                offset=moffset*dim3 + (d3-1)*msize
+                !$acc loop
+                DO n=1,msize
                   buffer_r(n+offset)=rval3d(value(n),d3)
                 ENDDO
-                offset=offset+send%size
               ENDDO
               CALL trace_end("copy_to_buffer")
@@ -1296 +1371 @@
               IF (is_omp_level_master) THEN
                 IF (mpi_threading_mode==MPI_THREAD_SERIALIZED) THEN
+                  CALL abort_acc("mpi_threading_mode==MPI_THREAD_SERIALIZED")
                   !$OMP CRITICAL   
                   CALL MPI_ISEND(buffer_r,send%size*dim3,MPI_REAL8,send%rank,        &
@@ -1301 +1377 @@
                   !$OMP END CRITICAL
                 ELSE IF (mpi_threading_mode==MPI_THREAD_MULTIPLE) THEN
+                  CALL abort_acc("mpi_threading_mode==MPI_THREAD_MULTIPLE")
                   CALL MPI_ISEND(buffer_r,send%size*dim3,MPI_REAL8,send%rank,        &
                     send%tag+1000*message%number,comm_icosa, message%mpi_req(ireq),ierr)
@@ -1335 +1412 @@
               src_rval3d=>field(recv%domain)%rval3d
               sgn=>recv%sign
+              msize=recv%size
 
               CALL trace_start("copy_data")
-
+              !$acc parallel loop collapse(2) default(present) async if (field(ind)%ondevice)
               DO d3=lbegin,lend
-                DO n=1,recv%size
+                DO n=1,msize
                   rval3d(value(n),d3)=src_rval3d(src_value(n),d3)*sgn(n)
                 ENDDO
@@ -1351 +1429 @@
               IF (is_omp_level_master) THEN
                 IF (mpi_threading_mode==MPI_THREAD_SERIALIZED) THEN
+                  CALL abort_acc("mpi_threading_mode==MPI_THREAD_SERIALIZED")
                   !$OMP CRITICAL
                   CALL MPI_IRECV(buffer_r,recv%size*dim3,MPI_REAL8,recv%rank,        &
@@ -1356 +1435 @@
                   !$OMP END CRITICAL
                 ELSE IF (mpi_threading_mode==MPI_THREAD_MULTIPLE) THEN
+                  CALL abort_acc("mpi_threading_mode==MPI_THREAD_MULTIPLE")
                   CALL MPI_IRECV(buffer_r,recv%size*dim3,MPI_REAL8,recv%rank,        &
                     recv%tag+1000*message%number,comm_icosa, message%mpi_req(ireq),ierr)
@@ -1390 +1470 @@
               ireq=send%ireq
               buffer_r=>message%buffers(ireq)%r
+              msize=send%size
+              moffset=send%offset
 
               CALL trace_start("copy_to_buffer")
+              !$acc parallel loop default(present) collapse(2) async if (field(ind)%ondevice)
               DO d4=1,dim4
-                offset=send%offset*dim3*dim4 + dim3*send%size*(d4-1) +               &
-                  (lbegin-1)*send%size
-
                 DO d3=lbegin,lend
-                  DO n=1,send%size
+                  offset=moffset*dim3*dim4 + dim3*msize*(d4-1) +               &
+                    (d3-1)*msize
+                  !$acc loop
+                  DO n=1,msize
                     buffer_r(n+offset)=rval4d(value(n),d3,d4)
                   ENDDO
-                  offset=offset+send%size
                 ENDDO
               ENDDO
@@ -1411 +1493 @@
               IF (is_omp_level_master) THEN
                 IF (mpi_threading_mode==MPI_THREAD_SERIALIZED) THEN
+                  CALL abort_acc("mpi_threading_mode==MPI_THREAD_SERIALIZED")
                   !$OMP CRITICAL
                   CALL MPI_ISEND(buffer_r,send%size*dim3*dim4,MPI_REAL8,send%rank,   &
@@ -1416 +1499 @@
                   !$OMP END CRITICAL
                 ELSE IF (mpi_threading_mode==MPI_THREAD_MULTIPLE) THEN
+                  CALL abort_acc("mpi_threading_mode==MPI_THREAD_MULTIPLE")
                   CALL MPI_ISEND(buffer_r,send%size*dim3*dim4,MPI_REAL8,send%rank,   &
                     send%tag+1000*message%number,comm_icosa, message%mpi_req(ireq),ierr)
@@ -1451 +1535 @@
               src_rval4d=>field(recv%domain)%rval4d
               sgn=>recv%sign
-
+              msize=recv%size
               CALL trace_start("copy_data")
+
+              !$acc parallel loop collapse(3) default(present) async if (field(ind)%ondevice)
               DO d4=1,dim4
                 DO d3=lbegin,lend
-                  DO n=1,recv%size
+                  DO n=1,msize
                     rval4d(value(n),d3,d4)=src_rval4d(src_value(n),d3,d4)*sgn(n)
                   ENDDO
@@ -1469 +1555 @@
               IF (is_omp_level_master) THEN
                 IF (mpi_threading_mode==MPI_THREAD_SERIALIZED) THEN
+                  CALL abort_acc("mpi_threading_mode==MPI_THREAD_SERIALIZED")
                  !$OMP CRITICAL           
                   CALL MPI_IRECV(buffer_r,recv%size*dim3*dim4,MPI_REAL8,recv%rank,   &
@@ -1474 +1561 @@
                   !$OMP END CRITICAL
                 ELSE IF (mpi_threading_mode==MPI_THREAD_MULTIPLE) THEN
+                  CALL abort_acc("mpi_threading_mode==MPI_THREAD_MULTIPLE")
                   CALL MPI_IRECV(buffer_r,recv%size*dim3*dim4,MPI_REAL8,recv%rank,   &
                     recv%tag+1000*message%number,comm_icosa, message%mpi_req(ireq),ierr) 
@@ -1485 +1573 @@
 
       IF (mpi_threading_mode==MPI_THREAD_FUNNELED .OR. mpi_threading_mode==MPI_THREAD_SINGLE) THEN
+!$acc wait
 !$OMP BARRIER
 !$OMP MASTER        
@@ -1492 +1581 @@
           msize=message%buffers(ireq)%size
           rank=message%buffers(ireq)%rank
-          CALL MPI_ISEND(buffer_r,msize,MPI_REAL8,rank,1000*message%number,          &
-            comm_icosa, message%mpi_req(ireq),ierr)
+          ! Using the "if" clause on the "host_data" directive seems to cause a crash at compilation
+          IF (message%ondevice) THEN
+            !$acc host_data use_device(buffer_r) ! if (message%ondevice)
+            CALL MPI_ISEND(buffer_r,msize,MPI_REAL8,rank,1000*message%number,          &
+              comm_icosa, message%mpi_req(ireq),ierr)
+            !$acc end host_data
+          ELSE
+            CALL MPI_ISEND(buffer_r,msize,MPI_REAL8,rank,1000*message%number,          &
+              comm_icosa, message%mpi_req(ireq),ierr)
+          ENDIF
         ENDDO
 
@@ -1500 +1597 @@
           msize=message%buffers(ireq)%size
           rank=message%buffers(ireq)%rank
-          CALL MPI_IRECV(buffer_r,msize,MPI_REAL8,rank,1000*message%number,          &
-            comm_icosa, message%mpi_req(ireq),ierr)
+          ! Using the "if" clause on the "host_data" directive seems to cause a crash at compilation
+          IF (message%ondevice) THEN
+            !$acc host_data use_device(buffer_r) ! if (message%ondevice)
+            CALL MPI_IRECV(buffer_r,msize,MPI_REAL8,rank,1000*message%number,          &
+              comm_icosa, message%mpi_req(ireq),ierr)
+            !$acc end host_data
+          ELSE
+            CALL MPI_IRECV(buffer_r,msize,MPI_REAL8,rank,1000*message%number,          &
+              comm_icosa, message%mpi_req(ireq),ierr)
+          ENDIF
         ENDDO
 
@@ -1552 +1657 @@
     INTEGER :: ind,n
     INTEGER :: dim3,dim4,d3,d4,lbegin,lend
-    INTEGER :: offset
+    INTEGER :: offset, msize, moffset
 
     message%open=.FALSE.
@@ -1586 +1691 @@
               sgn=>recv%sign
               offset=recv%offset
-              DO n=1,recv%size
+              msize=recv%size
+              !$acc parallel loop default(present) async if (field(ind)%ondevice)
+              DO n=1,msize
                 rval2d(value(n))=buffer_r(n+offset)*sgn(n)  
               ENDDO
@@ -1621 +1728 @@
     
               CALL distrib_level(1,dim3, lbegin,lend)
-              offset=recv%offset*dim3 + (lbegin-1)*recv%size
+              msize=recv%size
+              moffset=recv%offset
               CALL trace_start("copy_from_buffer")
               
               IF (req%vector) THEN
+                !$acc parallel loop default(present) async if (field(ind)%ondevice)
                 DO d3=lbegin,lend
-                  DO n=1,recv%size
+                  offset=moffset*dim3 + (d3-1)*msize
+                  !$acc loop
+                  DO n=1,msize
                     rval3d(value(n),d3)=buffer_r(n+offset)*sgn(n)  
                   ENDDO
-                  offset=offset+recv%size
                 ENDDO
               ELSE
+                !$acc parallel loop default(present) async if (field(ind)%ondevice)
                 DO d3=lbegin,lend
-                  DO n=1,recv%size
+                  offset=moffset*dim3 + (d3-1)*msize
+                  !$acc loop
+                  DO n=1,msize
                     rval3d(value(n),d3)=buffer_r(n+offset)  
                   ENDDO
-                  offset=offset+recv%size
                 ENDDO
               ENDIF
@@ -1670 +1782 @@
               CALL distrib_level(1,dim3, lbegin,lend)
               dim4=size(rval4d,3)
+              msize=recv%size
+              moffset=recv%offset
               CALL trace_start("copy_from_buffer")
+              !$acc parallel loop default(present) collapse(2) async if (field(ind)%ondevice)
               DO d4=1,dim4
-                offset=recv%offset*dim3*dim4 + dim3*recv%size*(d4-1) +               &
-                  (lbegin-1)*recv%size
                 DO d3=lbegin,lend
-                  DO n=1,recv%size
+                  offset=moffset*dim3*dim4 + dim3*msize*(d4-1) +               &
+                    (d3-1)*msize
+                  !$acc loop
+                  DO n=1,msize
                     rval4d(value(n),d3,d4)=buffer_r(n+offset)*sgn(n) 
                   ENDDO
-                  offset=offset+recv%size
                 ENDDO
               ENDDO
@@ -1697 +1812 @@
 !    CALL trace_end("wait_message_mpi")
 !$OMP BARRIER
-    
+  
     CALL exit_profile(id_mpi)
 

codes/icosagcm/trunk/src/physics/physics.f90

@@ -27 +27 @@
   SUBROUTINE init_physics
     USE mpipara
+    USE abort_mod
     USE etat0_mod
     USE physics_dcmip_mod, ONLY : init_physics_dcmip=>init_physics
@@ -118 +119 @@
     
 !$OMP END PARALLEL
+    IF (phys_type /= phys_none) THEN
+       CALL abort_acc("physics /= 'none'")
+    END IF
   END SUBROUTINE init_physics
 
   SUBROUTINE zero_du_phys()
+    USE abort_mod
     REAL(rstd), DIMENSION(:,:), POINTER :: du
     INTEGER :: ind
@@ -154 +159 @@
     USE etat0_heldsz_mod
     USE etat0_venus_mod, ONLY : phys_venus => physics
+    USE abort_mod
     INTEGER, INTENT(IN)   :: it
     TYPE(t_field),POINTER :: f_phis(:)

codes/icosagcm/trunk/src/sphere/geometry.f90

@@ -266 +266 @@
     INTEGER :: nb_it=0
     TYPE(t_domain),POINTER :: d
-    INTEGER :: ind,it,i,j,n,k
-    REAL(rstd) :: x1(3),x2(3)
+    INTEGER :: ind,it,i,j,n
     REAL(rstd) :: vect(3,6)
     REAL(rstd) :: centr(3)
@@ -690 +689 @@
 !    CALL surf_triangle(d%xyz(:,ii_begin,jj_begin),d%xyz(:,ii_begin,jj_end),d%xyz(:,ii_end,jj_begin),S)
  
+    CALL update_device_field(geom%Ai)
+
+    CALL update_device_field(geom%xyz_i)
+    CALL update_device_field(geom%lon_i)
+    CALL update_device_field(geom%lat_i)
+    CALL update_device_field(geom%elon_i)
+    CALL update_device_field(geom%elat_i)
+    CALL update_device_field(geom%centroid)
+
+    CALL update_device_field(geom%xyz_e)
+    CALL update_device_field(geom%lon_e)
+    CALL update_device_field(geom%lat_e)
+    CALL update_device_field(geom%elon_e)
+    CALL update_device_field(geom%elat_e)
+    CALL update_device_field(geom%ep_e)
+    CALL update_device_field(geom%et_e)
+
+    CALL update_device_field(geom%xyz_v)
+    CALL update_device_field(geom%de)
+    CALL update_device_field(geom%le)
+    CALL update_device_field(geom%le_de)
+    CALL update_device_field(geom%bi)
+    CALL update_device_field(geom%Av)
+    CALL update_device_field(geom%S1)
+    CALL update_device_field(geom%S2)
+    CALL update_device_field(geom%Riv)
+    CALL update_device_field(geom%Riv2)
+    CALL update_device_field(geom%ne)
+    CALL update_device_field(geom%Wee)
+    CALL update_device_field(geom%bi)
+    CALL update_device_field(geom%fv)
+
   END SUBROUTINE set_geometry
   

codes/icosagcm/trunk/src/time/euler_scheme.F90

@@ -1 +1 @@
 MODULE euler_scheme_mod
   USE field_mod
+  USE abort_mod
   IMPLICIT NONE
   PRIVATE
@@ -45 +46 @@
 
        IF(with_dps) THEN ! update ps/mass
+          CALL abort_acc("Euler_scheme/with_dps")
           IF(caldyn_eta==eta_mass) THEN ! update ps
              ps=f_ps(ind) ; dps=f_dps(ind) ;              
@@ -71 +73 @@
        du=f_du(ind) ; dtheta_rhodz=f_dtheta_rhodz(ind)
 
-       DO l=ll_begin,ll_end
+       CALL compute_euler_scheme(u, du, theta_rhodz, dtheta_rhodz)
+    ENDDO
+
+    CALL trace_end("Euler_scheme")
+
+    CONTAINS
+
+    SUBROUTINE compute_euler_scheme(u, du, theta_rhodz, dtheta_rhodz)
+       REAL(rstd),INTENT(INOUT) :: u(iim*3*jjm,llm)
+       REAL(rstd),INTENT(IN)    :: du(iim*3*jjm,llm)
+       REAL(rstd),INTENT(INOUT) :: theta_rhodz(iim*jjm,llm,nqdyn)
+       REAL(rstd),INTENT(IN)    :: dtheta_rhodz(iim*jjm,llm,nqdyn)
+
+       !$acc data present(theta_rhodz(:,:,:), u(:,:),du(:,:), dtheta_rhodz(:,:,:)) async
+
+       !$acc parallel loop async
+       DO l=ll_begin,ll_end
+          !$acc loop
           !DIR$ SIMD
           DO ij=ij_begin,ij_end
@@ -80 +99 @@
           ENDDO
        ENDDO
-    ENDDO
-
-    CALL trace_end("Euler_scheme")  
+       
+       !$acc end data
+    END SUBROUTINE compute_euler_scheme
 
   END SUBROUTINE Euler_scheme
@@ -97 +116 @@
     INTEGER :: l,ij
 
+    !$acc data present(hflux(:,:), wflux(:,:), hfluxt(:,:), wfluxt(:,:)) async
+
     IF(fluxt_zero) THEN
 
        fluxt_zero=.FALSE.
-
-       DO l=ll_begin,ll_end
+       !$acc parallel loop async
+       DO l=ll_begin,ll_end
+          !$acc loop
           !DIR$ SIMD
           DO ij=ij_begin_ext,ij_end_ext
@@ -111 +133 @@
 
        IF(caldyn_eta==eta_mass) THEN ! no need for vertical fluxes if eta_lag
+          !$acc parallel loop async
           DO l=ll_begin,ll_endp1
+             !$acc loop
              !DIR$ SIMD
              DO ij=ij_begin,ij_end
@@ -120 +144 @@
 
     ELSE
-
-       DO l=ll_begin,ll_end
+       !$acc parallel loop async
+       DO l=ll_begin,ll_end
+          !$acc loop
           !DIR$ SIMD
           DO ij=ij_begin_ext,ij_end_ext
@@ -131 +156 @@
 
        IF(caldyn_eta==eta_mass) THEN ! no need for vertical fluxes if eta_lag
+          !$acc parallel loop async
           DO l=ll_begin,ll_endp1
+             !$acc loop
              !DIR$ SIMD
              DO ij=ij_begin,ij_end
@@ -138 +165 @@
           ENDDO
        END IF
-
     END IF
-
+    !$acc end data
   END SUBROUTINE accumulate_fluxes
 
   SUBROUTINE legacy_to_DEC(f_ps, f_u)
+    USE icosa
+    USE disvert_mod
+    USE omp_para
+    USE trace
+    TYPE(t_field),POINTER :: f_ps(:), f_u(:)
+    REAL(rstd), POINTER :: ps(:), u(:,:)
+    INTEGER :: ind,ij,l
+
+    CALL trace_start("legacy_to_DEC")
+
+    DO ind=1,ndomain
+       IF (.NOT. assigned_domain(ind)) CYCLE
+       CALL swap_dimensions(ind)
+       CALL swap_geometry(ind)
+
+       IF(caldyn_eta==eta_mass .AND. is_omp_first_level) THEN ! update ps
+          ps=f_ps(ind)
+          !$acc parallel loop async default(present)
+          !DIR$ SIMD
+          DO ij=ij_begin,ij_end
+             ps(ij)=(ps(ij)-ptop)/g ! convert ps to column-integrated mass
+          ENDDO
+       END IF
+
+       u=f_u(ind)
+       !$acc parallel loop async default(present)
+       DO l=ll_begin,ll_end
+          !$acc loop
+          !DIR$ SIMD
+          DO ij=ij_begin,ij_end
+             u(ij+u_right,l)=u(ij+u_right,l)*de(ij+u_right)
+             u(ij+u_lup,l)=u(ij+u_lup,l)*de(ij+u_lup)
+             u(ij+u_ldown,l)=u(ij+u_ldown,l)*de(ij+u_ldown)
+          ENDDO
+       ENDDO
+    ENDDO
+
+    CALL trace_end("legacy_to_DEC")  
+
+  END SUBROUTINE Legacy_to_DEC
+
+  SUBROUTINE DEC_to_legacy(f_ps, f_u)
     USE icosa
     USE disvert_mod
@@ -158 +226 @@
        CALL swap_geometry(ind)
 
-       IF(caldyn_eta==eta_mass .AND. is_omp_first_level) THEN ! update ps
-          ps=f_ps(ind)
-          !DIR$ SIMD
-          DO ij=ij_begin,ij_end
-             ps(ij)=(ps(ij)-ptop)/g ! convert ps to column-integrated mass
-          ENDDO
-       END IF
-
-       u=f_u(ind)
-       DO l=ll_begin,ll_end
-          !DIR$ SIMD
-          DO ij=ij_begin,ij_end
-             u(ij+u_right,l)=u(ij+u_right,l)*de(ij+u_right)
-             u(ij+u_lup,l)=u(ij+u_lup,l)*de(ij+u_lup)
-             u(ij+u_ldown,l)=u(ij+u_ldown,l)*de(ij+u_ldown)
-          ENDDO
-       ENDDO
-    ENDDO
-
-    CALL trace_end("legacy_to_DEC")  
-  END SUBROUTINE Legacy_to_DEC
-
-  SUBROUTINE DEC_to_legacy(f_ps, f_u)
-    USE icosa
-    USE disvert_mod
-    USE omp_para
-    USE trace
-    TYPE(t_field),POINTER :: f_ps(:), f_u(:)
-    REAL(rstd), POINTER :: ps(:), u(:,:)
-    INTEGER :: ind,ij,l
-    CALL trace_start("legacy_to_DEC")
-
-    DO ind=1,ndomain
-       IF (.NOT. assigned_domain(ind)) CYCLE
-       CALL swap_dimensions(ind)
-       CALL swap_geometry(ind)
-
        IF(caldyn_eta==eta_mass .AND. is_omp_first_level) THEN
           ps=f_ps(ind)
+          !$acc parallel loop async default(present)
           !DIR$ SIMD
           DO ij=ij_begin,ij_end
@@ -204 +236 @@
        
        u=f_u(ind)
-       DO l=ll_begin,ll_end
+       !$acc parallel loop async default(present)
+       DO l=ll_begin,ll_end
+          !$acc loop
           !DIR$ SIMD
           DO ij=ij_begin,ij_end

codes/icosagcm/trunk/src/time/hevi_scheme.F90

@@ -4 +4 @@
   USE field_mod
   USE euler_scheme_mod
+  USE abort_mod
   IMPLICIT NONE
   PRIVATE
@@ -13 +14 @@
 
 CONTAINS
-
   SUBROUTINE set_coefs_ark23(dt)
     ! ARK2 scheme by Giraldo, Kelly, Constantinescu 2013
@@ -60 +60 @@
 
     CALL legacy_to_DEC(f_ps, f_u)
+
     DO j=1,nb_stage
        CALL caldyn_hevi((j==1) .AND. (MOD(it,itau_out)==0), taujj(j), &
@@ -76 +77 @@
           wflux=f_wflux(ind);     wfluxt=f_wfluxt(ind)
           CALL accumulate_fluxes(hflux,wflux, hfluxt,wfluxt, wj(j), fluxt_zero(ind))
+
        END DO
        ! update model state
@@ -88 +90 @@
           CALL update_3D(cjl(l,j), f_u, f_du_fast(:,l))
           IF(.NOT. hydrostatic) THEN
+             CALL abort_acc("HEVI_scheme/!hydrostatic")
              CALL update_3D(bjl(l,j), f_W, f_dW_slow(:,l))
              CALL update_3D(cjl(l,j), f_W, f_dW_fast(:,l))
@@ -93 +96 @@
              CALL update_3D(cjl(l,j), f_geopot, f_dPhi_fast(:,l))
           END IF
+
 !$OMP BARRIER       
        END DO
@@ -142 +146 @@
     INTENT(IN) :: dy
     INTEGER :: l
+    !$acc kernels async default(present)
     DO l=ll_begin,ll_end
        y(:,l)=y(:,l)+w*dy(:,l)
     ENDDO
+    !$acc end kernels
   END SUBROUTINE compute_update_3D
   
@@ -164 +170 @@
     INTENT(INOUT) :: y
     INTENT(IN) :: dy
+    !$acc kernels async default(present)
     y(:)=y(:)+w*dy(:)
+    !$acc end kernels
   END SUBROUTINE compute_update_2D
   

codes/icosagcm/trunk/src/time/timeloop_gcm.F90

@@ -94 +94 @@
     END SELECT
     
+    IF (scheme_family /= hevi) THEN
+       CALL abort_acc("scheme_family /= hevi")
+    END IF
+
     ! Time-independant orography
     CALL allocate_field(f_phis,field_t,type_real,name='phis')
@@ -103 +107 @@
     CALL allocate_field(f_u,field_u,type_real,llm,name='u')
     CALL allocate_field(f_geopot,field_t,type_real,llm+1,name='geopot')
-    CALL allocate_field(f_W,field_t,type_real,llm+1,name='W')
+    CALL allocate_field(f_W,field_t,type_real,llm+1,name='W') ! used only if .not. hydrostatic
     CALL allocate_field(f_q,field_t,type_real,llm,nqtot,'q')
     ! Mass fluxes
-    CALL allocate_field(f_hflux,field_u,type_real,llm)    ! instantaneous mass fluxes
-    CALL allocate_field(f_hfluxt,field_u,type_real,llm)   ! mass "fluxes" accumulated in time
+    CALL allocate_field(f_hflux,field_u,type_real,llm, ondevice=.TRUE.)    ! instantaneous mass fluxes
+    CALL allocate_field(f_hfluxt,field_u,type_real,llm,ondevice=.TRUE.)   ! mass "fluxes" accumulated in time
     CALL allocate_field(f_wflux,field_t,type_real,llm+1)  ! vertical mass fluxes
-    CALL allocate_field(f_wfluxt,field_t,type_real,llm+1,name='wfluxt')
+    CALL allocate_field(f_wfluxt,field_t,type_real,llm+1,name='wfluxt',ondevice=.TRUE.)
     
     SELECT CASE(scheme_family)
@@ -125 +129 @@
     CASE(hevi)
        ! Trends
-       CALL allocate_fields(nb_stage,f_dps_slow, field_t,type_real,name='dps_slow')
-       CALL allocate_fields(nb_stage,f_dmass_slow, field_t,type_real,llm, name='dmass_slow')
-       CALL allocate_fields(nb_stage,f_dtheta_rhodz_slow, field_t,type_real,llm,nqdyn,name='dtheta_rhodz_fast')
-       CALL allocate_fields(nb_stage,f_du_slow, field_u,type_real,llm,name='du_slow')
-       CALL allocate_fields(nb_stage,f_du_fast, field_u,type_real,llm,name='du_fast')
+       CALL allocate_fields(nb_stage,f_dps_slow, field_t,type_real,name='dps_slow', ondevice=.TRUE.)
+       CALL allocate_fields(nb_stage,f_dmass_slow, field_t,type_real,llm, name='dmass_slow', ondevice=.TRUE.)
+       CALL allocate_fields(nb_stage,f_dtheta_rhodz_slow, field_t,type_real,llm,nqdyn,name='dtheta_rhodz_fast', ondevice=.TRUE.)
+       CALL allocate_fields(nb_stage,f_du_slow, field_u,type_real,llm,name='du_slow', ondevice=.TRUE.)
+       CALL allocate_fields(nb_stage,f_du_fast, field_u,type_real,llm,name='du_fast', ondevice=.TRUE.)
        CALL allocate_fields(nb_stage,f_dW_slow, field_t,type_real,llm+1,name='dW_slow')
        CALL allocate_fields(nb_stage,f_dW_fast, field_t,type_real,llm+1,name='dW_fast')
@@ -172 +176 @@
   
   SUBROUTINE timeloop
+    USE abort_mod
     USE dissip_gcm_mod
     USE sponge_mod
@@ -211 +216 @@
        rhodz=f_rhodz(ind); mass=f_mass(ind); ps=f_ps(ind)
        IF(caldyn_eta==eta_mass) THEN
-          CALL compute_rhodz(.TRUE., ps, rhodz) ! save rhodz for transport scheme before dynamics update ps
+          CALL compute_rhodz(.TRUE., ps, rhodz, ondevice=.FALSE.) ! save rhodz for transport scheme before dynamics update ps
        ELSE
           DO l=ll_begin,ll_end
@@ -244 +249 @@
     CALL SYSTEM_CLOCK(start_clock, rate_clock)
     !$OMP END MASTER   
+    call update_device_field(f_ps)
+    call update_device_field(f_mass)
+    CALL update_device_field(f_theta_rhodz)
+    CALL update_device_field(f_u)
+    CALL update_device_field(f_q)
+    CALL update_device_field(f_geopot)
+    CALL update_device_field(f_wflux)
+    CALL update_device_field(f_rhodz)
+
 
     DO it=itau0+1,itau0+itaumax
@@ -263 +277 @@
           CALL wait_message(req_mass0)
           CALL send_message(f_theta_rhodz,req_theta_rhodz0) 
-          CALL wait_message(req_theta_rhodz0) 
+          CALL wait_message(req_theta_rhodz0)
           CALL send_message(f_u,req_u0)
           CALL wait_message(req_u0)
@@ -281 +295 @@
        SELECT CASE(scheme_family)
        CASE(explicit)
+          CALL abort_acc("explicit_scheme")
           CALL explicit_scheme(it, fluxt_zero)
        CASE(hevi)
@@ -298 +313 @@
                 CALL swap_geometry(ind)
                 mass=f_mass(ind); ps=f_ps(ind);
-                CALL compute_rhodz(.TRUE., ps, mass)
+                CALL compute_rhodz(.TRUE., ps, mass, ondevice=.TRUE.)
              END DO
           ENDIF
@@ -311 +326 @@
           CALL euler_scheme(.FALSE.)  ! update only u, theta
           IF (iflag_sponge > 0) THEN
+             CALL abort_acc("iflag_sponge>0")
              CALL sponge(f_u,f_du,f_theta_rhodz,f_dtheta_rhodz)
              CALL euler_scheme(.FALSE.)  ! update only u, theta
@@ -321 +337 @@
        END IF
        CALL exit_profile(id_dissip)
-       
+
        CALL enter_profile(id_adv)
        IF(MOD(it,itau_adv)==0) THEN
@@ -329 +345 @@
           ! At this point advect_tracer has obtained the halos of u and rhodz,
           ! needed for correct computation of kinetic energy
+          IF(diagflux_on) CALL abort_acc("diagflux_on")
           IF(diagflux_on) CALL diagflux_energy(adv_over_out, f_phis,f_rhodz,f_theta_rhodz,f_u, f_geopot,f_theta,f_buf_i, f_hfluxt)
 
@@ -340 +357 @@
              END DO
           ENDIF
+          IF(positive_theta) CALL abort_acc("positive_theta")
           IF(positive_theta) CALL copy_q_to_theta(f_theta_rhodz,f_rhodz,f_q)
        END IF
        CALL exit_profile(id_adv)
-       
+
        CALL enter_profile(id_diags)
 !       IF (MOD(it,itau_physics)==0) THEN
@@ -360 +378 @@
 
        IF (MOD(it,itau_check_conserv)==0) THEN
+          CALL update_host_field(f_ps)
+          CALL update_host_field(f_theta_rhodz)
+          CALL update_host_field(f_u)
+          CALL update_host_field(f_dps)
+          CALL update_host_field(f_q)
           CALL check_conserve_detailed(it, AAM_dyn, &
                f_ps,f_dps,f_u,f_theta_rhodz,f_phis)
@@ -367 +390 @@
        IF (mod(it,itau_out)==0 ) THEN
           CALL transfert_request(f_u,req_e1_vect)
+          CALL update_host_field(f_ps)              
+          CALL update_host_field(f_mass)
+          CALL update_host_field(f_theta_rhodz)
+          CALL update_host_field(f_geopot)
+          CALL update_host_field(f_u)
+          CALL update_host_field(f_q)
           CALL write_output_fields_basic(.FALSE.,f_phis, f_ps, f_mass, f_geopot, f_theta_rhodz, f_u, f_W, f_q)
        ENDIF
@@ -374 +403 @@
     END DO
     
+    CALL update_host_field(f_ps)
+    CALL update_host_field(f_theta_rhodz)
+    CALL update_host_field(f_u)
+    CALL update_host_field(f_q)
+    CALL update_host_field(f_geopot)
+
 !    CALL write_etat0(itau0+itaumax,f_ps, f_phis,f_theta_rhodz,f_u,f_q) 
     CALL write_etat0(itau0+itaumax,f_ps, f_phis,f_theta_rhodz,f_u,f_q, f_geopot, f_W) 
-    
+
+    CALL update_host_field(f_dps)    
     CALL check_conserve_detailed(it, AAM_dyn, &
          f_ps,f_dps,f_u,f_theta_rhodz,f_phis)

codes/icosagcm/trunk/src/transport/advect.F90

@@ -49 +49 @@
   !=======================================================================================
 
-  SUBROUTINE compute_gradq3d(qi,sqrt_leng,gradq3d,xyz_i,xyz_v)
+  SUBROUTINE compute_gradq3d(qi,sqrt_leng,gradq3d,xyz_i, xyz_v)
     USE trace
     USE omp_para
@@ -68 +68 @@
     REAL(rstd) :: x1,x2,x3
     REAL(rstd) :: dq(3)
-
+    
     CALL trace_start("compute_gradq3d1")
 
@@ -86 +86 @@
 !    END DO
 
-     DO l = ll_begin,ll_end 
+     !$acc data create(gradtri(:,:,:), arr(:), ar(:)) present(sqrt_leng(:), xyz_i(:,:), xyz_v(:,:), qi(:,:), gradq3d(:,:,:)) async
+
+     !$acc parallel loop collapse(2) async private(A, dq)
+     DO l = ll_begin,ll_end
 !DIR$ SIMD
       DO ij=ij_begin_ext,ij_end_ext
@@ -151 +154 @@
         
       ENDDO
+     ENDDO
       
+     !$acc parallel loop collapse(2) async private(A, dq)
+     DO l = ll_begin,ll_end
       DO ij=ij_begin_ext,ij_end_ext
 
@@ -219 +225 @@
 
 !DIR$ SIMD
+    !$acc parallel loop async
     DO ij=ij_begin,ij_end
        ar(ij) = arr(ij+z_up)+arr(ij+z_lup)+arr(ij+z_ldown)+arr(ij+z_down)+arr(ij+z_rdown)+arr(ij+z_rup)+1.e-50
@@ -225 +232 @@
     CALL trace_start2("compute_gradq3d2")
       
+    !$acc parallel loop collapse(3) async
     DO k=1,3
       DO l =ll_begin,ll_end
@@ -239 +247 @@
 
     !============================================================================================= LIMITING 
+    !$acc parallel loop collapse(2) async
     DO l =ll_begin,ll_end
 !DIR$ SIMD
@@ -251 +260 @@
              minq = min(qi(ij,l),qi(ij+t_right,l),qi(ij+t_lup,l),qi(ij+t_rup,l),qi(ij+t_left,l), &
                   qi(ij+t_rdown,l),qi(ij+t_ldown,l))
-             alphamx = (maxq - qi(ij,l)) ; alphamx = alphamx/(maxq_c - qi(ij,l) )
-             alphamx = max(alphamx,0.0)
-             alphami = (minq - qi(ij,l)); alphami = alphami/(minq_c - qi(ij,l))
-             alphami = max(alphami,0.0) 
+             IF ((maxq_c - qi(ij,l)) /= 0.0) THEN
+               alphamx = (maxq - qi(ij,l)) ; alphamx = alphamx/(maxq_c - qi(ij,l) )
+               alphamx = max(alphamx,0.0)
+             ELSE
+               alphamx = 0.0
+             ENDIF
+             IF ((minq_c - qi(ij,l)) /= 0.0) THEN
+               alphami = (minq - qi(ij,l)); alphami = alphami/(minq_c - qi(ij,l))
+               alphami = max(alphami,0.0)
+             ELSE
+               alphami = 0.0
+             ENDIF
              alpha   = min(alphamx,alphami,1.0)
 !             gradq3d(ij,l,:) = alpha*gradq3d(ij,l,:)
@@ -264 +281 @@
 
   CALL trace_end("compute_gradq3d3")
+
+  !$acc end data
   
   CONTAINS
@@ -329 +348 @@
 
   ! Backward trajectories, for use with Miura approach
-  SUBROUTINE compute_backward_traj(normal,tangent,ue,tau, cc)
+  SUBROUTINE compute_backward_traj(normal,tangent,ue,tau,cc)
+    USE geometry, ONLY : xyz_e, de, wee, le
     USE trace
     USE omp_para
@@ -345 +365 @@
 
     ! TODO : compute normal displacement ue*tau as hfluxt / mass(upwind) then reconstruct tangential displacement
-    
+
+    !$acc data present(ue(:,:), cc(:,:,:), normal(:,:), tangent(:,:), xyz_e(:,:), de(:), wee(:,:,:), le(:)) async
+
     ! reconstruct tangential wind then 3D wind at edge then cc = edge midpoint - u*tau
+    !$acc parallel loop private(up_e, v_e) collapse(2) gang vector async
     DO l = ll_begin,ll_end
 !DIR$ SIMD
@@ -397 +420 @@
        ENDDO
     END DO
-
+    !$acc end data
     CALL trace_end("compute_backward_traj")
 
@@ -404 +427 @@
   ! Horizontal transport (S. Dubey, T. Dubos)
   ! Slope-limited van Leer approach with hexagons
-  SUBROUTINE compute_advect_horiz(update_mass,diagflux_on, hfluxt,cc,gradq3d, mass,qi,qfluxt)
+  SUBROUTINE compute_advect_horiz(update_mass,diagflux_on, hfluxt,cc,gradq3d, mass, qi, qfluxt)
     USE trace
     USE omp_para
+    USE abort_mod
+    USE geometry, only : Ai, xyz_i
     IMPLICIT NONE
     LOGICAL, INTENT(IN)       :: update_mass, diagflux_on
@@ -415 +440 @@
     REAL(rstd), INTENT(INOUT) :: qi(iim*jjm,llm)
     REAL(rstd), INTENT(INOUT) :: qfluxt(3*iim*jjm,MERGE(llm,1,diagflux_on)) ! time-integrated tracer flux
-
-    REAL(rstd) :: dq,dmass,qe, newmass
+! metrics terms
+    
+    REAL(rstd) :: dq,dmass,qe,newmass
     REAL(rstd) :: qflux(3*iim*jjm,llm)
-    INTEGER :: ij,l
+    INTEGER :: ij,l,ij_tmp
+
+    IF(diagflux_on) CALL abort_acc("compute_advect_horiz : diagflux_on")
 
     CALL trace_start("compute_advect_horiz")
 #include "../kernels/advect_horiz.k90"
     CALL trace_end("compute_advect_horiz")
+
   END SUBROUTINE compute_advect_horiz
 

codes/icosagcm/trunk/src/transport/advect_tracer.F90

@@ -32 +32 @@
     INTEGER :: ind
 
-    CALL allocate_field(f_normal,field_u,type_real,3, name='normal')
-    CALL allocate_field(f_tangent,field_u,type_real,3, name='tangent')
-    CALL allocate_field(f_gradq3d,field_t,type_real,llm,3, name='gradq3d')
-    CALL allocate_field(f_cc,field_u,type_real,llm,3, name='cc')
-    CALL allocate_field(f_sqrt_leng,field_t,type_real, name='sqrt_leng')
-    CALL allocate_field(f_dzqw, field_t, type_real, llm, name='dzqw')
-    CALL allocate_field(f_adzqw, field_t, type_real, llm, name='adzqw')
-    CALL allocate_field(f_dzq, field_t, type_real, llm, name='dzq')
-    CALL allocate_field(f_wq, field_t, type_real, llm+1, name='wq')
+    CALL allocate_field(f_normal,field_u,type_real,3, name='normal',ondevice=.TRUE.)
+    CALL allocate_field(f_tangent,field_u,type_real,3, name='tangent',ondevice=.TRUE.)
+    CALL allocate_field(f_gradq3d,field_t,type_real,llm,3, name='gradq3d',ondevice=.TRUE.)
+    CALL allocate_field(f_cc,field_u,type_real,llm,3, name='cc',ondevice=.TRUE.)
+    CALL allocate_field(f_sqrt_leng,field_t,type_real, name='sqrt_leng',ondevice=.TRUE.)
+    CALL allocate_field(f_dzqw, field_t, type_real, llm, name='dzqw',ondevice=.TRUE.)
+    CALL allocate_field(f_adzqw, field_t, type_real, llm, name='adzqw',ondevice=.TRUE.)
+    CALL allocate_field(f_dzq, field_t, type_real, llm, name='dzq',ondevice=.TRUE.)
+    CALL allocate_field(f_wq, field_t, type_real, llm+1, name='wq',ondevice=.TRUE.)
     
     DO ind=1,ndomain
@@ -52 +52 @@
     END DO
 
+    CALL update_device_field(f_tangent)
+    CALL update_device_field(f_normal)
+    CALL update_device_field(f_sqrt_leng)
+
   END SUBROUTINE init_advect_tracer
 
@@ -60 +64 @@
     USE write_field_mod
     USE tracer_mod
+    USE abort_mod
     TYPE(t_field),POINTER :: f_hfluxt(:)   ! time-integrated horizontal mass flux
     TYPE(t_field),POINTER :: f_wfluxt(:)   ! time-integrated vertical mass flux
@@ -138 +143 @@
        END DO
 
-       CALL compute_backward_traj(tangent,normal,u,0.5*dt*itau_adv, cc) 
-
+       CALL compute_backward_traj(tangent,normal,u,0.5*dt*itau_adv, cc)
     END DO
 
@@ -174 +178 @@
 
           IF(frac>0.) THEN ! accumulate mass, mass flux and tracer mass
+             CALL abort_acc("frac>0.")
              qmasst  = f_qmasst(ind) 
+             !$acc kernels default(present) async
              qmasst(:,ll_begin:ll_end,k) = qmasst(:,ll_begin:ll_end,k) + &
                   frac*rhodz(:,ll_begin:ll_end)*q(:,ll_begin:ll_end,k)
+             !$acc end kernels
              IF(k==nq_last) THEN
                 masst  = f_masst(ind)
                 massfluxt  = f_massfluxt(ind)
+                !$acc kernels default(present) async
                 masst(:,ll_begin:ll_end) = masst(:,ll_begin:ll_end)+frac*rhodz(:,ll_begin:ll_end)
                 massfluxt(:,ll_begin:ll_end) = massfluxt(:,ll_begin:ll_end)+hfluxt(:,ll_begin:ll_end)
+                !$acc end kernels
              END IF
           END IF
@@ -189 +198 @@
       ENDIF
     END DO 
-    
+
     ! 1/2 vertical transport
 !!$OMP BARRIER
@@ -208 +217 @@
          IF (advection_scheme(k)==advect_vanleer) CALL vlz(k==nq_last, 0.5,wfluxt,rhodz, q(:,:,k),0, dzqw, adzqw, dzq, wq)
        END DO
-
     END DO
 
@@ -251 +259 @@
     ! finite difference of q
 
+     !$acc data present(q(:,:), mass(:,:), wfluxt(:,:), dzqw(:,:), adzqw(:,:), dzq(:,:), wq(:,:)) async
+
+     !$acc parallel loop async collapse(2)
      DO l=ll_beginp1,ll_end
 !DIR$ SIMD
@@ -265 +276 @@
     ! dzq = slope*dz, i.e. the reconstructed q varies by dzq inside level l
 
+     !$acc parallel loop async collapse(2)
      DO l=ll_beginp1,ll_endm1
 !DIR$ SIMD
@@ -281 +293 @@
     ! 0 slope in top and bottom layers
     IF (is_omp_first_level) THEN
+      !$acc parallel loop async
       DO ij=ijb,ije
            dzq(ij,1)=0.
@@ -287 +300 @@
       
     IF (is_omp_last_level) THEN
+      !$acc parallel loop async
       DO ij=ijb,ije
           dzq(ij,llm)=0.
@@ -297 +311 @@
     ! sigw = fraction of mass that leaves level l/l+1
     ! then amount of q leaving level l/l+1 = wq = w * qq
+     !$acc parallel loop collapse(2) async
      DO l=ll_beginp1,ll_end
 !DIR$ SIMD
@@ -313 +328 @@
     ! wq = 0 at top and bottom
     IF (is_omp_first_level) THEN
+      !$acc parallel loop async
        DO ij=ijb,ije
             wq(ij,1)=0.
@@ -319 +335 @@
     
     IF (is_omp_last_level) THEN
+      !$acc parallel loop async
       DO ij=ijb,ije
             wq(ij,llm+1)=0.
@@ -329 +346 @@
 
     ! update q, mass is updated only after all q's have been updated
+    !$acc parallel loop collapse(2) async
     DO l=ll_begin,ll_end
 !DIR$ SIMD
@@ -337 +355 @@
        ENDDO
     END DO
-
+    !$acc end data
     CALL trace_end("vlz")
 

codes/icosagcm/trunk/src/vertical/disvert.f90

@@ -1 +1 @@
 MODULE disvert_mod
   USE icosa
+  USE abort_mod
   REAL(rstd), SAVE, POINTER :: ap(:)
 !$OMP THREADPRIVATE(ap) 
@@ -146 +147 @@
     END IF
 
+    !$acc enter data copyin(ap(:), bp(:), mass_dak(:), mass_dbk(:)) async
+    
   END SUBROUTINE init_disvert  
   
-  SUBROUTINE compute_rhodz(comp, ps, rhodz)
+  SUBROUTINE compute_rhodz(comp, ps, rhodz, ondevice)
     USE icosa
     USE omp_para
@@ -154 +157 @@
     REAL(rstd), INTENT(IN) :: ps(iim*jjm)
     REAL(rstd), INTENT(INOUT) :: rhodz(iim*jjm,llm)
+    LOGICAL, INTENT(IN), OPTIONAL :: ondevice ! .TRUE. compute on device, .FALSE. stay on host
     REAL(rstd) :: m, err
     INTEGER :: l,i,j,ij,dd
+    LOGICAL :: ondevice_
+
     err=0.
+
+    ! by default, do not compute on device
+    ondevice_ = .FALSE.
+    IF (PRESENT(ondevice)) ondevice_ = ondevice
 
     IF(comp) THEN
@@ -164 +174 @@
     END IF
 
+    !$acc data present(ps(:), rhodz(:,:), ap(:), bp(:)) async if (ondevice_)
+
     IF(comp) THEN
-    DO l = ll_begin, ll_end
-       DO j=jj_begin-dd,jj_end+dd
-          DO i=ii_begin-dd,ii_end+dd
-             ij=(j-1)*iim+i
-             m = ( ap(l) - ap(l+1) + (bp(l)-bp(l+1))*ps(ij) )/g 
-                rhodz(ij,l) = m
-          ENDDO
-        ENDDO
-      ENDDO
-             ELSE
+      !$acc parallel loop collapse(3) async if (ondevice_)
       DO l = ll_begin, ll_end
         DO j=jj_begin-dd,jj_end+dd
@@ -180 +183 @@
             ij=(j-1)*iim+i
             m = ( ap(l) - ap(l+1) + (bp(l)-bp(l+1))*ps(ij) )/g
-                err = MAX(err,abs(m-rhodz(ij,l)))
+            rhodz(ij,l) = m
           ENDDO
-       ENDDO
-    ENDDO
-
-       IF(err>1e-10) THEN
+        ENDDO
+      ENDDO
+    ELSE
+      !$acc parallel loop reduction(max:err) collapse(3) async if (ondevice_)
+      DO l = ll_begin, ll_end
+        DO j=jj_begin-dd,jj_end+dd
+          DO i=ii_begin-dd,ii_end+dd
+            ij=(j-1)*iim+i
+            m = ( ap(l) - ap(l+1) + (bp(l)-bp(l+1))*ps(ij) )/g
+            err = MAX(err,abs(m-rhodz(ij,l)))
+          ENDDO
+        ENDDO
+      ENDDO
+
+      !$acc wait
+      IF(err>1e-10) THEN
           PRINT *, 'Discrepancy between ps and rhodz detected', err
           STOP
@@ -191 +206 @@
     ENDIF
 
+    !$acc end data
   END SUBROUTINE compute_rhodz
 
- 
+
   SUBROUTINE write_apbp
   USE icosa