Changeset 940 for codes


Ignore:
Timestamp:
07/05/19 15:13:09 (5 years ago)
Author:
dubos
Message:

devel : DySL for enstrophy-conserving scheme

Location:
codes/icosagcm/devel
Files:
1 deleted
7 edited

Legend:

Unmodified
Added
Removed
  • codes/icosagcm/devel/Python/src/unstructured/macros.jin

    r919 r940  
    396396#define AV Av(ij) 
    397397#define FV fv(ij) 
    398 #define WEE wee(itrisk,edge) 
     398#define WEE wee(itrisk,edge,1) 
    399399#define edge_ne(iedge,ij) 1. 
    400400#endif 
  • codes/icosagcm/devel/src/dynamics/compute_caldyn_Coriolis.F90

    r939 r940  
    11MODULE compute_caldyn_Coriolis_mod 
    22  USE prec, ONLY : rstd 
     3  USE caldyn_vars_mod, ONLY : caldyn_conserv, conserv_energy, conserv_enstrophy, conserv_gassmann 
    34  USE grid_param 
    45  USE earth_const 
     
    4647    END_BLOCK 
    4748  END_BLOCK 
    48 ! 
    49   FORALL_CELLS() 
    50     ON_EDGES 
    51       du_trisk=0. 
    52       FORALL_TRISK 
    53         du_trisk = du_trisk + WEE*hflux(EDGE_TRISK)*(qu(EDGE)+qu(EDGE_TRISK)) 
    54       END_BLOCK 
    55       du(EDGE) = du(EDGE) + .5*du_trisk 
    56     END_BLOCK 
    57   END_BLOCK 
    58  
     49  ! 
     50  SELECT CASE(caldyn_conserv) 
     51  CASE(conserv_energy) ! energy-conserving TRiSK 
     52     FORALL_CELLS() 
     53       ON_EDGES 
     54         du_trisk=0. 
     55         FORALL_TRISK 
     56           du_trisk = du_trisk + WEE*hflux(EDGE_TRISK)*(qu(EDGE)+qu(EDGE_TRISK)) 
     57         END_BLOCK 
     58         du(EDGE) = du(EDGE) + .5*du_trisk 
     59       END_BLOCK 
     60     END_BLOCK 
     61  CASE(conserv_enstrophy) ! enstrophy-conserving TRiSK 
     62     FORALL_CELLS() 
     63       ON_EDGES 
     64         du_trisk=0. 
     65         FORALL_TRISK 
     66           du_trisk = du_trisk + WEE*hflux(EDGE_TRISK) 
     67         END_BLOCK 
     68         du(EDGE) = du(EDGE) + du_trisk*qu(EDGE) 
     69       END_BLOCK 
     70     END_BLOCK 
     71  END SELECT 
    5972END_BLOCK 
    6073 
     
    6578    USE data_unstructured_mod, ONLY : enter_trace, exit_trace, & 
    6679         id_coriolis, left, right, primal_deg, primal_edge, primal_ne, & 
    67          trisk_deg, trisk, wee ! FIXME wee 
     80         trisk_deg, trisk 
    6881          
    6982    FIELD_U     :: hflux, Ftheta, qu, du ! IN, BUF, IN, INOUT 
     
    8093  SUBROUTINE compute_caldyn_Coriolis_hex(hflux,theta,qu, Ftheta, convm,dtheta_rhodz,du) 
    8194    USE icosa 
    82     USE caldyn_vars_mod 
    8395    REAL(rstd),INTENT(IN)    :: hflux(3*iim*jjm,llm)  ! hflux in kg/s 
    8496    REAL(rstd),INTENT(IN)    :: theta(iim*jjm,llm,nqdyn) ! active scalars 
  • codes/icosagcm/devel/src/kernels_hex/coriolis.k90

    r563 r940  
    4040   END DO 
    4141   ! 
    42    DO l = ll_begin, ll_end 
    43       !DIR$ SIMD 
    44       DO ij=ij_begin, ij_end 
    45          du_trisk=0. 
    46          du_trisk = du_trisk + wee(ij+u_right,1,1)*hflux(ij+u_rup,l)*(qu(ij+u_right,l)+qu(ij+u_rup,l)) 
    47          du_trisk = du_trisk + wee(ij+u_right,2,1)*hflux(ij+u_lup,l)*(qu(ij+u_right,l)+qu(ij+u_lup,l)) 
    48          du_trisk = du_trisk + wee(ij+u_right,3,1)*hflux(ij+u_left,l)*(qu(ij+u_right,l)+qu(ij+u_left,l)) 
    49          du_trisk = du_trisk + wee(ij+u_right,4,1)*hflux(ij+u_ldown,l)*(qu(ij+u_right,l)+qu(ij+u_ldown,l)) 
    50          du_trisk = du_trisk + wee(ij+u_right,5,1)*hflux(ij+u_rdown,l)*(qu(ij+u_right,l)+qu(ij+u_rdown,l)) 
    51          du_trisk = du_trisk + wee(ij+u_right,1,2)*hflux(ij+t_right+u_ldown,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_ldown,l)) 
    52          du_trisk = du_trisk + wee(ij+u_right,2,2)*hflux(ij+t_right+u_rdown,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_rdown,l)) 
    53          du_trisk = du_trisk + wee(ij+u_right,3,2)*hflux(ij+t_right+u_right,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_right,l)) 
    54          du_trisk = du_trisk + wee(ij+u_right,4,2)*hflux(ij+t_right+u_rup,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_rup,l)) 
    55          du_trisk = du_trisk + wee(ij+u_right,5,2)*hflux(ij+t_right+u_lup,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_lup,l)) 
    56          du(ij+u_right,l) = du(ij+u_right,l) + .5*du_trisk 
    57          du_trisk=0. 
    58          du_trisk = du_trisk + wee(ij+u_lup,1,1)*hflux(ij+u_left,l)*(qu(ij+u_lup,l)+qu(ij+u_left,l)) 
    59          du_trisk = du_trisk + wee(ij+u_lup,2,1)*hflux(ij+u_ldown,l)*(qu(ij+u_lup,l)+qu(ij+u_ldown,l)) 
    60          du_trisk = du_trisk + wee(ij+u_lup,3,1)*hflux(ij+u_rdown,l)*(qu(ij+u_lup,l)+qu(ij+u_rdown,l)) 
    61          du_trisk = du_trisk + wee(ij+u_lup,4,1)*hflux(ij+u_right,l)*(qu(ij+u_lup,l)+qu(ij+u_right,l)) 
    62          du_trisk = du_trisk + wee(ij+u_lup,5,1)*hflux(ij+u_rup,l)*(qu(ij+u_lup,l)+qu(ij+u_rup,l)) 
    63          du_trisk = du_trisk + wee(ij+u_lup,1,2)*hflux(ij+t_lup+u_right,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_right,l)) 
    64          du_trisk = du_trisk + wee(ij+u_lup,2,2)*hflux(ij+t_lup+u_rup,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_rup,l)) 
    65          du_trisk = du_trisk + wee(ij+u_lup,3,2)*hflux(ij+t_lup+u_lup,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_lup,l)) 
    66          du_trisk = du_trisk + wee(ij+u_lup,4,2)*hflux(ij+t_lup+u_left,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_left,l)) 
    67          du_trisk = du_trisk + wee(ij+u_lup,5,2)*hflux(ij+t_lup+u_ldown,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_ldown,l)) 
    68          du(ij+u_lup,l) = du(ij+u_lup,l) + .5*du_trisk 
    69          du_trisk=0. 
    70          du_trisk = du_trisk + wee(ij+u_ldown,1,1)*hflux(ij+u_rdown,l)*(qu(ij+u_ldown,l)+qu(ij+u_rdown,l)) 
    71          du_trisk = du_trisk + wee(ij+u_ldown,2,1)*hflux(ij+u_right,l)*(qu(ij+u_ldown,l)+qu(ij+u_right,l)) 
    72          du_trisk = du_trisk + wee(ij+u_ldown,3,1)*hflux(ij+u_rup,l)*(qu(ij+u_ldown,l)+qu(ij+u_rup,l)) 
    73          du_trisk = du_trisk + wee(ij+u_ldown,4,1)*hflux(ij+u_lup,l)*(qu(ij+u_ldown,l)+qu(ij+u_lup,l)) 
    74          du_trisk = du_trisk + wee(ij+u_ldown,5,1)*hflux(ij+u_left,l)*(qu(ij+u_ldown,l)+qu(ij+u_left,l)) 
    75          du_trisk = du_trisk + wee(ij+u_ldown,1,2)*hflux(ij+t_ldown+u_lup,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_lup,l)) 
    76          du_trisk = du_trisk + wee(ij+u_ldown,2,2)*hflux(ij+t_ldown+u_left,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_left,l)) 
    77          du_trisk = du_trisk + wee(ij+u_ldown,3,2)*hflux(ij+t_ldown+u_ldown,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_ldown,l)) 
    78          du_trisk = du_trisk + wee(ij+u_ldown,4,2)*hflux(ij+t_ldown+u_rdown,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_rdown,l)) 
    79          du_trisk = du_trisk + wee(ij+u_ldown,5,2)*hflux(ij+t_ldown+u_right,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_right,l)) 
    80          du(ij+u_ldown,l) = du(ij+u_ldown,l) + .5*du_trisk 
     42   SELECT CASE(caldyn_conserv) 
     43   CASE(conserv_energy) ! energy-conserving TRiSK 
     44      DO l = ll_begin, ll_end 
     45         !DIR$ SIMD 
     46         DO ij=ij_begin, ij_end 
     47            du_trisk=0. 
     48            du_trisk = du_trisk + wee(ij+u_right,1,1)*hflux(ij+u_rup,l)*(qu(ij+u_right,l)+qu(ij+u_rup,l)) 
     49            du_trisk = du_trisk + wee(ij+u_right,2,1)*hflux(ij+u_lup,l)*(qu(ij+u_right,l)+qu(ij+u_lup,l)) 
     50            du_trisk = du_trisk + wee(ij+u_right,3,1)*hflux(ij+u_left,l)*(qu(ij+u_right,l)+qu(ij+u_left,l)) 
     51            du_trisk = du_trisk + wee(ij+u_right,4,1)*hflux(ij+u_ldown,l)*(qu(ij+u_right,l)+qu(ij+u_ldown,l)) 
     52            du_trisk = du_trisk + wee(ij+u_right,5,1)*hflux(ij+u_rdown,l)*(qu(ij+u_right,l)+qu(ij+u_rdown,l)) 
     53            du_trisk = du_trisk + wee(ij+u_right,1,2)*hflux(ij+t_right+u_ldown,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_ldown,l)) 
     54            du_trisk = du_trisk + wee(ij+u_right,2,2)*hflux(ij+t_right+u_rdown,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_rdown,l)) 
     55            du_trisk = du_trisk + wee(ij+u_right,3,2)*hflux(ij+t_right+u_right,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_right,l)) 
     56            du_trisk = du_trisk + wee(ij+u_right,4,2)*hflux(ij+t_right+u_rup,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_rup,l)) 
     57            du_trisk = du_trisk + wee(ij+u_right,5,2)*hflux(ij+t_right+u_lup,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_lup,l)) 
     58            du(ij+u_right,l) = du(ij+u_right,l) + .5*du_trisk 
     59            du_trisk=0. 
     60            du_trisk = du_trisk + wee(ij+u_lup,1,1)*hflux(ij+u_left,l)*(qu(ij+u_lup,l)+qu(ij+u_left,l)) 
     61            du_trisk = du_trisk + wee(ij+u_lup,2,1)*hflux(ij+u_ldown,l)*(qu(ij+u_lup,l)+qu(ij+u_ldown,l)) 
     62            du_trisk = du_trisk + wee(ij+u_lup,3,1)*hflux(ij+u_rdown,l)*(qu(ij+u_lup,l)+qu(ij+u_rdown,l)) 
     63            du_trisk = du_trisk + wee(ij+u_lup,4,1)*hflux(ij+u_right,l)*(qu(ij+u_lup,l)+qu(ij+u_right,l)) 
     64            du_trisk = du_trisk + wee(ij+u_lup,5,1)*hflux(ij+u_rup,l)*(qu(ij+u_lup,l)+qu(ij+u_rup,l)) 
     65            du_trisk = du_trisk + wee(ij+u_lup,1,2)*hflux(ij+t_lup+u_right,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_right,l)) 
     66            du_trisk = du_trisk + wee(ij+u_lup,2,2)*hflux(ij+t_lup+u_rup,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_rup,l)) 
     67            du_trisk = du_trisk + wee(ij+u_lup,3,2)*hflux(ij+t_lup+u_lup,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_lup,l)) 
     68            du_trisk = du_trisk + wee(ij+u_lup,4,2)*hflux(ij+t_lup+u_left,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_left,l)) 
     69            du_trisk = du_trisk + wee(ij+u_lup,5,2)*hflux(ij+t_lup+u_ldown,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_ldown,l)) 
     70            du(ij+u_lup,l) = du(ij+u_lup,l) + .5*du_trisk 
     71            du_trisk=0. 
     72            du_trisk = du_trisk + wee(ij+u_ldown,1,1)*hflux(ij+u_rdown,l)*(qu(ij+u_ldown,l)+qu(ij+u_rdown,l)) 
     73            du_trisk = du_trisk + wee(ij+u_ldown,2,1)*hflux(ij+u_right,l)*(qu(ij+u_ldown,l)+qu(ij+u_right,l)) 
     74            du_trisk = du_trisk + wee(ij+u_ldown,3,1)*hflux(ij+u_rup,l)*(qu(ij+u_ldown,l)+qu(ij+u_rup,l)) 
     75            du_trisk = du_trisk + wee(ij+u_ldown,4,1)*hflux(ij+u_lup,l)*(qu(ij+u_ldown,l)+qu(ij+u_lup,l)) 
     76            du_trisk = du_trisk + wee(ij+u_ldown,5,1)*hflux(ij+u_left,l)*(qu(ij+u_ldown,l)+qu(ij+u_left,l)) 
     77            du_trisk = du_trisk + wee(ij+u_ldown,1,2)*hflux(ij+t_ldown+u_lup,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_lup,l)) 
     78            du_trisk = du_trisk + wee(ij+u_ldown,2,2)*hflux(ij+t_ldown+u_left,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_left,l)) 
     79            du_trisk = du_trisk + wee(ij+u_ldown,3,2)*hflux(ij+t_ldown+u_ldown,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_ldown,l)) 
     80            du_trisk = du_trisk + wee(ij+u_ldown,4,2)*hflux(ij+t_ldown+u_rdown,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_rdown,l)) 
     81            du_trisk = du_trisk + wee(ij+u_ldown,5,2)*hflux(ij+t_ldown+u_right,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_right,l)) 
     82            du(ij+u_ldown,l) = du(ij+u_ldown,l) + .5*du_trisk 
     83         END DO 
    8184      END DO 
    82    END DO 
     85   CASE(conserv_enstrophy) ! enstrophy-conserving TRiSK 
     86      DO l = ll_begin, ll_end 
     87         !DIR$ SIMD 
     88         DO ij=ij_begin, ij_end 
     89            du_trisk=0. 
     90            du_trisk = du_trisk + wee(ij+u_right,1,1)*hflux(ij+u_rup,l) 
     91            du_trisk = du_trisk + wee(ij+u_right,2,1)*hflux(ij+u_lup,l) 
     92            du_trisk = du_trisk + wee(ij+u_right,3,1)*hflux(ij+u_left,l) 
     93            du_trisk = du_trisk + wee(ij+u_right,4,1)*hflux(ij+u_ldown,l) 
     94            du_trisk = du_trisk + wee(ij+u_right,5,1)*hflux(ij+u_rdown,l) 
     95            du_trisk = du_trisk + wee(ij+u_right,1,2)*hflux(ij+t_right+u_ldown,l) 
     96            du_trisk = du_trisk + wee(ij+u_right,2,2)*hflux(ij+t_right+u_rdown,l) 
     97            du_trisk = du_trisk + wee(ij+u_right,3,2)*hflux(ij+t_right+u_right,l) 
     98            du_trisk = du_trisk + wee(ij+u_right,4,2)*hflux(ij+t_right+u_rup,l) 
     99            du_trisk = du_trisk + wee(ij+u_right,5,2)*hflux(ij+t_right+u_lup,l) 
     100            du(ij+u_right,l) = du(ij+u_right,l) + du_trisk*qu(ij+u_right,l) 
     101            du_trisk=0. 
     102            du_trisk = du_trisk + wee(ij+u_lup,1,1)*hflux(ij+u_left,l) 
     103            du_trisk = du_trisk + wee(ij+u_lup,2,1)*hflux(ij+u_ldown,l) 
     104            du_trisk = du_trisk + wee(ij+u_lup,3,1)*hflux(ij+u_rdown,l) 
     105            du_trisk = du_trisk + wee(ij+u_lup,4,1)*hflux(ij+u_right,l) 
     106            du_trisk = du_trisk + wee(ij+u_lup,5,1)*hflux(ij+u_rup,l) 
     107            du_trisk = du_trisk + wee(ij+u_lup,1,2)*hflux(ij+t_lup+u_right,l) 
     108            du_trisk = du_trisk + wee(ij+u_lup,2,2)*hflux(ij+t_lup+u_rup,l) 
     109            du_trisk = du_trisk + wee(ij+u_lup,3,2)*hflux(ij+t_lup+u_lup,l) 
     110            du_trisk = du_trisk + wee(ij+u_lup,4,2)*hflux(ij+t_lup+u_left,l) 
     111            du_trisk = du_trisk + wee(ij+u_lup,5,2)*hflux(ij+t_lup+u_ldown,l) 
     112            du(ij+u_lup,l) = du(ij+u_lup,l) + du_trisk*qu(ij+u_lup,l) 
     113            du_trisk=0. 
     114            du_trisk = du_trisk + wee(ij+u_ldown,1,1)*hflux(ij+u_rdown,l) 
     115            du_trisk = du_trisk + wee(ij+u_ldown,2,1)*hflux(ij+u_right,l) 
     116            du_trisk = du_trisk + wee(ij+u_ldown,3,1)*hflux(ij+u_rup,l) 
     117            du_trisk = du_trisk + wee(ij+u_ldown,4,1)*hflux(ij+u_lup,l) 
     118            du_trisk = du_trisk + wee(ij+u_ldown,5,1)*hflux(ij+u_left,l) 
     119            du_trisk = du_trisk + wee(ij+u_ldown,1,2)*hflux(ij+t_ldown+u_lup,l) 
     120            du_trisk = du_trisk + wee(ij+u_ldown,2,2)*hflux(ij+t_ldown+u_left,l) 
     121            du_trisk = du_trisk + wee(ij+u_ldown,3,2)*hflux(ij+t_ldown+u_ldown,l) 
     122            du_trisk = du_trisk + wee(ij+u_ldown,4,2)*hflux(ij+t_ldown+u_rdown,l) 
     123            du_trisk = du_trisk + wee(ij+u_ldown,5,2)*hflux(ij+t_ldown+u_right,l) 
     124            du(ij+u_ldown,l) = du(ij+u_ldown,l) + du_trisk*qu(ij+u_ldown,l) 
     125         END DO 
     126      END DO 
     127   END SELECT 
    83128   !---------------------------- coriolis ---------------------------------- 
    84129   !-------------------------------------------------------------------------- 
  • codes/icosagcm/devel/src/kernels_unst/coriolis.k90

    r683 r940  
    176176   !$OMP END DO 
    177177   ! 
    178    !$OMP DO SCHEDULE(STATIC) 
    179    DO edge = 1, edge_num 
    180       ! this VLOOP iterates over the TRISK stencil 
    181       SELECT CASE(trisk_deg(edge)) 
    182       CASE(4) 
    183          !DIR$ SIMD 
    184          DO l = 1, llm 
    185             du_trisk=0. 
    186             itrisk = 1 
    187             edge_trisk = trisk(1,edge) 
    188             du_trisk = du_trisk + wee(itrisk,edge)*hflux(l,edge_trisk)*(qu(l,edge)+qu(l,edge_trisk)) 
    189             itrisk = 2 
    190             edge_trisk = trisk(2,edge) 
    191             du_trisk = du_trisk + wee(itrisk,edge)*hflux(l,edge_trisk)*(qu(l,edge)+qu(l,edge_trisk)) 
    192             itrisk = 3 
    193             edge_trisk = trisk(3,edge) 
    194             du_trisk = du_trisk + wee(itrisk,edge)*hflux(l,edge_trisk)*(qu(l,edge)+qu(l,edge_trisk)) 
    195             itrisk = 4 
    196             edge_trisk = trisk(4,edge) 
    197             du_trisk = du_trisk + wee(itrisk,edge)*hflux(l,edge_trisk)*(qu(l,edge)+qu(l,edge_trisk)) 
    198             du(l,edge) = du(l,edge) + .5*du_trisk 
    199          END DO 
    200       CASE(10) 
    201          !DIR$ SIMD 
    202          DO l = 1, llm 
    203             du_trisk=0. 
    204             itrisk = 1 
    205             edge_trisk = trisk(1,edge) 
    206             du_trisk = du_trisk + wee(itrisk,edge)*hflux(l,edge_trisk)*(qu(l,edge)+qu(l,edge_trisk)) 
    207             itrisk = 2 
    208             edge_trisk = trisk(2,edge) 
    209             du_trisk = du_trisk + wee(itrisk,edge)*hflux(l,edge_trisk)*(qu(l,edge)+qu(l,edge_trisk)) 
    210             itrisk = 3 
    211             edge_trisk = trisk(3,edge) 
    212             du_trisk = du_trisk + wee(itrisk,edge)*hflux(l,edge_trisk)*(qu(l,edge)+qu(l,edge_trisk)) 
    213             itrisk = 4 
    214             edge_trisk = trisk(4,edge) 
    215             du_trisk = du_trisk + wee(itrisk,edge)*hflux(l,edge_trisk)*(qu(l,edge)+qu(l,edge_trisk)) 
    216             itrisk = 5 
    217             edge_trisk = trisk(5,edge) 
    218             du_trisk = du_trisk + wee(itrisk,edge)*hflux(l,edge_trisk)*(qu(l,edge)+qu(l,edge_trisk)) 
    219             itrisk = 6 
    220             edge_trisk = trisk(6,edge) 
    221             du_trisk = du_trisk + wee(itrisk,edge)*hflux(l,edge_trisk)*(qu(l,edge)+qu(l,edge_trisk)) 
    222             itrisk = 7 
    223             edge_trisk = trisk(7,edge) 
    224             du_trisk = du_trisk + wee(itrisk,edge)*hflux(l,edge_trisk)*(qu(l,edge)+qu(l,edge_trisk)) 
    225             itrisk = 8 
    226             edge_trisk = trisk(8,edge) 
    227             du_trisk = du_trisk + wee(itrisk,edge)*hflux(l,edge_trisk)*(qu(l,edge)+qu(l,edge_trisk)) 
    228             itrisk = 9 
    229             edge_trisk = trisk(9,edge) 
    230             du_trisk = du_trisk + wee(itrisk,edge)*hflux(l,edge_trisk)*(qu(l,edge)+qu(l,edge_trisk)) 
    231             itrisk = 10 
    232             edge_trisk = trisk(10,edge) 
    233             du_trisk = du_trisk + wee(itrisk,edge)*hflux(l,edge_trisk)*(qu(l,edge)+qu(l,edge_trisk)) 
    234             du(l,edge) = du(l,edge) + .5*du_trisk 
    235          END DO 
    236       CASE DEFAULT 
    237          !DIR$ SIMD 
    238          DO l = 1, llm 
    239             du_trisk=0. 
    240             DO itrisk = 1, trisk_deg(edge) 
    241                edge_trisk = trisk(itrisk,edge) 
    242                du_trisk = du_trisk + wee(itrisk,edge)*hflux(l,edge_trisk)*(qu(l,edge)+qu(l,edge_trisk)) 
    243             END DO 
    244             du(l,edge) = du(l,edge) + .5*du_trisk 
    245          END DO 
    246       END SELECT 
    247    END DO 
    248    !$OMP END DO 
     178   SELECT CASE(caldyn_conserv) 
     179   CASE(conserv_energy) ! energy-conserving TRiSK 
     180      !$OMP DO SCHEDULE(STATIC) 
     181      DO edge = 1, edge_num 
     182         ! this VLOOP iterates over the TRISK stencil 
     183         SELECT CASE(trisk_deg(edge)) 
     184         CASE(4) 
     185            !DIR$ SIMD 
     186            DO l = 1, llm 
     187               du_trisk=0. 
     188               itrisk = 1 
     189               edge_trisk = trisk(1,edge) 
     190               du_trisk = du_trisk + wee(itrisk,edge,1)*hflux(l,edge_trisk)*(qu(l,edge)+qu(l,edge_trisk)) 
     191               itrisk = 2 
     192               edge_trisk = trisk(2,edge) 
     193               du_trisk = du_trisk + wee(itrisk,edge,1)*hflux(l,edge_trisk)*(qu(l,edge)+qu(l,edge_trisk)) 
     194               itrisk = 3 
     195               edge_trisk = trisk(3,edge) 
     196               du_trisk = du_trisk + wee(itrisk,edge,1)*hflux(l,edge_trisk)*(qu(l,edge)+qu(l,edge_trisk)) 
     197               itrisk = 4 
     198               edge_trisk = trisk(4,edge) 
     199               du_trisk = du_trisk + wee(itrisk,edge,1)*hflux(l,edge_trisk)*(qu(l,edge)+qu(l,edge_trisk)) 
     200               du(l,edge) = du(l,edge) + .5*du_trisk 
     201            END DO 
     202         CASE(10) 
     203            !DIR$ SIMD 
     204            DO l = 1, llm 
     205               du_trisk=0. 
     206               itrisk = 1 
     207               edge_trisk = trisk(1,edge) 
     208               du_trisk = du_trisk + wee(itrisk,edge,1)*hflux(l,edge_trisk)*(qu(l,edge)+qu(l,edge_trisk)) 
     209               itrisk = 2 
     210               edge_trisk = trisk(2,edge) 
     211               du_trisk = du_trisk + wee(itrisk,edge,1)*hflux(l,edge_trisk)*(qu(l,edge)+qu(l,edge_trisk)) 
     212               itrisk = 3 
     213               edge_trisk = trisk(3,edge) 
     214               du_trisk = du_trisk + wee(itrisk,edge,1)*hflux(l,edge_trisk)*(qu(l,edge)+qu(l,edge_trisk)) 
     215               itrisk = 4 
     216               edge_trisk = trisk(4,edge) 
     217               du_trisk = du_trisk + wee(itrisk,edge,1)*hflux(l,edge_trisk)*(qu(l,edge)+qu(l,edge_trisk)) 
     218               itrisk = 5 
     219               edge_trisk = trisk(5,edge) 
     220               du_trisk = du_trisk + wee(itrisk,edge,1)*hflux(l,edge_trisk)*(qu(l,edge)+qu(l,edge_trisk)) 
     221               itrisk = 6 
     222               edge_trisk = trisk(6,edge) 
     223               du_trisk = du_trisk + wee(itrisk,edge,1)*hflux(l,edge_trisk)*(qu(l,edge)+qu(l,edge_trisk)) 
     224               itrisk = 7 
     225               edge_trisk = trisk(7,edge) 
     226               du_trisk = du_trisk + wee(itrisk,edge,1)*hflux(l,edge_trisk)*(qu(l,edge)+qu(l,edge_trisk)) 
     227               itrisk = 8 
     228               edge_trisk = trisk(8,edge) 
     229               du_trisk = du_trisk + wee(itrisk,edge,1)*hflux(l,edge_trisk)*(qu(l,edge)+qu(l,edge_trisk)) 
     230               itrisk = 9 
     231               edge_trisk = trisk(9,edge) 
     232               du_trisk = du_trisk + wee(itrisk,edge,1)*hflux(l,edge_trisk)*(qu(l,edge)+qu(l,edge_trisk)) 
     233               itrisk = 10 
     234               edge_trisk = trisk(10,edge) 
     235               du_trisk = du_trisk + wee(itrisk,edge,1)*hflux(l,edge_trisk)*(qu(l,edge)+qu(l,edge_trisk)) 
     236               du(l,edge) = du(l,edge) + .5*du_trisk 
     237            END DO 
     238         CASE DEFAULT 
     239            !DIR$ SIMD 
     240            DO l = 1, llm 
     241               du_trisk=0. 
     242               DO itrisk = 1, trisk_deg(edge) 
     243                  edge_trisk = trisk(itrisk,edge) 
     244                  du_trisk = du_trisk + wee(itrisk,edge,1)*hflux(l,edge_trisk)*(qu(l,edge)+qu(l,edge_trisk)) 
     245               END DO 
     246               du(l,edge) = du(l,edge) + .5*du_trisk 
     247            END DO 
     248         END SELECT 
     249      END DO 
     250      !$OMP END DO 
     251   CASE(conserv_enstrophy) ! enstrophy-conserving TRiSK 
     252      !$OMP DO SCHEDULE(STATIC) 
     253      DO edge = 1, edge_num 
     254         ! this VLOOP iterates over the TRISK stencil 
     255         SELECT CASE(trisk_deg(edge)) 
     256         CASE(4) 
     257            !DIR$ SIMD 
     258            DO l = 1, llm 
     259               du_trisk=0. 
     260               itrisk = 1 
     261               edge_trisk = trisk(1,edge) 
     262               du_trisk = du_trisk + wee(itrisk,edge,1)*hflux(l,edge_trisk) 
     263               itrisk = 2 
     264               edge_trisk = trisk(2,edge) 
     265               du_trisk = du_trisk + wee(itrisk,edge,1)*hflux(l,edge_trisk) 
     266               itrisk = 3 
     267               edge_trisk = trisk(3,edge) 
     268               du_trisk = du_trisk + wee(itrisk,edge,1)*hflux(l,edge_trisk) 
     269               itrisk = 4 
     270               edge_trisk = trisk(4,edge) 
     271               du_trisk = du_trisk + wee(itrisk,edge,1)*hflux(l,edge_trisk) 
     272               du(l,edge) = du(l,edge) + du_trisk*qu(l,edge) 
     273            END DO 
     274         CASE(10) 
     275            !DIR$ SIMD 
     276            DO l = 1, llm 
     277               du_trisk=0. 
     278               itrisk = 1 
     279               edge_trisk = trisk(1,edge) 
     280               du_trisk = du_trisk + wee(itrisk,edge,1)*hflux(l,edge_trisk) 
     281               itrisk = 2 
     282               edge_trisk = trisk(2,edge) 
     283               du_trisk = du_trisk + wee(itrisk,edge,1)*hflux(l,edge_trisk) 
     284               itrisk = 3 
     285               edge_trisk = trisk(3,edge) 
     286               du_trisk = du_trisk + wee(itrisk,edge,1)*hflux(l,edge_trisk) 
     287               itrisk = 4 
     288               edge_trisk = trisk(4,edge) 
     289               du_trisk = du_trisk + wee(itrisk,edge,1)*hflux(l,edge_trisk) 
     290               itrisk = 5 
     291               edge_trisk = trisk(5,edge) 
     292               du_trisk = du_trisk + wee(itrisk,edge,1)*hflux(l,edge_trisk) 
     293               itrisk = 6 
     294               edge_trisk = trisk(6,edge) 
     295               du_trisk = du_trisk + wee(itrisk,edge,1)*hflux(l,edge_trisk) 
     296               itrisk = 7 
     297               edge_trisk = trisk(7,edge) 
     298               du_trisk = du_trisk + wee(itrisk,edge,1)*hflux(l,edge_trisk) 
     299               itrisk = 8 
     300               edge_trisk = trisk(8,edge) 
     301               du_trisk = du_trisk + wee(itrisk,edge,1)*hflux(l,edge_trisk) 
     302               itrisk = 9 
     303               edge_trisk = trisk(9,edge) 
     304               du_trisk = du_trisk + wee(itrisk,edge,1)*hflux(l,edge_trisk) 
     305               itrisk = 10 
     306               edge_trisk = trisk(10,edge) 
     307               du_trisk = du_trisk + wee(itrisk,edge,1)*hflux(l,edge_trisk) 
     308               du(l,edge) = du(l,edge) + du_trisk*qu(l,edge) 
     309            END DO 
     310         CASE DEFAULT 
     311            !DIR$ SIMD 
     312            DO l = 1, llm 
     313               du_trisk=0. 
     314               DO itrisk = 1, trisk_deg(edge) 
     315                  edge_trisk = trisk(itrisk,edge) 
     316                  du_trisk = du_trisk + wee(itrisk,edge,1)*hflux(l,edge_trisk) 
     317               END DO 
     318               du(l,edge) = du(l,edge) + du_trisk*qu(l,edge) 
     319            END DO 
     320         END SELECT 
     321      END DO 
     322      !$OMP END DO 
     323   END SELECT 
    249324   !---------------------------- coriolis ---------------------------------- 
    250325   !-------------------------------------------------------------------------- 
  • codes/icosagcm/devel/src/unstructured/caldyn_unstructured.F90

    r935 r940  
    163163  FIELD_MASS  :: convm 
    164164  FIELD_THETA :: theta, dtheta_rhodz 
     165  INTEGER, PARAMETER :: conserv_energy=1, conserv_enstrophy=2, caldyn_conserv = conserv_enstrophy ! FIXME 
    165166  DECLARE_INDICES 
    166167  DECLARE_EDGES 
     
    211212!----------------------------- Unused ----------------------------- 
    212213 
     214#ifdef BEGIN_DYSL 
     215KERNEL(gradient) 
     216  FORALL_CELLS_EXT() 
     217    ON_EDGES 
     218      grad(EDGE) = SIGN*(b(CELL2)-b(CELL1)) 
     219    END_BLOCK 
     220  END_BLOCK 
     221END_BLOCK 
     222 
     223KERNEL(div) 
     224  FORALL_CELLS_EXT() 
     225    ON_PRIMAL 
     226      div_ij=0. 
     227      FORALL_EDGES 
     228        div_ij = div_ij + SIGN*LE_DE*u(EDGE) 
     229      END_BLOCK 
     230      divu(CELL) = div_ij / AI 
     231    END_BLOCK 
     232  END_BLOCK 
     233END_BLOCK 
     234#endif END_DYSL 
     235 
    213236SUBROUTINE gradient(b,grad) BINDC(gradient) 
    214237  FIELD_MASS :: b 
  • codes/icosagcm/devel/src/unstructured/data_unstructured.F90

    r936 r940  
    4343  NUM1(max_nb_stage), BIND(C)              :: tauj       ! diagonal of fast Butcher tableau 
    4444  NUM2(max_nb_stage,max_nb_stage), BIND(C) :: cslj, cflj ! slow and fast modified Butcher tableaus 
    45   NUM2(:,:), POINTER          :: centroid, xyz_v, Riv2, wee, ap,bp, mass_bl, mass_dak, mass_dbk 
    46  
     45  NUM2(:,:), POINTER          :: centroid, xyz_v, Riv2, ap,bp, mass_bl, mass_dak, mass_dbk 
     46  NUM3(:,:,:), POINTER        :: wee 
    4747  INTEGER(C_INT), BIND(C) :: comm_icosa 
    4848 
     
    147147#define ALLOC1(v,n1) IF(ASSOCIATED(v)) DEALLOCATE(v) ; ALLOCATE(v(n1)) 
    148148#define ALLOC2(v,n1,n2) IF(ASSOCIATED(v)) DEALLOCATE(v) ; ALLOCATE(v(n1,n2)) 
     149#define ALLOC3(v,n1,n2,n3) IF(ASSOCIATED(v)) DEALLOCATE(v) ; ALLOCATE(v(n1,n2,n3)) 
    149150 
    150151  SUBROUTINE init_mesh( &  
     
    225226    ALLOC1(le_de,edge_num) 
    226227    ALLOC2(Riv2, max_dual_deg, dual_num) 
    227     ALLOC2(wee, max_trisk_deg, edge_num) 
     228    ALLOC3(wee, max_trisk_deg, edge_num, 1) 
    228229    Ai(:) = Ai_(:) 
    229230    Av(:) = Av_(:) 
     
    231232    le_de(:) = le_de_(:) 
    232233    Riv2(:,:)=Riv2_(:,:) 
    233     wee(:,:) = wee_(:,:) 
     234    wee(:,:,1) = wee_(:,:) 
    234235    IF(is_mpi_master) THEN 
    235236       PRINT *, 'Max Ai : ',    MAXVAL(ABS(Ai)) 
  • codes/icosagcm/devel/src/unstructured/init_unstructured.f90

    r883 r940  
    282282    ALLOCATE(Riv2, source = Ddata_read2) 
    283283    CALL read_from_gridfile(id_nc, 'float', 'wee') 
    284     ALLOCATE(wee, source = Ddata_read2) 
     284    ALLOCATE(wee(SIZE(Ddata_read2,1), SIZE(Ddata_read2,2), 1)) 
     285    wee(:,:,1) = Ddata_read2(:,:) 
    285286 
    286287    ! read cell centers and bounds 
Note: See TracChangeset for help on using the changeset viewer.