Ignore:
Timestamp:
01/09/14 09:56:11 (10 years ago)
Author:
ymipsl
Message:

Add new openMP parallelism based on distribution of domains on threads. There is no more limitation of number of threads by MPI process.

YM

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r174 r186  
    4949  !======================================================================================= 
    5050 
    51   SUBROUTINE compute_gradq3d(qi,one_over_sqrt_leng,gradq3d) 
     51  SUBROUTINE compute_gradq3d(qi_in,one_over_sqrt_leng_in,gradq3d_out,xyz_i,xyz_v) 
    5252    USE trace 
    5353    USE omp_para 
    5454    IMPLICIT NONE 
    55     REAL(rstd),INTENT(IN)  :: qi(iim*jjm,llm) 
    56     REAL(rstd),INTENT(IN)  :: one_over_sqrt_leng(iim*jjm) 
    57     REAL(rstd),INTENT(OUT) :: gradq3d(iim*jjm,llm,3)  
     55    REAL(rstd),INTENT(IN)  :: qi_in(iim*jjm,llm) 
     56    REAL(rstd),INTENT(IN)  :: one_over_sqrt_leng_in(iim*jjm) 
     57    REAL(rstd),INTENT(IN)  :: xyz_i(iim*jjm,3) 
     58    REAL(rstd),INTENT(IN)  :: xyz_v(2*iim*jjm,3) 
     59    REAL(rstd),INTENT(OUT) :: gradq3d_out(iim*jjm,llm,3)  
    5860    REAL(rstd) :: maxq,minq,minq_c,maxq_c  
    5961    REAL(rstd) :: alphamx,alphami,alpha ,maggrd 
    … …  
    6365    REAL(rstd) :: gradtri(2*iim*jjm,llm,3)  
    6466    INTEGER :: ij,k,ind,l 
    65  
    66     CALL trace_start("compute_gradq3d") 
     67    REAL(rstd)  :: qi(iim*jjm,llm) 
     68    REAL(rstd)  :: one_over_sqrt_leng(iim*jjm) 
     69    REAL(rstd) :: gradq3d(iim*jjm,llm,3)  
     70    REAL(rstd) :: detx,dety,detz,det 
     71    REAL(rstd) :: A(3,3), a11,a12,a13,a21,a22,a23,a31,a32,a33 
     72    REAL(rstd) :: x1,x2,x3 
     73    REAL(rstd) :: dq(3) 
     74 
     75    qi=qi_in 
     76    one_over_sqrt_leng=one_over_sqrt_leng_in 
     77     
     78    CALL trace_start("compute_gradq3d1") 
    6779 
    6880    ! TODO : precompute ar, drop arr as output argument of gradq ? 
    … …  
    7183    ! Compute gradient at triangles solving a linear system 
    7284    ! arr = area of triangle joining centroids of hexagons 
     85!     DO l = ll_begin,ll_end  
     86!!$SIMD 
     87!      DO ij=ij_begin_ext,ij_end_ext 
     88!!             CALL gradq(ij,l,ij+t_rup,ij+t_lup,ij+z_up,qi,gradtri(ij+z_up,l,:),arr(ij+z_up)) 
     89!!             CALL gradq(ij,l,ij+t_ldown,ij+t_rdown,ij+z_down,qi,gradtri(ij+z_down,l,:),arr(ij+z_down)) 
     90!             CALL gradq(ij,l,ij+t_rup,ij+t_lup,ij+z_up,qi,gradtri(ij+z_up,l,1),gradtri(ij+z_up,l,2),gradtri(ij+z_up,l,3),arr(ij+z_up)) 
     91!             CALL gradq(ij,l,ij+t_ldown,ij+t_rdown,ij+z_down,qi,gradtri(ij+z_down,l,1),gradtri(ij+z_down,l,2),gradtri(ij+z_down,l,3),arr(ij+z_down)) 
     92!       END DO 
     93!    END DO 
     94 
    7395     DO l = ll_begin,ll_end  
    7496!$SIMD 
    7597      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)) 
    78        END DO 
     98!       CALL gradq(ij,l,ij+t_rup,ij+t_lup,ij+z_up,qi,gradtri(ij+z_up,l,1),gradtri(ij+z_up,l,2),gradtri(ij+z_up,l,3),arr(ij+z_up)) 
     99 
     100         
     101        A(1,1)=xyz_i(ij+t_rup,1)-xyz_i(ij,1);  A(1,2)=xyz_i(ij+t_rup,2)-xyz_i(ij,2); A(1,3)=xyz_i(ij+t_rup,3)-xyz_i(ij,3)  
     102        A(2,1)=xyz_i(ij+t_lup,1)-xyz_i(ij,1);  A(2,2)=xyz_i(ij+t_lup,2)-xyz_i(ij,2); A(2,3)=xyz_i(ij+t_lup,3)-xyz_i(ij,3)  
     103        A(3,1)=xyz_v(ij+z_up,1);               A(3,2)= xyz_v(ij+z_up,2);             A(3,3)=xyz_v(ij+z_up,3) 
     104     
     105        dq(1) = qi(ij+t_rup,l)-qi(ij,l) 
     106        dq(2) = qi(ij+t_lup,l)-qi(ij,l) 
     107        dq(3) = 0.0 
     108 
     109 
     110!        CALL determinant(A(1,1),A(2,1),A(3,1),A(1,2),A(2,2),A(3,2),A(1,3),A(2,3),A(3,3),det) 
     111 
     112         a11=A(1,1) ; a12=A(2,1) ; a13=A(3,1) 
     113         a21=A(1,2) ; a22=A(2,2) ; a23=A(3,2) 
     114         a31=A(1,3) ; a32=A(2,3) ; a33=A(3,3) 
     115 
     116         x1 =  a11 * (a22 * a33 - a23 * a32) 
     117         x2 =  a12 * (a21 * a33 - a23 * a31) 
     118         x3 =  a13 * (a21 * a32 - a22 * a31) 
     119         det =  x1 - x2 + x3 
     120                  
     121!        CALL determinant(dq(1),dq(2),dq(3),A(1,2),A(2,2),A(3,2),A(1,3),A(2,3),A(3,3),detx) 
     122 
     123         a11=dq(1)  ; a12=dq(2)  ; a13=dq(3) 
     124         a21=A(1,2) ; a22=A(2,2) ; a23=A(3,2) 
     125         a31=A(1,3) ; a32=A(2,3) ; a33=A(3,3) 
     126 
     127         x1 =  a11 * (a22 * a33 - a23 * a32) 
     128         x2 =  a12 * (a21 * a33 - a23 * a31) 
     129         x3 =  a13 * (a21 * a32 - a22 * a31) 
     130         detx =  x1 - x2 + x3 
     131         
     132!        CALL determinant(A(1,1),A(2,1),A(3,1),dq(1),dq(2),dq(3),A(1,3),A(2,3),A(3,3),dety) 
     133 
     134         a11=A(1,1) ; a12=A(2,1) ; a13=A(3,1) 
     135         a21=dq(1)  ; a22=dq(2)  ; a23=dq(3) 
     136         a31=A(1,3) ; a32=A(2,3) ; a33=A(3,3) 
     137 
     138         x1 =  a11 * (a22 * a33 - a23 * a32) 
     139         x2 =  a12 * (a21 * a33 - a23 * a31) 
     140         x3 =  a13 * (a21 * a32 - a22 * a31) 
     141         dety =  x1 - x2 + x3 
     142 
     143!        CALL determinant(A(1,1),A(2,1),A(3,1),A(1,2),A(2,2),A(3,2),dq(1),dq(2),dq(3),detz) 
     144 
     145         a11=A(1,1) ; a12=A(2,1) ; a13=A(3,1) 
     146         a21=A(1,2) ; a22=A(2,2) ; a23=A(3,2) 
     147         a31=dq(1)  ; a32=dq(2)  ; a33=dq(3) 
     148 
     149         x1 =  a11 * (a22 * a33 - a23 * a32) 
     150         x2 =  a12 * (a21 * a33 - a23 * a31) 
     151         x3 =  a13 * (a21 * a32 - a22 * a31) 
     152         detz =  x1 - x2 + x3 
     153 
     154        gradtri(ij+z_up,l,1) = detx 
     155        gradtri(ij+z_up,l,2) = dety 
     156        gradtri(ij+z_up,l,3) = detz 
     157        arr(ij+z_up) = det 
     158         
     159      ENDDO 
     160       
     161      DO ij=ij_begin_ext,ij_end_ext 
     162 
     163 
     164!        CALL gradq(ij,l,ij+t_ldown,ij+t_rdown,ij+z_down,qi,gradtri(ij+z_down,l,1),gradtri(ij+z_down,l,2),gradtri(ij+z_down,l,3),arr(ij+z_down)) 
     165 
     166        A(1,1)=xyz_i(ij+t_ldown,1)-xyz_i(ij,1);  A(1,2)=xyz_i(ij+t_ldown,2)-xyz_i(ij,2); A(1,3)=xyz_i(ij+t_ldown,3)-xyz_i(ij,3)  
     167        A(2,1)=xyz_i(ij+t_rdown,1)-xyz_i(ij,1);  A(2,2)=xyz_i(ij+t_rdown,2)-xyz_i(ij,2); A(2,3)=xyz_i(ij+t_rdown,3)-xyz_i(ij,3)  
     168        A(3,1)=xyz_v(ij+z_down,1);               A(3,2)= xyz_v(ij+z_down,2);             A(3,3)=xyz_v(ij+z_down,3) 
     169  
     170        dq(1) = qi(ij+t_ldown,l)-qi(ij,l) 
     171        dq(2) = qi(ij+t_rdown,l)-qi(ij,l) 
     172        dq(3) = 0.0 
     173 
     174 
     175!        CALL determinant(A(1,1),A(2,1),A(3,1),A(1,2),A(2,2),A(3,2),A(1,3),A(2,3),A(3,3),det) 
     176 
     177         a11=A(1,1) ; a12=A(2,1) ; a13=A(3,1) 
     178         a21=A(1,2) ; a22=A(2,2) ; a23=A(3,2) 
     179         a31=A(1,3) ; a32=A(2,3) ; a33=A(3,3) 
     180 
     181         x1 =  a11 * (a22 * a33 - a23 * a32) 
     182         x2 =  a12 * (a21 * a33 - a23 * a31) 
     183         x3 =  a13 * (a21 * a32 - a22 * a31) 
     184         det =  x1 - x2 + x3 
     185                  
     186!        CALL determinant(dq(1),dq(2),dq(3),A(1,2),A(2,2),A(3,2),A(1,3),A(2,3),A(3,3),detx) 
     187 
     188         a11=dq(1)  ; a12=dq(2)  ; a13=dq(3) 
     189         a21=A(1,2) ; a22=A(2,2) ; a23=A(3,2) 
     190         a31=A(1,3) ; a32=A(2,3) ; a33=A(3,3) 
     191 
     192         x1 =  a11 * (a22 * a33 - a23 * a32) 
     193         x2 =  a12 * (a21 * a33 - a23 * a31) 
     194         x3 =  a13 * (a21 * a32 - a22 * a31) 
     195         detx =  x1 - x2 + x3 
     196         
     197!        CALL determinant(A(1,1),A(2,1),A(3,1),dq(1),dq(2),dq(3),A(1,3),A(2,3),A(3,3),dety) 
     198 
     199         a11=A(1,1) ; a12=A(2,1) ; a13=A(3,1) 
     200         a21=dq(1)  ; a22=dq(2)  ; a23=dq(3) 
     201         a31=A(1,3) ; a32=A(2,3) ; a33=A(3,3) 
     202 
     203         x1 =  a11 * (a22 * a33 - a23 * a32) 
     204         x2 =  a12 * (a21 * a33 - a23 * a31) 
     205         x3 =  a13 * (a21 * a32 - a22 * a31) 
     206         dety =  x1 - x2 + x3 
     207 
     208!        CALL determinant(A(1,1),A(2,1),A(3,1),A(1,2),A(2,2),A(3,2),dq(1),dq(2),dq(3),detz) 
     209 
     210         a11=A(1,1) ; a12=A(2,1) ; a13=A(3,1) 
     211         a21=A(1,2) ; a22=A(2,2) ; a23=A(3,2) 
     212         a31=dq(1)  ; a32=dq(2)  ; a33=dq(3) 
     213 
     214         x1 =  a11 * (a22 * a33 - a23 * a32) 
     215         x2 =  a12 * (a21 * a33 - a23 * a31) 
     216         x3 =  a13 * (a21 * a32 - a22 * a31) 
     217         detz =  x1 - x2 + x3 
     218 
     219         gradtri(ij+z_down,l,1) = detx 
     220         gradtri(ij+z_down,l,2) = dety 
     221         gradtri(ij+z_down,l,3) = detz 
     222         arr(ij+z_down) = det 
     223 
     224      END DO 
    79225    END DO 
    80226 
    81227!$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 
     228    DO ij=ij_begin,ij_end 
     229       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 
    84230    ENDDO 
     231    CALL trace_end("compute_gradq3d1") 
     232    CALL trace_start2("compute_gradq3d2") 
    85233       
    86234    DO k=1,3 
    … …  
    94242      END DO 
    95243    ENDDO 
     244    CALL trace_end2("compute_gradq3d2") 
     245    CALL trace_start("compute_gradq3d3") 
    96246 
    97247    !============================================================================================= LIMITING  
    … …  
    99249!$SIMD 
    100250      DO ij=ij_begin,ij_end 
    101              maggrd =  dot_product(gradq3d(ij,l,:),gradq3d(ij,l,:)) 
     251!             maggrd =  dot_product(gradq3d(ij,l,:),gradq3d(ij,l,:)) 
     252             maggrd = gradq3d(ij,l,1)*gradq3d(ij,l,1) + gradq3d(ij,l,2)*gradq3d(ij,l,2) + gradq3d(ij,l,3)*gradq3d(ij,l,3)  
    102253             maggrd = sqrt(maggrd)  
    103254             maxq_c = qi(ij,l) + maggrd*one_over_sqrt_leng(ij) 
    … …  
    112263             alphami = max(alphami,0.0)  
    113264             alpha   = min(alphamx,alphami,1.0) 
    114              gradq3d(ij,l,:) = alpha*gradq3d(ij,l,:) 
     265!             gradq3d(ij,l,:) = alpha*gradq3d(ij,l,:) 
     266             gradq3d(ij,l,1) = alpha*gradq3d(ij,l,1) 
     267             gradq3d(ij,l,2) = alpha*gradq3d(ij,l,2) 
     268             gradq3d(ij,l,3) = alpha*gradq3d(ij,l,3) 
    115269       END DO 
    116270    END DO 
    117271 
    118   CALL trace_end("compute_gradq3d") 
     272  CALL trace_end("compute_gradq3d3") 
     273   
     274  gradq3d_out=gradq3d 
     275   
     276  CONTAINS 
     277 
     278    SUBROUTINE gradq(n0,l,n1,n2,n3,q,dq1,dq2,dq3,det) 
     279    IMPLICIT NONE 
     280    INTEGER, INTENT(IN) :: n0,l,n1,n2,n3 
     281    REAL(rstd), INTENT(IN)     :: q(iim*jjm,llm) 
     282!    REAL(rstd), INTENT(OUT)    :: dq(3), det 
     283    REAL(rstd), INTENT(OUT)    :: dq1,dq2,dq3,det 
     284    REAL(rstd)    :: dq(3) 
     285     
     286    REAL(rstd)  ::detx,dety,detz 
     287    INTEGER    :: info 
     288    INTEGER    :: IPIV(3) 
     289 
     290    REAL(rstd) :: A(3,3) 
     291    REAL(rstd) :: B(3) 
     292     
     293    ! TODO : replace A by A1,A2,A3 
     294    A(1,1)=xyz_i(n1,1)-xyz_i(n0,1);  A(1,2)=xyz_i(n1,2)-xyz_i(n0,2); A(1,3)=xyz_i(n1,3)-xyz_i(n0,3)  
     295    A(2,1)=xyz_i(n2,1)-xyz_i(n0,1);  A(2,2)=xyz_i(n2,2)-xyz_i(n0,2); A(2,3)=xyz_i(n2,3)-xyz_i(n0,3)  
     296    A(3,1)=xyz_v(n3,1);              A(3,2)= xyz_v(n3,2);            A(3,3)=xyz_v(n3,3) 
     297 
     298    dq(1) = q(n1,l)-q(n0,l) 
     299    dq(2) = q(n2,l)-q(n0,l) 
     300    dq(3) = 0.0 
     301 
     302    !    CALL DGESV(3,1,A,3,IPIV,dq(:),3,info) 
     303!    CALL determinant(A(:,1),A(:,2),A(:,3),det) 
     304!    CALL determinant(dq,A(:,2),A(:,3),detx) 
     305!    CALL determinant(A(:,1),dq,A(:,3),dety) 
     306!    CALL determinant(A(:,1),A(:,2),dq,detz) 
     307!    dq(1) = detx 
     308!    dq(2) = dety 
     309!    dq(3) = detz  
     310 
     311    CALL determinant(A(1,1),A(2,1),A(3,1),A(1,2),A(2,2),A(3,2),A(1,3),A(2,3),A(3,3),det) 
     312    CALL determinant(dq(1),dq(2),dq(3),A(1,2),A(2,2),A(3,2),A(1,3),A(2,3),A(3,3),dq1) 
     313    CALL determinant(A(1,1),A(2,1),A(3,1),dq(1),dq(2),dq(3),A(1,3),A(2,3),A(3,3),dq2) 
     314    CALL determinant(A(1,1),A(2,1),A(3,1),A(1,2),A(2,2),A(3,2),dq(1),dq(2),dq(3),dq3) 
     315 
     316  END SUBROUTINE gradq 
     317 
     318  !========================================================================== 
     319!  PURE SUBROUTINE determinant(a1,a2,a3,det) 
     320!    IMPLICIT NONE 
     321!    REAL(rstd), DIMENSION(3), INTENT(IN) :: a1,a2,a3 
     322!    REAL(rstd), INTENT(OUT) :: det 
     323!    REAL(rstd) ::  x1,x2,x3 
     324!    x1 =  a1(1) * (a2(2) * a3(3) - a2(3) * a3(2)) 
     325!    x2 =  a1(2) * (a2(1) * a3(3) - a2(3) * a3(1)) 
     326!    x3 =  a1(3) * (a2(1) * a3(2) - a2(2) * a3(1)) 
     327!    det =  x1 - x2 + x3 
     328!  END SUBROUTINE determinant 
     329 
     330   SUBROUTINE determinant(a11,a12,a13,a21,a22,a23,a31,a32,a33,det) 
     331    IMPLICIT NONE 
     332    REAL(rstd), INTENT(IN) :: a11,a12,a13,a21,a22,a23,a31,a32,a33 
     333    REAL(rstd), INTENT(OUT) :: det 
     334    REAL(rstd) ::  x1,x2,x3 
     335    x1 =  a11 * (a22 * a33 - a23 * a32) 
     336    x2 =  a12 * (a21 * a33 - a23 * a31) 
     337    x3 =  a13 * (a21 * a32 - a22 * a31) 
     338    det =  x1 - x2 + x3 
     339  END SUBROUTINE determinant 
     340 
     341 
    119342  END SUBROUTINE compute_gradq3d 
    120343 
    … …  
    222445             IF (ne(ij,right)*hfluxt(ij+u_right,l)>0) THEN  
    223446                ed = cc(ij+u_right,l,:) - xyz_i(ij,:) 
    224                 qe = qi(ij,l)+sum2(gradq3d(ij,l,:),ed)  
     447!                qe = qi(ij,l)+sum2(gradq3d(ij,l,:),ed)  
     448                 qe = qi(ij,l)+gradq3d(ij,l,1)*ed(1)+gradq3d(ij,l,2)*ed(2)+gradq3d(ij,l,3)*ed(3)  
    225449             ELSE 
    226450                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)  
     451!                qe = qi(ij+t_right,l)+sum2(gradq3d(ij+t_right,l,:),ed)  
     452                qe = qi(ij+t_right,l) + gradq3d(ij+t_right,l,1)*ed(1)+gradq3d(ij+t_right,l,2)*ed(2)+gradq3d(ij+t_right,l,3)*ed(3)  
    228453             ENDIF 
    229454             qflux(ij+u_right,l) = hfluxt(ij+u_right,l)*qe 
    … …  
    231456             IF (ne(ij,lup)*hfluxt(ij+u_lup,l)>0) THEN  
    232457                ed = cc(ij+u_lup,l,:) - xyz_i(ij,:) 
    233                 qe = qi(ij,l)+sum2(gradq3d(ij,l,:),ed) 
     458!                qe = qi(ij,l)+sum2(gradq3d(ij,l,:),ed) 
     459                qe = qi(ij,l)  + gradq3d(ij,l,1)*ed(1)+gradq3d(ij,l,2)*ed(2)+gradq3d(ij,l,3)*ed(3)  
    234460             ELSE 
    235461                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)  
     462!                qe = qi(ij+t_lup,l)+sum2(gradq3d(ij+t_lup,l,:),ed)  
     463                qe = qi(ij+t_lup,l) + gradq3d(ij+t_lup,l,1)*ed(1)+gradq3d(ij+t_lup,l,2)*ed(2)+gradq3d(ij+t_lup,l,3)*ed(3)  
    237464             ENDIF 
    238465             qflux(ij+u_lup,l) = hfluxt(ij+u_lup,l)*qe  
    … …  
    240467             IF (ne(ij,ldown)*hfluxt(ij+u_ldown,l)>0) THEN  
    241468                ed = cc(ij+u_ldown,l,:) - xyz_i(ij,:) 
    242                 qe = qi(ij,l)+sum2(gradq3d(ij,l,:),ed)  
     469!                qe = qi(ij,l)+sum2(gradq3d(ij,l,:),ed)  
     470                qe = qi(ij,l) + gradq3d(ij,l,1)*ed(1)+gradq3d(ij,l,2)*ed(2)+gradq3d(ij,l,3)*ed(3)  
    243471             ELSE 
    244472                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)  
     473!                qe = qi(ij+t_ldown,l)+sum2(gradq3d(ij+t_ldown,l,:),ed)  
     474                qe = qi(ij+t_ldown,l)+gradq3d(ij+t_ldown,l,1)*ed(1)+gradq3d(ij+t_ldown,l,2)*ed(2)+gradq3d(ij+t_ldown,l,3)*ed(3)  
    246475             ENDIF 
    247476             qflux(ij+u_ldown,l) = hfluxt(ij+u_ldown,l)*qe 
    … …  
    289518  END SUBROUTINE compute_advect_horiz 
    290519 
    291   !========================================================================== 
    292   PURE SUBROUTINE gradq(n0,l,n1,n2,n3,q,dq,det) 
    293     IMPLICIT NONE 
    294     INTEGER, INTENT(IN) :: n0,l,n1,n2,n3 
    295     REAL(rstd), INTENT(IN)     :: q(iim*jjm,llm) 
    296     REAL(rstd), INTENT(OUT)    :: dq(3), det 
    297      
    298     REAL(rstd)  ::detx,dety,detz 
    299     INTEGER    :: info 
    300     INTEGER    :: IPIV(3) 
    301  
    302     REAL(rstd) :: A(3,3) 
    303     REAL(rstd) :: B(3) 
    304  
    305     ! TODO : replace A by A1,A2,A3 
    306     A(1,1)=xyz_i(n1,1)-xyz_i(n0,1);  A(1,2)=xyz_i(n1,2)-xyz_i(n0,2); A(1,3)=xyz_i(n1,3)-xyz_i(n0,3)  
    307     A(2,1)=xyz_i(n2,1)-xyz_i(n0,1);  A(2,2)=xyz_i(n2,2)-xyz_i(n0,2); A(2,3)=xyz_i(n2,3)-xyz_i(n0,3)  
    308     A(3,1)=xyz_v(n3,1);              A(3,2)= xyz_v(n3,2);            A(3,3)=xyz_v(n3,3) 
    309  
    310     dq(1) = q(n1,l)-q(n0,l) 
    311     dq(2) = q(n2,l)-q(n0,l) 
    312     dq(3) = 0.0 
    313     !    CALL DGESV(3,1,A,3,IPIV,dq(:),3,info) 
    314     CALL determinant(A(:,1),A(:,2),A(:,3),det) 
    315     CALL determinant(dq,A(:,2),A(:,3),detx) 
    316     CALL determinant(A(:,1),dq,A(:,3),dety) 
    317     CALL determinant(A(:,1),A(:,2),dq,detz) 
    318     dq(1) = detx 
    319     dq(2) = dety 
    320     dq(3) = detz  
    321   END SUBROUTINE gradq 
    322  
    323   !========================================================================== 
    324   PURE SUBROUTINE determinant(a1,a2,a3,det) 
    325     IMPLICIT NONE 
    326     REAL(rstd), DIMENSION(3), INTENT(IN) :: a1,a2,a3 
    327     REAL(rstd), INTENT(OUT) :: det 
    328     REAL(rstd) ::  x1,x2,x3 
    329     x1 =  a1(1) * (a2(2) * a3(3) - a2(3) * a3(2)) 
    330     x2 =  a1(2) * (a2(1) * a3(3) - a2(3) * a3(1)) 
    331     x3 =  a1(3) * (a2(1) * a3(2) - a2(2) * a3(1)) 
    332     det =  x1 - x2 + x3 
    333   END SUBROUTINE determinant 
    334520 
    335521END MODULE advect_mod 
Note: See TracChangeset for help on using the changeset viewer.