Changeset 174


Ignore:
Timestamp:
10/07/13 17:48:24 (11 years ago)
Author:
ymipsl
Message:

Transform 2 loops on i and j in one loop ij for efficiency, vectorization and future GPU programing

YM

Location:
codes/icosagcm/trunk/src
Files:
8 edited

Legend:

Unmodified
Added
Removed

  • codes/icosagcm/trunk/src/advect.f90

    r151 r174  
    1717    REAL(rstd),INTENT(OUT) :: one_over_sqrt_leng(iim*jjm) 
    1818 
    19     INTEGER :: i,j,n 
    20  
    21     DO j=jj_begin-1,jj_end+1 
    22        DO i=ii_begin-1,ii_end+1 
    23           n=(j-1)*iim+i 
    24  
    25           CALL cross_product2(xyz_v(n+z_rdown,:),xyz_v(n+z_rup,:),normal(n+u_right,:)) 
    26           normal(n+u_right,:)=normal(n+u_right,:)/sqrt(sum(normal(n+u_right,:)**2)+1e-50)*ne(n,right) 
    27  
    28           CALL cross_product2(xyz_v(n+z_up,:),xyz_v(n+z_lup,:),normal(n+u_lup,:)) 
    29           normal(n+u_lup,:)=normal(n+u_lup,:)/sqrt(sum(normal(n+u_lup,:)**2)+1e-50)*ne(n,lup) 
    30  
    31           CALL cross_product2(xyz_v(n+z_ldown,:),xyz_v(n+z_down,:),normal(n+u_ldown,:))   
    32           normal(n+u_ldown,:)=normal(n+u_ldown,:)/sqrt(sum(normal(n+u_ldown,:)**2)+1e-50)*ne(n,ldown) 
    33  
    34           tangent(n+u_right,:)=xyz_v(n+z_rup,:)-xyz_v(n+z_rdown,:)  
    35           tangent(n+u_right,:)=tangent(n+u_right,:)/sqrt(sum(tangent(n+u_right,:)**2)+1e-50) 
    36  
    37           tangent(n+u_lup,:)=xyz_v(n+z_lup,:)-xyz_v(n+z_up,:)  
    38           tangent(n+u_lup,:)=tangent(n+u_lup,:)/sqrt(sum(tangent(n+u_lup,:)**2)+1e-50) 
    39  
    40           tangent(n+u_ldown,:)=xyz_v(n+z_down,:)-xyz_v(n+z_ldown,:) 
    41           tangent(n+u_ldown,:)=tangent(n+u_ldown,:)/sqrt(sum(tangent(n+u_ldown,:)**2)+1e-50) 
    42  
    43           one_over_sqrt_leng(n) = 1./sqrt(max(sum((xyz_v(n+z_up,:) - xyz_i(n,:))**2),sum((xyz_v(n+z_down,:) - xyz_i(n,:))**2), & 
    44                                           sum((xyz_v(n+z_rup,:) - xyz_i(n,:))**2),sum((xyz_v(n+z_rdown,:) - xyz_i(n,:))**2),  & 
    45                                           sum((xyz_v(n+z_lup,:) - xyz_i(n,:))**2),sum((xyz_v(n+z_ldown,:) - xyz_i(n,:))**2)) ) 
    46        ENDDO 
     19    INTEGER :: ij 
     20 
     21!$SIMD 
     22      DO ij=ij_begin,ij_end 
     23 
     24          CALL cross_product2(xyz_v(ij+z_rdown,:),xyz_v(ij+z_rup,:),normal(ij+u_right,:)) 
     25          normal(ij+u_right,:)=normal(ij+u_right,:)/sqrt(sum(normal(ij+u_right,:)**2)+1e-50)*ne(ij,right) 
     26 
     27          CALL cross_product2(xyz_v(ij+z_up,:),xyz_v(ij+z_lup,:),normal(ij+u_lup,:)) 
     28          normal(ij+u_lup,:)=normal(ij+u_lup,:)/sqrt(sum(normal(ij+u_lup,:)**2)+1e-50)*ne(ij,lup) 
     29 
     30          CALL cross_product2(xyz_v(ij+z_ldown,:),xyz_v(ij+z_down,:),normal(ij+u_ldown,:))        
     31          normal(ij+u_ldown,:)=normal(ij+u_ldown,:)/sqrt(sum(normal(ij+u_ldown,:)**2)+1e-50)*ne(ij,ldown) 
     32 
     33          tangent(ij+u_right,:)=xyz_v(ij+z_rup,:)-xyz_v(ij+z_rdown,:)  
     34          tangent(ij+u_right,:)=tangent(ij+u_right,:)/sqrt(sum(tangent(ij+u_right,:)**2)+1e-50) 
     35 
     36          tangent(ij+u_lup,:)=xyz_v(ij+z_lup,:)-xyz_v(ij+z_up,:)  
     37          tangent(ij+u_lup,:)=tangent(ij+u_lup,:)/sqrt(sum(tangent(ij+u_lup,:)**2)+1e-50) 
     38 
     39          tangent(ij+u_ldown,:)=xyz_v(ij+z_down,:)-xyz_v(ij+z_ldown,:) 
     40          tangent(ij+u_ldown,:)=tangent(ij+u_ldown,:)/sqrt(sum(tangent(ij+u_ldown,:)**2)+1e-50) 
     41 
     42          one_over_sqrt_leng(ij) = 1./sqrt(max(sum((xyz_v(ij+z_up,:) - xyz_i(ij,:))**2),sum((xyz_v(ij+z_down,:) - xyz_i(ij,:))**2), & 
     43                                          sum((xyz_v(ij+z_rup,:) - xyz_i(ij,:))**2),sum((xyz_v(ij+z_rdown,:) - xyz_i(ij,:))**2),  & 
     44                                          sum((xyz_v(ij+z_lup,:) - xyz_i(ij,:))**2),sum((xyz_v(ij+z_ldown,:) - xyz_i(ij,:))**2)) ) 
    4745    ENDDO 
    4846 
    … …  
    6462    REAL(rstd) :: ar(iim*jjm) 
    6563    REAL(rstd) :: gradtri(2*iim*jjm,llm,3)  
    66     INTEGER :: i,j,n,k,ind,l 
     64    INTEGER :: ij,k,ind,l 
    6765 
    6866    CALL trace_start("compute_gradq3d") 
    … …  
    7472    ! arr = area of triangle joining centroids of hexagons 
    7573     DO l = ll_begin,ll_end  
    76        DO j=jj_begin-1,jj_end+1 
    77           DO i=ii_begin-1,ii_end+1 
    78              n=(j-1)*iim+i    
    79              CALL gradq(n,l,n+t_rup,n+t_lup,n+z_up,qi,gradtri(n+z_up,l,:),arr(n+z_up)) 
    80              CALL gradq(n,l,n+t_ldown,n+t_rdown,n+z_down,qi,gradtri(n+z_down,l,:),arr(n+z_down)) 
    81           END DO 
     74!$SIMD 
     75      DO ij=ij_begin_ext,ij_end_ext 
     76             CALL gradq(ij,l,ij+t_rup,ij+t_lup,ij+z_up,qi,gradtri(ij+z_up,l,:),arr(ij+z_up)) 
     77             CALL gradq(ij,l,ij+t_ldown,ij+t_rdown,ij+z_down,qi,gradtri(ij+z_down,l,:),arr(ij+z_down)) 
    8278       END DO 
    8379    END DO 
    8480 
    85     DO j=jj_begin-1,jj_end+1 
    86       DO i=ii_begin-1,ii_end+1 
    87          n=(j-1)*iim+i    
    88          ar(n) = arr(n+z_up)+arr(n+z_lup)+arr(n+z_ldown)+arr(n+z_down)+arr(n+z_rdown)+arr(n+z_rup)+1.e-50 
    89       ENDDO 
     81!$SIMD 
     82      DO ij=ij_begin,ij_end 
     83         ar(ij) = arr(ij+z_up)+arr(ij+z_lup)+arr(ij+z_ldown)+arr(ij+z_down)+arr(ij+z_rdown)+arr(ij+z_rup)+1.e-50 
    9084    ENDDO 
    9185       
    9286    DO k=1,3 
    9387      DO l =ll_begin,ll_end 
    94         DO j=jj_begin,jj_end 
    95            DO i=ii_begin,ii_end 
    96               n=(j-1)*iim+i 
    97               gradq3d(n,l,k) = ( gradtri(n+z_up,l,k) + gradtri(n+z_down,l,k)   +   & 
    98                                  gradtri(n+z_rup,l,k) +  gradtri(n+z_ldown,l,k)   +   &  
    99                                  gradtri(n+z_lup,l,k)+    gradtri(n+z_rdown,l,k) ) / ar(n)   
    100            END DO 
     88!$SIMD 
     89      DO ij=ij_begin,ij_end 
     90              gradq3d(ij,l,k) = ( gradtri(ij+z_up,l,k) + gradtri(ij+z_down,l,k)   +   & 
     91                                 gradtri(ij+z_rup,l,k) +  gradtri(ij+z_ldown,l,k)   +   &  
     92                                 gradtri(ij+z_lup,l,k)+    gradtri(ij+z_rdown,l,k) ) / ar(ij)   
    10193        END DO 
    10294      END DO 
    … …  
    10597    !============================================================================================= LIMITING  
    10698    DO l =ll_begin,ll_end 
    107        DO j=jj_begin,jj_end 
    108           DO i=ii_begin,ii_end 
    109              n=(j-1)*iim+i 
    110              maggrd =  dot_product(gradq3d(n,l,:),gradq3d(n,l,:)) 
     99!$SIMD 
     100      DO ij=ij_begin,ij_end 
     101             maggrd =  dot_product(gradq3d(ij,l,:),gradq3d(ij,l,:)) 
    111102             maggrd = sqrt(maggrd)  
    112              maxq_c = qi(n,l) + maggrd*one_over_sqrt_leng(n) 
    113              minq_c = qi(n,l) - maggrd*one_over_sqrt_leng(n) 
    114              maxq = max(qi(n,l),qi(n+t_right,l),qi(n+t_lup,l),qi(n+t_rup,l),qi(n+t_left,l), & 
    115                   qi(n+t_rdown,l),qi(n+t_ldown,l)) 
    116              minq = min(qi(n,l),qi(n+t_right,l),qi(n+t_lup,l),qi(n+t_rup,l),qi(n+t_left,l), & 
    117                   qi(n+t_rdown,l),qi(n+t_ldown,l)) 
    118              alphamx = (maxq - qi(n,l)) ; alphamx = alphamx/(maxq_c - qi(n,l) ) 
     103             maxq_c = qi(ij,l) + maggrd*one_over_sqrt_leng(ij) 
     104             minq_c = qi(ij,l) - maggrd*one_over_sqrt_leng(ij) 
     105             maxq = max(qi(ij,l),qi(ij+t_right,l),qi(ij+t_lup,l),qi(ij+t_rup,l),qi(ij+t_left,l), & 
     106                  qi(ij+t_rdown,l),qi(ij+t_ldown,l)) 
     107             minq = min(qi(ij,l),qi(ij+t_right,l),qi(ij+t_lup,l),qi(ij+t_rup,l),qi(ij+t_left,l), & 
     108                  qi(ij+t_rdown,l),qi(ij+t_ldown,l)) 
     109             alphamx = (maxq - qi(ij,l)) ; alphamx = alphamx/(maxq_c - qi(ij,l) ) 
    119110             alphamx = max(alphamx,0.0) 
    120              alphami = (minq - qi(n,l)); alphami = alphami/(minq_c - qi(n,l)) 
     111             alphami = (minq - qi(ij,l)); alphami = alphami/(minq_c - qi(ij,l)) 
    121112             alphami = max(alphami,0.0)  
    122113             alpha   = min(alphamx,alphami,1.0) 
    123              gradq3d(n,l,:) = alpha*gradq3d(n,l,:) 
    124           END DO 
     114             gradq3d(ij,l,:) = alpha*gradq3d(ij,l,:) 
    125115       END DO 
    126116    END DO 
    … …  
    141131 
    142132    REAL(rstd) :: v_e(3), up_e, qe, ed(3)     
    143     INTEGER :: i,j,n,l 
     133    INTEGER :: ij,l 
    144134 
    145135    CALL trace_start("compute_backward_traj") 
    … …  
    149139    ! reconstruct tangential wind then 3D wind at edge then cc = edge midpoint - u*tau 
    150140    DO l = ll_begin,ll_end 
    151        DO j=jj_begin-1,jj_end+1 
    152           DO i=ii_begin-1,ii_end+1 
    153              n=(j-1)*iim+i 
    154  
    155              up_e =1/de(n+u_right)*(                                                   & 
    156                   wee(n+u_right,1,1)*le(n+u_rup)*ue(n+u_rup,l)+                        & 
    157                   wee(n+u_right,2,1)*le(n+u_lup)*ue(n+u_lup,l)+                        & 
    158                   wee(n+u_right,3,1)*le(n+u_left)*ue(n+u_left,l)+                      & 
    159                   wee(n+u_right,4,1)*le(n+u_ldown)*ue(n+u_ldown,l)+                    & 
    160                   wee(n+u_right,5,1)*le(n+u_rdown)*ue(n+u_rdown,l)+                    &  
    161                   wee(n+u_right,1,2)*le(n+t_right+u_ldown)*ue(n+t_right+u_ldown,l)+    & 
    162                   wee(n+u_right,2,2)*le(n+t_right+u_rdown)*ue(n+t_right+u_rdown,l)+    & 
    163                   wee(n+u_right,3,2)*le(n+t_right+u_right)*ue(n+t_right+u_right,l)+    & 
    164                   wee(n+u_right,4,2)*le(n+t_right+u_rup)*ue(n+t_right+u_rup,l)+        & 
    165                   wee(n+u_right,5,2)*le(n+t_right+u_lup)*ue(n+t_right+u_lup,l)         & 
     141!$SIMD 
     142      DO ij=ij_begin,ij_end 
     143             up_e =1/de(ij+u_right)*(                                                   & 
     144                  wee(ij+u_right,1,1)*le(ij+u_rup)*ue(ij+u_rup,l)+                        & 
     145                  wee(ij+u_right,2,1)*le(ij+u_lup)*ue(ij+u_lup,l)+                        & 
     146                  wee(ij+u_right,3,1)*le(ij+u_left)*ue(ij+u_left,l)+                      & 
     147                  wee(ij+u_right,4,1)*le(ij+u_ldown)*ue(ij+u_ldown,l)+                    & 
     148                  wee(ij+u_right,5,1)*le(ij+u_rdown)*ue(ij+u_rdown,l)+                    &  
     149                  wee(ij+u_right,1,2)*le(ij+t_right+u_ldown)*ue(ij+t_right+u_ldown,l)+    & 
     150                  wee(ij+u_right,2,2)*le(ij+t_right+u_rdown)*ue(ij+t_right+u_rdown,l)+    & 
     151                  wee(ij+u_right,3,2)*le(ij+t_right+u_right)*ue(ij+t_right+u_right,l)+    & 
     152                  wee(ij+u_right,4,2)*le(ij+t_right+u_rup)*ue(ij+t_right+u_rup,l)+        & 
     153                  wee(ij+u_right,5,2)*le(ij+t_right+u_lup)*ue(ij+t_right+u_lup,l)         & 
    166154                  ) 
    167              v_e = ue(n+u_right,l)*normal(n+u_right,:) + up_e*tangent(n+u_right,:) 
    168              cc(n+u_right,l,:) = xyz_e(n+u_right,:) - v_e*tau 
    169  
    170              up_e=1/de(n+u_lup)*(                                                      & 
    171                   wee(n+u_lup,1,1)*le(n+u_left)*ue(n+u_left,l)+                        & 
    172                   wee(n+u_lup,2,1)*le(n+u_ldown)*ue(n+u_ldown,l)+                      & 
    173                   wee(n+u_lup,3,1)*le(n+u_rdown)*ue(n+u_rdown,l)+                      & 
    174                   wee(n+u_lup,4,1)*le(n+u_right)*ue(n+u_right,l)+                      & 
    175                   wee(n+u_lup,5,1)*le(n+u_rup)*ue(n+u_rup,l)+                          &  
    176                   wee(n+u_lup,1,2)*le(n+t_lup+u_right)*ue(n+t_lup+u_right,l)+          & 
    177                   wee(n+u_lup,2,2)*le(n+t_lup+u_rup)*ue(n+t_lup+u_rup,l)+              & 
    178                   wee(n+u_lup,3,2)*le(n+t_lup+u_lup)*ue(n+t_lup+u_lup,l)+              & 
    179                   wee(n+u_lup,4,2)*le(n+t_lup+u_left)*ue(n+t_lup+u_left,l)+            & 
    180                   wee(n+u_lup,5,2)*le(n+t_lup+u_ldown)*ue(n+t_lup+u_ldown,l)           & 
     155             v_e = ue(ij+u_right,l)*normal(ij+u_right,:) + up_e*tangent(ij+u_right,:) 
     156             cc(ij+u_right,l,:) = xyz_e(ij+u_right,:) - v_e*tau 
     157 
     158             up_e=1/de(ij+u_lup)*(                                                      & 
     159                  wee(ij+u_lup,1,1)*le(ij+u_left)*ue(ij+u_left,l)+                        & 
     160                  wee(ij+u_lup,2,1)*le(ij+u_ldown)*ue(ij+u_ldown,l)+                      & 
     161                  wee(ij+u_lup,3,1)*le(ij+u_rdown)*ue(ij+u_rdown,l)+                      & 
     162                  wee(ij+u_lup,4,1)*le(ij+u_right)*ue(ij+u_right,l)+                      & 
     163                  wee(ij+u_lup,5,1)*le(ij+u_rup)*ue(ij+u_rup,l)+                          &  
     164                  wee(ij+u_lup,1,2)*le(ij+t_lup+u_right)*ue(ij+t_lup+u_right,l)+          & 
     165                  wee(ij+u_lup,2,2)*le(ij+t_lup+u_rup)*ue(ij+t_lup+u_rup,l)+              & 
     166                  wee(ij+u_lup,3,2)*le(ij+t_lup+u_lup)*ue(ij+t_lup+u_lup,l)+              & 
     167                  wee(ij+u_lup,4,2)*le(ij+t_lup+u_left)*ue(ij+t_lup+u_left,l)+            & 
     168                  wee(ij+u_lup,5,2)*le(ij+t_lup+u_ldown)*ue(ij+t_lup+u_ldown,l)           & 
    181169                  ) 
    182              v_e = ue(n+u_lup,l)*normal(n+u_lup,:) + up_e*tangent(n+u_lup,:) 
    183              cc(n+u_lup,l,:) = xyz_e(n+u_lup,:) - v_e*tau 
    184  
    185              up_e=1/de(n+u_ldown)*(                                                    & 
    186                   wee(n+u_ldown,1,1)*le(n+u_rdown)*ue(n+u_rdown,l)+                    & 
    187                   wee(n+u_ldown,2,1)*le(n+u_right)*ue(n+u_right,l)+                    & 
    188                   wee(n+u_ldown,3,1)*le(n+u_rup)*ue(n+u_rup,l)+                        & 
    189                   wee(n+u_ldown,4,1)*le(n+u_lup)*ue(n+u_lup,l)+                        & 
    190                   wee(n+u_ldown,5,1)*le(n+u_left)*ue(n+u_left,l)+                      &  
    191                   wee(n+u_ldown,1,2)*le(n+t_ldown+u_lup)*ue(n+t_ldown+u_lup,l)+        & 
    192                   wee(n+u_ldown,2,2)*le(n+t_ldown+u_left)*ue(n+t_ldown+u_left,l)+      & 
    193                   wee(n+u_ldown,3,2)*le(n+t_ldown+u_ldown)*ue(n+t_ldown+u_ldown,l)+    & 
    194                   wee(n+u_ldown,4,2)*le(n+t_ldown+u_rdown)*ue(n+t_ldown+u_rdown,l)+    & 
    195                   wee(n+u_ldown,5,2)*le(n+t_ldown+u_right)*ue(n+t_ldown+u_right,l)     & 
     170             v_e = ue(ij+u_lup,l)*normal(ij+u_lup,:) + up_e*tangent(ij+u_lup,:) 
     171             cc(ij+u_lup,l,:) = xyz_e(ij+u_lup,:) - v_e*tau 
     172             
     173 
     174             up_e=1/de(ij+u_ldown)*(                                                    & 
     175                  wee(ij+u_ldown,1,1)*le(ij+u_rdown)*ue(ij+u_rdown,l)+                    & 
     176                  wee(ij+u_ldown,2,1)*le(ij+u_right)*ue(ij+u_right,l)+                    & 
     177                  wee(ij+u_ldown,3,1)*le(ij+u_rup)*ue(ij+u_rup,l)+                        & 
     178                  wee(ij+u_ldown,4,1)*le(ij+u_lup)*ue(ij+u_lup,l)+                        & 
     179                  wee(ij+u_ldown,5,1)*le(ij+u_left)*ue(ij+u_left,l)+                      &  
     180                  wee(ij+u_ldown,1,2)*le(ij+t_ldown+u_lup)*ue(ij+t_ldown+u_lup,l)+        & 
     181                  wee(ij+u_ldown,2,2)*le(ij+t_ldown+u_left)*ue(ij+t_ldown+u_left,l)+      & 
     182                  wee(ij+u_ldown,3,2)*le(ij+t_ldown+u_ldown)*ue(ij+t_ldown+u_ldown,l)+    & 
     183                  wee(ij+u_ldown,4,2)*le(ij+t_ldown+u_rdown)*ue(ij+t_ldown+u_rdown,l)+    & 
     184                  wee(ij+u_ldown,5,2)*le(ij+t_ldown+u_right)*ue(ij+t_ldown+u_right,l)     & 
    196185                  ) 
    197              v_e = ue(n+u_ldown,l)*normal(n+u_ldown,:) + up_e*tangent(n+u_ldown,:) 
    198              cc(n+u_ldown,l,:) = xyz_e(n+u_ldown,:) - v_e*tau 
    199           ENDDO 
     186             v_e = ue(ij+u_ldown,l)*normal(ij+u_ldown,:) + up_e*tangent(ij+u_ldown,:) 
     187             cc(ij+u_ldown,l,:) = xyz_e(ij+u_ldown,:) - v_e*tau 
    200188       ENDDO 
    201189    END DO 
    … …  
    220208    REAL(rstd) :: dq,dmass,qe,ed(3), newmass 
    221209    REAL(rstd) :: qflux(3*iim*jjm,llm) 
    222     INTEGER :: i,j,n,k,l 
     210    INTEGER :: ij,k,l 
    223211 
    224212    CALL trace_start("compute_advect_horiz") 
    … …  
    228216    ! ne*hfluxt>0 iff outgoing 
    229217    DO l = ll_begin,ll_end 
    230        DO j=jj_begin-1,jj_end+1 
    231           DO i=ii_begin-1,ii_end+1 
    232              n=(j-1)*iim+i 
    233              if (ne(n,right)*hfluxt(n+u_right,l)>0) then  
    234                 ed = cc(n+u_right,l,:) - xyz_i(n,:) 
    235                 qe = qi(n,l)+sum2(gradq3d(n,l,:),ed)  
    236              else 
    237                 ed = cc(n+u_right,l,:) - xyz_i(n+t_right,:) 
    238                 qe = qi(n+t_right,l)+sum2(gradq3d(n+t_right,l,:),ed)  
    239              endif 
    240              qflux(n+u_right,l) = hfluxt(n+u_right,l)*qe 
     218 
     219!$SIMD 
     220      DO ij=ij_begin_ext,ij_end_ext 
     221 
     222             IF (ne(ij,right)*hfluxt(ij+u_right,l)>0) THEN  
     223                ed = cc(ij+u_right,l,:) - xyz_i(ij,:) 
     224                qe = qi(ij,l)+sum2(gradq3d(ij,l,:),ed)  
     225             ELSE 
     226                ed = cc(ij+u_right,l,:) - xyz_i(ij+t_right,:) 
     227                qe = qi(ij+t_right,l)+sum2(gradq3d(ij+t_right,l,:),ed)  
     228             ENDIF 
     229             qflux(ij+u_right,l) = hfluxt(ij+u_right,l)*qe 
    241230              
    242              if (ne(n,lup)*hfluxt(n+u_lup,l)>0) then  
    243                 ed = cc(n+u_lup,l,:) - xyz_i(n,:) 
    244                 qe = qi(n,l)+sum2(gradq3d(n,l,:),ed) 
    245              else 
    246                 ed = cc(n+u_lup,l,:) - xyz_i(n+t_lup,:) 
    247                 qe = qi(n+t_lup,l)+sum2(gradq3d(n+t_lup,l,:),ed)  
    248              endif 
    249              qflux(n+u_lup,l) = hfluxt(n+u_lup,l)*qe  
    250               
    251              if (ne(n,ldown)*hfluxt(n+u_ldown,l)>0) then  
    252                 ed = cc(n+u_ldown,l,:) - xyz_i(n,:) 
    253                 qe = qi(n,l)+sum2(gradq3d(n,l,:),ed)  
    254              else 
    255                 ed = cc(n+u_ldown,l,:) - xyz_i(n+t_ldown,:) 
    256                 qe = qi(n+t_ldown,l)+sum2(gradq3d(n+t_ldown,l,:),ed)  
    257              endif 
    258              qflux(n+u_ldown,l) = hfluxt(n+u_ldown,l)*qe 
    259           end do 
    260        end do 
     231             IF (ne(ij,lup)*hfluxt(ij+u_lup,l)>0) THEN  
     232                ed = cc(ij+u_lup,l,:) - xyz_i(ij,:) 
     233                qe = qi(ij,l)+sum2(gradq3d(ij,l,:),ed) 
     234             ELSE 
     235                ed = cc(ij+u_lup,l,:) - xyz_i(ij+t_lup,:) 
     236                qe = qi(ij+t_lup,l)+sum2(gradq3d(ij+t_lup,l,:),ed)  
     237             ENDIF 
     238             qflux(ij+u_lup,l) = hfluxt(ij+u_lup,l)*qe  
     239 
     240             IF (ne(ij,ldown)*hfluxt(ij+u_ldown,l)>0) THEN  
     241                ed = cc(ij+u_ldown,l,:) - xyz_i(ij,:) 
     242                qe = qi(ij,l)+sum2(gradq3d(ij,l,:),ed)  
     243             ELSE 
     244                ed = cc(ij+u_ldown,l,:) - xyz_i(ij+t_ldown,:) 
     245                qe = qi(ij+t_ldown,l)+sum2(gradq3d(ij+t_ldown,l,:),ed)  
     246             ENDIF 
     247             qflux(ij+u_ldown,l) = hfluxt(ij+u_ldown,l)*qe 
     248       END DO 
    261249    END DO 
    262250 
    263251    ! update q and, if update_mass, update mass 
    264252    DO l = ll_begin,ll_end 
    265        DO j=jj_begin,jj_end 
    266           DO i=ii_begin,ii_end 
    267              n=(j-1)*iim+i   
     253!$SIMD 
     254      DO ij=ij_begin,ij_end 
    268255             ! sign convention as Ringler et al. (2010) eq. 21 
    269              dmass =   hfluxt(n+u_right,l)*ne(n,right)   & 
    270                      + hfluxt(n+u_lup,l)  *ne(n,lup)     & 
    271                      + hfluxt(n+u_ldown,l)*ne(n,ldown)   & 
    272                      + hfluxt(n+u_rup,l)  *ne(n,rup)     & 
    273                      + hfluxt(n+u_left,l) *ne(n,left)    & 
    274                      + hfluxt(n+u_rdown,l)*ne(n,rdown) 
    275  
    276              dq    =   qflux(n+u_right,l) *ne(n,right)   & 
    277                      + qflux(n+u_lup,l)   *ne(n,lup)     & 
    278                      + qflux(n+u_ldown,l) *ne(n,ldown)   & 
    279                      + qflux(n+u_rup,l)   *ne(n,rup)     & 
    280                      + qflux(n+u_left,l)  *ne(n,left)    & 
    281                      + qflux(n+u_rdown,l) *ne(n,rdown) 
    282  
    283              newmass = mass(n,l) - dmass/Ai(n) 
    284              qi(n,l) = (qi(n,l)*mass(n,l) - dq/Ai(n) ) / newmass 
    285              IF(update_mass) mass(n,l) = newmass 
    286           END DO 
     256             dmass =   hfluxt(ij+u_right,l)*ne(ij,right)   & 
     257                     + hfluxt(ij+u_lup,l)  *ne(ij,lup)     & 
     258                     + hfluxt(ij+u_ldown,l)*ne(ij,ldown)   & 
     259                     + hfluxt(ij+u_rup,l)  *ne(ij,rup)     & 
     260                     + hfluxt(ij+u_left,l) *ne(ij,left)    & 
     261                     + hfluxt(ij+u_rdown,l)*ne(ij,rdown) 
     262 
     263             dq    =   qflux(ij+u_right,l) *ne(ij,right)   & 
     264                     + qflux(ij+u_lup,l)   *ne(ij,lup)     & 
     265                     + qflux(ij+u_ldown,l) *ne(ij,ldown)   & 
     266                     + qflux(ij+u_rup,l)   *ne(ij,rup)     & 
     267                     + qflux(ij+u_left,l)  *ne(ij,left)    & 
     268                     + qflux(ij+u_rdown,l) *ne(ij,rdown) 
     269 
     270              
     271             newmass = mass(ij,l) - dmass/Ai(ij) 
     272             qi(ij,l) = (qi(ij,l)*mass(ij,l) - dq/Ai(ij) ) / newmass 
     273             IF(update_mass) mass(ij,l) = newmass 
     274 
    287275       END DO 
    288276    END DO 

  • codes/icosagcm/trunk/src/advect_tracer.f90

    r151 r174  
    1010  TYPE(t_field),POINTER :: f_one_over_sqrt_leng(:) 
    1111 
    12   TYPE(t_message) :: req_u, req_wfluxt, req_q, req_rhodz, req_gradq3d 
     12  TYPE(t_message) :: req_u, req_cc, req_wfluxt, req_q, req_rhodz, req_gradq3d 
    1313 
    1414  REAL(rstd), PARAMETER :: pente_max=2.0 ! for vlz 
    … …  
    8181      first=.FALSE. 
    8282      CALL init_message(f_u,req_e1_vect,req_u) 
     83      CALL init_message(f_cc,req_e1_scal,req_cc) 
    8384      CALL init_message(f_wfluxt,req_i1,req_wfluxt) 
    8485      CALL init_message(f_q,req_i1,req_q) 
    … …  
    124125    END DO 
    125126 
     127    CALL send_message(f_cc,req_cc) 
     128 
    126129 
    127130    ! horizontal transport - split in two to place transfer of gradq3d 
    … …  
    138141 
    139142       CALL send_message(f_gradq3d,req_gradq3d) 
     143       CALL wait_message(req_cc) 
    140144       CALL wait_message(req_gradq3d) 
    141145 
    … …  
    205209 
    206210    REAL(rstd) :: dzqmax, newmass, sigw, qq, w 
    207     INTEGER :: i,ij,l,j 
     211    INTEGER :: i,ij,l,j,ijb,ije 
    208212 
    209213    CALL trace_start("vlz") 
     214      
     215     ijb=((jj_begin-halo)-1)*iim+ii_begin-halo 
     216     ije = ((jj_end+halo)-1)*iim+ii_end+halo 
    210217 
    211218    ! finite difference of q 
    212219 
    213220     DO l=ll_beginp1,ll_end 
    214        DO j=jj_begin-halo,jj_end+halo 
    215           DO i=ii_begin-halo,ii_end+halo 
    216              ij=(j-1)*iim+i 
    217              dzqw(ij,l)=q(ij,l)-q(ij,l-1) 
    218              adzqw(ij,l)=abs(dzqw(ij,l)) 
    219           ENDDO 
     221!$SIMD 
     222       DO ij=ijb,ije 
     223         dzqw(ij,l)=q(ij,l)-q(ij,l-1) 
     224         adzqw(ij,l)=abs(dzqw(ij,l)) 
    220225       ENDDO 
    221226    ENDDO 
    … …  
    228233 
    229234     DO l=ll_beginp1,ll_endm1 
    230        DO j=jj_begin-halo,jj_end+halo 
    231           DO i=ii_begin-halo,ii_end+halo 
    232              ij=(j-1)*iim+i 
    233              IF(dzqw(ij,l)*dzqw(ij,l+1).gt.0.) THEN 
    234                 dzq(ij,l) = 0.5*( dzqw(ij,l)+dzqw(ij,l+1) ) 
    235                 dzqmax    = pente_max * min( adzqw(ij,l),adzqw(ij,l+1) ) 
    236                 dzq(ij,l) = sign( min(abs(dzq(ij,l)),dzqmax) , dzq(ij,l) )  ! NB : sign(a,b)=a*sign(b) 
    237              ELSE 
    238                 dzq(ij,l)=0. 
    239              ENDIF 
    240           ENDDO 
     235!$SIMD 
     236       DO ij=ijb,ije  
     237         IF(dzqw(ij,l)*dzqw(ij,l+1).gt.0.) THEN 
     238             dzq(ij,l) = 0.5*( dzqw(ij,l)+dzqw(ij,l+1) ) 
     239             dzqmax    = pente_max * min( adzqw(ij,l),adzqw(ij,l+1) ) 
     240             dzq(ij,l) = sign( min(abs(dzq(ij,l)),dzqmax) , dzq(ij,l) )  ! NB : sign(a,b)=a*sign(b) 
     241          ELSE 
     242             dzq(ij,l)=0. 
     243          ENDIF 
    241244       ENDDO 
    242245    ENDDO 
    … …  
    245248    ! 0 slope in top and bottom layers 
    246249    IF (omp_first) THEN 
    247       DO j=jj_begin-halo,jj_end+halo 
    248          DO i=ii_begin-halo,ii_end+halo 
    249            ij=(j-1)*iim+i 
     250      DO ij=ijb,ije 
    250251           dzq(ij,1)=0. 
    251          ENDDO 
    252252      ENDDO 
    253253    ENDIF 
    254254       
    255255    IF (omp_last) THEN 
    256       DO j=jj_begin-halo,jj_end+halo 
    257          DO i=ii_begin-halo,ii_end+halo 
     256      DO ij=ijb,ije 
    258257          dzq(ij,llm)=0. 
    259          ENDDO 
    260258      ENDDO 
    261259    ENDIF 
    … …  
    267265    ! then amount of q leaving level l/l+1 = wq = w * qq 
    268266     DO l=ll_beginp1,ll_end 
    269        DO j=jj_begin-halo,jj_end+halo 
    270           DO i=ii_begin-halo,ii_end+halo 
    271              ij=(j-1)*iim+i 
     267!$SIMD 
     268       DO ij=ijb,ije 
    272269             w = fac*wfluxt(ij,l) 
    273270             IF(w>0.) THEN  ! upward transport, upwind side is at level l  
    … …  
    279276             ENDIF 
    280277             wq(ij,l) = w*qq 
    281           ENDDO 
    282278       ENDDO 
    283279    END DO 
    284280    ! wq = 0 at top and bottom 
    285281    IF (omp_first) THEN 
    286       DO j=jj_begin-halo,jj_end+halo 
    287          DO i=ii_begin-halo,ii_end+halo 
    288             ij=(j-1)*iim+i 
     282       DO ij=ijb,ije 
    289283            wq(ij,1)=0. 
    290          ENDDO 
    291284      END DO 
    292285    ENDIF 
    293286     
    294287    IF (omp_last) THEN 
    295       DO j=jj_begin-halo,jj_end+halo 
    296          DO i=ii_begin-halo,ii_end+halo 
    297             ij=(j-1)*iim+i 
     288      DO ij=ijb,ije 
    298289            wq(ij,llm+1)=0. 
    299          ENDDO 
    300290      END DO 
    301291    ENDIF 
    … …  
    307297    ! update q, mass is updated only after all q's have been updated 
    308298    DO l=ll_begin,ll_end 
    309        DO j=jj_begin-halo,jj_end+halo 
    310           DO i=ii_begin-halo,ii_end+halo 
    311              ij=(j-1)*iim+i 
     299!$SIMD 
     300       DO ij=ijb,ije 
    312301             newmass = mass(ij,l) + fac*(wfluxt(ij,l)-wfluxt(ij,l+1)) 
    313302             q(ij,l) = ( q(ij,l)*mass(ij,l) + wq(ij,l)-wq(ij,l+1) ) / newmass 
    314303             IF(update_mass) mass(ij,l)=newmass 
    315           ENDDO 
    316304       ENDDO 
    317305    END DO 

  • codes/icosagcm/trunk/src/caldyn_gcm.f90

    r171 r174  
    9797          wwuu=f_wwuu(ind) 
    9898           
    99           DO j=jj_begin-1,jj_end+1 
    100              DO i=ii_begin-1,ii_end+1 
    101                 ij=(j-1)*iim+i 
    102                 ! lower BCs : geopot=phis, wflux=0, wwuu=0 
    103                 geopot(ij,1) = phis(ij) 
    104                 wflux(ij,1) = 0. 
    105                 wwuu(ij+u_right,1)=0    
    106                 wwuu(ij+u_lup,1)=0    
    107                 wwuu(ij+u_ldown,1)=0 
    108                 ! top BCs : wflux=0, wwuu=0 
    109                 wflux(ij,llm+1)  = 0. 
    110                 wwuu(ij+u_right,llm+1)=0    
    111                 wwuu(ij+u_lup,llm+1)=0    
    112                 wwuu(ij+u_ldown,llm+1)=0 
    113              ENDDO 
     99          DO ij=ij_begin_ext,ij_end_ext 
     100              ! lower BCs : geopot=phis, wflux=0, wwuu=0 
     101              geopot(ij,1) = phis(ij) 
     102              wflux(ij,1) = 0. 
     103              wwuu(ij+u_right,1)=0    
     104              wwuu(ij+u_lup,1)=0    
     105              wwuu(ij+u_ldown,1)=0 
     106              ! top BCs : wflux=0, wwuu=0 
     107              wflux(ij,llm+1)  = 0. 
     108              wwuu(ij+u_right,llm+1)=0    
     109              wwuu(ij+u_lup,llm+1)=0    
     110              wwuu(ij+u_ldown,llm+1)=0 
    114111          ENDDO 
    115112       END DO 
    … …  
    303300!  REAL(rstd) :: qv(2*iim*jjm,llm)     ! potential velocity   
    304301   
    305  
    306302  CALL trace_start("compute_pvort")   
    307303 
    … …  
    316312     DO l = ll_begin,ll_end 
    317313        CALL test_message(req_u)  
    318         DO j=jj_begin-1,jj_end+1 
    319            DO i=ii_begin-1,ii_end+1 
    320               ij=(j-1)*iim+i 
    321               m = ( mass_dak(l)+ps(ij)*mass_dbk(l) )/g 
    322               rhodz(ij,l) = m 
    323               theta(ij,l) = theta_rhodz(ij,l)/rhodz(ij,l) 
    324            ENDDO 
     314!DIR$ SIMD 
     315        DO ij=ij_begin_ext,ij_end_ext 
     316           m = ( mass_dak(l)+ps(ij)*mass_dbk(l) )/g 
     317           rhodz(ij,l) = m 
     318           theta(ij,l) = theta_rhodz(ij,l)/rhodz(ij,l) 
    325319        ENDDO 
    326320     ENDDO 
    … …  
    328322     DO l = ll_begin,ll_end 
    329323        CALL test_message(req_u)  
    330         DO j=jj_begin-1,jj_end+1 
    331            DO i=ii_begin-1,ii_end+1 
    332               ij=(j-1)*iim+i 
    333               theta(ij,l) = theta_rhodz(ij,l)/rhodz(ij,l) 
    334            ENDDO 
    335         ENDDO 
     324!DIR$ SIMD 
     325       DO ij=ij_begin_ext,ij_end_ext 
     326         theta(ij,l) = theta_rhodz(ij,l)/rhodz(ij,l) 
     327       ENDDO 
    336328     ENDDO 
    337329  END IF 
    … …  
    341333!!! Compute shallow-water potential vorticity 
    342334  DO l = ll_begin,ll_end 
    343  
    344     DO j=jj_begin-1,jj_end+1 
    345        DO i=ii_begin-1,ii_end+1 
    346           ij=(j-1)*iim+i 
    347            
     335!DIR$ SIMD 
     336     DO ij=ij_begin_ext,ij_end_ext          
    348337          etav= 1./Av(ij+z_up)*(  ne_rup        * u(ij+u_rup,l)        * de(ij+u_rup)         & 
    349338                                + ne_left * u(ij+t_rup+u_left,l) * de(ij+t_rup+u_left)  & 
    … …  
    366355          qv(ij+z_down,l) =( etav+fv(ij+z_down) )/hv 
    367356 
    368         ENDDO 
    369357      ENDDO 
    370358 
    371       DO j=jj_begin,jj_end 
    372          DO i=ii_begin,ii_end 
    373             ij=(j-1)*iim+i 
    374             qu(ij+u_right,l) = 0.5*(qv(ij+z_rdown,l)+qv(ij+z_rup,l))   
    375             qu(ij+u_lup,l) = 0.5*(qv(ij+z_up,l)+qv(ij+z_lup,l))   
    376             qu(ij+u_ldown,l) = 0.5*(qv(ij+z_ldown,l)+qv(ij+z_down,l))   
    377          END DO 
     359!DIR$ SIMD 
     360      DO ij=ij_begin,ij_end 
     361         qu(ij+u_right,l) = 0.5*(qv(ij+z_rdown,l)+qv(ij+z_rup,l))   
     362         qu(ij+u_lup,l) = 0.5*(qv(ij+z_up,l)+qv(ij+z_lup,l))   
     363         qu(ij+u_ldown,l) = 0.5*(qv(ij+z_ldown,l)+qv(ij+z_down,l))   
    378364      END DO 
    379365       
    … …  
    407393       DO l = 1,llm 
    408394          !$OMP DO SCHEDULE(STATIC)  
    409           DO j=jj_begin-1,jj_end+1 
    410              DO i=ii_begin-1,ii_end+1 
    411                 ij=(j-1)*iim+i 
    412                 p_ik = ptop + mass_ak(l) + mass_bk(l)*ps(ij) ! FIXME : leave ps for the moment ; change ps to Ms later 
    413                 !         p_ik = ptop + g*(mass_ak(l)+ mass_bk(l)*ps(i,j)) 
    414                 exner_ik = cpp * (p_ik/preff) ** kappa 
    415                 pk(ij,l) = exner_ik 
    416                 ! specific volume v = kappa*theta*pi/p = dphi/g/rhodz 
    417                 geopot(ij,l+1) = geopot(ij,l) + (g*kappa)*rhodz(ij,l)*theta(ij,l)*exner_ik/p_ik 
    418              ENDDO 
     395!DIR$ SIMD 
     396         DO ij=ij_begin_ext,ij_end_ext          
     397            p_ik = ptop + mass_ak(l) + mass_bk(l)*ps(ij) ! FIXME : leave ps for the moment ; change ps to Ms later 
     398            !         p_ik = ptop + g*(mass_ak(l)+ mass_bk(l)*ps(i,j)) 
     399            exner_ik = cpp * (p_ik/preff) ** kappa 
     400            pk(ij,l) = exner_ik 
     401            ! specific volume v = kappa*theta*pi/p = dphi/g/rhodz 
     402            geopot(ij,l+1) = geopot(ij,l) + (g*kappa)*rhodz(ij,l)*theta(ij,l)*exner_ik/p_ik 
    419403          ENDDO 
    420404       ENDDO 
    … …  
    427411 
    428412       ! uppermost layer 
    429        DO j=jj_begin-1,jj_end+1 
    430           DO i=ii_begin-1,ii_end+1 
    431              ij=(j-1)*iim+i 
    432              pk(ij,llm) = ptop + (.5*g)*rhodz(ij,llm) 
    433           END DO 
     413!DIR$ SIMD 
     414       DO ij=ij_begin_ext,ij_end_ext          
     415         pk(ij,llm) = ptop + (.5*g)*rhodz(ij,llm) 
    434416       END DO 
    435417       ! other layers 
    436418       DO l = llm-1, 1, -1 
    437419          !$OMP DO SCHEDULE(STATIC)  
    438           DO j=jj_begin-1,jj_end+1 
    439              DO i=ii_begin-1,ii_end+1 
    440                 ij=(j-1)*iim+i 
    441                 pk(ij,l) = pk(ij,l+1) + (.5*g)*(rhodz(ij,l)+rhodz(ij,l+1)) 
    442              END DO 
     420!DIR$ SIMD 
     421          DO ij=ij_begin_ext,ij_end_ext          
     422             pk(ij,l) = pk(ij,l+1) + (.5*g)*(rhodz(ij,l)+rhodz(ij,l+1)) 
    443423          END DO 
    444424       END DO 
    445425       ! surface pressure (for diagnostics) 
    446        DO j=jj_begin-1,jj_end+1 
    447           DO i=ii_begin-1,ii_end+1 
    448              ij=(j-1)*iim+i 
    449              ps(ij) = pk(ij,1) + (.5*g)*rhodz(ij,1) 
    450           END DO 
     426       DO ij=ij_begin_ext,ij_end_ext          
     427         ps(ij) = pk(ij,1) + (.5*g)*rhodz(ij,1) 
    451428       END DO 
    452429 
    … …  
    454431          DO l = 1,llm 
    455432             !$OMP DO SCHEDULE(STATIC)  
    456              DO j=jj_begin-1,jj_end+1 
    457                 DO i=ii_begin-1,ii_end+1 
    458                    ij=(j-1)*iim+i 
    459                    geopot(ij,l+1) = geopot(ij,l) + g*rhodz(ij,l) 
    460                 ENDDO 
     433!DIR$ SIMD 
     434             DO ij=ij_begin_ext,ij_end_ext          
     435               geopot(ij,l+1) = geopot(ij,l) + g*rhodz(ij,l) 
    461436             ENDDO 
    462437          ENDDO 
    … …  
    464439          DO l = 1,llm 
    465440             !$OMP DO SCHEDULE(STATIC)  
    466              DO j=jj_begin-1,jj_end+1 
    467                 DO i=ii_begin-1,ii_end+1 
    468                    ij=(j-1)*iim+i 
    469                    p_ik = pk(ij,l) 
    470                    exner_ik = cpp * (p_ik/preff) ** kappa 
    471                    geopot(ij,l+1) = geopot(ij,l) + (g*kappa)*rhodz(ij,l)*theta(ij,l)*exner_ik/p_ik  
    472                    pk(ij,l) = exner_ik 
    473                 ENDDO 
     441!DIR$ SIMD 
     442              DO ij=ij_begin_ext,ij_end_ext          
     443                p_ik = pk(ij,l) 
     444                exner_ik = cpp * (p_ik/preff) ** kappa 
     445                geopot(ij,l+1) = geopot(ij,l) + (g*kappa)*rhodz(ij,l)*theta(ij,l)*exner_ik/p_ik  
     446                pk(ij,l) = exner_ik 
    474447             ENDDO 
    475448          ENDDO 
    … …  
    514487!!!  Compute mass and theta fluxes 
    515488    IF (caldyn_conserv==energy) CALL test_message(req_qu)  
    516     DO j=jj_begin-1,jj_end+1 
    517       DO i=ii_begin-1,ii_end+1 
    518         ij=(j-1)*iim+i 
     489!DIR$ SIMD 
     490    DO ij=ij_begin_ext,ij_end_ext          
    519491        hflux(ij+u_right,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_right,l))*u(ij+u_right,l)*le(ij+u_right) 
    520492        hflux(ij+u_lup,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_lup,l))*u(ij+u_lup,l)*le(ij+u_lup) 
    … …  
    524496        Ftheta(ij+u_lup,l)=0.5*(theta(ij,l)+theta(ij+t_lup,l))*hflux(ij+u_lup,l) 
    525497        Ftheta(ij+u_ldown,l)=0.5*(theta(ij,l)+theta(ij+t_ldown,l))*hflux(ij+u_ldown,l) 
    526       ENDDO 
    527498    ENDDO 
    528499 
    529500!!! compute horizontal divergence of fluxes 
    530     DO j=jj_begin,jj_end 
    531       DO i=ii_begin,ii_end 
    532         ij=(j-1)*iim+i   
     501!DIR$ SIMD 
     502    DO ij=ij_begin,ij_end          
    533503        ! convm = -div(mass flux), sign convention as in Ringler et al. 2012, eq. 21 
    534504        convm(ij,l)= -1./Ai(ij)*(ne_right*hflux(ij+u_right,l)   +  & 
    … …  
    547517                                 ne_ldown*Ftheta(ij+u_ldown,l)    +  &   
    548518                                 ne_rdown*Ftheta(ij+u_rdown,l)) 
    549       ENDDO 
    550519    ENDDO 
    551520 
    … …  
    560529 
    561530        DO l=ll_begin,ll_end 
    562           DO j=jj_begin,jj_end 
    563              DO i=ii_begin,ii_end 
    564                 ij=(j-1)*iim+i 
    565  
     531!DIR$ SIMD 
     532          DO ij=ij_begin,ij_end          
     533     
    566534                uu = wee(ij+u_right,1,1)*hflux(ij+u_rup,l)*(qu(ij+u_right,l)+qu(ij+u_rup,l))+                             & 
    567535                     wee(ij+u_right,2,1)*hflux(ij+u_lup,l)*(qu(ij+u_right,l)+qu(ij+u_lup,l))+                             & 
    … …  
    601569                du(ij+u_ldown,l) = .5*uu/de(ij+u_ldown) 
    602570                 
    603              ENDDO 
    604571          ENDDO 
    605572       ENDDO 
    … …  
    608575   
    609576        DO l=ll_begin,ll_end 
    610           DO j=jj_begin,jj_end 
    611              DO i=ii_begin,ii_end 
    612                 ij=(j-1)*iim+i 
     577!DIR$ SIMD 
     578          DO ij=ij_begin,ij_end          
    613579 
    614580                uu = wee(ij+u_right,1,1)*hflux(ij+u_rup,l)+                           & 
    … …  
    649615                du(ij+u_ldown,l) = qu(ij+u_ldown,l)*uu/de(ij+u_ldown) 
    650616      
    651              ENDDO 
    652617          ENDDO 
    653618       ENDDO 
    … …  
    661626 
    662627       DO l=ll_begin,ll_end 
    663           DO j=jj_begin,jj_end 
    664              DO i=ii_begin,ii_end 
    665                 ij=(j-1)*iim+i 
     628!DIR$ SIMD 
     629          DO ij=ij_begin,ij_end          
    666630                 
    667631                berni(ij,l) = pk(ij,l) + & 
    … …  
    673637                                     le(ij+u_rdown)*de(ij+u_rdown)*u(ij+u_rdown,l)**2 )   
    674638                pk(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1)) 
    675              ENDDO 
    676639          ENDDO 
    677640       ENDDO 
    … …  
    680643 
    681644       DO l=ll_begin,ll_end 
    682           DO j=jj_begin,jj_end 
    683              DO i=ii_begin,ii_end 
    684                 ij=(j-1)*iim+i 
     645!DIR$ SIMD 
     646            DO ij=ij_begin,ij_end          
    685647                 
    686648                berni(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1)) & 
    … …  
    691653                                     le(ij+u_ldown)*de(ij+u_ldown)*u(ij+u_ldown,l)**2 +    & 
    692654                                     le(ij+u_rdown)*de(ij+u_rdown)*u(ij+u_rdown,l)**2 )   
    693         
    694              ENDDO 
    695655          ENDDO 
    696656       ENDDO 
    … …  
    700660!!! Add gradients of Bernoulli and Exner functions to du 
    701661    DO l=ll_begin,ll_end 
    702        DO j=jj_begin,jj_end 
    703           DO i=ii_begin,ii_end 
    704              ij=(j-1)*iim+i 
     662!DIR$ SIMD 
     663       DO ij=ij_begin,ij_end          
    705664         
    706665             du(ij+u_right,l) = du(ij+u_right,l) + 1/de(ij+u_right) * (             & 
    … …  
    720679                                   + ne_ldown*berni(ij,l)+ne_rup*berni(ij+t_ldown,l) ) 
    721680 
    722           ENDDO 
    723681       ENDDO 
    724682    ENDDO 
    … …  
    750708    REAL(rstd) :: p_ik, exner_ik 
    751709 
     710 
    752711  CALL trace_start("compute_caldyn_vert") 
    753712 
    … …  
    757716    IF (caldyn_conserv==energy) CALL test_message(req_qu)  
    758717!$OMP DO SCHEDULE(STATIC)  
    759     DO j=jj_begin,jj_end 
    760       DO i=ii_begin,ii_end 
    761         ij=(j-1)*iim+i 
     718!DIR$ SIMD 
     719    DO ij=ij_begin,ij_end          
    762720        convm(ij,l) = convm(ij,l) + convm(ij,l+1) 
    763       ENDDO 
    764721    ENDDO 
    765722  ENDDO 
    … …  
    770727! compute dps 
    771728  IF (omp_first) THEN 
    772     DO j=jj_begin,jj_end 
    773       DO i=ii_begin,ii_end 
    774         ij=(j-1)*iim+i 
     729!DIR$ SIMD 
     730     DO ij=ij_begin,ij_end          
    775731        ! dps/dt = -int(div flux)dz 
    776732        dps(ij) = convm(ij,1) * g  
    777       ENDDO 
    778733    ENDDO 
    779734  ENDIF 
    … …  
    782737  DO l=ll_beginp1,ll_end 
    783738    IF (caldyn_conserv==energy) CALL test_message(req_qu)  
    784     DO j=jj_begin,jj_end 
    785       DO i=ii_begin,ii_end 
    786         ij=(j-1)*iim+i 
     739!DIR$ SIMD 
     740      DO ij=ij_begin,ij_end          
    787741        ! w = int(z,ztop,div(flux)dz) + B(eta)dps/dt 
    788742        ! => w>0 for upward transport 
    789743        wflux( ij, l ) = bp(l) * convm( ij, 1 ) - convm( ij, l )  
    790       ENDDO 
    791744    ENDDO 
    792745  ENDDO 
    793746 
    794747  DO l=ll_begin,ll_endm1 
    795     DO j=jj_begin,jj_end 
    796       DO i=ii_begin,ii_end 
    797         ij=(j-1)*iim+i 
     748!DIR$ SIMD 
     749     DO ij=ij_begin,ij_end          
    798750        dtheta_rhodz(ij, l ) = dtheta_rhodz(ij, l )  -   0.5 * (  wflux(ij,l+1) * (theta(ij,l) + theta(ij,l+1))) 
    799       ENDDO 
    800751    ENDDO 
    801752  ENDDO 
    802753   
    803754  DO l=ll_beginp1,ll_end 
    804     DO j=jj_begin,jj_end 
    805       DO i=ii_begin,ii_end 
    806         ij=(j-1)*iim+i 
     755!DIR$ SIMD 
     756      DO ij=ij_begin,ij_end          
    807757        dtheta_rhodz(ij, l ) = dtheta_rhodz(ij, l )  +   0.5 * ( wflux(ij,l  ) * (theta(ij,l-1) + theta(ij,l) ) ) 
    808       ENDDO 
    809758    ENDDO 
    810759  ENDDO 
    … …  
    812761! Compute vertical transport 
    813762  DO l=ll_beginp1,ll_end 
    814     DO j=jj_begin,jj_end 
    815       DO i=ii_begin,ii_end 
    816         ij=(j-1)*iim+i 
     763!DIR$ SIMD 
     764      DO ij=ij_begin,ij_end          
    817765        wwuu(ij+u_right,l) = 0.5*( wflux(ij,l) + wflux(ij+t_right,l)) * (u(ij+u_right,l) - u(ij+u_right,l-1)) 
    818766        wwuu(ij+u_lup,l) = 0.5* ( wflux(ij,l) + wflux(ij+t_lup,l)) * (u(ij+u_lup,l) - u(ij+u_lup,l-1)) 
    819767        wwuu(ij+u_ldown,l) = 0.5*( wflux(ij,l) + wflux(ij+t_ldown,l)) * (u(ij+u_ldown,l) - u(ij+u_ldown,l-1)) 
    820        ENDDO 
    821768     ENDDO 
    822769   ENDDO 
    … …  
    827774! Add vertical transport to du 
    828775  DO l=ll_begin,ll_end 
    829     DO j=jj_begin,jj_end 
    830       DO i=ii_begin,ii_end 
    831         ij=(j-1)*iim+i 
     776!DIR$ SIMD 
     777     DO ij=ij_begin,ij_end          
    832778        du(ij+u_right, l )   = du(ij+u_right,l)  - (wwuu(ij+u_right,l+1)+ wwuu(ij+u_right,l)) / (rhodz(ij,l)+rhodz(ij+t_right,l)) 
    833779        du(ij+u_lup, l )     = du(ij+u_lup,l)    - (wwuu(ij+u_lup,l+1)  + wwuu(ij+u_lup,l))   / (rhodz(ij,l)+rhodz(ij+t_lup,l)) 
    834780        du(ij+u_ldown, l )   = du(ij+u_ldown,l)  - (wwuu(ij+u_ldown,l+1)+ wwuu(ij+u_ldown,l)) / (rhodz(ij,l)+rhodz(ij+t_ldown,l)) 
    835       ENDDO 
    836781    ENDDO 
    837782  ENDDO       

  • codes/icosagcm/trunk/src/dimensions.f90

    r151 r174  
    99    INTEGER :: ii_nb 
    1010    INTEGER :: jj_nb 
    11      
     11    INTEGER :: ij_begin 
     12    INTEGER :: ij_end 
     13    INTEGER :: ij_begin_ext 
     14    INTEGER :: ij_end_ext 
     15         
    1216    INTEGER :: t_right 
    1317    INTEGER :: t_rup 
    … …  
    5660     ii_nb=d%ii_nb 
    5761     jj_nb=d%jj_nb 
    58       
     62     
     63     ij_begin=(jj_begin-1)*iim+ii_begin 
     64     ij_end=(jj_end-1)*iim+ii_end  
     65     ij_begin_ext=((jj_begin-1)-1)*iim+ii_begin-1 
     66     ij_end_ext = ((jj_end+1)-1)*iim+ii_end+1 
     67         
    5968     t_right=d%t_right 
    6069     t_rup=d%t_rup 

  • codes/icosagcm/trunk/src/dissip_gcm.f90

    r167 r174  
    7777  CHARACTER(len=255)    :: rayleigh_friction_key 
    7878  REAL(rstd)            :: mintau 
     79  INTEGER               :: seed_size 
     80  INTEGER,ALLOCATABLE   :: seed(:) 
    7981   
    8082             
    81   INTEGER :: i,j,n,ind,it,iter 
     83  INTEGER :: i,j,ij,ind,it,iter 
    8284 
    8385   rayleigh_friction_key='none' 
    … …  
    142144    cdivgrad=-1 
    143145    cgradrot=-1 
    144    
    145     CALL RANDOM_SEED() 
    146146     
    147147    DO ind=1,ndomain 
    … …  
    152152      DO j=jj_begin,jj_end 
    153153        DO i=ii_begin,ii_end 
    154           n=(j-1)*iim+i 
     154          ij=(j-1)*iim+i    
    155155          CALL RANDOM_NUMBER(r) 
    156           u(n+u_right)=r-0.5 
     156          u(ij+u_right)=r-0.5 
    157157          CALL RANDOM_NUMBER(r) 
    158           u(n+u_lup)=r-0.5 
     158          u(ij+u_lup)=r-0.5 
    159159          CALL RANDOM_NUMBER(r) 
    160           u(n+u_ldown)=r-0.5 
    161        ENDDO 
    162       ENDDO 
    163          
     160          u(ij+u_ldown)=r-0.5 
     161        ENDDO 
     162      ENDDO         
    164163    ENDDO 
    165164   
    … …  
    191190        DO j=jj_begin,jj_end 
    192191          DO i=ii_begin,ii_end 
    193             n=(j-1)*iim+i 
    194             if (le(n+u_right)>1e-100) dumax=MAX(dumax,ABS(du(n+u_right))) 
    195             if (le(n+u_lup)>1e-100) dumax=MAX(dumax,ABS(du(n+u_lup))) 
    196             if (le(n+u_ldown)>1e-100) dumax=MAX(dumax,ABS(du(n+u_ldown))) 
     192            ij=(j-1)*iim+i    
     193            if (le(ij+u_right)>1e-100) dumax=MAX(dumax,ABS(du(ij+u_right))) 
     194            if (le(ij+u_lup)>1e-100)   dumax=MAX(dumax,ABS(du(ij+u_lup))) 
     195            if (le(ij+u_ldown)>1e-100) dumax=MAX(dumax,ABS(du(ij+u_ldown))) 
    197196          ENDDO 
    198197        ENDDO 
    … …  
    228227      u=f_u(ind) 
    229228 
    230      DO j=jj_begin,jj_end 
     229      DO j=jj_begin,jj_end 
    231230        DO i=ii_begin,ii_end 
    232           n=(j-1)*iim+i 
     231          ij=(j-1)*iim+i    
    233232          CALL RANDOM_NUMBER(r) 
    234           u(n+u_right)=r-0.5 
     233          u(ij+u_right)=r-0.5 
    235234          CALL RANDOM_NUMBER(r) 
    236           u(n+u_lup)=r-0.5 
     235          u(ij+u_lup)=r-0.5 
    237236          CALL RANDOM_NUMBER(r) 
    238           u(n+u_ldown)=r-0.5 
    239        ENDDO 
    240       ENDDO 
    241          
     237          u(ij+u_ldown)=r-0.5 
     238        ENDDO 
     239      ENDDO         
    242240    ENDDO 
    243241 
    … …  
    268266        DO j=jj_begin,jj_end 
    269267          DO i=ii_begin,ii_end 
    270             n=(j-1)*iim+i 
    271             if (le(n+u_right)>1e-100) dumax=MAX(dumax,ABS(du(n+u_right))) 
    272             if (le(n+u_lup)>1e-100) dumax=MAX(dumax,ABS(du(n+u_lup))) 
    273             if (le(n+u_ldown)>1e-100) dumax=MAX(dumax,ABS(du(n+u_ldown))) 
     268            ij=(j-1)*iim+i    
     269            if (le(ij+u_right)>1e-100) dumax=MAX(dumax,ABS(du(ij+u_right))) 
     270            if (le(ij+u_lup)>1e-100) dumax=MAX(dumax,ABS(du(ij+u_lup))) 
     271            if (le(ij+u_ldown)>1e-100) dumax=MAX(dumax,ABS(du(ij+u_ldown))) 
    274272          ENDDO 
    275273        ENDDO 
    … …  
    304302      theta=f_theta(ind) 
    305303  
    306        DO j=jj_begin,jj_end 
    307           DO i=ii_begin,ii_end 
    308             n=(j-1)*iim+i 
    309             CALL RANDOM_NUMBER(r) 
    310             theta(n)=r-0.5 
    311          ENDDO 
    312         ENDDO 
    313          
     304      DO j=jj_begin,jj_end 
     305        DO i=ii_begin,ii_end 
     306          ij=(j-1)*iim+i    
     307          CALL RANDOM_NUMBER(r) 
     308          theta(ij)=r-0.5 
     309        ENDDO 
     310      ENDDO   
    314311    ENDDO 
    315312 
    … …  
    339336        DO j=jj_begin,jj_end 
    340337          DO i=ii_begin,ii_end 
    341             n=(j-1)*iim+i 
    342             dthetamax=MAX(dthetamax,ABS(dtheta(n))) 
     338            ij=(j-1)*iim+i    
     339            dthetamax=MAX(dthetamax,ABS(dtheta(ij))) 
    343340          ENDDO 
    344341        ENDDO 
    345342      ENDDO 
     343       
    346344      IF (using_mpi) THEN 
    347345        CALL MPI_ALLREDUCE(dthetamax,dthetamax1,1,MPI_REAL8,MPI_MAX,comm_icosa,ierr) 
    348346        dthetamax=dthetamax1 
    349347      ENDIF 
     348       
    350349      IF (is_mpi_root) PRINT *,"divgrad : it :",it ,": dthetamax",dthetamax 
    351350 
    … …  
    417416 
    418417    INTEGER :: ind, shear 
    419     INTEGER :: l,i,j,n 
     418    INTEGER :: l,ij 
    420419 
    421420!$OMP BARRIER 
    … …  
    444443 
    445444      DO l=ll_begin,ll_end 
    446         DO j=jj_begin,jj_end 
    447           DO i=ii_begin,ii_end 
    448             n=(j-1)*iim+i 
    449  
    450             due(n+u_right,l) = -0.5*( due_diss1(n+u_right,l)/tau_graddiv(l) + due_diss2(n+u_right,l)/tau_gradrot(l))*itau_dissip  
    451             due(n+u_lup,l)   = -0.5*( due_diss1(n+u_lup,l)  /tau_graddiv(l) + due_diss2(n+u_lup,l)  /tau_gradrot(l))*itau_dissip 
    452             due(n+u_ldown,l) = -0.5*( due_diss1(n+u_ldown,l)/tau_graddiv(l) + due_diss2(n+u_ldown,l)/tau_gradrot(l))*itau_dissip 
    453  
    454             dtheta_rhodz(n,l) = -0.5*dtheta_rhodz_diss(n,l)/tau_divgrad(l)*itau_dissip 
    455           ENDDO 
     445!$SIMD 
     446        DO ij=ij_begin,ij_end 
     447 
     448            due(ij+u_right,l) = -0.5*( due_diss1(ij+u_right,l)/tau_graddiv(l) + due_diss2(ij+u_right,l)/tau_gradrot(l))*itau_dissip  
     449            due(ij+u_lup,l)   = -0.5*( due_diss1(ij+u_lup,l)  /tau_graddiv(l) + due_diss2(ij+u_lup,l)  /tau_gradrot(l))*itau_dissip 
     450            due(ij+u_ldown,l) = -0.5*( due_diss1(ij+u_ldown,l)/tau_graddiv(l) + due_diss2(ij+u_ldown,l)/tau_gradrot(l))*itau_dissip 
     451 
     452            dtheta_rhodz(ij,l) = -0.5*dtheta_rhodz_diss(ij,l)/tau_divgrad(l)*itau_dissip 
    456453        ENDDO 
    457454      ENDDO 
    … …  
    490487        LOGICAL :: hybrid_eta 
    491488        INTEGER :: shift_u, shift_t, zcoords, nn 
    492         z = (phi(n,l)+phi(n+shift_t,l))/(2.*g) 
     489        z = (phi(ij,l)+phi(ij+shift_t,l))/(2.*g) 
    493490        IF(z>zh) THEN  ! relax only in the sponge layer z>zh 
    494491 
    495            nn = n+shift_u 
     492           nn = ij+shift_u 
    496493           CALL xyz2lonlat(xyz_e(nn,:),lon,lat) 
    497494           zcoords = 1 ; hybrid_eta = .FALSE. ! use z instead of p or hyam/hybm 
    … …  
    522519    REAL(rstd),POINTER    :: due(:,:) 
    523520    INTEGER :: ind 
    524     INTEGER :: it,i,j,l,ij 
     521    INTEGER :: it,l,ij 
    525522        
    526523    CALL trace_start("gradiv") 
    … …  
    532529      due=f_due(ind)  
    533530      DO  l = ll_begin, ll_end 
    534         DO j=jj_begin,jj_end 
    535           DO i=ii_begin,ii_end 
    536             ij=(j-1)*iim+i 
     531!$SIMD 
     532        DO ij=ij_begin,ij_end 
    537533             due(ij+u_right,l)=ue(ij+u_right,l) 
    538534             due(ij+u_lup,l)=ue(ij+u_lup,l) 
    539535             due(ij+u_ldown,l)=ue(ij+u_ldown,l) 
    540           ENDDO 
    541536        ENDDO 
    542537      ENDDO 
    … …  
    571566    REAL(rstd),POINTER    :: due(:,:) 
    572567    INTEGER :: ind 
    573     INTEGER :: it,i,j,l,ij 
     568    INTEGER :: it,l,ij 
    574569        
    575570    CALL trace_start("gradrot") 
    … …  
    581576      due=f_due(ind)  
    582577      DO  l = ll_begin, ll_end 
    583         DO j=jj_begin,jj_end 
    584           DO i=ii_begin,ii_end 
    585             ij=(j-1)*iim+i 
     578!$SIMD 
     579        DO ij=ij_begin,ij_end 
    586580             due(ij+u_right,l)=ue(ij+u_right,l) 
    587581             due(ij+u_lup,l)=ue(ij+u_lup,l) 
    588582             due(ij+u_ldown,l)=ue(ij+u_ldown,l) 
    589           ENDDO 
    590583        ENDDO 
    591584      ENDDO 
    … …  
    663656 
    664657    INTEGER :: ind 
    665     INTEGER :: it,i,j,l,ij 
     658    INTEGER :: it,l,ij 
    666659 
    667660    CALL trace_start("divgrad") 
    … …  
    674667      dtheta_rhodz=f_dtheta_rhodz(ind)  
    675668      DO  l = ll_begin, ll_end 
    676         DO j=jj_begin,jj_end 
    677           DO i=ii_begin,ii_end 
    678             ij=(j-1)*iim+i 
     669!$SIMD 
     670        DO ij=ij_begin,ij_end 
    679671            dtheta_rhodz(ij,l) = theta_rhodz(ij,l) / mass(ij,l) 
    680           ENDDO 
    681672        ENDDO 
    682673      ENDDO 
    … …  
    702693     
    703694      DO  l = ll_begin, ll_end 
    704         DO j=jj_begin,jj_end 
    705           DO i=ii_begin,ii_end 
    706             ij=(j-1)*iim+i 
     695!$SIMD 
     696        DO ij=ij_begin,ij_end 
    707697            dtheta_rhodz(ij,l) = dtheta_rhodz(ij,l) * mass(ij,l) 
    708           ENDDO 
    709698        ENDDO 
    710699      ENDDO 
    … …  
    725714    REAL(rstd) :: divu_i(iim*jjm,llb:lle) 
    726715     
    727     INTEGER :: i,j,n,l 
     716    INTEGER :: ij,l 
    728717 
    729718    DO l=llb,lle 
    730       DO j=jj_begin,jj_end 
    731         DO i=ii_begin,ii_end 
    732           n=(j-1)*iim+i 
    733           divu_i(n,l)=1./Ai(n)*(ne(n,right)*ue(n+u_right,l)*le(n+u_right)  +  & 
    734                              ne(n,rup)*ue(n+u_rup,l)*le(n+u_rup)        +  &   
    735                              ne(n,lup)*ue(n+u_lup,l)*le(n+u_lup)        +  &   
    736                              ne(n,left)*ue(n+u_left,l)*le(n+u_left)     +  &   
    737                              ne(n,ldown)*ue(n+u_ldown,l)*le(n+u_ldown)  +  &   
    738                              ne(n,rdown)*ue(n+u_rdown,l)*le(n+u_rdown)) 
    739         ENDDO 
     719!$SIMD 
     720      DO ij=ij_begin,ij_end 
     721          divu_i(ij,l)=1./Ai(ij)*(ne(ij,right)*ue(ij+u_right,l)*le(ij+u_right)  +  & 
     722                             ne(ij,rup)*ue(ij+u_rup,l)*le(ij+u_rup)        +  &   
     723                             ne(ij,lup)*ue(ij+u_lup,l)*le(ij+u_lup)        +  &   
     724                             ne(ij,left)*ue(ij+u_left,l)*le(ij+u_left)     +  &   
     725                             ne(ij,ldown)*ue(ij+u_ldown,l)*le(ij+u_ldown)  +  &   
     726                             ne(ij,rdown)*ue(ij+u_rdown,l)*le(ij+u_rdown)) 
    740727      ENDDO 
    741728    ENDDO 
    742729     
    743730    DO l=llb,lle 
    744       DO j=jj_begin,jj_end 
    745         DO i=ii_begin,ii_end 
     731!$SIMD 
     732      DO ij=ij_begin,ij_end 
    746733  
    747           n=(j-1)*iim+i 
    748   
    749           gradivu_e(n+u_right,l)=-1/de(n+u_right)*(ne(n,right)*divu_i(n,l)+ ne(n+t_right,left)*divu_i(n+t_right,l) )        
    750  
    751           gradivu_e(n+u_lup,l)=-1/de(n+u_lup)*(ne(n,lup)*divu_i(n,l)+ ne(n+t_lup,rdown)*divu_i(n+t_lup,l))        
    752      
    753           gradivu_e(n+u_ldown,l)=-1/de(n+u_ldown)*(ne(n,ldown)*divu_i(n,l)+ne(n+t_ldown,rup)*divu_i(n+t_ldown,l) ) 
    754  
    755         ENDDO 
     734          gradivu_e(ij+u_right,l)=-1/de(ij+u_right)*(ne(ij,right)*divu_i(ij,l)+ ne(ij+t_right,left)*divu_i(ij+t_right,l) )        
     735 
     736          gradivu_e(ij+u_lup,l)=-1/de(ij+u_lup)*(ne(ij,lup)*divu_i(ij,l)+ ne(ij+t_lup,rdown)*divu_i(ij+t_lup,l))        
     737     
     738          gradivu_e(ij+u_ldown,l)=-1/de(ij+u_ldown)*(ne(ij,ldown)*divu_i(ij,l)+ne(ij+t_ldown,rup)*divu_i(ij+t_ldown,l) ) 
     739 
    756740      ENDDO 
    757741    ENDDO 
    758742 
    759743    DO l=llb,lle 
    760       DO j=jj_begin,jj_end 
    761         DO i=ii_begin,ii_end 
    762           n=(j-1)*iim+i 
    763           gradivu_e(n+u_right,l)=-gradivu_e(n+u_right,l)*cgraddiv 
    764           gradivu_e(n+u_lup,l)=-gradivu_e(n+u_lup,l)*cgraddiv 
    765           gradivu_e(n+u_ldown,l)=-gradivu_e(n+u_ldown,l)*cgraddiv 
    766         ENDDO 
     744!$SIMD 
     745      DO ij=ij_begin,ij_end 
     746          gradivu_e(ij+u_right,l)=-gradivu_e(ij+u_right,l)*cgraddiv 
     747          gradivu_e(ij+u_lup,l)=-gradivu_e(ij+u_lup,l)*cgraddiv 
     748          gradivu_e(ij+u_ldown,l)=-gradivu_e(ij+u_ldown,l)*cgraddiv 
    767749      ENDDO 
    768750    ENDDO 
    … …  
    780762    REAL(rstd)  :: grad_e(3*iim*jjm,llb:lle) 
    781763 
    782     INTEGER :: i,j,n,l 
     764    INTEGER :: ij,l 
    783765 
    784766        
    785767    DO l=llb,lle 
    786       DO j=jj_begin-1,jj_end+1 
    787         DO i=ii_begin-1,ii_end+1 
    788     
    789           n=(j-1)*iim+i 
     768!$SIMD 
     769      DO ij=ij_begin_ext,ij_end_ext 
    790770  
    791           grad_e(n+u_right,l)=-1/de(n+u_right)*(ne(n,right)*theta(n,l)+ ne(n+t_right,left)*theta(n+t_right,l) )         
    792  
    793           grad_e(n+u_lup,l)=-1/de(n+u_lup)*(ne(n,lup)*theta(n,l)+ ne(n+t_lup,rdown)*theta(n+t_lup,l ))        
    794      
    795           grad_e(n+u_ldown,l)=-1/de(n+u_ldown)*(ne(n,ldown)*theta(n,l)+ne(n+t_ldown,rup)*theta(n+t_ldown,l) ) 
    796  
    797         ENDDO 
     771          grad_e(ij+u_right,l)=-1/de(ij+u_right)*(ne(ij,right)*theta(ij,l)+ ne(ij+t_right,left)*theta(ij+t_right,l) )         
     772 
     773          grad_e(ij+u_lup,l)=-1/de(ij+u_lup)*(ne(ij,lup)*theta(ij,l)+ ne(ij+t_lup,rdown)*theta(ij+t_lup,l ))        
     774     
     775          grad_e(ij+u_ldown,l)=-1/de(ij+u_ldown)*(ne(ij,ldown)*theta(ij,l)+ne(ij+t_ldown,rup)*theta(ij+t_ldown,l) ) 
     776 
    798777      ENDDO 
    799778    ENDDO 
    … …  
    801780     
    802781    DO l=llb,lle 
    803       DO j=jj_begin,jj_end 
    804         DO i=ii_begin,ii_end 
    805           n=(j-1)*iim+i 
    806           divgrad_i(n,l)=1./Ai(n)*(ne(n,right)*grad_e(n+u_right,l)*le(n+u_right)  +  & 
    807                              ne(n,rup)*grad_e(n+u_rup,l)*le(n+u_rup)              +  &   
    808                              ne(n,lup)*grad_e(n+u_lup,l)*le(n+u_lup)              +  &   
    809                              ne(n,left)*grad_e(n+u_left,l)*le(n+u_left)           +  &   
    810                              ne(n,ldown)*grad_e(n+u_ldown,l)*le(n+u_ldown)        +  &   
    811                              ne(n,rdown)*grad_e(n+u_rdown,l)*le(n+u_rdown)) 
    812         ENDDO 
     782!$SIMD 
     783      DO ij=jj_begin,ij_end 
     784 
     785          divgrad_i(ij,l)=1./Ai(ij)*(ne(ij,right)*grad_e(ij+u_right,l)*le(ij+u_right)  +  & 
     786                             ne(ij,rup)*grad_e(ij+u_rup,l)*le(ij+u_rup)              +  &   
     787                             ne(ij,lup)*grad_e(ij+u_lup,l)*le(ij+u_lup)              +  &   
     788                             ne(ij,left)*grad_e(ij+u_left,l)*le(ij+u_left)           +  &   
     789                             ne(ij,ldown)*grad_e(ij+u_ldown,l)*le(ij+u_ldown)        +  &   
     790                             ne(ij,rdown)*grad_e(ij+u_rdown,l)*le(ij+u_rdown)) 
    813791      ENDDO 
    814792    ENDDO 
    815793    
    816794    DO l=llb,lle 
    817       DO j=jj_begin,jj_end 
    818         DO i=ii_begin,ii_end 
    819           n=(j-1)*iim+i 
    820           divgrad_i(n,l)=-divgrad_i(n,l)*cdivgrad 
    821         ENDDO 
     795      DO ij=ij_begin,ij_end 
     796          divgrad_i(ij,l)=-divgrad_i(ij,l)*cdivgrad 
    822797      ENDDO 
    823798    ENDDO 
    … …  
    835810    REAL(rstd) :: rot_v(2*iim*jjm,llb:lle) 
    836811 
    837     INTEGER :: i,j,n,l 
     812    INTEGER :: ij,l 
    838813      
    839814    DO l=llb,lle 
    840       DO j=jj_begin-1,jj_end+1 
    841         DO i=ii_begin-1,ii_end+1 
    842           n=(j-1)*iim+i 
    843          
    844           rot_v(n+z_up,l)= 1./Av(n+z_up)*(  ne(n,rup)*ue(n+u_rup,l)*de(n+u_rup)                     & 
    845                                 + ne(n+t_rup,left)*ue(n+t_rup+u_left,l)*de(n+t_rup+u_left)          & 
    846                                 - ne(n,lup)*ue(n+u_lup,l)*de(n+u_lup) )  
     815!$SIMD 
     816      DO ij=ij_begin_ext,ij_end_ext 
     817         
     818          rot_v(ij+z_up,l)= 1./Av(ij+z_up)*(  ne(ij,rup)*ue(ij+u_rup,l)*de(ij+u_rup)                     & 
     819                                + ne(ij+t_rup,left)*ue(ij+t_rup+u_left,l)*de(ij+t_rup+u_left)          & 
     820                                - ne(ij,lup)*ue(ij+u_lup,l)*de(ij+u_lup) )  
    847821                               
    848           rot_v(n+z_down,l) = 1./Av(n+z_down)*( ne(n,ldown)*ue(n+u_ldown,l)*de(n+u_ldown)                 & 
    849                                      + ne(n+t_ldown,right)*ue(n+t_ldown+u_right,l)*de(n+t_ldown+u_right)  & 
    850                                      - ne(n,rdown)*ue(n+u_rdown,l)*de(n+u_rdown) ) 
    851    
    852         ENDDO 
     822          rot_v(ij+z_down,l) = 1./Av(ij+z_down)*( ne(ij,ldown)*ue(ij+u_ldown,l)*de(ij+u_ldown)                 & 
     823                                     + ne(ij+t_ldown,right)*ue(ij+t_ldown+u_right,l)*de(ij+t_ldown+u_right)  & 
     824                                     - ne(ij,rdown)*ue(ij+u_rdown,l)*de(ij+u_rdown) ) 
     825   
    853826      ENDDO 
    854827    ENDDO                               
    855828   
    856829    DO l=llb,lle 
    857       DO j=jj_begin,jj_end 
    858         DO i=ii_begin,ii_end 
    859           n=(j-1)*iim+i 
    860    
    861           gradrot_e(n+u_right,l)=1/le(n+u_right)*ne(n,right)*(rot_v(n+z_rdown,l)-rot_v(n+z_rup,l))  
     830!$SIMD 
     831      DO ij=ij_begin,ij_end 
     832   
     833          gradrot_e(ij+u_right,l)=1/le(ij+u_right)*ne(ij,right)*(rot_v(ij+z_rdown,l)-rot_v(ij+z_rup,l))  
    862834          
    863           gradrot_e(n+u_lup,l)=1/le(n+u_lup)*ne(n,lup)*(rot_v(n+z_up,l)-rot_v(n+z_lup,l))  
    864          
    865           gradrot_e(n+u_ldown,l)=1/le(n+u_ldown)*ne(n,ldown)*(rot_v(n+z_ldown,l)-rot_v(n+z_down,l)) 
    866          
    867         ENDDO 
     835          gradrot_e(ij+u_lup,l)=1/le(ij+u_lup)*ne(ij,lup)*(rot_v(ij+z_up,l)-rot_v(ij+z_lup,l))  
     836         
     837          gradrot_e(ij+u_ldown,l)=1/le(ij+u_ldown)*ne(ij,ldown)*(rot_v(ij+z_ldown,l)-rot_v(ij+z_down,l)) 
     838         
    868839      ENDDO 
    869840    ENDDO 
    870841 
    871842    DO l=llb,lle 
    872       DO j=jj_begin,jj_end 
    873         DO i=ii_begin,ii_end 
    874           n=(j-1)*iim+i 
    875      
    876           gradrot_e(n+u_right,l)=-gradrot_e(n+u_right,l)*cgradrot        
    877           gradrot_e(n+u_lup,l)=-gradrot_e(n+u_lup,l)*cgradrot 
    878           gradrot_e(n+u_ldown,l)=-gradrot_e(n+u_ldown,l)*cgradrot 
    879          
    880         ENDDO 
     843!$SIMD 
     844      DO ij=ij_begin,ij_end 
     845     
     846          gradrot_e(ij+u_right,l)=-gradrot_e(ij+u_right,l)*cgradrot        
     847          gradrot_e(ij+u_lup,l)=-gradrot_e(ij+u_lup,l)*cgradrot 
     848          gradrot_e(ij+u_ldown,l)=-gradrot_e(ij+u_ldown,l)*cgradrot 
     849         
    881850      ENDDO 
    882851    ENDDO   

  • codes/icosagcm/trunk/src/domain.f90

    r151 r174  
    7272  USE grid_param 
    7373  USE mpipara 
     74  USE IOIPSL 
    7475  IMPLICIT NONE 
    7576  INTEGER :: ind,nf,ni,nj,i,j 
    7677  INTEGER :: quotient, rest 
     78  INTEGER :: halo_i,halo_j 
     79   
    7780  TYPE(t_domain),POINTER :: d 
    7881  
    … …  
    8184    ALLOCATE(domglo_rank(ndomain_glo)) 
    8285    ALLOCATE(domglo_loc_ind(ndomain_glo)) 
     86    
     87    halo_i=0 
     88    CALL getin("halo_i",halo_i) 
     89    halo_j=1 
     90    CALL getin("halo_j",halo_j) 
    8391 
    8492    ind=0 
    … …  
    122130 
    123131 
    124           d%iim=d%ii_nb+2*halo 
    125           d%jjm=d%jj_nb+2*halo 
    126           d%ii_begin=halo+1 
    127           d%jj_begin=halo+1 
     132          d%iim=d%ii_nb+2*halo  + halo_i*2 
     133          d%jjm=d%jj_nb+2*halo  + halo_j*2 
     134          d%ii_begin=halo+1  + halo_i 
     135          d%jj_begin=halo+1  + halo_j 
    128136          d%ii_end=d%ii_begin+d%ii_nb-1 
    129137          d%jj_end=d%jj_begin+d%jj_nb-1 

  • codes/icosagcm/trunk/src/domain_param.f90

    r95 r174  
    33INTEGER :: nsplit_i 
    44INTEGER :: nsplit_j 
    5 INTEGER :: halo=3 
     5INTEGER :: halo=1 
    66 
    77INTEGER, PARAMETER :: default_nsplit_i=1 

  • codes/icosagcm/trunk/src/timeloop_gcm.f90

    r171 r174  
    336336             ps=f_ps(ind) ; dps=f_dps(ind) ;               
    337337             IF (omp_first) THEN 
    338                 DO j=jj_begin,jj_end 
    339                    DO i=ii_begin,ii_end 
    340                       ij=(j-1)*iim+i 
    341                       ps(ij)=ps(ij)+dt*dps(ij) 
    342                    ENDDO 
     338!$SIMD 
     339                DO ij=ij_begin,ij_end 
     340                    ps(ij)=ps(ij)+dt*dps(ij) 
    343341                ENDDO 
    344342             ENDIF  
    … …  
    346344             mass=f_mass(ind) ; dmass=f_dmass(ind) ;               
    347345             DO l=1,llm 
    348                 DO j=jj_begin,jj_end 
    349                    DO i=ii_begin,ii_end 
    350                       ij=(j-1)*iim+i 
    351                       mass(ij,l)=mass(ij,l)+dt*dmass(ij,l) 
    352                    ENDDO 
     346!$SIMD 
     347                DO ij=ij_begin,ij_end 
     348                    mass(ij,l)=mass(ij,l)+dt*dmass(ij,l) 
    353349                ENDDO 
    354350             END DO 
    … …  
    364360 
    365361       DO l=ll_begin,ll_end 
    366          DO j=jj_begin,jj_end 
    367            DO i=ii_begin,ii_end 
    368              ij=(j-1)*iim+i 
     362!$SIMD 
     363         DO ij=ij_begin,ij_end 
    369364             u(ij+u_right,l)=u(ij+u_right,l)+dt*du(ij+u_right,l) 
    370365             u(ij+u_lup,l)=u(ij+u_lup,l)+dt*du(ij+u_lup,l) 
    371366             u(ij+u_ldown,l)=u(ij+u_ldown,l)+dt*du(ij+u_ldown,l) 
    372367             theta_rhodz(ij,l)=theta_rhodz(ij,l)+dt*dtheta_rhodz(ij,l) 
    373            ENDDO 
    374368         ENDDO 
    375369       ENDDO 
    … …  
    401395                
    402396               IF (stage==1) THEN ! first stage : save model state in XXm1 
    403                   DO j=jj_begin,jj_end 
    404                      DO i=ii_begin,ii_end 
    405                         ij=(j-1)*iim+i 
    406                         psm1(ij)=ps(ij) 
    407                      ENDDO 
    408                   ENDDO 
     397!$SIMD 
     398                 DO ij=ij_begin,ij_end 
     399                   psm1(ij)=ps(ij) 
     400                 ENDDO 
    409401               ENDIF 
    410402             
    411403               ! updates are of the form x1 := x0 + tau*f(x1) 
    412                DO j=jj_begin,jj_end 
    413                   DO i=ii_begin,ii_end 
    414                      ij=(j-1)*iim+i 
    415                      ps(ij)=psm1(ij)+tau*dps(ij) 
    416                   ENDDO 
     404!$SIMD 
     405               DO ij=ij_begin,ij_end 
     406                   ps(ij)=psm1(ij)+tau*dps(ij) 
    417407               ENDDO 
    418408            ENDDO 
    … …  
    428418            IF (stage==1) THEN ! first stage : save model state in XXm1 
    429419               DO l=ll_begin,ll_end 
    430                   DO j=jj_begin,jj_end 
    431                      DO i=ii_begin,ii_end 
    432                         ij=(j-1)*iim+i 
    433                         massm1(ij,l)=mass(ij,l) 
    434                      ENDDO 
    435                   ENDDO 
     420!$SIMD 
     421                 DO ij=ij_begin,ij_end 
     422                    massm1(ij,l)=mass(ij,l) 
     423                 ENDDO 
    436424               ENDDO 
    437425            END IF 
    … …  
    439427            ! updates are of the form x1 := x0 + tau*f(x1) 
    440428            DO l=ll_begin,ll_end 
    441                DO j=jj_begin,jj_end 
    442                   DO i=ii_begin,ii_end 
    443                      ij=(j-1)*iim+i 
    444                      mass(ij,l)=massm1(ij,l)+tau*dmass(ij,l) 
    445                   ENDDO 
     429!$SIMD 
     430               DO ij=ij_begin,ij_end 
     431                   mass(ij,l)=massm1(ij,l)+tau*dmass(ij,l) 
    446432               ENDDO 
    447433            ENDDO 
    … …  
    462448         IF (stage==1) THEN ! first stage : save model state in XXm1 
    463449           DO l=ll_begin,ll_end 
    464              DO j=jj_begin,jj_end 
    465                DO i=ii_begin,ii_end 
    466                  ij=(j-1)*iim+i 
     450!$SIMD 
     451                DO ij=ij_begin,ij_end 
    467452                 um1(ij+u_right,l)=u(ij+u_right,l) 
    468453                 um1(ij+u_lup,l)=u(ij+u_lup,l) 
    469454                 um1(ij+u_ldown,l)=u(ij+u_ldown,l) 
    470455                 theta_rhodzm1(ij,l)=theta_rhodz(ij,l) 
    471                ENDDO 
    472456             ENDDO 
    473457           ENDDO 
    … …  
    475459 
    476460         DO l=ll_begin,ll_end 
    477            DO j=jj_begin,jj_end 
    478              DO i=ii_begin,ii_end 
    479                ij=(j-1)*iim+i 
     461!$SIMD 
     462           DO ij=ij_begin,ij_end 
    480463               u(ij+u_right,l)=um1(ij+u_right,l)+tau*du(ij+u_right,l) 
    481464               u(ij+u_lup,l)=um1(ij+u_lup,l)+tau*du(ij+u_lup,l) 
    482465               u(ij+u_ldown,l)=um1(ij+u_ldown,l)+tau*du(ij+u_ldown,l) 
    483466               theta_rhodz(ij,l)=theta_rhodzm1(ij,l)+tau*dtheta_rhodz(ij,l) 
    484              ENDDO 
    485467           ENDDO 
    486468         ENDDO 
    … …  
    564546 
    565547       DO l=ll_begin,ll_end 
    566          DO j=jj_begin-1,jj_end+1 
    567            DO i=ii_begin-1,ii_end+1 
    568              ij=(j-1)*iim+i 
     548!$SIMD 
     549         DO ij=ij_begin_ext,ij_end_ext 
    569550             hfluxt(ij+u_right,l) = tau*hflux(ij+u_right,l) 
    570551             hfluxt(ij+u_lup,l) = tau*hflux(ij+u_lup,l) 
    571552             hfluxt(ij+u_ldown,l) = tau*hflux(ij+u_ldown,l) 
    572            ENDDO 
    573553         ENDDO 
    574554       ENDDO 
    … …  
    576556       IF(caldyn_eta==eta_mass) THEN ! no need for vertical fluxes if eta_lag 
    577557          DO l=ll_begin,ll_endp1 
    578              DO j=jj_begin,jj_end 
    579                 DO i=ii_begin,ii_end 
    580                    ij=(j-1)*iim+i 
    581                    wfluxt(ij,l) = tau*wflux(ij,l) 
    582                 ENDDO 
     558!$SIMD 
     559             DO ij=ij_begin,ij_end 
     560                wfluxt(ij,l) = tau*wflux(ij,l) 
    583561             ENDDO 
    584562          ENDDO 
    … …  
    588566 
    589567       DO l=ll_begin,ll_end 
    590          DO j=jj_begin-1,jj_end+1 
    591            DO i=ii_begin-1,ii_end+1 
    592              ij=(j-1)*iim+i 
     568!$SIMD 
     569         DO ij=ij_begin_ext,ij_end_ext 
    593570             hfluxt(ij+u_right,l) = hfluxt(ij+u_right,l)+tau*hflux(ij+u_right,l) 
    594571             hfluxt(ij+u_lup,l) = hfluxt(ij+u_lup,l)+tau*hflux(ij+u_lup,l) 
    595572             hfluxt(ij+u_ldown,l) = hfluxt(ij+u_ldown,l)+tau*hflux(ij+u_ldown,l) 
    596            ENDDO 
    597573         ENDDO 
    598574       ENDDO 
    … …  
    600576       IF(caldyn_eta==eta_mass) THEN ! no need for vertical fluxes if eta_lag 
    601577          DO l=ll_begin,ll_endp1 
    602              DO j=jj_begin,jj_end 
    603                 DO i=ii_begin,ii_end 
    604                    ij=(j-1)*iim+i 
     578!$SIMD 
     579             DO ij=ij_begin,ij_end 
    605580                   wfluxt(ij,l) = wfluxt(ij,l)+tau*wflux(ij,l) 
    606                 ENDDO 
    607581             ENDDO 
    608582          ENDDO 
    … …  
    613587  END SUBROUTINE accumulate_fluxes 
    614588   
    615   FUNCTION maxval_i(p) 
    616     USE icosa 
    617     IMPLICIT NONE 
    618     REAL(rstd), DIMENSION(iim*jjm) :: p 
    619     REAL(rstd) :: maxval_i 
    620     INTEGER :: j, ij 
    621      
    622     maxval_i=p((jj_begin-1)*iim+ii_begin) 
    623      
    624     DO j=jj_begin-1,jj_end+1 
    625        ij=(j-1)*iim 
    626        maxval_i = MAX(maxval_i, MAXVAL(p(ij+ii_begin:ij+ii_end))) 
    627     END DO 
    628   END FUNCTION maxval_i 
    629  
    630   FUNCTION maxval_ik(p) 
    631     USE icosa 
    632     IMPLICIT NONE 
    633     REAL(rstd) :: p(iim*jjm, llm) 
    634     REAL(rstd) :: maxval_ik(llm) 
    635     INTEGER :: l,j, ij 
    636      
    637     DO l=1,llm 
    638        maxval_ik(l)=p((jj_begin-1)*iim+ii_begin,l) 
    639        DO j=jj_begin-1,jj_end+1 
    640           ij=(j-1)*iim 
    641           maxval_ik(l) = MAX(maxval_ik(l), MAXVAL(p(ij+ii_begin:ij+ii_end,l))) 
    642        END DO 
    643     END DO 
    644   END FUNCTION maxval_ik 
     589!  FUNCTION maxval_i(p) 
     590!    USE icosa 
     591!    IMPLICIT NONE 
     592!    REAL(rstd), DIMENSION(iim*jjm) :: p 
     593!    REAL(rstd) :: maxval_i 
     594!    INTEGER :: j, ij 
     595!     
     596!    maxval_i=p((jj_begin-1)*iim+ii_begin) 
     597!     
     598!    DO j=jj_begin-1,jj_end+1 
     599!       ij=(j-1)*iim 
     600!       maxval_i = MAX(maxval_i, MAXVAL(p(ij+ii_begin:ij+ii_end))) 
     601!    END DO 
     602!  END FUNCTION maxval_i 
     603 
     604!  FUNCTION maxval_ik(p) 
     605!    USE icosa 
     606!    IMPLICIT NONE 
     607!    REAL(rstd) :: p(iim*jjm, llm) 
     608!    REAL(rstd) :: maxval_ik(llm) 
     609!    INTEGER :: l,j, ij 
     610!     
     611!    DO l=1,llm 
     612!       maxval_ik(l)=p((jj_begin-1)*iim+ii_begin,l) 
     613!       DO j=jj_begin-1,jj_end+1 
     614!          ij=(j-1)*iim 
     615!          maxval_ik(l) = MAX(maxval_ik(l), MAXVAL(p(ij+ii_begin:ij+ii_end,l))) 
     616!       END DO 
     617!    END DO 
     618!  END FUNCTION maxval_ik 
    645619 
    646620END MODULE timeloop_gcm_mod 
Note: See TracChangeset for help on using the changeset viewer.