Changeset 658


Ignore:
Timestamp:
12/30/17 02:00:38 (6 years ago)
Author:
dubos
Message:

devel/unstructured : updated kernels

Location:
codes/icosagcm/devel/src
Files:
1 added
16 edited

Legend:

Unmodified
Added
Removed

  • codes/icosagcm/devel/src/kernels_unst/caldyn_dmass.k90

    r624 r658  
    33   !$OMP DO SCHEDULE(STATIC) 
    44   DO ij = 1, primal_num 
     5      !DIR$ SIMD 
    56      DO l = 1, llm 
    67         convm(l,ij) = mass_dbk(l,ij) * dmass_col(ij) 

  • codes/icosagcm/devel/src/kernels_unst/caldyn_fast.k90

    r614 r658  
    66      !$OMP DO SCHEDULE(STATIC) 
    77      DO ij = 1, primal_num 
     8         !DIR$ SIMD 
    89         DO l = 1, llm 
    910            berni(l,ij) = pk(l,ij) 
    … …  
    1617      !$OMP DO SCHEDULE(STATIC) 
    1718      DO ij = 1, primal_num 
     19         !DIR$ SIMD 
    1820         DO l = 1, llm 
    1921            berni(l,ij) = .5*(geopot(l,ij)+geopot(l+1,ij)) 
    … …  
    2426      !$OMP DO SCHEDULE(STATIC) 
    2527      DO ij = 1, primal_num 
     28         !DIR$ SIMD 
    2629         DO l = 1, llm 
    2730            berni(l,ij) = .5*(geopot(l,ij)+geopot(l+1,ij)) 
    … …  
    3942      ij_left = left(edge) 
    4043      ij_right = right(edge) 
     44      !DIR$ SIMD 
    4145      DO l = 1, llm 
    4246         due = .5*(theta(l,ij_left,1)+theta(l,ij_right,1))*(pk(l,ij_right)-pk(l,ij_left)) + berni(l,ij_right)-berni(l,ij_left) 

  • codes/icosagcm/devel/src/kernels_unst/caldyn_slow_NH.k90

    r614 r658  
    77      l=1 
    88      w_il(l,ij) = 2.*W(l,ij)/rhodz(kup,ij) 
     9      !DIR$ SIMD 
    910      DO l = 2, llm 
    1011         kdown = l-1 
    … …  
    3435      v_el(l,edge) = .5*le_de(edge)*(u(kup,edge)+u(kdown,edge)) ! checked 
    3536      G_el(l,edge) = v_el(l,edge)*W_el - DePhil(l,edge)*W2_el 
     37      !DIR$ SIMD 
    3638      DO l = 2, llm 
    3739         kdown = l-1 
    … …  
    6466   !$OMP DO SCHEDULE(STATIC) 
    6567   DO ij = 1, primal_num 
    66       DO l = 1, llm+1 
    67          gPhi2=0. 
    68          dP=0. 
    69          divG=0 
    70          DO iedge = 1, primal_deg(ij) 
    71             edge = primal_edge(iedge,ij) 
    72             gPhi2 = gPhi2 + le_de(edge)*DePhil(l,edge)**2 
    73             dP = dP + le_de(edge)*DePhil(l,edge)*v_el(l,edge) 
    74             divG = divG + primal_ne(iedge,ij)*G_el(l,edge) ! -div(G_el), G_el already has le_de 
    75          END DO 
    76          gradPhi2(l,ij) = 1./(2.*Ai(ij)) * gPhi2 
    77          dPhi(l,ij) = gradPhi2(l,ij)*w_il(l,ij) - 1./(2.*Ai(ij))*dP 
    78          dW(l,ij) = (-1./Ai(ij))*divG 
    79       END DO 
     68      ! this VLOOP iterates over primal cell edges 
     69      SELECT CASE(primal_deg(ij)) 
     70      CASE(4) 
     71         edge1 = primal_edge(1,ij) 
     72         edge2 = primal_edge(2,ij) 
     73         edge3 = primal_edge(3,ij) 
     74         edge4 = primal_edge(4,ij) 
     75         le_de1 = le_de(edge1) 
     76         le_de2 = le_de(edge2) 
     77         le_de3 = le_de(edge3) 
     78         le_de4 = le_de(edge4) 
     79         sign1 = primal_ne(1,ij) 
     80         sign2 = primal_ne(2,ij) 
     81         sign3 = primal_ne(3,ij) 
     82         sign4 = primal_ne(4,ij) 
     83         !DIR$ SIMD 
     84         DO l = 1, llm+1 
     85            gPhi2=0. 
     86            dP=0. 
     87            divG=0 
     88            gPhi2 = gPhi2 + le_de1*DePhil(l,edge1)**2 
     89            dP = dP + le_de1*DePhil(l,edge1)*v_el(l,edge1) 
     90            divG = divG + sign1*G_el(l,edge1) ! -div(G_el), G_el already has le_de 
     91            gPhi2 = gPhi2 + le_de2*DePhil(l,edge2)**2 
     92            dP = dP + le_de2*DePhil(l,edge2)*v_el(l,edge2) 
     93            divG = divG + sign2*G_el(l,edge2) ! -div(G_el), G_el already has le_de 
     94            gPhi2 = gPhi2 + le_de3*DePhil(l,edge3)**2 
     95            dP = dP + le_de3*DePhil(l,edge3)*v_el(l,edge3) 
     96            divG = divG + sign3*G_el(l,edge3) ! -div(G_el), G_el already has le_de 
     97            gPhi2 = gPhi2 + le_de4*DePhil(l,edge4)**2 
     98            dP = dP + le_de4*DePhil(l,edge4)*v_el(l,edge4) 
     99            divG = divG + sign4*G_el(l,edge4) ! -div(G_el), G_el already has le_de 
     100            gradPhi2(l,ij) = 1./(2.*Ai(ij)) * gPhi2 
     101            dPhi(l,ij) = gradPhi2(l,ij)*w_il(l,ij) - 1./(2.*Ai(ij))*dP 
     102            dW(l,ij) = (-1./Ai(ij))*divG 
     103         END DO 
     104      CASE(5) 
     105         edge1 = primal_edge(1,ij) 
     106         edge2 = primal_edge(2,ij) 
     107         edge3 = primal_edge(3,ij) 
     108         edge4 = primal_edge(4,ij) 
     109         edge5 = primal_edge(5,ij) 
     110         le_de1 = le_de(edge1) 
     111         le_de2 = le_de(edge2) 
     112         le_de3 = le_de(edge3) 
     113         le_de4 = le_de(edge4) 
     114         le_de5 = le_de(edge5) 
     115         sign1 = primal_ne(1,ij) 
     116         sign2 = primal_ne(2,ij) 
     117         sign3 = primal_ne(3,ij) 
     118         sign4 = primal_ne(4,ij) 
     119         sign5 = primal_ne(5,ij) 
     120         !DIR$ SIMD 
     121         DO l = 1, llm+1 
     122            gPhi2=0. 
     123            dP=0. 
     124            divG=0 
     125            gPhi2 = gPhi2 + le_de1*DePhil(l,edge1)**2 
     126            dP = dP + le_de1*DePhil(l,edge1)*v_el(l,edge1) 
     127            divG = divG + sign1*G_el(l,edge1) ! -div(G_el), G_el already has le_de 
     128            gPhi2 = gPhi2 + le_de2*DePhil(l,edge2)**2 
     129            dP = dP + le_de2*DePhil(l,edge2)*v_el(l,edge2) 
     130            divG = divG + sign2*G_el(l,edge2) ! -div(G_el), G_el already has le_de 
     131            gPhi2 = gPhi2 + le_de3*DePhil(l,edge3)**2 
     132            dP = dP + le_de3*DePhil(l,edge3)*v_el(l,edge3) 
     133            divG = divG + sign3*G_el(l,edge3) ! -div(G_el), G_el already has le_de 
     134            gPhi2 = gPhi2 + le_de4*DePhil(l,edge4)**2 
     135            dP = dP + le_de4*DePhil(l,edge4)*v_el(l,edge4) 
     136            divG = divG + sign4*G_el(l,edge4) ! -div(G_el), G_el already has le_de 
     137            gPhi2 = gPhi2 + le_de5*DePhil(l,edge5)**2 
     138            dP = dP + le_de5*DePhil(l,edge5)*v_el(l,edge5) 
     139            divG = divG + sign5*G_el(l,edge5) ! -div(G_el), G_el already has le_de 
     140            gradPhi2(l,ij) = 1./(2.*Ai(ij)) * gPhi2 
     141            dPhi(l,ij) = gradPhi2(l,ij)*w_il(l,ij) - 1./(2.*Ai(ij))*dP 
     142            dW(l,ij) = (-1./Ai(ij))*divG 
     143         END DO 
     144      CASE(6) 
     145         edge1 = primal_edge(1,ij) 
     146         edge2 = primal_edge(2,ij) 
     147         edge3 = primal_edge(3,ij) 
     148         edge4 = primal_edge(4,ij) 
     149         edge5 = primal_edge(5,ij) 
     150         edge6 = primal_edge(6,ij) 
     151         le_de1 = le_de(edge1) 
     152         le_de2 = le_de(edge2) 
     153         le_de3 = le_de(edge3) 
     154         le_de4 = le_de(edge4) 
     155         le_de5 = le_de(edge5) 
     156         le_de6 = le_de(edge6) 
     157         sign1 = primal_ne(1,ij) 
     158         sign2 = primal_ne(2,ij) 
     159         sign3 = primal_ne(3,ij) 
     160         sign4 = primal_ne(4,ij) 
     161         sign5 = primal_ne(5,ij) 
     162         sign6 = primal_ne(6,ij) 
     163         !DIR$ SIMD 
     164         DO l = 1, llm+1 
     165            gPhi2=0. 
     166            dP=0. 
     167            divG=0 
     168            gPhi2 = gPhi2 + le_de1*DePhil(l,edge1)**2 
     169            dP = dP + le_de1*DePhil(l,edge1)*v_el(l,edge1) 
     170            divG = divG + sign1*G_el(l,edge1) ! -div(G_el), G_el already has le_de 
     171            gPhi2 = gPhi2 + le_de2*DePhil(l,edge2)**2 
     172            dP = dP + le_de2*DePhil(l,edge2)*v_el(l,edge2) 
     173            divG = divG + sign2*G_el(l,edge2) ! -div(G_el), G_el already has le_de 
     174            gPhi2 = gPhi2 + le_de3*DePhil(l,edge3)**2 
     175            dP = dP + le_de3*DePhil(l,edge3)*v_el(l,edge3) 
     176            divG = divG + sign3*G_el(l,edge3) ! -div(G_el), G_el already has le_de 
     177            gPhi2 = gPhi2 + le_de4*DePhil(l,edge4)**2 
     178            dP = dP + le_de4*DePhil(l,edge4)*v_el(l,edge4) 
     179            divG = divG + sign4*G_el(l,edge4) ! -div(G_el), G_el already has le_de 
     180            gPhi2 = gPhi2 + le_de5*DePhil(l,edge5)**2 
     181            dP = dP + le_de5*DePhil(l,edge5)*v_el(l,edge5) 
     182            divG = divG + sign5*G_el(l,edge5) ! -div(G_el), G_el already has le_de 
     183            gPhi2 = gPhi2 + le_de6*DePhil(l,edge6)**2 
     184            dP = dP + le_de6*DePhil(l,edge6)*v_el(l,edge6) 
     185            divG = divG + sign6*G_el(l,edge6) ! -div(G_el), G_el already has le_de 
     186            gradPhi2(l,ij) = 1./(2.*Ai(ij)) * gPhi2 
     187            dPhi(l,ij) = gradPhi2(l,ij)*w_il(l,ij) - 1./(2.*Ai(ij))*dP 
     188            dW(l,ij) = (-1./Ai(ij))*divG 
     189         END DO 
     190      CASE DEFAULT 
     191         !DIR$ SIMD 
     192         DO l = 1, llm+1 
     193            gPhi2=0. 
     194            dP=0. 
     195            divG=0 
     196            DO iedge = 1, primal_deg(ij) 
     197               edge = primal_edge(iedge,ij) 
     198               gPhi2 = gPhi2 + le_de(edge)*DePhil(l,edge)**2 
     199               dP = dP + le_de(edge)*DePhil(l,edge)*v_el(l,edge) 
     200               divG = divG + primal_ne(iedge,ij)*G_el(l,edge) ! -div(G_el), G_el already has le_de 
     201            END DO 
     202            gradPhi2(l,ij) = 1./(2.*Ai(ij)) * gPhi2 
     203            dPhi(l,ij) = gradPhi2(l,ij)*w_il(l,ij) - 1./(2.*Ai(ij))*dP 
     204            dW(l,ij) = (-1./Ai(ij))*divG 
     205         END DO 
     206      END SELECT 
    80207   END DO 
    81208   !$OMP END DO 
    … …  
    85212   !$OMP DO SCHEDULE(STATIC) 
    86213   DO ij = 1, primal_num 
    87       DO l = 1, llm 
    88          u2=0. 
    89          DO iedge = 1, primal_deg(ij) 
    90             edge = primal_edge(iedge,ij) 
    91             u2 = u2 + le_de(edge)*u(l,edge)**2 
    92          END DO 
    93          berni(l,ij) = 1./(4.*Ai(ij)) * u2 - .25*( gradPhi2(l,ij)*w_il(l,ij)**2 + gradPhi2(l+1,ij)*w_il(l+1,ij)**2 ) 
    94       END DO 
     214      ! this VLOOP iterates over primal cell edges 
     215      SELECT CASE(primal_deg(ij)) 
     216      CASE(4) 
     217         edge1 = primal_edge(1,ij) 
     218         edge2 = primal_edge(2,ij) 
     219         edge3 = primal_edge(3,ij) 
     220         edge4 = primal_edge(4,ij) 
     221         le_de1 = le_de(edge1) 
     222         le_de2 = le_de(edge2) 
     223         le_de3 = le_de(edge3) 
     224         le_de4 = le_de(edge4) 
     225         !DIR$ SIMD 
     226         DO l = 1, llm 
     227            u2=0. 
     228            u2 = u2 + le_de1*u(l,edge1)**2 
     229            u2 = u2 + le_de2*u(l,edge2)**2 
     230            u2 = u2 + le_de3*u(l,edge3)**2 
     231            u2 = u2 + le_de4*u(l,edge4)**2 
     232            berni(l,ij) = 1./(4.*Ai(ij)) * u2 - .25*( gradPhi2(l,ij)*w_il(l,ij)**2 + gradPhi2(l+1,ij)*w_il(l+1,ij)**2 ) 
     233         END DO 
     234      CASE(5) 
     235         edge1 = primal_edge(1,ij) 
     236         edge2 = primal_edge(2,ij) 
     237         edge3 = primal_edge(3,ij) 
     238         edge4 = primal_edge(4,ij) 
     239         edge5 = primal_edge(5,ij) 
     240         le_de1 = le_de(edge1) 
     241         le_de2 = le_de(edge2) 
     242         le_de3 = le_de(edge3) 
     243         le_de4 = le_de(edge4) 
     244         le_de5 = le_de(edge5) 
     245         !DIR$ SIMD 
     246         DO l = 1, llm 
     247            u2=0. 
     248            u2 = u2 + le_de1*u(l,edge1)**2 
     249            u2 = u2 + le_de2*u(l,edge2)**2 
     250            u2 = u2 + le_de3*u(l,edge3)**2 
     251            u2 = u2 + le_de4*u(l,edge4)**2 
     252            u2 = u2 + le_de5*u(l,edge5)**2 
     253            berni(l,ij) = 1./(4.*Ai(ij)) * u2 - .25*( gradPhi2(l,ij)*w_il(l,ij)**2 + gradPhi2(l+1,ij)*w_il(l+1,ij)**2 ) 
     254         END DO 
     255      CASE(6) 
     256         edge1 = primal_edge(1,ij) 
     257         edge2 = primal_edge(2,ij) 
     258         edge3 = primal_edge(3,ij) 
     259         edge4 = primal_edge(4,ij) 
     260         edge5 = primal_edge(5,ij) 
     261         edge6 = primal_edge(6,ij) 
     262         le_de1 = le_de(edge1) 
     263         le_de2 = le_de(edge2) 
     264         le_de3 = le_de(edge3) 
     265         le_de4 = le_de(edge4) 
     266         le_de5 = le_de(edge5) 
     267         le_de6 = le_de(edge6) 
     268         !DIR$ SIMD 
     269         DO l = 1, llm 
     270            u2=0. 
     271            u2 = u2 + le_de1*u(l,edge1)**2 
     272            u2 = u2 + le_de2*u(l,edge2)**2 
     273            u2 = u2 + le_de3*u(l,edge3)**2 
     274            u2 = u2 + le_de4*u(l,edge4)**2 
     275            u2 = u2 + le_de5*u(l,edge5)**2 
     276            u2 = u2 + le_de6*u(l,edge6)**2 
     277            berni(l,ij) = 1./(4.*Ai(ij)) * u2 - .25*( gradPhi2(l,ij)*w_il(l,ij)**2 + gradPhi2(l+1,ij)*w_il(l+1,ij)**2 ) 
     278         END DO 
     279      CASE DEFAULT 
     280         !DIR$ SIMD 
     281         DO l = 1, llm 
     282            u2=0. 
     283            DO iedge = 1, primal_deg(ij) 
     284               edge = primal_edge(iedge,ij) 
     285               u2 = u2 + le_de(edge)*u(l,edge)**2 
     286            END DO 
     287            berni(l,ij) = 1./(4.*Ai(ij)) * u2 - .25*( gradPhi2(l,ij)*w_il(l,ij)**2 + gradPhi2(l+1,ij)*w_il(l+1,ij)**2 ) 
     288         END DO 
     289      END SELECT 
    95290   END DO 
    96291   !$OMP END DO 
    … …  
    99294      ij_left = left(edge) 
    100295      ij_right = right(edge) 
     296      !DIR$ SIMD 
    101297      DO l = 1, llm 
    102298         ! Compute mass flux and grad(berni) 

  • codes/icosagcm/devel/src/kernels_unst/caldyn_slow_hydro.k90

    r614 r658  
    55      ij_left = left(edge) 
    66      ij_right = right(edge) 
     7      !DIR$ SIMD 
    78      DO l = 1, llm 
    89         uu = .5*(rhodz(l,ij_left)+rhodz(l,ij_right))*u(l,edge) 
    … …  
    1314   !$OMP DO SCHEDULE(STATIC) 
    1415   DO ij = 1, primal_num 
    15       DO l = 1, llm 
    16          ke=0.d0 
    17          DO iedge = 1, primal_deg(ij) 
    18             edge = primal_edge(iedge,ij) 
    19             ke = ke + le_de(edge)*u(l,edge)**2 
     16      ! this VLOOP iterates over primal cell edges 
     17      SELECT CASE(primal_deg(ij)) 
     18      CASE(4) 
     19         edge1 = primal_edge(1,ij) 
     20         edge2 = primal_edge(2,ij) 
     21         edge3 = primal_edge(3,ij) 
     22         edge4 = primal_edge(4,ij) 
     23         le_de1 = le_de(edge1) 
     24         le_de2 = le_de(edge2) 
     25         le_de3 = le_de(edge3) 
     26         le_de4 = le_de(edge4) 
     27         !DIR$ SIMD 
     28         DO l = 1, llm 
     29            ke=0.d0 
     30            ke = ke + le_de1*u(l,edge1)**2 
     31            ke = ke + le_de2*u(l,edge2)**2 
     32            ke = ke + le_de3*u(l,edge3)**2 
     33            ke = ke + le_de4*u(l,edge4)**2 
     34            BERNI(l,ij)=ke*(.25/Ai(ij)) 
    2035         END DO 
    21          BERNI(l,ij)=ke*(.25/Ai(ij)) 
    22       END DO 
     36      CASE(5) 
     37         edge1 = primal_edge(1,ij) 
     38         edge2 = primal_edge(2,ij) 
     39         edge3 = primal_edge(3,ij) 
     40         edge4 = primal_edge(4,ij) 
     41         edge5 = primal_edge(5,ij) 
     42         le_de1 = le_de(edge1) 
     43         le_de2 = le_de(edge2) 
     44         le_de3 = le_de(edge3) 
     45         le_de4 = le_de(edge4) 
     46         le_de5 = le_de(edge5) 
     47         !DIR$ SIMD 
     48         DO l = 1, llm 
     49            ke=0.d0 
     50            ke = ke + le_de1*u(l,edge1)**2 
     51            ke = ke + le_de2*u(l,edge2)**2 
     52            ke = ke + le_de3*u(l,edge3)**2 
     53            ke = ke + le_de4*u(l,edge4)**2 
     54            ke = ke + le_de5*u(l,edge5)**2 
     55            BERNI(l,ij)=ke*(.25/Ai(ij)) 
     56         END DO 
     57      CASE(6) 
     58         edge1 = primal_edge(1,ij) 
     59         edge2 = primal_edge(2,ij) 
     60         edge3 = primal_edge(3,ij) 
     61         edge4 = primal_edge(4,ij) 
     62         edge5 = primal_edge(5,ij) 
     63         edge6 = primal_edge(6,ij) 
     64         le_de1 = le_de(edge1) 
     65         le_de2 = le_de(edge2) 
     66         le_de3 = le_de(edge3) 
     67         le_de4 = le_de(edge4) 
     68         le_de5 = le_de(edge5) 
     69         le_de6 = le_de(edge6) 
     70         !DIR$ SIMD 
     71         DO l = 1, llm 
     72            ke=0.d0 
     73            ke = ke + le_de1*u(l,edge1)**2 
     74            ke = ke + le_de2*u(l,edge2)**2 
     75            ke = ke + le_de3*u(l,edge3)**2 
     76            ke = ke + le_de4*u(l,edge4)**2 
     77            ke = ke + le_de5*u(l,edge5)**2 
     78            ke = ke + le_de6*u(l,edge6)**2 
     79            BERNI(l,ij)=ke*(.25/Ai(ij)) 
     80         END DO 
     81      CASE DEFAULT 
     82         !DIR$ SIMD 
     83         DO l = 1, llm 
     84            ke=0.d0 
     85            DO iedge = 1, primal_deg(ij) 
     86               edge = primal_edge(iedge,ij) 
     87               ke = ke + le_de(edge)*u(l,edge)**2 
     88            END DO 
     89            BERNI(l,ij)=ke*(.25/Ai(ij)) 
     90         END DO 
     91      END SELECT 
    2392   END DO 
    2493   !$OMP END DO 
    … …  
    2897         ij_left = left(edge) 
    2998         ij_right = right(edge) 
     99         !DIR$ SIMD 
    30100         DO l = 1, llm 
    31101            du(l,edge) = 1.*(berni(l,ij_left)-berni(l,ij_right)) ! minus gradient 
    … …  
    38108         ij_left = left(edge) 
    39109         ij_right = right(edge) 
     110         !DIR$ SIMD 
    40111         DO l = 1, llm 
    41112            du(l,edge) = du(l,edge) + 1.*(berni(l,ij_left)-berni(l,ij_right)) ! minus gradient 

  • codes/icosagcm/devel/src/kernels_unst/caldyn_solver.k90

    r614 r658  
    11   !-------------------------------------------------------------------------- 
    22   !---------------------------- caldyn_solver ---------------------------------- 
    3    !$OMP DO SCHEDULE(STATIC) 
    4    DO ij = 1, primal_num 
    5       l=1 
    6       m_il(l,ij) = .5*rhodz(l,ij) 
    7       DO l = 2, llm 
    8          m_il(l,ij) = .5*(rhodz(l,ij)+rhodz(l-1,ij)) 
    9       END DO 
    10       l=llm+1 
    11       m_il(l,ij) = .5*rhodz(l-1,ij) 
    12    END DO 
    13    !$OMP END DO 
    14    ! 
    15    IF(tau>0) THEN ! solve implicit problem for geopotential 
    16       CALL compute_NH_geopot(tau,PHI_BOT_VAR, rhodz, m_il, theta, W, geopot) 
    17    END IF 
    183   ! 
    194   ! Compute pressure (pres) and Exner function (pk) 
    … …  
    238   gamma = 1./(1.-kappa) 
    249   vreff = Rd*Treff/preff ! reference specific volume 
    25    Cvd = cpp-Rd 
     10   Cvd = 1./(cpp-Rd) 
     11   Rd_preff = kappa*cpp/preff 
    2612   SELECT CASE(caldyn_thermo) 
    2713   CASE(thermo_theta) 
    2814      !$OMP DO SCHEDULE(STATIC) 
    2915      DO ij = 1, primal_num 
     16         !DIR$ SIMD 
    3017         DO l = 1, llm 
    31             rho_ij = g*rhodz(l,ij)/(geopot(l+1,ij)-geopot(l,ij)) 
    32             X_ij = (cpp/preff)*kappa*theta(l,ij,1)*rho_ij 
     18            rho_ij = 1./(geopot(l+1,ij)-geopot(l,ij)) 
     19            rho_ij = rho_ij*g*rhodz(l,ij) 
     20            X_ij = Rd_preff*theta(l,ij,1)*rho_ij 
    3321            ! kappa.theta.rho = p/exner 
    3422            ! => X = (p/p0)/(exner/Cp) 
    … …  
    4432      !$OMP DO SCHEDULE(STATIC) 
    4533      DO ij = 1, primal_num 
     34         !DIR$ SIMD 
    4635         DO l = 1, llm 
    47             rho_ij = g*rhodz(l,ij)/(geopot(l+1,ij)-geopot(l,ij)) 
    48             T_ij = Treff*exp( (theta(l,ij,1)+Rd*log(vreff*rho_ij))/Cvd ) 
     36            rho_ij = 1./(geopot(l+1,ij)-geopot(l,ij)) 
     37            rho_ij = rho_ij*g*rhodz(l,ij) 
     38            T_ij = Treff*exp( (theta(l,ij,1)+Rd*log(vreff*rho_ij))*Cvd ) 
    4939            pres(l,ij) = rho_ij*Rd*T_ij 
    5040            pk(l,ij) = T_ij 
    … …  
    6454      W(l,ij) = W(l,ij)+tau*dW(l,ij) ! update W 
    6555      dPhi(l,ij) = g*g*W(l,ij)/m_il(l,ij) 
     56      !DIR$ SIMD 
    6657      DO l = 2, llm 
    6758         dW(l,ij) = (1./g)*(pres(l-1,ij)-pres(l,ij)) - m_il(l,ij) 
    … …  
    8071   !$OMP DO SCHEDULE(STATIC) 
    8172   DO ij = 1, primal_num 
     73      !DIR$ SIMD 
    8274      DO l = 1, llm 
    8375         ! compute du = -0.5*g^2.grad(w^2) 
    … …  
    9082      ij_left = left(edge) 
    9183      ij_right = right(edge) 
     84      !DIR$ SIMD 
    9285      DO l = 1, llm 
    9386         du(l,edge) = 1.*(berni(l,ij_left)-berni(l,ij_right)) 

  • codes/icosagcm/devel/src/kernels_unst/caldyn_vert.k90

    r624 r658  
    44      !$OMP DO SCHEDULE(STATIC) 
    55      DO ij = 1, primal_num 
     6         !DIR$ SIMD 
    67         DO l = 2, llm 
    78            dtheta_rhodz(l,ij,iq) = dtheta_rhodz(l,ij,iq) + 0.5*(theta(l,ij,iq)+theta(l-1,ij,iq))*wflux(l,ij) 
    … …  
    1112      !$OMP DO SCHEDULE(STATIC) 
    1213      DO ij = 1, primal_num 
     14         !DIR$ SIMD 
    1315         DO l = 1, llm-1 
    1416            dtheta_rhodz(l,ij,iq) = dtheta_rhodz(l,ij,iq) - 0.5*(theta(l,ij,iq)+theta(l+1,ij,iq))*wflux(l+1,ij) 
    … …  
    2325         ij_left = left(edge) 
    2426         ij_right = right(edge) 
     27         !DIR$ SIMD 
    2528         DO l = 2, llm 
    2629            wwuu(l,edge) = .25*(wflux(l,ij_left)+wflux(l,ij_right))*(u(l,edge)+u(l-1,edge)) ! Fu 
    … …  
    3437         ij_left = left(edge) 
    3538         ij_right = right(edge) 
     39         !DIR$ SIMD 
    3640         DO l = 1, llm 
    3741            dFu_deta = wwuu(l+1,edge)-wwuu(l,edge) ! d/deta (F*u) 
    … …  
    4650         ij_left = left(edge) 
    4751         ij_right = right(edge) 
     52         !DIR$ SIMD 
    4853         DO l = 2, llm 
    4954            wwuu(l,edge) = .5*(wflux(l,ij_left)+wflux(l,ij_right))*(u(l,edge)-u(l-1,edge)) 
    … …  
    5762         ij_left = left(edge) 
    5863         ij_right = right(edge) 
     64         !DIR$ SIMD 
    5965         DO l = 1, llm 
    6066            du(l,edge) = du(l,edge) - (wwuu(l,edge)+wwuu(l+1,edge)) / (rhodz(l,ij_left)+rhodz(l,ij_right)) 

  • codes/icosagcm/devel/src/kernels_unst/caldyn_vert_NH.k90

    r614 r658  
    99      eta_dot(l,ij) = wflux_ij / mass(l,ij) 
    1010      wcov(l,ij) = w_ij*(geopot(l+1,ij)-geopot(l,ij)) 
     11      !DIR$ SIMD 
    1112      DO l = 2, llm-1 
    1213         w_ij = .5*( W(l,ij)+W(l+1,ij) )/mass(l,ij) 
    … …  
    2930      ij_left = left(edge) 
    3031      ij_right = right(edge) 
     32      !DIR$ SIMD 
    3133      DO l = 1, llm 
    3234         du(l,edge) = du(l,edge) - .5*(wcov(l,ij_left)+wcov(l,ij_right))*1.*(eta_dot(l,ij_right)-eta_dot(l,ij_left)) 
    … …  
    3840   !$OMP DO SCHEDULE(STATIC) 
    3941   DO ij = 1, primal_num 
     42      !DIR$ SIMD 
    4043      DO l = 2, llm 
    4144         dPhi(l,ij)=dPhi(l,ij)-wflux(l,ij)*(geopot(l+1,ij)-geopot(l-1,ij))/(mass(l-1,ij)+mass(l,ij)) 
    … …  
    4952      l=1 
    5053      dW(l,ij) = dW(l,ij) - W_etadot(l,ij) 
     54      !DIR$ SIMD 
    5155      DO l = 2, llm 
    5256         dW(l,ij) = dW(l,ij) + W_etadot(l-1,ij) - W_etadot(l,ij) 

  • codes/icosagcm/devel/src/kernels_unst/compute_NH_geopot.k90

    r614 r658  
    44   g2=g*g 
    55   gm2 = 1./g2 
     6   vreff = Treff*cpp/preff*kappa 
    67   gamma = 1./(1.-kappa) 
    78   !$OMP BARRIER 
    … …  
    3536            DO l = 1,llm 
    3637               rho_ij = (g*m_ik(l,ij))/(Phi_il(l+1,ij)-Phi_il(l,ij)) 
    37                X_ij = Treff*exp(theta(l,ij)/cpp) ! theta = Tref.exp(s/Cp) 
    38                X_ij = (cpp/preff)*kappa*X_ij*rho_ij 
    39                p_ik(l,ij) = preff*(X_ij**gamma) 
     38               X_ij = log(vreff*rho_ij) + theta(l,ij)/cpp 
     39               p_ik(l,ij) = preff*exp(X_ij*gamma) 
    4040               c2_mik = gamma*p_ik(l,ij)/(rho_ij*m_ik(l,ij)) ! c^2 = gamma*R*T = gamma*p/rho 
    4141               A_ik(l,ij) = c2_mik*(tau/g*rho_ij)**2 

  • codes/icosagcm/devel/src/kernels_unst/coriolis.k90

    r614 r658  
    77         ij_left = left(edge) 
    88         ij_right = right(edge) 
     9         !DIR$ SIMD 
    910         DO l = 1, llm 
    1011            Ftheta(l,edge) = .5*(theta(l,ij_left,iq)+theta(l,ij_right,iq))*hflux(l,edge) 
    … …  
    1415      !$OMP DO SCHEDULE(STATIC) 
    1516      DO ij = 1, primal_num 
    16          DO l = 1, llm 
    17             divF=0. 
    18             DO iedge = 1, primal_deg(ij) 
    19                edge = primal_edge(iedge,ij) 
    20                divF = divF + Ftheta(l,edge)*primal_ne(iedge,ij) 
    21             END DO 
    22             dtheta_rhodz(l,ij,iq) = -divF / Ai(ij) 
    23          END DO 
     17         ! this VLOOP iterates over primal cell edges 
     18         SELECT CASE(primal_deg(ij)) 
     19         CASE(4) 
     20            edge1 = primal_edge(1,ij) 
     21            edge2 = primal_edge(2,ij) 
     22            edge3 = primal_edge(3,ij) 
     23            edge4 = primal_edge(4,ij) 
     24            sign1 = primal_ne(1,ij) 
     25            sign2 = primal_ne(2,ij) 
     26            sign3 = primal_ne(3,ij) 
     27            sign4 = primal_ne(4,ij) 
     28            !DIR$ SIMD 
     29            DO l = 1, llm 
     30               divF=0. 
     31               divF = divF + Ftheta(l,edge1)*sign1 
     32               divF = divF + Ftheta(l,edge2)*sign2 
     33               divF = divF + Ftheta(l,edge3)*sign3 
     34               divF = divF + Ftheta(l,edge4)*sign4 
     35               dtheta_rhodz(l,ij,iq) = -divF / Ai(ij) 
     36            END DO 
     37         CASE(5) 
     38            edge1 = primal_edge(1,ij) 
     39            edge2 = primal_edge(2,ij) 
     40            edge3 = primal_edge(3,ij) 
     41            edge4 = primal_edge(4,ij) 
     42            edge5 = primal_edge(5,ij) 
     43            sign1 = primal_ne(1,ij) 
     44            sign2 = primal_ne(2,ij) 
     45            sign3 = primal_ne(3,ij) 
     46            sign4 = primal_ne(4,ij) 
     47            sign5 = primal_ne(5,ij) 
     48            !DIR$ SIMD 
     49            DO l = 1, llm 
     50               divF=0. 
     51               divF = divF + Ftheta(l,edge1)*sign1 
     52               divF = divF + Ftheta(l,edge2)*sign2 
     53               divF = divF + Ftheta(l,edge3)*sign3 
     54               divF = divF + Ftheta(l,edge4)*sign4 
     55               divF = divF + Ftheta(l,edge5)*sign5 
     56               dtheta_rhodz(l,ij,iq) = -divF / Ai(ij) 
     57            END DO 
     58         CASE(6) 
     59            edge1 = primal_edge(1,ij) 
     60            edge2 = primal_edge(2,ij) 
     61            edge3 = primal_edge(3,ij) 
     62            edge4 = primal_edge(4,ij) 
     63            edge5 = primal_edge(5,ij) 
     64            edge6 = primal_edge(6,ij) 
     65            sign1 = primal_ne(1,ij) 
     66            sign2 = primal_ne(2,ij) 
     67            sign3 = primal_ne(3,ij) 
     68            sign4 = primal_ne(4,ij) 
     69            sign5 = primal_ne(5,ij) 
     70            sign6 = primal_ne(6,ij) 
     71            !DIR$ SIMD 
     72            DO l = 1, llm 
     73               divF=0. 
     74               divF = divF + Ftheta(l,edge1)*sign1 
     75               divF = divF + Ftheta(l,edge2)*sign2 
     76               divF = divF + Ftheta(l,edge3)*sign3 
     77               divF = divF + Ftheta(l,edge4)*sign4 
     78               divF = divF + Ftheta(l,edge5)*sign5 
     79               divF = divF + Ftheta(l,edge6)*sign6 
     80               dtheta_rhodz(l,ij,iq) = -divF / Ai(ij) 
     81            END DO 
     82         CASE DEFAULT 
     83            !DIR$ SIMD 
     84            DO l = 1, llm 
     85               divF=0. 
     86               DO iedge = 1, primal_deg(ij) 
     87                  edge = primal_edge(iedge,ij) 
     88                  divF = divF + Ftheta(l,edge)*primal_ne(iedge,ij) 
     89               END DO 
     90               dtheta_rhodz(l,ij,iq) = -divF / Ai(ij) 
     91            END DO 
     92         END SELECT 
    2493      END DO 
    2594      !$OMP END DO 
    … …  
    2897   !$OMP DO SCHEDULE(STATIC) 
    2998   DO ij = 1, primal_num 
    30       DO l = 1, llm 
    31          divF=0. 
    32          DO iedge = 1, primal_deg(ij) 
    33             edge = primal_edge(iedge,ij) 
    34             divF = divF + hflux(l,edge)*primal_ne(iedge,ij) 
    35          END DO 
    36          convm(l,ij) = -divF / Ai(ij) 
    37       END DO 
     99      ! this VLOOP iterates over primal cell edges 
     100      SELECT CASE(primal_deg(ij)) 
     101      CASE(4) 
     102         edge1 = primal_edge(1,ij) 
     103         edge2 = primal_edge(2,ij) 
     104         edge3 = primal_edge(3,ij) 
     105         edge4 = primal_edge(4,ij) 
     106         sign1 = primal_ne(1,ij) 
     107         sign2 = primal_ne(2,ij) 
     108         sign3 = primal_ne(3,ij) 
     109         sign4 = primal_ne(4,ij) 
     110         !DIR$ SIMD 
     111         DO l = 1, llm 
     112            divF=0. 
     113            divF = divF + hflux(l,edge1)*sign1 
     114            divF = divF + hflux(l,edge2)*sign2 
     115            divF = divF + hflux(l,edge3)*sign3 
     116            divF = divF + hflux(l,edge4)*sign4 
     117            convm(l,ij) = -divF / Ai(ij) 
     118         END DO 
     119      CASE(5) 
     120         edge1 = primal_edge(1,ij) 
     121         edge2 = primal_edge(2,ij) 
     122         edge3 = primal_edge(3,ij) 
     123         edge4 = primal_edge(4,ij) 
     124         edge5 = primal_edge(5,ij) 
     125         sign1 = primal_ne(1,ij) 
     126         sign2 = primal_ne(2,ij) 
     127         sign3 = primal_ne(3,ij) 
     128         sign4 = primal_ne(4,ij) 
     129         sign5 = primal_ne(5,ij) 
     130         !DIR$ SIMD 
     131         DO l = 1, llm 
     132            divF=0. 
     133            divF = divF + hflux(l,edge1)*sign1 
     134            divF = divF + hflux(l,edge2)*sign2 
     135            divF = divF + hflux(l,edge3)*sign3 
     136            divF = divF + hflux(l,edge4)*sign4 
     137            divF = divF + hflux(l,edge5)*sign5 
     138            convm(l,ij) = -divF / Ai(ij) 
     139         END DO 
     140      CASE(6) 
     141         edge1 = primal_edge(1,ij) 
     142         edge2 = primal_edge(2,ij) 
     143         edge3 = primal_edge(3,ij) 
     144         edge4 = primal_edge(4,ij) 
     145         edge5 = primal_edge(5,ij) 
     146         edge6 = primal_edge(6,ij) 
     147         sign1 = primal_ne(1,ij) 
     148         sign2 = primal_ne(2,ij) 
     149         sign3 = primal_ne(3,ij) 
     150         sign4 = primal_ne(4,ij) 
     151         sign5 = primal_ne(5,ij) 
     152         sign6 = primal_ne(6,ij) 
     153         !DIR$ SIMD 
     154         DO l = 1, llm 
     155            divF=0. 
     156            divF = divF + hflux(l,edge1)*sign1 
     157            divF = divF + hflux(l,edge2)*sign2 
     158            divF = divF + hflux(l,edge3)*sign3 
     159            divF = divF + hflux(l,edge4)*sign4 
     160            divF = divF + hflux(l,edge5)*sign5 
     161            divF = divF + hflux(l,edge6)*sign6 
     162            convm(l,ij) = -divF / Ai(ij) 
     163         END DO 
     164      CASE DEFAULT 
     165         !DIR$ SIMD 
     166         DO l = 1, llm 
     167            divF=0. 
     168            DO iedge = 1, primal_deg(ij) 
     169               edge = primal_edge(iedge,ij) 
     170               divF = divF + hflux(l,edge)*primal_ne(iedge,ij) 
     171            END DO 
     172            convm(l,ij) = -divF / Ai(ij) 
     173         END DO 
     174      END SELECT 
    38175   END DO 
    39176   !$OMP END DO 
    … …  
    41178   !$OMP DO SCHEDULE(STATIC) 
    42179   DO edge = 1, edge_num 
    43       DO l = 1, llm 
    44          du_trisk=0. 
    45          DO itrisk = 1, trisk_deg(edge) 
    46             edge_trisk = trisk(itrisk,edge) 
    47             du_trisk = du_trisk + wee(itrisk,edge)*hflux(l,edge_trisk)*(qu(l,edge)+qu(l,edge_trisk)) 
    48          END DO 
    49          du(l,edge) = du(l,edge) + .5*du_trisk 
    50       END DO 
     180      ! this VLOOP iterates over the TRISK stencil 
     181      SELECT CASE(trisk_deg(edge)) 
     182      CASE(10) 
     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            itrisk = 5 
     199            edge_trisk = trisk(5,edge) 
     200            du_trisk = du_trisk + wee(itrisk,edge)*hflux(l,edge_trisk)*(qu(l,edge)+qu(l,edge_trisk)) 
     201            itrisk = 6 
     202            edge_trisk = trisk(6,edge) 
     203            du_trisk = du_trisk + wee(itrisk,edge)*hflux(l,edge_trisk)*(qu(l,edge)+qu(l,edge_trisk)) 
     204            itrisk = 7 
     205            edge_trisk = trisk(7,edge) 
     206            du_trisk = du_trisk + wee(itrisk,edge)*hflux(l,edge_trisk)*(qu(l,edge)+qu(l,edge_trisk)) 
     207            itrisk = 8 
     208            edge_trisk = trisk(8,edge) 
     209            du_trisk = du_trisk + wee(itrisk,edge)*hflux(l,edge_trisk)*(qu(l,edge)+qu(l,edge_trisk)) 
     210            itrisk = 9 
     211            edge_trisk = trisk(9,edge) 
     212            du_trisk = du_trisk + wee(itrisk,edge)*hflux(l,edge_trisk)*(qu(l,edge)+qu(l,edge_trisk)) 
     213            itrisk = 10 
     214            edge_trisk = trisk(10,edge) 
     215            du_trisk = du_trisk + wee(itrisk,edge)*hflux(l,edge_trisk)*(qu(l,edge)+qu(l,edge_trisk)) 
     216            du(l,edge) = du(l,edge) + .5*du_trisk 
     217         END DO 
     218      CASE(4) 
     219         !DIR$ SIMD 
     220         DO l = 1, llm 
     221            du_trisk=0. 
     222            itrisk = 1 
     223            edge_trisk = trisk(1,edge) 
     224            du_trisk = du_trisk + wee(itrisk,edge)*hflux(l,edge_trisk)*(qu(l,edge)+qu(l,edge_trisk)) 
     225            itrisk = 2 
     226            edge_trisk = trisk(2,edge) 
     227            du_trisk = du_trisk + wee(itrisk,edge)*hflux(l,edge_trisk)*(qu(l,edge)+qu(l,edge_trisk)) 
     228            itrisk = 3 
     229            edge_trisk = trisk(3,edge) 
     230            du_trisk = du_trisk + wee(itrisk,edge)*hflux(l,edge_trisk)*(qu(l,edge)+qu(l,edge_trisk)) 
     231            itrisk = 4 
     232            edge_trisk = trisk(4,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 
    51247   END DO 
    52248   !$OMP END DO 

  • codes/icosagcm/devel/src/kernels_unst/div.k90

    r614 r658  
    33   !$OMP DO SCHEDULE(STATIC) 
    44   DO ij = 1, primal_num 
    5       DO l = 1, llm 
    6          div_ij=0. 
    7          DO iedge = 1, primal_deg(ij) 
    8             edge = primal_edge(iedge,ij) 
    9             div_ij = div_ij + primal_ne(iedge,ij)*le_de(edge)*u(l,edge) 
     5      ! this VLOOP iterates over primal cell edges 
     6      SELECT CASE(primal_deg(ij)) 
     7      CASE(4) 
     8         edge1 = primal_edge(1,ij) 
     9         edge2 = primal_edge(2,ij) 
     10         edge3 = primal_edge(3,ij) 
     11         edge4 = primal_edge(4,ij) 
     12         le_de1 = le_de(edge1) 
     13         le_de2 = le_de(edge2) 
     14         le_de3 = le_de(edge3) 
     15         le_de4 = le_de(edge4) 
     16         sign1 = primal_ne(1,ij) 
     17         sign2 = primal_ne(2,ij) 
     18         sign3 = primal_ne(3,ij) 
     19         sign4 = primal_ne(4,ij) 
     20         !DIR$ SIMD 
     21         DO l = 1, llm 
     22            div_ij=0. 
     23            div_ij = div_ij + sign1*le_de1*u(l,edge1) 
     24            div_ij = div_ij + sign2*le_de2*u(l,edge2) 
     25            div_ij = div_ij + sign3*le_de3*u(l,edge3) 
     26            div_ij = div_ij + sign4*le_de4*u(l,edge4) 
     27            divu(l,ij) = div_ij / Ai(ij) 
    1028         END DO 
    11          divu(l,ij) = div_ij / Ai(ij) 
    12       END DO 
     29      CASE(5) 
     30         edge1 = primal_edge(1,ij) 
     31         edge2 = primal_edge(2,ij) 
     32         edge3 = primal_edge(3,ij) 
     33         edge4 = primal_edge(4,ij) 
     34         edge5 = primal_edge(5,ij) 
     35         le_de1 = le_de(edge1) 
     36         le_de2 = le_de(edge2) 
     37         le_de3 = le_de(edge3) 
     38         le_de4 = le_de(edge4) 
     39         le_de5 = le_de(edge5) 
     40         sign1 = primal_ne(1,ij) 
     41         sign2 = primal_ne(2,ij) 
     42         sign3 = primal_ne(3,ij) 
     43         sign4 = primal_ne(4,ij) 
     44         sign5 = primal_ne(5,ij) 
     45         !DIR$ SIMD 
     46         DO l = 1, llm 
     47            div_ij=0. 
     48            div_ij = div_ij + sign1*le_de1*u(l,edge1) 
     49            div_ij = div_ij + sign2*le_de2*u(l,edge2) 
     50            div_ij = div_ij + sign3*le_de3*u(l,edge3) 
     51            div_ij = div_ij + sign4*le_de4*u(l,edge4) 
     52            div_ij = div_ij + sign5*le_de5*u(l,edge5) 
     53            divu(l,ij) = div_ij / Ai(ij) 
     54         END DO 
     55      CASE(6) 
     56         edge1 = primal_edge(1,ij) 
     57         edge2 = primal_edge(2,ij) 
     58         edge3 = primal_edge(3,ij) 
     59         edge4 = primal_edge(4,ij) 
     60         edge5 = primal_edge(5,ij) 
     61         edge6 = primal_edge(6,ij) 
     62         le_de1 = le_de(edge1) 
     63         le_de2 = le_de(edge2) 
     64         le_de3 = le_de(edge3) 
     65         le_de4 = le_de(edge4) 
     66         le_de5 = le_de(edge5) 
     67         le_de6 = le_de(edge6) 
     68         sign1 = primal_ne(1,ij) 
     69         sign2 = primal_ne(2,ij) 
     70         sign3 = primal_ne(3,ij) 
     71         sign4 = primal_ne(4,ij) 
     72         sign5 = primal_ne(5,ij) 
     73         sign6 = primal_ne(6,ij) 
     74         !DIR$ SIMD 
     75         DO l = 1, llm 
     76            div_ij=0. 
     77            div_ij = div_ij + sign1*le_de1*u(l,edge1) 
     78            div_ij = div_ij + sign2*le_de2*u(l,edge2) 
     79            div_ij = div_ij + sign3*le_de3*u(l,edge3) 
     80            div_ij = div_ij + sign4*le_de4*u(l,edge4) 
     81            div_ij = div_ij + sign5*le_de5*u(l,edge5) 
     82            div_ij = div_ij + sign6*le_de6*u(l,edge6) 
     83            divu(l,ij) = div_ij / Ai(ij) 
     84         END DO 
     85      CASE DEFAULT 
     86         !DIR$ SIMD 
     87         DO l = 1, llm 
     88            div_ij=0. 
     89            DO iedge = 1, primal_deg(ij) 
     90               edge = primal_edge(iedge,ij) 
     91               div_ij = div_ij + primal_ne(iedge,ij)*le_de(edge)*u(l,edge) 
     92            END DO 
     93            divu(l,ij) = div_ij / Ai(ij) 
     94         END DO 
     95      END SELECT 
    1396   END DO 
    1497   !$OMP END DO 

  • codes/icosagcm/devel/src/kernels_unst/gradient.k90

    r614 r658  
    55      ij_left = left(edge) 
    66      ij_right = right(edge) 
     7      !DIR$ SIMD 
    78      DO l = 1, llm 
    89         grad(l,edge) = 1.*(b(l,ij_right)-b(l,ij_left)) 

  • codes/icosagcm/devel/src/kernels_unst/pvort_only.k90

    r614 r658  
    33   !$OMP DO SCHEDULE(STATIC) 
    44   DO ij = 1, dual_num 
    5       DO l = 1, llm 
    6          etav = 0.d0 
    7          DO iedge = 1, dual_deg(ij) 
    8             edge = dual_edge(iedge,ij) 
    9             etav = etav + dual_ne(iedge,ij)*u(l,edge) 
     5      ! this VLOOP iterates over dual cell edges 
     6      SELECT CASE(dual_deg(ij)) 
     7      CASE(3) 
     8         edge1 = dual_edge(1,ij) 
     9         edge2 = dual_edge(2,ij) 
     10         edge3 = dual_edge(3,ij) 
     11         sign1 = dual_ne(1,ij) 
     12         sign2 = dual_ne(2,ij) 
     13         sign3 = dual_ne(3,ij) 
     14         vertex1 = dual_vertex(1,ij) 
     15         vertex2 = dual_vertex(2,ij) 
     16         vertex3 = dual_vertex(3,ij) 
     17         !DIR$ SIMD 
     18         DO l = 1, llm 
     19            etav = 0.d0 
     20            etav = etav + sign1*u(l,edge1) 
     21            etav = etav + sign2*u(l,edge2) 
     22            etav = etav + sign3*u(l,edge3) 
     23            hv=0. 
     24            hv = hv + Riv2(1,ij)*rhodz(l,vertex1) 
     25            hv = hv + Riv2(2,ij)*rhodz(l,vertex2) 
     26            hv = hv + Riv2(3,ij)*rhodz(l,vertex3) 
     27            qv(l,ij) = (etav + fv(ij)*Av(ij) )/(hv*Av(ij)) 
    1028         END DO 
    11          hv=0. 
    12          DO ivertex = 1, dual_deg(ij) 
    13             vertex = dual_vertex(ivertex,ij) 
    14             hv = hv + Riv2(ivertex,ij)*rhodz(l,vertex) 
     29      CASE(4) 
     30         edge1 = dual_edge(1,ij) 
     31         edge2 = dual_edge(2,ij) 
     32         edge3 = dual_edge(3,ij) 
     33         edge4 = dual_edge(4,ij) 
     34         sign1 = dual_ne(1,ij) 
     35         sign2 = dual_ne(2,ij) 
     36         sign3 = dual_ne(3,ij) 
     37         sign4 = dual_ne(4,ij) 
     38         vertex1 = dual_vertex(1,ij) 
     39         vertex2 = dual_vertex(2,ij) 
     40         vertex3 = dual_vertex(3,ij) 
     41         vertex4 = dual_vertex(4,ij) 
     42         !DIR$ SIMD 
     43         DO l = 1, llm 
     44            etav = 0.d0 
     45            etav = etav + sign1*u(l,edge1) 
     46            etav = etav + sign2*u(l,edge2) 
     47            etav = etav + sign3*u(l,edge3) 
     48            etav = etav + sign4*u(l,edge4) 
     49            hv=0. 
     50            hv = hv + Riv2(1,ij)*rhodz(l,vertex1) 
     51            hv = hv + Riv2(2,ij)*rhodz(l,vertex2) 
     52            hv = hv + Riv2(3,ij)*rhodz(l,vertex3) 
     53            hv = hv + Riv2(4,ij)*rhodz(l,vertex4) 
     54            qv(l,ij) = (etav + fv(ij)*Av(ij) )/(hv*Av(ij)) 
    1555         END DO 
    16          qv(l,ij) = (etav + fv(ij)*Av(ij) )/(hv*Av(ij)) 
    17       END DO 
     56      CASE DEFAULT 
     57         !DIR$ SIMD 
     58         DO l = 1, llm 
     59            etav = 0.d0 
     60            DO iedge = 1, dual_deg(ij) 
     61               edge = dual_edge(iedge,ij) 
     62               etav = etav + dual_ne(iedge,ij)*u(l,edge) 
     63            END DO 
     64            hv=0. 
     65            DO ivertex = 1, dual_deg(ij) 
     66               vertex = dual_vertex(ivertex,ij) 
     67               hv = hv + Riv2(ivertex,ij)*rhodz(l,vertex) 
     68            END DO 
     69            qv(l,ij) = (etav + fv(ij)*Av(ij) )/(hv*Av(ij)) 
     70         END DO 
     71      END SELECT 
    1872   END DO 
    1973   !$OMP END DO 
    … …  
    2276      ij_up = up(edge) 
    2377      ij_down = down(edge) 
     78      !DIR$ SIMD 
    2479      DO l = 1, llm 
    2580         qu(l,edge)=0.5d0*(qv(l,ij_down)+qv(l,ij_up)) 

  • codes/icosagcm/devel/src/kernels_unst/theta.k90

    r614 r658  
    1616      !$OMP DO SCHEDULE(STATIC) 
    1717      DO ij = 1, primal_num 
     18         !DIR$ SIMD 
    1819         DO l = 1, llm 
    1920            ! FIXME : formula below (used in DYNAMICO) is for dak, dbk based on pressure rather than mass 
    … …  
    2829      !$OMP DO SCHEDULE(STATIC) 
    2930      DO ij = 1, primal_num 
     31         !DIR$ SIMD 
    3032         DO l = 1, llm 
    3133            theta(l,ij,iq) = theta_rhodz(l,ij,iq)/rhodz(l,ij) 

  • codes/icosagcm/devel/src/unstructured/caldyn_unstructured.F90

    r650 r658  
    2121#define DECLARE_VERTICES INTEGER VERTICES 
    2222#define PHI_BOT(ij) Phi_bot 
    23 #define PHI_BOT_VAR 0. 
    2423 
    2524#define BINDC_(thename) BIND(C, name=#thename) 
    … …  
    4039#define HASNAN(field) (ANY(.NOT.ABS(field)<1e20)) 
    4140 
     41#define START_TRACE(id,nprimal,ndual,nedge) CALL enter_trace(id, 8*llm*((nprimal)*primal_num+(ndual)*dual_num+(nedge)*edge_num) ) 
     42#define STOP_TRACE CALL exit_trace() 
     43 
    4244!----------------------------- Non-Hydrostatic ----------------------------- 
    4345 
    44 SUBROUTINE compute_NH_geopot(tau,dummy, m_ik, m_il, theta, W_il, Phi_il) 
    45   FIELD_MASS   :: m_ik, theta, p_ik, A_ik, C_ik   ! IN*2,LOCAL*3 
    46   FIELD_GEOPOT :: m_il, W_il, Phi_il, Phi_star_il, R_il, x_il, B_il, D_il  ! IN,INOUT*2, LOCAL*5  
    47   DBL :: tau, dummy, gamma, rho_ij, X_ij, Y_ij, wil, tau2_g, g2, gm2, ml_g2, c2_mik 
    48   DECLARE_INDICES 
     46SUBROUTINE compute_NH_geopot(tau, m_ik, m_il, theta, W_il, Phi_il) 
     47  FIELD_MASS   :: m_ik, theta   ! IN*2 
     48  FIELD_GEOPOT :: m_il, W_il, Phi_il, Phi_star_il  ! IN,INOUT*2, LOCAL*5  
     49  DBL :: tau, gamma, tau2_g, tau2_g2, g2, gm2, vreff, Rd_preff 
    4950  INTEGER :: iter 
     51  DECLARE_INDICES 
     52  DBL :: rho_ij, X_ij, Y_ij, wil, rho_c2_mik, c2_mik, ml_g2 
     53#define COLUMN 0 
     54#if COLUMN  
     55  DOUBLE1(llm)  :: pk, Ak, Ck 
     56  DOUBLE1(llm+1):: Rl, Bl, Dl, xl 
     57#define p_ik(l,ij) pk(l) 
     58#define A_ik(l,ij) Ak(l) 
     59#define C_ik(l,ij) Ck(l) 
     60#define R_il(l,ij) Rl(l) 
     61#define B_il(l,ij) Bl(l) 
     62#define D_il(l,ij) Dl(l) 
     63#define x_il(l,ij) xl(l) 
     64#else 
     65  FIELD_MASS :: p_ik, A_ik, C_ik 
     66  FIELD_GEOPOT :: R_il, B_il, D_il, x_il 
     67#endif 
     68 
     69  START_TRACE(id_NH_geopot, 7,0,0) 
    5070#include "../kernels_unst/compute_NH_geopot.k90" 
     71  STOP_TRACE   
     72 
     73#if COLUMN 
     74#undef p_ik 
     75#undef A_ik 
     76#undef C_ik 
     77#undef R_il 
     78#undef B_il 
     79#undef D_il 
     80#undef x_il 
     81#endif 
     82#undef COLUMN 
    5183END SUBROUTINE compute_NH_geopot 
    5284 
    … …  
    5991  DECLARE_EDGES 
    6092  DBL :: W_el, W2_el, gPhi2, dP, divG, u2, uu 
     93  START_TRACE(id_slow_NH, 5,0,3) 
    6194#include "../kernels_unst/caldyn_slow_NH.k90" 
     95  STOP_TRACE 
    6296END SUBROUTINE compute_caldyn_slow_NH 
    6397 
    … …  
    69103  FIELD_U      :: du                     ! OUT 
    70104  DECLARE_INDICES 
    71   DBL :: X_ij, rho_ij, T_ij, gamma, Cvd, vreff 
     105  DBL :: X_ij, rho_ij, T_ij, gamma, Cvd, vreff, Rd_preff 
     106#include "../kernels_unst/caldyn_mil.k90" 
     107  IF(tau>0) THEN ! solve implicit problem for geopotential 
     108    CALL compute_NH_geopot(tau, rhodz, m_il, theta, W, geopot) 
     109  END IF 
     110  START_TRACE(id_solver, 7,0,1) 
    72111#include "../kernels_unst/caldyn_solver.k90" 
     112  STOP_TRACE 
    73113END SUBROUTINE compute_caldyn_solver 
    74114 
    … …  
    79119  DECLARE_INDICES 
    80120  DBL :: w_ij, wflux_ij 
     121  START_TRACE(id_vert_NH, 6,0,1) 
    81122#include "../kernels_unst/caldyn_vert_NH.k90" 
     123  STOP_TRACE 
    82124END SUBROUTINE compute_caldyn_vert_NH 
    83125 
    … …  
    91133  DECLARE_INDICES 
    92134  DBL :: gdz, ke, uu, chi, gv, exner_ik, temp_ik, p_ik, qv, Rmix 
     135  START_TRACE(id_geopot, 3,0,3) 
    93136#include "../kernels_unst/compute_geopot.k90" 
     137  STOP_TRACE 
    94138END SUBROUTINE compute_geopot 
    95139! 
    … …  
    102146  LOGICAL, PARAMETER :: zero=.TRUE. 
    103147  DBL :: ke, uu 
     148  START_TRACE(id_slow_hydro, 3,0,3) 
    104149#include "../kernels_unst/caldyn_slow_hydro.k90" 
     150  STOP_TRACE 
    105151END SUBROUTINE compute_caldyn_slow_hydro 
    106152 
    … …  
    118164  wwuu=0. 
    119165  !$OMP BARRIER 
     166  START_TRACE(id_vert, 5,0,3) 
    120167#include "../kernels_unst/caldyn_wflux.k90" 
    121168#include "../kernels_unst/caldyn_dmass.k90" 
    122169#include "../kernels_unst/caldyn_vert.k90" 
     170  STOP_TRACE 
    123171END SUBROUTINE caldyn_vert 
    124172 
    … …  
    130178  DECLARE_EDGES 
    131179  DBL :: divF, du_trisk 
     180  START_TRACE(id_coriolis, 3,4,0) ! primal, dual, edge 
    132181#include "../kernels_unst/coriolis.k90" 
     182  STOP_TRACE 
    133183END SUBROUTINE 
    134184 
    … …  
    139189  DECLARE_INDICES 
    140190  DBL :: m 
     191  START_TRACE(id_theta, 3,0,0) ! primal, dual, edge 
    141192#include "../kernels_unst/theta.k90" 
     193  STOP_TRACE 
    142194END SUBROUTINE 
    143195 
    … …  
    150202  DECLARE_VERTICES 
    151203  DBL :: etav, hv 
     204  START_TRACE(id_pvort_only, 1,1,2) ! primal, dual, edge 
    152205#include "../kernels_unst/pvort_only.k90" 
     206  STOP_TRACE 
    153207END SUBROUTINE compute_pvort_only 
    154208 
    … …  
    162216  DECLARE_EDGES 
    163217  DBL          :: due 
    164  
     218  START_TRACE(id_fast, 4,0,2) ! primal, dual, edge 
    165219#include "../kernels_unst/caldyn_fast.k90" 
    166  
     220  STOP_TRACE 
    167221END SUBROUTINE compute_caldyn_fast 
    168222 

  • codes/icosagcm/devel/src/unstructured/data_unstructured.F90

    r651 r658  
    4444  INTEGER(C_INT), BIND(C) :: comm_icosa 
    4545 
    46   INTEGER, PARAMETER :: id_pvort_only=1, id_slow_hydro=2, id_fast=3, id_coriolis=4, id_theta=5, & 
    47        id_vert_NH=6, id_solver=7, id_slow_NH=8, id_NH_geopot=9, id_vert=10, id_NH_Phi_star=11, nb_routines=11 
     46  INTEGER, PARAMETER :: id_dev1=1, id_dev2=2, & 
     47       id_pvort_only=3, id_slow_hydro=4, id_fast=5, id_coriolis=6, id_theta=7, id_geopot=8, id_vert=9, & 
     48       id_solver=10, id_slow_NH=11, id_NH_geopot=12, id_vert_NH=13, id_update=14, nb_routines=14  
    4849  DBL, PRIVATE :: start_time, time_spent(nb_routines) ! time spent in each kernel 
    4950  INTEGER, PRIVATE :: current_id, nb_calls(nb_routines), bytes(nb_routines) ! bytes read or written be each kernel 
    5051  CHARACTER(len = 10) :: id_name(nb_routines) = & 
    51        (/'pvort_only', 'slow_hydro', 'fast      ', 'coriolis  ', 'theta     ', & 
    52        'vert_NH   ', 'solver    ', 'slow_NH   ', 'NH_geopot ', 'vert      ', & 
    53        'Phi_star  '/) 
     52       (/'dev1      ', 'dev2      ', & 
     53       'pvort_only', 'slow_hydro', 'fast      ', 'coriolis  ', 'theta     ', 'geopot    ', 'vert      ', & 
     54       'solver    ', 'slow_NH   ', 'NH_geopot ', 'vert_NH   ',  'update    ' /) 
    5455 
    5556CONTAINS 

  • codes/icosagcm/devel/src/unstructured/timestep_unstructured.F90

    r652 r658  
    219219    ! 
    220220    INTEGER :: l, ij 
     221    CALL enter_trace(id_update, 8*sz*(j+1) ) 
    221222    !  PRINT *, clj(1:j,j) 
    222223    SELECT CASE(j) 
    … …  
    263264       END DO 
    264265    END SELECT 
     266    CALL exit_trace() 
    265267  END SUBROUTINE update 
    266268   
Note: See TracChangeset for help on using the changeset viewer.