Ignore:
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/src/dynamics
Files:
3 edited
1 moved

Legend:

Unmodified
Added
Removed

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

    r899 r953  
    1515   
    1616  SUBROUTINE init_caldyn 
     17    USE abort_mod 
    1718    CHARACTER(len=255) :: def 
    1819    INTEGER            :: ind 
    … …  
    2122    hydrostatic=.TRUE. 
    2223    CALL getin("hydrostatic",hydrostatic) 
    23    
     24 
     25    IF (.NOT. hydrostatic) THEN 
     26       CALL abort_acc("hydrostatic /= .TRUE.") 
     27    END IF 
     28 
    2429    dysl=.NOT.hydrostatic ! dysl defaults to .FALSE. if hydrostatic, .TRUE. if NH                                               
    2530    CALL getin("dysl",dysl) 
    … …  
    3439    dysl_caldyn_coriolis=dysl 
    3540    CALL getin("dysl_caldyn_coriolis",dysl_caldyn_coriolis) 
     41 
     42    ! dysl is not supported with openACC 
     43    IF (dysl_geopot .OR. dysl_caldyn_fast .OR. dysl_slow_hydro .OR. dysl_pvort_only .OR. dysl_caldyn_coriolis) THEN 
     44       CALL abort_acc("dysl not supported") 
     45    END IF 
    3646 
    3747    def='energy' 
    … …  
    122132  SUBROUTINE allocate_caldyn 
    123133    CALL allocate_field(f_out_u,field_u,type_real,llm)  
    124     CALL allocate_field(f_qu,field_u,type_real,llm)  
    125     CALL allocate_field(f_qv,field_z,type_real,llm)  
    126     CALL allocate_field(f_pk,    field_t,type_real,llm,  name='pk') 
    127     CALL allocate_field(f_wwuu,  field_u,type_real,llm+1,name='wwuu')     
     134    CALL allocate_field(f_qu,field_u,type_real,llm, ondevice=.TRUE.)  
     135    CALL allocate_field(f_qv,field_z,type_real,llm, ondevice=.TRUE.)  
     136    CALL allocate_field(f_pk,    field_t,type_real,llm,  name='pk', ondevice=.TRUE.) 
     137    CALL allocate_field(f_wwuu,  field_u,type_real,llm+1,name='wwuu', ondevice=.TRUE.)     
    128138    CALL allocate_field(f_planetvel,  field_u,type_real, name='planetvel') ! planetary velocity at r=a 
    129139    IF(.NOT.hydrostatic) THEN 

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

    r933 r953  
    2626    USE output_field_mod 
    2727    USE checksum_mod 
     28    USE abort_mod 
    2829    IMPLICIT NONE 
    2930    LOGICAL,INTENT(IN)    :: write_out 
    … …  
    8687       CALL wait_message(req_ps) ! COM00 
    8788    ELSE 
     89       CALL abort_acc("HEVI_scheme/!eta_mass") 
    8890       CALL send_message(f_mass,req_mass) ! COM00 
    8991       CALL wait_message(req_mass) ! COM00 
    … …  
    9395 
    9496    IF(.NOT.hydrostatic) THEN 
     97       CALL abort_acc("HEVI_scheme/!hydrostatic") 
    9598       CALL send_message(f_geopot,req_geopot) ! COM03 
    9699       CALL wait_message(req_geopot) ! COM03 
    … …  
    112115       du=f_du_fast(ind) 
    113116       IF(hydrostatic) THEN 
    114           du(:,:)=0. 
     117          !$acc kernels present(du) async 
     118          du(:,:)=0.0d0 
     119          !$acc end kernels 
    115120          CALL compute_geopot(mass,theta, ps,pk,geopot) 
    116121       ELSE 
     122          CALL abort_acc("HEVI_scheme/!hydrostatic") 
    117123          phis = f_phis(ind) 
    118124          W = f_W(ind) 
    … …  
    161167          CALL compute_caldyn_slow_hydro(u,mass,hflux,du, .TRUE.) 
    162168       ELSE 
     169          CALL abort_acc("HEVI_scheme/!hydrostatic") 
    163170          W = f_W(ind) 
    164171          dW = f_dW_slow(ind) 
    … …  
    170177          CALL compute_caldyn_slow_NH(u,mass,geopot,W, F_el,gradPhi2,w_il, hflux,du,dPhi,dW) 
    171178       END IF 
    172        CALL compute_caldyn_Coriolis(hflux,theta,qu, convm,dtheta_rhodz,du) 
     179       CALL compute_caldyn_Coriolis(hflux,theta,qu,convm,dtheta_rhodz,du) 
     180        
    173181       IF(caldyn_eta==eta_mass) THEN 
    174182          wflux=f_wflux(ind) 
    … …  
    177185          CALL compute_caldyn_vert(u,theta,mass,convm, wflux,wwuu, dps, dtheta_rhodz, du) 
    178186          IF(.NOT.hydrostatic) THEN 
     187             CALL abort_acc("HEVI_scheme/!hydrostatic") 
    179188             W_etadot=f_Wetadot(ind) 
    180189             CALL compute_caldyn_vert_NH(mass,geopot,W,wflux, W_etadot, du,dPhi,dW) 

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

    r899 r953  
    55  USE omp_para 
    66  USE trace 
     7  USE abort_mod 
    78  IMPLICIT NONE 
    89  PRIVATE 
    … …  
    2526 
    2627  !**************************** Geopotential ***************************** 
    27  
    2828  SUBROUTINE compute_geopot(rhodz,theta, ps,pk,geopot)  
    2929    REAL(rstd),INTENT(IN)    :: rhodz(iim*jjm,llm) 
    … …  
    4646 
    4747    IF(dysl_geopot) THEN 
     48      CALL abort_acc("HEVI_scheme/!dysl_geopot") 
    4849#include "../kernels/compute_geopot.k90" 
    4950    ELSE 
    … …  
    5253    ! Works also when caldyn_eta=eta_mass           
    5354 
     55    !$acc data present(rhodz(:,:), ps(:), pk(:,:), geopot(:,:), theta(:,:,:)) async 
     56 
    5457    IF(boussinesq) THEN ! compute geopotential and pk=Lagrange multiplier 
     58       CALL abort_acc("HEVI_scheme/boussinesq") 
    5559       ! specific volume 1 = dphi/g/rhodz 
    5660       !         IF (is_omp_level_master) THEN ! no openMP on vertical due to dependency 
    5761       DO l = 1,llm 
     62          !$acc parallel loop 
    5863          !DIR$ SIMD 
    5964          DO ij=ij_omp_begin_ext,ij_omp_end_ext          
    6065             geopot(ij,l+1) = geopot(ij,l) + g*rhodz(ij,l) 
    6166          ENDDO 
     67          !$acc end parallel loop 
    6268       ENDDO 
    6369       ! use hydrostatic balance with theta*rhodz to find pk (Lagrange multiplier=pressure)  
    6470       ! uppermost layer 
     71       !$acc parallel loop 
    6572       !DIR$ SIMD 
    6673       DO ij=ij_begin_ext,ij_end_ext          
    6774          pk(ij,llm) = ptop + (.5*g)*theta(ij,llm,1)*rhodz(ij,llm) 
    6875       END DO 
     76       !$acc end parallel loop 
    6977       ! other layers 
    7078       DO l = llm-1, 1, -1 
    7179          !          !$OMP DO SCHEDULE(STATIC)  
     80          !$acc parallel loop 
    7281          !DIR$ SIMD 
    7382          DO ij=ij_begin_ext,ij_end_ext          
    7483             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)) 
    7584          END DO 
     85          !$acc end parallel loop 
    7686       END DO 
    7787       ! now pk contains the Lagrange multiplier (pressure) 
    … …  
    8191       SELECT CASE(caldyn_thermo) 
    8292          CASE(thermo_theta, thermo_entropy) 
     93             !$acc parallel loop async 
    8394             !DIR$ SIMD 
    8495             DO ij=ij_omp_begin_ext,ij_omp_end_ext 
    8596                pk(ij,llm) = ptop + (.5*g)*rhodz(ij,llm) 
    8697             END DO 
     98 
    8799             ! other layers 
     100             ! We use kernels here instead of "loop seq" + "loop gang vector" 
     101             ! to be sure the code really abides the standard. In practice, 
     102             ! it seems like the compiler interchanges the loops. 
     103             !$acc kernels async 
    88104             DO l = llm-1, 1, -1 
    89105                !DIR$ SIMD 
    … …  
    92108                END DO 
    93109             END DO 
     110             !$acc end kernels 
    94111             ! surface pressure (for diagnostics) 
    95112             IF(caldyn_eta==eta_lag) THEN 
     113                !$acc parallel loop async 
    96114                DO ij=ij_omp_begin_ext,ij_omp_end_ext          
    97115                   ps(ij) = pk(ij,1) + (.5*g)*rhodz(ij,1) 
    … …  
    99117             END IF 
    100118          CASE(thermo_moist) ! theta(ij,l,2) = qv = mv/md 
     119             CALL abort_acc("HEVI_scheme/thermo_moist") 
     120             !$acc parallel loop 
    101121             !DIR$ SIMD 
    102122             DO ij=ij_omp_begin_ext,ij_omp_end_ext 
    103123                pk(ij,llm) = ptop + (.5*g)*rhodz(ij,llm)*(1.+theta(ij,l,2)) 
    104124             END DO 
     125             !$acc end parallel loop 
     126 
    105127             ! other layers 
    106128             DO l = llm-1, 1, -1 
     129                !$acc parallel loop 
    107130                !DIR$ SIMD 
    108131                DO ij=ij_omp_begin_ext,ij_omp_end_ext          
    … …  
    111134                        rhodz(ij,l+1)*(1.+theta(ij,l+1,2)) ) 
    112135                END DO 
     136                !$acc end parallel loop 
    113137             END DO 
    114138             ! surface pressure (for diagnostics) 
    115139             IF(caldyn_eta==eta_lag) THEN 
     140                !$acc parallel loop 
    116141                DO ij=ij_omp_begin_ext,ij_omp_end_ext          
    117142                   ps(ij) = pk(ij,1) + (.5*g)*rhodz(ij,1)*(1.+theta(ij,l,2)) 
    118143                END DO 
     144                !$acc end parallel loop 
    119145             END IF 
    120146          END SELECT 
    121147 
    122        DO l = 1,llm 
    123           SELECT CASE(caldyn_thermo) 
     148        SELECT CASE(caldyn_thermo) 
    124149          CASE(thermo_theta) 
    125              !DIR$ SIMD 
    126              DO ij=ij_omp_begin_ext,ij_omp_end_ext 
    127                 p_ik = pk(ij,l) 
    128                 exner_ik = cpp * (p_ik/preff) ** kappa 
    129                 pk(ij,l) = exner_ik 
    130                 ! specific volume v = kappa*theta*pi/p = dphi/g/rhodz 
    131                 geopot(ij,l+1) = geopot(ij,l) + (g*kappa)*rhodz(ij,l)*theta(ij,l,1)*exner_ik/p_ik 
    132              ENDDO 
     150            ! We use kernels here instead of "loop seq" + "loop gang vector" 
     151            ! to be sure the code really abides the standard. In practice, 
     152            ! it seems like the compiler interchanges the loops. 
     153            !$acc kernels async 
     154            DO l = 1,llm 
     155               !DIR$ SIMD 
     156               DO ij=ij_omp_begin_ext,ij_omp_end_ext 
     157                  p_ik = pk(ij,l) 
     158                  exner_ik = cpp * (p_ik/preff) ** kappa 
     159                  pk(ij,l) = exner_ik 
     160                  ! specific volume v = kappa*theta*pi/p = dphi/g/rhodz 
     161                  geopot(ij,l+1) = geopot(ij,l) + (g*kappa)*rhodz(ij,l)*theta(ij,l,1)*exner_ik/p_ik 
     162               ENDDO 
     163            ENDDO 
     164            !$acc end kernels 
     165 
    133166          CASE(thermo_entropy) ! theta is in fact entropy = cpp*log(theta/Treff) = cpp*log(T/Treff) - Rd*log(p/preff) 
    134              !DIR$ SIMD 
    135              DO ij=ij_omp_begin_ext,ij_omp_end_ext 
    136                 p_ik = pk(ij,l) 
    137                 temp_ik = Treff*exp((theta(ij,l,1) + Rd*log(p_ik/preff))/cpp) 
    138                 pk(ij,l) = temp_ik 
    139                 ! specific volume v = Rd*T/p = dphi/g/rhodz 
    140                 geopot(ij,l+1) = geopot(ij,l) + (g*Rd)*rhodz(ij,l)*temp_ik/p_ik 
     167             CALL abort_acc("HEVI_scheme/thermo_entropy") 
     168             DO l = 1,llm 
     169               !$acc parallel loop 
     170               !DIR$ SIMD 
     171               DO ij=ij_omp_begin_ext,ij_omp_end_ext 
     172                  p_ik = pk(ij,l) 
     173                  temp_ik = Treff*exp((theta(ij,l,1) + Rd*log(p_ik/preff))/cpp) 
     174                  pk(ij,l) = temp_ik 
     175                  ! specific volume v = Rd*T/p = dphi/g/rhodz 
     176                  geopot(ij,l+1) = geopot(ij,l) + (g*Rd)*rhodz(ij,l)*temp_ik/p_ik 
     177               ENDDO 
     178               !$acc end parallel loop 
    141179             ENDDO 
    142180          CASE(thermo_moist) ! theta is moist pseudo-entropy per dry air mass 
    143              DO ij=ij_omp_begin_ext,ij_omp_end_ext 
    144                 p_ik = pk(ij,l) 
    145                 qv = theta(ij,l,2) ! water vaper mixing ratio = mv/md 
    146                 Rmix = Rd+qv*Rv 
    147                 chi = ( theta(ij,l,1) + Rmix*log(p_ik/preff) ) / (cpp + qv*cppv) ! log(T/Treff) 
    148                 temp_ik = Treff*exp(chi) 
    149                 pk(ij,l) = temp_ik 
    150                 ! specific volume v = R*T/p = dphi/g/rhodz 
    151                 ! R = (Rd + qv.Rv)/(1+qv) 
    152                 geopot(ij,l+1) = geopot(ij,l) + g*Rmix*rhodz(ij,l)*temp_ik/(p_ik*(1+qv)) 
     181             CALL abort_acc("HEVI_scheme/thermo_moist") 
     182             DO l = 1,llm 
     183               !$acc parallel loop 
     184               DO ij=ij_omp_begin_ext,ij_omp_end_ext 
     185                  p_ik = pk(ij,l) 
     186                  qv = theta(ij,l,2) ! water vaper mixing ratio = mv/md 
     187                  Rmix = Rd+qv*Rv 
     188                  chi = ( theta(ij,l,1) + Rmix*log(p_ik/preff) ) / (cpp + qv*cppv) ! log(T/Treff) 
     189                  temp_ik = Treff*exp(chi) 
     190                  pk(ij,l) = temp_ik 
     191                  ! specific volume v = R*T/p = dphi/g/rhodz 
     192                  ! R = (Rd + qv.Rv)/(1+qv) 
     193                  geopot(ij,l+1) = geopot(ij,l) + g*Rmix*rhodz(ij,l)*temp_ik/(p_ik*(1+qv)) 
     194               ENDDO 
     195               !$acc end parallel loop 
    153196             ENDDO 
    154197          CASE DEFAULT 
    155198             STOP 
    156           END SELECT 
    157        ENDDO 
     199        END SELECT 
    158200    END IF 
    159  
     201    !$acc end data 
    160202    END IF ! dysl 
    161203 
    … …  
    167209  END SUBROUTINE compute_geopot 
    168210 
    169   SUBROUTINE compute_caldyn_vert(u,theta,rhodz,convm, wflux,wwuu, dps,dtheta_rhodz,du) 
     211  SUBROUTINE compute_caldyn_vert(u, theta, rhodz, convm, wflux, wwuu, dps, dtheta_rhodz, du) 
     212    USE disvert_mod, ONLY : bp 
    170213    REAL(rstd),INTENT(IN)  :: u(iim*3*jjm,llm) 
    171214    REAL(rstd),INTENT(IN)  :: theta(iim*jjm,llm,nqdyn) 
    … …  
    193236    CALL trace_start("compute_caldyn_vert") 
    194237 
     238    !$acc data async & 
     239    !$acc present(rhodz(:,:), u(:,:), wwuu(:,:), wflux(:,:), dps(:), convm(:,:), du(:,:), dtheta_rhodz(:,:,:), theta(:,:,:), bp(:)) 
     240 
     241 
    195242    !$OMP BARRIER    
    196243!!! cumulate mass flux convergence from top to bottom 
    197244    !  IF (is_omp_level_master) THEN 
     245    ! We use kernels here instead of "loop seq" + "loop gang vector" 
     246    ! to be sure the code really abides the standard. In practice, 
     247    ! it seems like the compiler interchanges the loops. 
     248    !$acc kernels async 
    198249    DO  l = llm-1, 1, -1 
    199250       !    IF (caldyn_conserv==energy) CALL test_message(req_qu)  
    … …  
    205256       ENDDO 
    206257    ENDDO 
     258    !$acc end kernels 
    207259    !  ENDIF 
    208260 
    … …  
    213265    ! compute dps 
    214266    IF (is_omp_first_level) THEN 
     267       !$acc parallel loop  async 
    215268       !DIR$ SIMD 
    216269       DO ij=ij_begin,ij_end          
    … …  
    221274 
    222275!!! Compute vertical mass flux (l=1,llm+1 done by caldyn_BC) 
     276    !$acc parallel loop collapse(2) async 
    223277    DO l=ll_beginp1,ll_end 
    224278       !    IF (caldyn_conserv==energy) CALL test_message(req_qu)  
    … …  
    233287    !--> flush wflux   
    234288    !$OMP BARRIER 
    235  
     289    !$acc parallel loop collapse(3) async 
    236290    DO iq=1,nqdyn 
    237291       DO l=ll_begin,ll_endm1 
    … …  
    242296          END DO 
    243297       END DO 
     298    END DO 
     299     
     300    !$acc parallel loop collapse(3) async 
     301    DO iq=1,nqdyn 
    244302       DO l=ll_beginp1,ll_end 
    245303          !DIR$ SIMD 
    … …  
    252310 
    253311    ! Compute vertical transport 
     312    !$acc parallel loop collapse(2) async 
    254313    DO l=ll_beginp1,ll_end 
    255314       !DIR$ SIMD 
    … …  
    265324 
    266325    ! Add vertical transport to du 
     326    !$acc parallel loop collapse(2) async 
    267327    DO l=ll_begin,ll_end 
    268328       !DIR$ SIMD 
    … …  
    283343    !   ENDDO 
    284344 
     345    !$acc end data 
    285346    CALL trace_end("compute_caldyn_vert") 
    286347 

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

    r899 r953  
    66  USE transfert_mod 
    77  USE caldyn_kernels_base_mod 
     8  USE abort_mod 
    89  IMPLICIT NONE 
    910  PRIVATE 
    … …  
    2021 
    2122  SUBROUTINE compute_theta(ps,theta_rhodz, rhodz,theta) 
     23    USE disvert_mod, ONLY : mass_dbk, mass_dak 
    2224    REAL(rstd),INTENT(IN)    :: ps(iim*jjm) 
    2325    REAL(rstd),INTENT(IN)    :: theta_rhodz(iim*jjm,llm,nqdyn) 
    … …  
    2628    INTEGER :: ij,l,iq 
    2729    REAL(rstd) :: m 
     30 
    2831    CALL trace_start("compute_theta") 
    2932 
     33    !$acc data async & 
     34    !$acc present(ps(:), theta_rhodz(:,:,:), rhodz(:,:), theta(:,:,:), mass_dak(:), mass_dbk(:)) 
     35 
    3036    IF(caldyn_eta==eta_mass) THEN ! Compute mass 
     37       !$acc parallel loop collapse(2) async 
    3138       DO l = ll_begin,ll_end 
    3239          !DIR$ SIMD 
    … …  
    3845    END IF 
    3946 
    40     DO l = ll_begin,ll_end 
    41        DO iq=1,nqdyn 
     47    !$acc parallel loop collapse(3) async 
     48    DO iq=1,nqdyn 
     49       DO l = ll_begin,ll_end 
    4250          !DIR$ SIMD 
    4351          DO ij=ij_begin_ext,ij_end_ext 
    … …  
    4755    END DO 
    4856 
     57    !$acc end data 
    4958    CALL trace_end("compute_theta") 
    5059  END SUBROUTINE compute_theta 
    5160 
    5261  SUBROUTINE compute_pvort_only(u,rhodz,qu,qv) 
     62    USE geometry, ONLY : Av, Riv2, fv 
    5363    REAL(rstd),INTENT(IN)  :: u(iim*3*jjm,llm) 
    5464    REAL(rstd),INTENT(INOUT) :: rhodz(iim*jjm,llm) 
    … …  
    6171    CALL trace_start("compute_pvort_only")   
    6272!!! Compute shallow-water potential vorticity 
     73 
    6374    IF(dysl_pvort_only) THEN 
     75      CALL abort_acc("HEVI_scheme/dysl_pvort_only") 
    6476#include "../kernels/pvort_only.k90" 
    6577    ELSE 
    6678 
     79    !$acc data async & 
     80    !$acc present(rhodz(:,:), u(:,:), qu(:,:), qv(:,:), Av(:), Riv2(:,:), fv(:)) 
     81     
    6782    radius_m2=radius**(-2) 
     83    !$acc parallel loop collapse(2) async 
    6884    DO l = ll_begin,ll_end 
    6985       !DIR$ SIMD 
    … …  
    85101          qv(ij+z_down,l) =( etav+fv(ij+z_down) )/hv 
    86102       ENDDO 
    87  
     103    END DO 
     104 
     105    !$acc parallel loop collapse(2) async 
     106    DO l = ll_begin,ll_end 
    88107       !DIR$ SIMD 
    89108       DO ij=ij_begin,ij_end 
    … …  
    94113 
    95114    ENDDO 
     115 
     116    !$acc end data 
    96117    
    97118    END IF ! dysl 
     119 
    98120    CALL trace_end("compute_pvort_only") 
    99121 
    … …  
    409431 
    410432    IF(dysl_caldyn_fast) THEN 
     433       CALL abort_acc("HEVI_scheme/dysl_caldyn_fast") 
    411434#include "../kernels/caldyn_fast.k90" 
    412435    ELSE 
    413436 
    414     ! Compute Bernoulli term 
    415     IF(boussinesq) THEN 
    416        DO l=ll_begin,ll_end 
    417           !DIR$ SIMD 
    418           DO ij=ij_begin,ij_end          
    419              berni(ij,l) = pk(ij,l) 
    420              ! from now on pk contains the vertically-averaged geopotential 
    421              pk(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1)) 
     437 
     438       ! Default case : ported to openacc 
     439       IF(.not. boussinesq .and. caldyn_thermo /= thermo_moist) THEN 
     440          !$acc data async & 
     441          !$acc present(rhodz(:,:), u(:,:),pk(:,:), geopot(:,:), theta(:,:,:), du(:,:)) 
     442 
     443#define BERNI(ij) .5*(geopot(ij,l)+geopot(ij,l+1)) 
     444          !$acc parallel loop collapse(2) async 
     445          DO l=ll_begin,ll_end 
     446             !DIR$ SIMD 
     447             DO ij=ij_begin,ij_end 
     448                due_right =  0.5*(theta(ij,l,1)+theta(ij+t_right,l,1))  & 
     449                     *(pk(ij+t_right,l)-pk(ij,l))        & 
     450                     +  BERNI(ij+t_right)-BERNI(ij) 
     451                due_lup = 0.5*(theta(ij,l,1)+theta(ij+t_lup,l,1))    & 
     452                     *(pk(ij+t_lup,l)-pk(ij,l))          & 
     453                     +  BERNI(ij+t_lup)-BERNI(ij) 
     454                due_ldown = 0.5*(theta(ij,l,1)+theta(ij+t_ldown,l,1)) & 
     455                     *(pk(ij+t_ldown,l)-pk(ij,l))      & 
     456                     +   BERNI(ij+t_ldown)-BERNI(ij) 
     457                du(ij+u_right,l) = du(ij+u_right,l) - ne_right*due_right 
     458                du(ij+u_lup,l)   = du(ij+u_lup,l)   - ne_lup*due_lup 
     459                du(ij+u_ldown,l) = du(ij+u_ldown,l) - ne_ldown*due_ldown 
     460                u(ij+u_right,l) = u(ij+u_right,l) + tau*du(ij+u_right,l) 
     461                u(ij+u_lup,l)   = u(ij+u_lup,l)   + tau*du(ij+u_lup,l) 
     462                u(ij+u_ldown,l) = u(ij+u_ldown,l) + tau*du(ij+u_ldown,l) 
     463             END DO 
    422464          END DO 
    423        END DO 
    424     ELSE ! compressible 
    425  
    426        DO l=ll_begin,ll_end 
    427           SELECT CASE(caldyn_thermo) 
    428           CASE(thermo_theta) ! vdp = theta.dpi => B = Phi 
    429              !DIR$ SIMD 
    430              DO ij=ij_begin,ij_end          
    431                 berni(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1)) 
    432              END DO 
    433           CASE(thermo_entropy) ! vdp = dG + sdT => B = Phi + G, G=h-Ts=T*(cpp-s)  
    434              !DIR$ SIMD 
    435              DO ij=ij_begin,ij_end          
    436                 berni(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1)) & 
    437                      + pk(ij,l)*(cpp-theta(ij,l,1)) ! pk=temperature, theta=entropy 
    438              END DO 
    439           CASE(thermo_moist)  
    440              !DIR$ SIMD 
    441              DO ij=ij_begin,ij_end          
    442                 ! du/dt = grad(Bd)+rv.grad(Bv)+s.grad(T) 
    443                 ! Bd = Phi + gibbs_d 
    444                 ! Bv = Phi + gibbs_v 
    445                 ! pk=temperature, theta=entropy 
    446                 qv = theta(ij,l,2) 
    447                 temp = pk(ij,l) 
    448                 chi = log(temp/Treff) 
    449                 nu = (chi*(cpp+qv*cppv)-theta(ij,l,1))/(Rd+qv*Rv) ! log(p/preff) 
    450                 berni(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1)) & 
    451                      + temp*(cpp*(1.-chi)+Rd*nu) 
    452                 berniv(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1)) & 
    453                      + temp*(cppv*(1.-chi)+Rv*nu) 
    454             END DO 
    455           END SELECT 
    456        END DO 
    457  
    458     END IF ! Boussinesq/compressible 
     465 
     466#undef BERNI 
     467 
     468          !$acc end data 
     469       ELSE ! Not default case : not ported to openacc 
     470          CALL abort_acc("HEVI_scheme/compute_caldyn_fast") 
     471          !!$acc parallel loop gang private(berni(:,:), berniv(:,:)) async 
     472          DO l=ll_begin,ll_end 
     473             ! Compute Bernoulliterm 
     474             IF(boussinesq) THEN 
     475                !!$acc loop vector 
     476                !DIR$ SIMD 
     477                DO ij=ij_begin,ij_end 
     478                   berni(ij,l) = pk(ij,l) 
     479                   ! from now on pk contains the vertically-averaged geopotential 
     480                   pk(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1)) 
     481                END DO 
     482             ELSE ! compressible 
     483                SELECT CASE(caldyn_thermo) 
     484                CASE(thermo_theta) ! vdp = theta.dpi => B = Phi 
     485                   CALL dynamico_abort("Case already dealt with") 
     486                CASE(thermo_entropy) ! vdp = dG + sdT => B = Phi + G, G=h-Ts=T*(cpp-s) 
     487                   !$acc loop vector 
     488                   !DIR$ SIMD 
     489                   DO ij=ij_begin,ij_end 
     490                      berni(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1)) & 
     491                           + pk(ij,l)*(cpp-theta(ij,l,1)) ! pk=temperature, theta=entropy 
     492                   END DO 
     493                CASE(thermo_moist) 
     494!!$acc loop vector 
     495                   !DIR$ SIMD 
     496                   DO ij=ij_begin,ij_end 
     497                      ! du/dt = grad(Bd)+rv.grad(Bv)+s.grad(T) 
     498                      ! Bd = Phi + gibbs_d 
     499                      ! Bv = Phi + gibbs_v 
     500                      ! pk=temperature, theta=entropy 
     501                      qv = theta(ij,l,2) 
     502                      temp = pk(ij,l) 
     503                      chi = log(temp/Treff) 
     504                      nu = (chi*(cpp+qv*cppv)-theta(ij,l,1))/(Rd+qv*Rv) ! log(p/preff) 
     505                      berni(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1)) & 
     506                           + temp*(cpp*(1.-chi)+Rd*nu) 
     507                      berniv(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1)) & 
     508                           + temp*(cppv*(1.-chi)+Rv*nu) 
     509                   END DO 
     510                END SELECT 
     511             END IF ! Boussinesq/compressible 
    459512 
    460513!!! u:=u+tau*du, du = -grad(B)-theta.grad(pi) 
    461     DO l=ll_begin,ll_end 
    462        IF(caldyn_thermo == thermo_moist) THEN 
    463           !DIR$ SIMD 
    464           DO ij=ij_begin,ij_end          
    465              due_right =  berni(ij+t_right,l)-berni(ij,l)       & 
    466                   + 0.5*(theta(ij,l,1)+theta(ij+t_right,l,1))   & 
    467                        *(pk(ij+t_right,l)-pk(ij,l))             & 
    468                   + 0.5*(theta(ij,l,2)+theta(ij+t_right,l,2))   & 
    469                        *(berniv(ij+t_right,l)-berniv(ij,l)) 
    470               
    471              due_lup = berni(ij+t_lup,l)-berni(ij,l)            & 
    472                   + 0.5*(theta(ij,l,1)+theta(ij+t_lup,l,1))     & 
    473                        *(pk(ij+t_lup,l)-pk(ij,l))               & 
    474                   + 0.5*(theta(ij,l,2)+theta(ij+t_lup,l,2))     & 
    475                        *(berniv(ij+t_lup,l)-berniv(ij,l)) 
    476               
    477              due_ldown = berni(ij+t_ldown,l)-berni(ij,l)        & 
    478                   + 0.5*(theta(ij,l,1)+theta(ij+t_ldown,l,1))   & 
    479                        *(pk(ij+t_ldown,l)-pk(ij,l))             & 
    480                   + 0.5*(theta(ij,l,2)+theta(ij+t_ldown,l,2))   & 
    481                        *(berniv(ij+t_ldown,l)-berniv(ij,l)) 
    482               
    483              du(ij+u_right,l) = du(ij+u_right,l) - ne_right*due_right 
    484              du(ij+u_lup,l)   = du(ij+u_lup,l)   - ne_lup*due_lup 
    485              du(ij+u_ldown,l) = du(ij+u_ldown,l) - ne_ldown*due_ldown 
    486              u(ij+u_right,l) = u(ij+u_right,l) + tau*du(ij+u_right,l) 
    487              u(ij+u_lup,l)   = u(ij+u_lup,l)   + tau*du(ij+u_lup,l) 
    488              u(ij+u_ldown,l) = u(ij+u_ldown,l) + tau*du(ij+u_ldown,l) 
    489           END DO 
    490        ELSE 
    491           !DIR$ SIMD 
    492           DO ij=ij_begin,ij_end          
    493              due_right =  0.5*(theta(ij,l,1)+theta(ij+t_right,l,1))  & 
    494                   *(pk(ij+t_right,l)-pk(ij,l))        & 
    495                   +  berni(ij+t_right,l)-berni(ij,l) 
    496              due_lup = 0.5*(theta(ij,l,1)+theta(ij+t_lup,l,1))    & 
    497                   *(pk(ij+t_lup,l)-pk(ij,l))          & 
    498                   +  berni(ij+t_lup,l)-berni(ij,l) 
    499              due_ldown = 0.5*(theta(ij,l,1)+theta(ij+t_ldown,l,1)) & 
    500                   *(pk(ij+t_ldown,l)-pk(ij,l))      & 
    501                   +   berni(ij+t_ldown,l)-berni(ij,l) 
    502              du(ij+u_right,l) = du(ij+u_right,l) - ne_right*due_right 
    503              du(ij+u_lup,l)   = du(ij+u_lup,l)   - ne_lup*due_lup 
    504              du(ij+u_ldown,l) = du(ij+u_ldown,l) - ne_ldown*due_ldown 
    505              u(ij+u_right,l) = u(ij+u_right,l) + tau*du(ij+u_right,l) 
    506              u(ij+u_lup,l)   = u(ij+u_lup,l)   + tau*du(ij+u_lup,l) 
    507              u(ij+u_ldown,l) = u(ij+u_ldown,l) + tau*du(ij+u_ldown,l) 
     514             IF(caldyn_thermo == thermo_moist) THEN 
     515!!$acc loop vector 
     516                !DIR$ SIMD 
     517                DO ij=ij_begin,ij_end 
     518                   due_right =  berni(ij+t_right,l)-berni(ij,l)       & 
     519                        + 0.5*(theta(ij,l,1)+theta(ij+t_right,l,1))   & 
     520                        *(pk(ij+t_right,l)-pk(ij,l))             & 
     521                        + 0.5*(theta(ij,l,2)+theta(ij+t_right,l,2))   & 
     522                        *(berniv(ij+t_right,l)-berniv(ij,l)) 
     523 
     524                   due_lup = berni(ij+t_lup,l)-berni(ij,l)            & 
     525                        + 0.5*(theta(ij,l,1)+theta(ij+t_lup,l,1))     & 
     526                        *(pk(ij+t_lup,l)-pk(ij,l))               & 
     527                        + 0.5*(theta(ij,l,2)+theta(ij+t_lup,l,2))     & 
     528                        *(berniv(ij+t_lup,l)-berniv(ij,l)) 
     529 
     530                   due_ldown = berni(ij+t_ldown,l)-berni(ij,l)        & 
     531                        + 0.5*(theta(ij,l,1)+theta(ij+t_ldown,l,1))   & 
     532                        *(pk(ij+t_ldown,l)-pk(ij,l))             & 
     533                        + 0.5*(theta(ij,l,2)+theta(ij+t_ldown,l,2))   & 
     534                        *(berniv(ij+t_ldown,l)-berniv(ij,l)) 
     535 
     536                   du(ij+u_right,l) = du(ij+u_right,l) - ne_right*due_right 
     537                   du(ij+u_lup,l)   = du(ij+u_lup,l)   - ne_lup*due_lup 
     538                   du(ij+u_ldown,l) = du(ij+u_ldown,l) - ne_ldown*due_ldown 
     539                   u(ij+u_right,l) = u(ij+u_right,l) + tau*du(ij+u_right,l) 
     540                   u(ij+u_lup,l)   = u(ij+u_lup,l)   + tau*du(ij+u_lup,l) 
     541                   u(ij+u_ldown,l) = u(ij+u_ldown,l) + tau*du(ij+u_ldown,l) 
     542                END DO 
     543             ELSE 
     544                CALL dynamico_abort("Case already dealt with") 
     545             END IF 
    508546          END DO 
    509547       END IF 
    510     END DO 
    511  
    512     END IF ! dysl 
     548    END IF 
    513549    CALL trace_end("compute_caldyn_fast") 
    514550 
    515551  END SUBROUTINE compute_caldyn_fast 
    516  
     552     
    517553  SUBROUTINE compute_caldyn_Coriolis(hflux,theta,qu, convm,dtheta_rhodz,du) 
     554    USE geometry, ONLY : Ai, wee 
    518555    REAL(rstd),INTENT(IN)  :: hflux(3*iim*jjm,llm)  ! hflux in kg/s 
    519556    REAL(rstd),INTENT(IN)  :: theta(iim*jjm,llm,nqdyn) ! active scalars 
    … …  
    534571 
    535572    ELSE 
    536 #define FTHETA(ij) Ftheta(ij,1) 
    537  
    538     DO l=ll_begin, ll_end 
    539        ! compute theta flux 
    540        DO iq=1,nqdyn 
    541        !DIR$ SIMD 
     573#define FTHETA(ij) Ftheta(ij,l) 
     574 
     575    !$acc data async & 
     576    !$acc create(Ftheta(3*iim*jjm,llm)) & 
     577    !$acc present(qu(:,:), hflux(:,:), convm(:,:), du(:,:), dtheta_rhodz(:,:,:), theta(:,:,:), Ai(:), wee(:,:,:)) 
     578 
     579    ! compute theta flux 
     580    DO iq=1,nqdyn 
     581       !$acc parallel loop collapse(2) present(Ftheta, theta) async 
     582       DO l=ll_begin, ll_end 
     583          !DIR$ SIMD 
    542584          DO ij=ij_begin_ext,ij_end_ext       
    543585             FTHETA(ij+u_right) = 0.5*(theta(ij,l,iq)+theta(ij+t_right,l,iq)) & 
    … …  
    548590                  * hflux(ij+u_ldown,l) 
    549591          END DO 
     592       END DO 
     593 
     594       !$acc parallel loop collapse(2) present(Ftheta) async 
     595       DO l=ll_begin, ll_end 
    550596          ! horizontal divergence of fluxes 
    551        !DIR$ SIMD 
     597          !DIR$ SIMD 
    552598          DO ij=ij_begin,ij_end          
    553599             ! dtheta_rhodz =  -div(flux.theta) 
    … …  
    561607          END DO 
    562608       END DO 
    563  
     609    END DO 
     610 
     611    !$acc parallel loop collapse(2) present(Ai) async 
     612    DO l=ll_begin, ll_end 
    564613       !DIR$ SIMD 
    565614       DO ij=ij_begin,ij_end          
    … …  
    578627 
    579628    CASE(energy) ! energy-conserving TRiSK 
    580  
     629       !$acc parallel loop collapse(2) present(qu, wee, du) async 
    581630       DO l=ll_begin,ll_end 
    582631          !DIR$ SIMD 
    … …  
    622671 
    623672    CASE(enstrophy) ! enstrophy-conserving TRiSK 
    624  
     673        
     674       !$acc parallel loop collapse(2) present(wee, du) async 
    625675       DO l=ll_begin,ll_end 
    626676          !DIR$ SIMD 
    … …  
    665715          END DO 
    666716       END DO 
     717       !$acc end parallel loop 
    667718    CASE DEFAULT 
    668719       STOP 
    669720    END SELECT 
     721 
     722    !$acc end data 
     723 
    670724#undef FTHETA 
    671725 
    … …  
    676730  END SUBROUTINE compute_caldyn_Coriolis 
    677731 
    678   SUBROUTINE compute_caldyn_slow_hydro(u,rhodz,hflux,du, zero) 
     732  SUBROUTINE compute_caldyn_slow_hydro(u,rhodz,hflux,du,zero) 
     733    USE geometry, ONLY : Ai, le_de 
    679734    LOGICAL, INTENT(IN) :: zero 
    680735    REAL(rstd),INTENT(IN)  :: u(3*iim*jjm,llm)    ! prognostic "velocity" 
    … …  
    684739     
    685740    REAL(rstd) :: berni(iim*jjm,llm)  ! Bernoulli function 
    686     REAL(rstd) :: berni1(iim*jjm)  ! Bernoulli function 
    687741    REAL(rstd) :: uu_right, uu_lup, uu_ldown, ke, uu 
    688742    INTEGER :: ij,l 
    … …  
    698752     ELSE 
    699753 
    700 #define BERNI(ij) berni1(ij) 
     754#define BERNI(ij) berni(ij,l) 
    701755 
    702756    DO l = ll_begin, ll_end 
    703        !  Compute mass fluxes 
    704757       IF (caldyn_conserv==energy) CALL test_message(req_qu)  
     758    END DO 
     759     
     760    !$acc data async & 
     761    !$acc create(berni(:,ll_begin:ll_end)) & 
     762    !$acc present(rhodz(:,ll_begin:ll_end), u(:,ll_begin:ll_end), hflux(:,:), du(:,ll_begin:ll_end), Ai(:), le_de(:)) 
     763 
     764    !  Compute mass fluxes 
     765    !$acc parallel loop collapse(2) async 
     766    DO l = ll_begin, ll_end 
    705767       !DIR$ SIMD 
    706768       DO ij=ij_begin_ext,ij_end_ext 
    … …  
    715777          hflux(ij+u_ldown,l)=uu_ldown 
    716778       ENDDO 
    717        ! Compute Bernoulli=kinetic energy 
     779    END DO 
     780 
     781    ! Compute Bernoulli=kinetic energy 
     782    !$acc parallel loop collapse(2) async 
     783    DO l = ll_begin, ll_end 
    718784       !DIR$ SIMD 
    719785       DO ij=ij_begin,ij_end          
    … …  
    726792               le_de(ij+u_rdown)*u(ij+u_rdown,l)**2 )   
    727793       ENDDO 
    728        ! Compute du=-grad(Bernoulli) 
    729        IF(zero) THEN 
     794    END DO 
     795 
     796    ! Compute du=-grad(Bernoulli) 
     797    IF(zero) THEN 
     798       !$acc parallel loop collapse(2) async 
     799       DO l = ll_begin, ll_end 
    730800          !DIR$ SIMD 
    731801          DO ij=ij_begin,ij_end 
    … …  
    734804             du(ij+u_ldown,l) = ne_ldown*(BERNI(ij)-BERNI(ij+t_ldown)) 
    735805          END DO 
    736        ELSE 
     806       END DO 
     807    ELSE 
     808       !$acc parallel loop collapse(2) async 
     809       DO l = ll_begin, ll_end 
    737810          !DIR$ SIMD 
    738811          DO ij=ij_begin,ij_end 
    … …  
    744817                  ne_ldown*(BERNI(ij)-BERNI(ij+t_ldown)) 
    745818          END DO 
    746        END IF 
    747     END DO 
     819       END DO 
     820    END IF 
     821 
     822    !$acc end data 
    748823 
    749824#undef BERNI 
    750825 
    751826    END IF ! dysl 
     827 
    752828    CALL trace_end("compute_caldyn_slow_hydro")     
    753829  END SUBROUTINE compute_caldyn_slow_hydro 
Note: See TracChangeset for help on using the changeset viewer.