Ignore:
Timestamp:
06/13/19 16:45:42 (5 years ago)
Author:
adurocher
Message:

trunk : Fixed aliasing issue in compute_divgrad/gradiv/gradrot in dissip

Added compute_*_inplace to avoid aliasing

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r899 r900  
    199199          u=f_u(ind) 
    200200          du=f_du(ind) 
    201           CALL compute_gradiv(u,du,1,1) 
    202           u=du 
     201          CALL compute_gradiv_inplace(u,1,1) 
     202          du=u 
    203203        ENDDO 
    204204      ENDDO 
    … …  
    292292          u=f_u(ind) 
    293293          du=f_du(ind) 
    294           CALL compute_gradrot(u,du,1,1) 
    295           u=du 
     294          CALL compute_gradrot_inplace(u,1,1) 
     295          du=u 
    296296        ENDDO 
    297297      ENDDO 
    … …  
    380380          theta=f_theta(ind) 
    381381          dtheta=f_dtheta(ind) 
    382           CALL compute_divgrad(theta,dtheta,1,1) 
    383           theta=dtheta 
     382          CALL compute_divgrad_inplace(theta,1,1) 
     383          dtheta=theta 
    384384        ENDDO 
    385385      ENDDO 
    … …  
    673673        CALL swap_geometry(ind) 
    674674        due=f_due(ind)  
    675         CALL compute_gradiv(due,due,ll_begin,ll_end) 
     675        CALL compute_gradiv_inplace(due,ll_begin,ll_end) 
    676676      ENDDO 
    677677    ENDDO 
    … …  
    722722        CALL swap_geometry(ind) 
    723723        due=f_due(ind)  
    724         CALL compute_gradrot(due,due,ll_begin,ll_end) 
     724        CALL compute_gradrot_inplace(due,ll_begin,ll_end) 
    725725      ENDDO 
    726726 
    … …  
    763763        CALL swap_geometry(ind) 
    764764        dtheta=f_dtheta(ind)  
    765         CALL compute_divgrad(dtheta,dtheta,ll_begin,ll_end) 
     765        CALL compute_divgrad_inplace(dtheta,ll_begin,ll_end) 
    766766      ENDDO 
    767767 
    … …  
    814814        CALL swap_geometry(ind) 
    815815        dtheta_rhodz=f_dtheta_rhodz(ind)  
    816         CALL compute_divgrad(dtheta_rhodz,dtheta_rhodz,ll_begin,ll_end) 
     816        CALL compute_divgrad_inplace(dtheta_rhodz,ll_begin,ll_end) 
    817817      ENDDO 
    818818 
    … …  
    846846    REAL(rstd),INTENT(IN)  :: ue(iim*3*jjm,llm) 
    847847    REAL(rstd),INTENT(OUT) :: gradivu_e(iim*3*jjm,llm) 
     848     
     849    gradivu_e = ue 
     850    CALL compute_gradiv_inplace(gradivu_e, llb, lle) 
     851  END SUBROUTINE compute_gradiv 
     852   
     853  SUBROUTINE compute_gradiv_inplace(ue_gradivu_e,llb,lle) 
     854  USE icosa 
     855  IMPLICIT NONE 
     856    INTEGER,INTENT(IN)     :: llb 
     857    INTEGER,INTENT(IN)     :: lle 
     858    REAL(rstd),INTENT(INOUT) :: ue_gradivu_e(iim*3*jjm,llm) 
    848859    REAL(rstd) :: divu_i(iim*jjm,llb:lle) 
    849860     
    850861    INTEGER :: ij,l 
    851862 
     863    DO l=llb,lle 
     864      !DIR$ SIMD 
     865      DO ij=ij_begin,ij_end 
     866          divu_i(ij,l)=1./Ai(ij)*(ne(ij,right)*ue_gradivu_e(ij+u_right,l)*le(ij+u_right)  +  & 
     867                             ne(ij,rup)*ue_gradivu_e(ij+u_rup,l)*le(ij+u_rup)        +  &   
     868                             ne(ij,lup)*ue_gradivu_e(ij+u_lup,l)*le(ij+u_lup)        +  &   
     869                             ne(ij,left)*ue_gradivu_e(ij+u_left,l)*le(ij+u_left)     +  &   
     870                             ne(ij,ldown)*ue_gradivu_e(ij+u_ldown,l)*le(ij+u_ldown)  +  &   
     871                             ne(ij,rdown)*ue_gradivu_e(ij+u_rdown,l)*le(ij+u_rdown)) 
     872      ENDDO 
     873    ENDDO 
     874     
     875    DO l=llb,lle 
     876      !DIR$ SIMD 
     877      DO ij=ij_begin,ij_end 
     878          ue_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) )  
     879          ue_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))    
     880          ue_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) ) 
     881      ENDDO 
     882    ENDDO 
     883 
     884    DO l=llb,lle 
     885      !DIR$ SIMD 
     886      DO ij=ij_begin,ij_end 
     887          ue_gradivu_e(ij+u_right,l)=-ue_gradivu_e(ij+u_right,l)*cgraddiv 
     888          ue_gradivu_e(ij+u_lup,l)=-ue_gradivu_e(ij+u_lup,l)*cgraddiv 
     889          ue_gradivu_e(ij+u_ldown,l)=-ue_gradivu_e(ij+u_ldown,l)*cgraddiv 
     890      ENDDO 
     891    ENDDO    
     892  END SUBROUTINE compute_gradiv_inplace 
     893   
     894  SUBROUTINE compute_divgrad(theta,divgrad_i,llb,lle) 
     895  USE icosa 
     896  IMPLICIT NONE 
     897    INTEGER,INTENT(IN)     :: llb 
     898    INTEGER,INTENT(IN)     :: lle 
     899    REAL(rstd),INTENT(IN) :: theta(iim*jjm,1:lle) 
     900    REAL(rstd),INTENT(OUT) :: divgrad_i(iim*jjm,1:lle) 
     901     
     902    divgrad_i = theta 
     903    CALL compute_divgrad_inplace(divgrad_i, llb, lle) 
     904  END SUBROUTINE 
     905   
     906  SUBROUTINE compute_divgrad_inplace(theta_divgrad_i,llb,lle) 
     907  USE icosa 
     908  IMPLICIT NONE 
     909    INTEGER,INTENT(IN)     :: llb 
     910    INTEGER,INTENT(IN)     :: lle 
     911    REAL(rstd),INTENT(INOUT) :: theta_divgrad_i(iim*jjm,1:lle) 
     912    REAL(rstd)  :: grad_e(3*iim*jjm,llb:lle) 
     913 
     914    INTEGER :: ij,l     
     915    DO l=llb,lle 
     916      !DIR$ SIMD 
     917      DO ij=ij_begin_ext,ij_end_ext  
     918          grad_e(ij+u_right,l)=-1/de(ij+u_right)*(ne(ij,right)*theta_divgrad_i(ij,l)+ ne(ij+t_right,left)*theta_divgrad_i(ij+t_right,l) ) 
     919          grad_e(ij+u_lup,l)=-1/de(ij+u_lup)*(ne(ij,lup)*theta_divgrad_i(ij,l)+ ne(ij+t_lup,rdown)*theta_divgrad_i(ij+t_lup,l )) 
     920          grad_e(ij+u_ldown,l)=-1/de(ij+u_ldown)*(ne(ij,ldown)*theta_divgrad_i(ij,l)+ne(ij+t_ldown,rup)*theta_divgrad_i(ij+t_ldown,l) ) 
     921      ENDDO 
     922    ENDDO     
     923     
     924    DO l=llb,lle 
     925      !DIR$ SIMD 
     926      DO ij=ij_begin,ij_end 
     927          theta_divgrad_i(ij,l)=1./Ai(ij)*(ne(ij,right)*grad_e(ij+u_right,l)*le(ij+u_right)  +  & 
     928                             ne(ij,rup)*grad_e(ij+u_rup,l)*le(ij+u_rup)              +  & 
     929                             ne(ij,lup)*grad_e(ij+u_lup,l)*le(ij+u_lup)              +  & 
     930                             ne(ij,left)*grad_e(ij+u_left,l)*le(ij+u_left)           +  & 
     931                             ne(ij,ldown)*grad_e(ij+u_ldown,l)*le(ij+u_ldown)        +  &  
     932                             ne(ij,rdown)*grad_e(ij+u_rdown,l)*le(ij+u_rdown)) 
     933      ENDDO 
     934    ENDDO 
     935    
     936    DO l=llb,lle 
     937      DO ij=ij_begin,ij_end 
     938          theta_divgrad_i(ij,l)=-theta_divgrad_i(ij,l)*cdivgrad 
     939      ENDDO 
     940    ENDDO 
     941  END SUBROUTINE compute_divgrad_inplace 
     942 
     943  SUBROUTINE compute_gradrot(ue,gradrot_e,llb,lle) 
     944    INTEGER,INTENT(IN)     :: llb 
     945    INTEGER,INTENT(IN)     :: lle 
     946    REAL(rstd),INTENT(IN) :: ue(iim*3*jjm,lle) 
     947    REAL(rstd),INTENT(OUT) :: gradrot_e(iim*3*jjm,lle) 
     948     
     949    gradrot_e = ue 
     950    CALL compute_gradrot_inplace(gradrot_e,llb,lle) 
     951  END SUBROUTINE 
     952     
     953  SUBROUTINE compute_gradrot_inplace(ue_gradrot_e,llb,lle) 
     954  USE icosa 
     955  IMPLICIT NONE 
     956    INTEGER,INTENT(IN)     :: llb 
     957    INTEGER,INTENT(IN)     :: lle 
     958    REAL(rstd),INTENT(INOUT) :: ue_gradrot_e(iim*3*jjm,lle) 
     959    REAL(rstd) :: rot_v(2*iim*jjm,llb:lle) 
     960 
     961    INTEGER :: ij,l 
     962 
     963    DO l=llb,lle 
     964!DIR$ SIMD 
     965      DO ij=ij_begin_ext,ij_end_ext         
     966          rot_v(ij+z_up,l)= 1./Av(ij+z_up)*(  ne(ij,rup)*ue_gradrot_e(ij+u_rup,l)*de(ij+u_rup) & 
     967                                + ne(ij+t_rup,left)*ue_gradrot_e(ij+t_rup+u_left,l)*de(ij+t_rup+u_left) & 
     968                                - ne(ij,lup)*ue_gradrot_e(ij+u_lup,l)*de(ij+u_lup) )                                
     969          rot_v(ij+z_down,l) = 1./Av(ij+z_down)*( ne(ij,ldown)*ue_gradrot_e(ij+u_ldown,l)*de(ij+u_ldown) & 
     970                                     + ne(ij+t_ldown,right)*ue_gradrot_e(ij+t_ldown+u_right,l)*de(ij+t_ldown+u_right) & 
     971                                     - ne(ij,rdown)*ue_gradrot_e(ij+u_rdown,l)*de(ij+u_rdown) ) 
     972      ENDDO 
     973    ENDDO                               
     974   
    852975    DO l=llb,lle 
    853976!DIR$ SIMD 
    854977      DO ij=ij_begin,ij_end 
    855           divu_i(ij,l)=1./Ai(ij)*(ne(ij,right)*ue(ij+u_right,l)*le(ij+u_right)  +  & 
    856                              ne(ij,rup)*ue(ij+u_rup,l)*le(ij+u_rup)        +  &   
    857                              ne(ij,lup)*ue(ij+u_lup,l)*le(ij+u_lup)        +  &   
    858                              ne(ij,left)*ue(ij+u_left,l)*le(ij+u_left)     +  &   
    859                              ne(ij,ldown)*ue(ij+u_ldown,l)*le(ij+u_ldown)  +  &   
    860                              ne(ij,rdown)*ue(ij+u_rdown,l)*le(ij+u_rdown)) 
    861       ENDDO 
    862     ENDDO 
    863      
     978          ue_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))        
     979          ue_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))         
     980          ue_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)) 
     981      ENDDO 
     982    ENDDO 
     983 
    864984    DO l=llb,lle 
    865985!DIR$ SIMD 
    866       DO ij=ij_begin,ij_end 
    867   
    868           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) )        
    869  
    870           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))        
    871      
    872           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) ) 
    873  
    874       ENDDO 
    875     ENDDO 
    876  
    877     DO l=llb,lle 
    878 !DIR$ SIMD 
    879       DO ij=ij_begin,ij_end 
    880           gradivu_e(ij+u_right,l)=-gradivu_e(ij+u_right,l)*cgraddiv 
    881           gradivu_e(ij+u_lup,l)=-gradivu_e(ij+u_lup,l)*cgraddiv 
    882           gradivu_e(ij+u_ldown,l)=-gradivu_e(ij+u_ldown,l)*cgraddiv 
    883       ENDDO 
    884     ENDDO 
    885  
    886     
    887   END SUBROUTINE compute_gradiv 
    888    
    889   SUBROUTINE compute_divgrad(theta,divgrad_i,llb,lle) 
    890   USE icosa 
    891   IMPLICIT NONE 
    892     INTEGER,INTENT(IN)     :: llb 
    893     INTEGER,INTENT(IN)     :: lle 
    894     REAL(rstd),INTENT(IN)  :: theta(iim*jjm,llm) 
    895     REAL(rstd),INTENT(OUT) :: divgrad_i(iim*jjm,llm) 
    896     REAL(rstd)  :: grad_e(3*iim*jjm,llb:lle) 
    897  
    898     INTEGER :: ij,l 
    899  
    900         
    901     DO l=llb,lle 
    902 !DIR$ SIMD 
    903       DO ij=ij_begin_ext,ij_end_ext 
    904   
    905           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) )         
    906  
    907           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 ))        
    908      
    909           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) ) 
    910  
    911       ENDDO 
    912     ENDDO 
    913      
    914      
    915     DO l=llb,lle 
    916 !DIR$ SIMD 
    917       DO ij=ij_begin,ij_end 
    918  
    919           divgrad_i(ij,l)=1./Ai(ij)*(ne(ij,right)*grad_e(ij+u_right,l)*le(ij+u_right)  +  & 
    920                              ne(ij,rup)*grad_e(ij+u_rup,l)*le(ij+u_rup)              +  &   
    921                              ne(ij,lup)*grad_e(ij+u_lup,l)*le(ij+u_lup)              +  &   
    922                              ne(ij,left)*grad_e(ij+u_left,l)*le(ij+u_left)           +  &   
    923                              ne(ij,ldown)*grad_e(ij+u_ldown,l)*le(ij+u_ldown)        +  &   
    924                              ne(ij,rdown)*grad_e(ij+u_rdown,l)*le(ij+u_rdown)) 
    925       ENDDO 
    926     ENDDO 
    927     
    928     DO l=llb,lle 
    929       DO ij=ij_begin,ij_end 
    930           divgrad_i(ij,l)=-divgrad_i(ij,l)*cdivgrad 
    931       ENDDO 
    932     ENDDO 
    933  
    934   END SUBROUTINE compute_divgrad 
    935  
    936      
    937   SUBROUTINE compute_gradrot(ue,gradrot_e,llb,lle) 
    938   USE icosa 
    939   IMPLICIT NONE 
    940     INTEGER,INTENT(IN)     :: llb 
    941     INTEGER,INTENT(IN)     :: lle 
    942     REAL(rstd),INTENT(IN)  :: ue(iim*3*jjm,llm) 
    943     REAL(rstd),INTENT(OUT) :: gradrot_e(iim*3*jjm,llm) 
    944     REAL(rstd) :: rot_v(2*iim*jjm,llb:lle) 
    945  
    946     INTEGER :: ij,l 
    947       
    948     DO l=llb,lle 
    949 !DIR$ SIMD 
    950       DO ij=ij_begin_ext,ij_end_ext 
    951          
    952           rot_v(ij+z_up,l)= 1./Av(ij+z_up)*(  ne(ij,rup)*ue(ij+u_rup,l)*de(ij+u_rup)                     & 
    953                                 + ne(ij+t_rup,left)*ue(ij+t_rup+u_left,l)*de(ij+t_rup+u_left)          & 
    954                                 - ne(ij,lup)*ue(ij+u_lup,l)*de(ij+u_lup) )  
    955                                
    956           rot_v(ij+z_down,l) = 1./Av(ij+z_down)*( ne(ij,ldown)*ue(ij+u_ldown,l)*de(ij+u_ldown)                 & 
    957                                      + ne(ij+t_ldown,right)*ue(ij+t_ldown+u_right,l)*de(ij+t_ldown+u_right)  & 
    958                                      - ne(ij,rdown)*ue(ij+u_rdown,l)*de(ij+u_rdown) ) 
    959    
    960       ENDDO 
    961     ENDDO                               
    962    
    963     DO l=llb,lle 
    964 !DIR$ SIMD 
    965       DO ij=ij_begin,ij_end 
    966    
    967           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))  
    968           
    969           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))  
    970          
    971           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)) 
    972          
    973       ENDDO 
    974     ENDDO 
    975  
    976     DO l=llb,lle 
    977 !DIR$ SIMD 
    978       DO ij=ij_begin,ij_end 
    979      
    980           gradrot_e(ij+u_right,l)=-gradrot_e(ij+u_right,l)*cgradrot        
    981           gradrot_e(ij+u_lup,l)=-gradrot_e(ij+u_lup,l)*cgradrot 
    982           gradrot_e(ij+u_ldown,l)=-gradrot_e(ij+u_ldown,l)*cgradrot 
    983          
     986      DO ij=ij_begin,ij_end     
     987          ue_gradrot_e(ij+u_right,l)=-ue_gradrot_e(ij+u_right,l)*cgradrot        
     988          ue_gradrot_e(ij+u_lup,l)=-ue_gradrot_e(ij+u_lup,l)*cgradrot 
     989          ue_gradrot_e(ij+u_ldown,l)=-ue_gradrot_e(ij+u_ldown,l)*cgradrot         
    984990      ENDDO 
    985991    ENDDO   
    986     
    987   END SUBROUTINE compute_gradrot 
     992  END SUBROUTINE compute_gradrot_inplace 
    988993 
    989994 
Note: See TracChangeset for help on using the changeset viewer.