Ignore:
Timestamp:
07/18/12 18:39:59 (12 years ago)
Author:
dubos
Message:

Review of 3D transport code
Fewer interations in dissip_gcm
This revision is certainly broken ; see next revision

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r19 r22  
    1       MODULE advect_mod 
    2    USE icosa 
    3         IMPLICIT NONE  
    4  
    5  
    6   TYPE(t_field),POINTER :: f_normal(:) 
    7   TYPE(t_field),POINTER :: f_tangent(:) 
    8   TYPE(t_field),POINTER ::  f_gradq3d(:) 
    9   REAL(rstd),POINTER :: gradq3d(:,:,:)  
    10   REAL(rstd),POINTER :: normal(:,:) 
    11   REAL(rstd),POINTER :: tangent(:,:) 
    12                          
     1MODULE advect_mod 
     2  USE icosa 
     3  IMPLICIT NONE  
     4 
    135CONTAINS 
    14           
    15          
    16   SUBROUTINE allocate_advect 
    17   USE field_mod 
    18   USE domain_mod 
    19   USE dimensions 
    20   USE geometry 
    21   USE metric 
    22   IMPLICIT NONE 
    23  
    24     CALL allocate_field(f_normal,field_u,type_real,3) 
    25     CALL allocate_field(f_tangent,field_u,type_real,3) 
    26     CALL allocate_field(f_gradq3d,field_t,type_real,llm,3) 
    27  
    28   END SUBROUTINE allocate_advect 
    29 !========================================================================== 
    30   SUBROUTINE swap_advect(ind) 
    31   USE domain_mod 
    32   USE dimensions 
    33   USE geometry 
    34   USE metric 
    35   IMPLICIT NONE 
    36   INTEGER,INTENT(IN) :: ind 
    37    
    38       normal=f_normal(ind) 
    39       tangent=f_tangent(ind) 
    40       gradq3d = f_gradq3d(ind)  
    41  
    42   END SUBROUTINE swap_advect 
    43 !========================================================================== 
    44      
    45   SUBROUTINE init_advect 
    46   USE domain_mod 
    47   USE dimensions 
    48   USE geometry 
    49   USE metric 
    50   USE vector 
    51   IMPLICIT NONE 
    52    INTEGER :: ind,i,j,n 
    53     
    54    CALL allocate_advect 
    55     
     6 
     7  !========================================================================== 
     8 
     9  SUBROUTINE init_advect(normal,tangent) 
     10    USE domain_mod 
     11    USE dimensions 
     12    USE geometry 
     13    USE metric 
     14    USE vector 
     15    IMPLICIT NONE 
     16    REAL(rstd),INTENT(OUT) :: normal(3*iim*jjm,3) 
     17    REAL(rstd),INTENT(OUT) :: tangent(3*iim*jjm,3) 
     18 
     19    INTEGER :: i,j,n 
     20 
    5621    DO j=jj_begin-1,jj_end+1 
    57       DO i=ii_begin-1,ii_end+1 
    58         n=(j-1)*iim+i 
    59              
    60         CALL cross_product2(xyz_v(n+z_rdown,:),xyz_v(n+z_rup,:),normal(n+u_right,:)) 
    61         normal(n+u_right,:)=normal(n+u_right,:)/sqrt(sum(normal(n+u_right,:)**2)+1e-50)*ne(n,right) 
    62           
    63         CALL cross_product2(xyz_v(n+z_up,:),xyz_v(n+z_lup,:),normal(n+u_lup,:)) 
    64          normal(n+u_lup,:)=normal(n+u_lup,:)/sqrt(sum(normal(n+u_lup,:)**2)+1e-50)*ne(n,lup) 
    65          
    66         CALL cross_product2(xyz_v(n+z_ldown,:),xyz_v(n+z_down,:),normal(n+u_ldown,:))     
    67         normal(n+u_ldown,:)=normal(n+u_ldown,:)/sqrt(sum(normal(n+u_ldown,:)**2)+1e-50)*ne(n,ldown) 
    68                        
    69         tangent(n+u_right,:)=xyz_v(n+z_rup,:)-xyz_v(n+z_rdown,:)  
    70         tangent(n+u_right,:)=tangent(n+u_right,:)/sqrt(sum(tangent(n+u_right,:)**2)+1e-50) 
    71               
    72         tangent(n+u_lup,:)=xyz_v(n+z_lup,:)-xyz_v(n+z_up,:)  
    73         tangent(n+u_lup,:)=tangent(n+u_lup,:)/sqrt(sum(tangent(n+u_lup,:)**2)+1e-50) 
    74       
    75         tangent(n+u_ldown,:)=xyz_v(n+z_down,:)-xyz_v(n+z_ldown,:) 
    76         tangent(n+u_ldown,:)=tangent(n+u_ldown,:)/sqrt(sum(tangent(n+u_ldown,:)**2)+1e-50) 
    77       END DO 
    78      ENDDO       
     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       END DO 
     43    ENDDO 
     44 
    7945  END SUBROUTINE init_advect 
    80 !======================================================================================= 
    81  
    82  
    83  
    84 !======================================================================================= 
    85   SUBROUTINE advect1(qi) 
    86   USE domain_mod 
    87   USE dimensions 
    88   USE geometry 
    89   USE metric 
    90   IMPLICIT NONE 
    91     REAL(rstd),INTENT(IN) :: qi(iim*jjm,llm) 
     46 
     47  !======================================================================================= 
     48 
     49  SUBROUTINE compute_gradq3d(qi,gradq3d) 
     50    USE domain_mod 
     51    USE dimensions 
     52    USE geometry 
     53    USE metric 
     54    IMPLICIT NONE 
     55    REAL(rstd),INTENT(IN)  :: qi(iim*jjm,llm) 
     56    REAL(rstd),INTENT(OUT) :: gradq3d(iim*jm,llm,3)  
    9257    REAL(rstd) :: maxq,minq,minq_c,maxq_c  
    9358    REAL(rstd) :: alphamx,alphami,alpha ,maggrd,leng 
    … …  
    9762    REAL(rstd) :: ar 
    9863    INTEGER :: i,j,n,k,ind,l 
    99 !========================================================================================== GRADIENT  
    100        Do l = 1,llm  
    101         DO j=jj_begin-1,jj_end+1 
    102          DO i=ii_begin-1,ii_end+1 
    103           n=(j-1)*iim+i    
    104           CALL gradq(n,n+t_rup,n+t_lup,n+z_up,qi(:,l),gradtri(n+z_up,l,:),arr(n+z_up)) 
    105           CALL gradq(n,n+t_ldown,n+t_rdown,n+z_down,qi(:,l),gradtri(n+z_down,l,:),arr(n+z_down)) 
    106          END DO 
    107         END DO 
    108        END DO 
    109  
    110 !      Do l =1,llm 
    111         DO j=jj_begin,jj_end 
     64    !========================================================================================== GRADIENT  
     65    Do l = 1,llm  
     66       DO j=jj_begin-1,jj_end+1 
     67          DO i=ii_begin-1,ii_end+1 
     68             n=(j-1)*iim+i    
     69             CALL gradq(n,n+t_rup,n+t_lup,n+z_up,qi(:,l),gradtri(n+z_up,l,:),arr(n+z_up)) 
     70             CALL gradq(n,n+t_ldown,n+t_rdown,n+z_down,qi(:,l),gradtri(n+z_down,l,:),arr(n+z_down)) 
     71          END DO 
     72       END DO 
     73    END DO 
     74 
     75    !      Do l =1,llm 
     76    DO j=jj_begin,jj_end 
     77       DO i=ii_begin,ii_end 
     78          n=(j-1)*iim+i 
     79          gradq3d(n,:,:) = gradtri(n+z_up,:,:) + gradtri(n+z_down,:,:)   +   & 
     80               gradtri(n+z_rup,:,:) +  gradtri(n+z_ldown,:,:)   +   &  
     81               gradtri(n+z_lup,:,:)+    gradtri(n+z_rdown,:,:)   
     82          ar = arr(n+z_up)+arr(n+z_lup)+arr(n+z_ldown)+arr(n+z_down)+arr(n+z_rdown)+arr(n+z_rup) 
     83          gradq3d(n,:,:) = gradq3d(n,:,:)/(ar+1.e-50)  
     84       END DO 
     85    END DO 
     86    !    END DO 
     87 
     88    !============================================================================================= LIMITING  
     89    ! GO TO 120  
     90    DO l =1,llm 
     91       DO j=jj_begin,jj_end 
    11292          DO i=ii_begin,ii_end 
    11393             n=(j-1)*iim+i 
    114              gradq3d(n,:,:) = gradtri(n+z_up,:,:) + gradtri(n+z_down,:,:)   +   & 
    115                            gradtri(n+z_rup,:,:) +  gradtri(n+z_ldown,:,:)   +   &  
    116                           gradtri(n+z_lup,:,:)+    gradtri(n+z_rdown,:,:)   
    117              ar = arr(n+z_up)+arr(n+z_lup)+arr(n+z_ldown)+arr(n+z_down)+arr(n+z_rdown)+arr(n+z_rup) 
    118              gradq3d(n,:,:) = gradq3d(n,:,:)/(ar+1.e-50)  
    119            END DO 
    120          END DO     
    121 !        END DO 
    122   
    123 !============================================================================================= LIMITING  
    124 ! GO TO 120  
    125   Do l =1,llm 
    126    DO j=jj_begin,jj_end 
    127      DO i=ii_begin,ii_end 
    128        n=(j-1)*iim+i 
    129        maggrd =  dot_product(gradq3d(n,l,:),gradq3d(n,l,:)) 
    130        maggrd = sqrt(maggrd)  
    131        leng = max(sum((xyz_v(n+z_up,:) - xyz_i(n,:))**2),sum((xyz_v(n+z_down,:) - xyz_i(n,:))**2), & 
    132        sum((xyz_v(n+z_rup,:) - xyz_i(n,:))**2),sum((xyz_v(n+z_rdown,:) - xyz_i(n,:))**2),  & 
    133        sum((xyz_v(n+z_lup,:) - xyz_i(n,:))**2),sum((xyz_v(n+z_ldown,:) - xyz_i(n,:))**2)) 
    134        maxq_c = qi(n,l) + maggrd*sqrt(leng) 
    135        minq_c = qi(n,l) - maggrd*sqrt(leng) 
    136        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), & 
    137              qi(n+t_rdown,l),qi(n+t_ldown,l)) 
    138        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), & 
    139              qi(n+t_rdown,l),qi(n+t_ldown,l)) 
    140        alphamx = (maxq - qi(n,l)) ; alphamx = alphamx/(maxq_c - qi(n,l) ) 
    141        alphamx = max(alphamx,0.0) 
    142        alphami = (minq - qi(n,l)); alphami = alphami/(minq_c - qi(n,l)) 
    143        alphami = max(alphami,0.0)  
    144        alpha   = min(alphamx,alphami,1.0) 
    145        gradq3d(n,l,:) = alpha*gradq3d(n,l,:) 
    146       END DO 
    147     END DO 
    148    END DO 
    149          END SUBROUTINE ADVECT1           
    150 !===================================================================================================      
    151          SUBROUTINE advect2(qi,him,ue,he,bigt) 
    152   USE domain_mod 
    153   USE dimensions 
    154   USE geometry 
    155   USE metric 
    156   IMPLICIT NONE 
     94             maggrd =  dot_product(gradq3d(n,l,:),gradq3d(n,l,:)) 
     95             maggrd = sqrt(maggrd)  
     96             leng = max(sum((xyz_v(n+z_up,:) - xyz_i(n,:))**2),sum((xyz_v(n+z_down,:) - xyz_i(n,:))**2), & 
     97                  sum((xyz_v(n+z_rup,:) - xyz_i(n,:))**2),sum((xyz_v(n+z_rdown,:) - xyz_i(n,:))**2),  & 
     98                  sum((xyz_v(n+z_lup,:) - xyz_i(n,:))**2),sum((xyz_v(n+z_ldown,:) - xyz_i(n,:))**2)) 
     99             maxq_c = qi(n,l) + maggrd*sqrt(leng) 
     100             minq_c = qi(n,l) - maggrd*sqrt(leng) 
     101             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), & 
     102                  qi(n+t_rdown,l),qi(n+t_ldown,l)) 
     103             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), & 
     104                  qi(n+t_rdown,l),qi(n+t_ldown,l)) 
     105             alphamx = (maxq - qi(n,l)) ; alphamx = alphamx/(maxq_c - qi(n,l) ) 
     106             alphamx = max(alphamx,0.0) 
     107             alphami = (minq - qi(n,l)); alphami = alphami/(minq_c - qi(n,l)) 
     108             alphami = max(alphami,0.0)  
     109             alpha   = min(alphamx,alphami,1.0) 
     110             gradq3d(n,l,:) = alpha*gradq3d(n,l,:) 
     111          END DO 
     112       END DO 
     113    END DO 
     114  END SUBROUTINE compute_gradq3d 
     115 
     116  !===================================================================================================    
     117  SUBROUTINE compute_advect_horiz(normal,tangent,qi,him,ue,he,bigt) 
     118    USE domain_mod 
     119    USE dimensions 
     120    USE geometry 
     121    USE metric 
     122    IMPLICIT NONE 
     123    REAL(rstd),INTENT(IN)    :: normal(3*iim*jjm,3) 
     124    REAL(rstd),INTENT(IN)    :: tangent(3*iim*jjm,3) 
    157125    REAL(rstd),INTENT(INOUT) :: qi(iim*jjm,llm) 
     126    REAL(rstd),INTENT(IN)    :: gradq3d(iim*jm,llm,3)  
    158127    REAL(rstd),INTENT(INOUT) :: him(iim*jjm,llm) 
    159     REAL(rstd),INTENT(IN) :: ue(iim*3*jjm,llm) 
    160     REAL(rstd),INTENT(IN) :: he(3*iim*jjm,llm) 
    161     REAL(rstd),INTENT(IN)::bigt 
     128    REAL(rstd),INTENT(IN)    :: ue(iim*3*jjm,llm) 
     129    REAL(rstd),INTENT(IN)    :: he(3*iim*jjm,llm) ! mass flux 
     130    REAL(rstd),INTENT(IN)    :: bigt 
     131 
    162132    REAL(rstd) :: dqi(iim*jjm,llm),dhi(iim*jjm,llm)  
    163133    REAL(rstd) :: cc(3*iim*jjm,llm,3)  
    164134    REAL(rstd) :: v_e(3*iim*jjm,llm,3) 
    165135    REAL(rstd) :: up_e 
    166      
     136 
    167137    REAL(rstd) :: qe(3*iim*jjm,llm) 
    168138    REAL(rstd) :: ed(3)     
    … …  
    172142 
    173143 
    174 !go to  123 
    175    DO l = 1,llm 
     144    !go to  123 
     145    DO l = 1,llm 
     146       DO j=jj_begin,jj_end 
     147          DO i=ii_begin,ii_end 
     148             n=(j-1)*iim+i 
     149 
     150             up_e =1/de(n+u_right)*(                                                   & 
     151                  wee(n+u_right,1,1)*le(n+u_rup)*ue(n+u_rup,l)+                        & 
     152                  wee(n+u_right,2,1)*le(n+u_lup)*ue(n+u_lup,l)+                        & 
     153                  wee(n+u_right,3,1)*le(n+u_left)*ue(n+u_left,l)+                      & 
     154                  wee(n+u_right,4,1)*le(n+u_ldown)*ue(n+u_ldown,l)+                    & 
     155                  wee(n+u_right,5,1)*le(n+u_rdown)*ue(n+u_rdown,l)+                    &  
     156                  wee(n+u_right,1,2)*le(n+t_right+u_ldown)*ue(n+t_right+u_ldown,l)+    & 
     157                  wee(n+u_right,2,2)*le(n+t_right+u_rdown)*ue(n+t_right+u_rdown,l)+    & 
     158                  wee(n+u_right,3,2)*le(n+t_right+u_right)*ue(n+t_right+u_right,l)+    & 
     159                  wee(n+u_right,4,2)*le(n+t_right+u_rup)*ue(n+t_right+u_rup,l)+        & 
     160                  wee(n+u_right,5,2)*le(n+t_right+u_lup)*ue(n+t_right+u_lup,l)         & 
     161                  ) 
     162             v_e(n+u_right,l,:)= ue(n+u_right,l)*normal(n+u_right,:) + up_e*tangent(n+u_right,:) 
     163 
     164             up_e=1/de(n+u_lup)*(                                                        & 
     165                  wee(n+u_lup,1,1)*le(n+u_left)*ue(n+u_left,l)+                        & 
     166                  wee(n+u_lup,2,1)*le(n+u_ldown)*ue(n+u_ldown,l)+                      & 
     167                  wee(n+u_lup,3,1)*le(n+u_rdown)*ue(n+u_rdown,l)+                      & 
     168                  wee(n+u_lup,4,1)*le(n+u_right)*ue(n+u_right,l)+                      & 
     169                  wee(n+u_lup,5,1)*le(n+u_rup)*ue(n+u_rup,l)+                          &  
     170                  wee(n+u_lup,1,2)*le(n+t_lup+u_right)*ue(n+t_lup+u_right,l)+          & 
     171                  wee(n+u_lup,2,2)*le(n+t_lup+u_rup)*ue(n+t_lup+u_rup,l)+              & 
     172                  wee(n+u_lup,3,2)*le(n+t_lup+u_lup)*ue(n+t_lup+u_lup,l)+              & 
     173                  wee(n+u_lup,4,2)*le(n+t_lup+u_left)*ue(n+t_lup+u_left,l)+            & 
     174                  wee(n+u_lup,5,2)*le(n+t_lup+u_ldown)*ue(n+t_lup+u_ldown,l)           & 
     175                  ) 
     176             v_e(n+u_lup,l,:)= ue(n+u_lup,l)*normal(n+u_lup,:) + up_e*tangent(n+u_lup,:) 
     177 
     178             up_e=1/de(n+u_ldown)*(                                                    & 
     179                  wee(n+u_ldown,1,1)*le(n+u_rdown)*ue(n+u_rdown,l)+                    & 
     180                  wee(n+u_ldown,2,1)*le(n+u_right)*ue(n+u_right,l)+                    & 
     181                  wee(n+u_ldown,3,1)*le(n+u_rup)*ue(n+u_rup,l)+                        & 
     182                  wee(n+u_ldown,4,1)*le(n+u_lup)*ue(n+u_lup,l)+                        & 
     183                  wee(n+u_ldown,5,1)*le(n+u_left)*ue(n+u_left,l)+                      &  
     184                  wee(n+u_ldown,1,2)*le(n+t_ldown+u_lup)*ue(n+t_ldown+u_lup,l)+        & 
     185                  wee(n+u_ldown,2,2)*le(n+t_ldown+u_left)*ue(n+t_ldown+u_left,l)+      & 
     186                  wee(n+u_ldown,3,2)*le(n+t_ldown+u_ldown)*ue(n+t_ldown+u_ldown,l)+    & 
     187                  wee(n+u_ldown,4,2)*le(n+t_ldown+u_rdown)*ue(n+t_ldown+u_rdown,l)+    & 
     188                  wee(n+u_ldown,5,2)*le(n+t_ldown+u_right)*ue(n+t_ldown+u_right,l)     & 
     189                  ) 
     190             v_e(n+u_ldown,l,:)= ue(n+u_ldown,l)*normal(n+u_ldown,:) + up_e*tangent(n+u_ldown,:) 
     191 
     192             cc(n+u_right,l,:) = xyz_e(n+u_right,:) - v_e(n+u_right,l,:)*0.5*bigt    !! redge is mid point of edge i  
     193             cc(n+u_lup,l,:)   = xyz_e(n+u_lup,:)   - v_e(n+u_lup,l,:)*0.5*bigt  
     194             cc(n+u_ldown,l,:) = xyz_e(n+u_ldown,:) - v_e(n+u_ldown,l,:)*0.5*bigt  
     195          ENDDO 
     196       ENDDO 
     197    END DO 
     198    !123         continue            
     199    !========================================================================================================== 
     200    DO l = 1,llm 
     201       DO j=jj_begin-1,jj_end+1 
     202          DO i=ii_begin-1,ii_end+1 
     203             n=(j-1)*iim+i 
     204             if (ne(n,right)*ue(n+u_right,l)>0) then  
     205                ed = cc(n+u_right,l,:) - xyz_i(n,:) 
     206                qe(n+u_right,l)=qi(n,l)+sum2(gradq3d(n,l,:),ed)  
     207             else 
     208                ed = cc(n+u_right,l,:) - xyz_i(n+t_right,:) 
     209                qe(n+u_right,l)=qi(n+t_right,l)+sum2(gradq3d(n+t_right,l,:),ed)  
     210             endif 
     211             if (ne(n,lup)*ue(n+u_lup,l)>0) then  
     212                ed = cc(n+u_lup,l,:) - xyz_i(n,:) 
     213                qe(n+u_lup,l)=qi(n,l)+sum2(gradq3d(n,l,:),ed) 
     214             else 
     215                ed = cc(n+u_lup,l,:) - xyz_i(n+t_lup,:) 
     216                qe(n+u_lup,l)=qi(n+t_lup,l)+sum2(gradq3d(n+t_lup,l,:),ed)  
     217             endif 
     218             if (ne(n,ldown)*ue(n+u_ldown,l)>0) then  
     219                ed = cc(n+u_ldown,l,:) - xyz_i(n,:) 
     220                qe(n+u_ldown,l)=qi(n,l)+ sum2(gradq3d(n,l,:),ed)  
     221             else 
     222                ed = cc(n+u_ldown,l,:) - xyz_i(n+t_ldown,:) 
     223                qe(n+u_ldown,l)=qi(n+t_ldown,l)+sum2(gradq3d(n+t_ldown,l,:),ed)  
     224             endif 
     225          end do 
     226       end do 
     227    END DO 
     228 
     229 
     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          flxx(n+u_right,:) = he(n+u_right,:)*qe(n+u_right,:)*ne(n,right) 
     234          flxx(n+u_lup,:) = he(n+u_lup,:)*qe(n+u_lup,:)*ne(n,lup)  
     235          flxx(n+u_ldown,:) = he(n+u_ldown,:)*qe(n+u_ldown,:)*ne(n,ldown) 
     236       ENDDO 
     237    ENDDO 
     238 
    176239    DO j=jj_begin,jj_end 
    177240       DO i=ii_begin,ii_end 
    178         n=(j-1)*iim+i 
    179  
    180                     up_e =1/de(n+u_right)*(                                                   & 
    181                          wee(n+u_right,1,1)*le(n+u_rup)*ue(n+u_rup,l)+                        & 
    182                          wee(n+u_right,2,1)*le(n+u_lup)*ue(n+u_lup,l)+                        & 
    183                          wee(n+u_right,3,1)*le(n+u_left)*ue(n+u_left,l)+                      & 
    184                          wee(n+u_right,4,1)*le(n+u_ldown)*ue(n+u_ldown,l)+                    & 
    185                          wee(n+u_right,5,1)*le(n+u_rdown)*ue(n+u_rdown,l)+                    &  
    186                          wee(n+u_right,1,2)*le(n+t_right+u_ldown)*ue(n+t_right+u_ldown,l)+    & 
    187                          wee(n+u_right,2,2)*le(n+t_right+u_rdown)*ue(n+t_right+u_rdown,l)+    & 
    188                          wee(n+u_right,3,2)*le(n+t_right+u_right)*ue(n+t_right+u_right,l)+    & 
    189                          wee(n+u_right,4,2)*le(n+t_right+u_rup)*ue(n+t_right+u_rup,l)+        & 
    190                          wee(n+u_right,5,2)*le(n+t_right+u_lup)*ue(n+t_right+u_lup,l)         & 
    191                                         ) 
    192        v_e(n+u_right,l,:)= ue(n+u_right,l)*normal(n+u_right,:) + up_e*tangent(n+u_right,:) 
    193    
    194                 up_e=1/de(n+u_lup)*(                                                        & 
    195                          wee(n+u_lup,1,1)*le(n+u_left)*ue(n+u_left,l)+                        & 
    196                          wee(n+u_lup,2,1)*le(n+u_ldown)*ue(n+u_ldown,l)+                      & 
    197                          wee(n+u_lup,3,1)*le(n+u_rdown)*ue(n+u_rdown,l)+                      & 
    198                          wee(n+u_lup,4,1)*le(n+u_right)*ue(n+u_right,l)+                      & 
    199                          wee(n+u_lup,5,1)*le(n+u_rup)*ue(n+u_rup,l)+                          &  
    200                          wee(n+u_lup,1,2)*le(n+t_lup+u_right)*ue(n+t_lup+u_right,l)+          & 
    201                          wee(n+u_lup,2,2)*le(n+t_lup+u_rup)*ue(n+t_lup+u_rup,l)+              & 
    202                          wee(n+u_lup,3,2)*le(n+t_lup+u_lup)*ue(n+t_lup+u_lup,l)+              & 
    203                          wee(n+u_lup,4,2)*le(n+t_lup+u_left)*ue(n+t_lup+u_left,l)+            & 
    204                          wee(n+u_lup,5,2)*le(n+t_lup+u_ldown)*ue(n+t_lup+u_ldown,l)           & 
    205                                   ) 
    206        v_e(n+u_lup,l,:)= ue(n+u_lup,l)*normal(n+u_lup,:) + up_e*tangent(n+u_lup,:) 
    207      
    208                   up_e=1/de(n+u_ldown)*(                                                    & 
    209                          wee(n+u_ldown,1,1)*le(n+u_rdown)*ue(n+u_rdown,l)+                    & 
    210                          wee(n+u_ldown,2,1)*le(n+u_right)*ue(n+u_right,l)+                    & 
    211                          wee(n+u_ldown,3,1)*le(n+u_rup)*ue(n+u_rup,l)+                        & 
    212                          wee(n+u_ldown,4,1)*le(n+u_lup)*ue(n+u_lup,l)+                        & 
    213                          wee(n+u_ldown,5,1)*le(n+u_left)*ue(n+u_left,l)+                      &  
    214                          wee(n+u_ldown,1,2)*le(n+t_ldown+u_lup)*ue(n+t_ldown+u_lup,l)+        & 
    215                          wee(n+u_ldown,2,2)*le(n+t_ldown+u_left)*ue(n+t_ldown+u_left,l)+      & 
    216                          wee(n+u_ldown,3,2)*le(n+t_ldown+u_ldown)*ue(n+t_ldown+u_ldown,l)+    & 
    217                          wee(n+u_ldown,4,2)*le(n+t_ldown+u_rdown)*ue(n+t_ldown+u_rdown,l)+    & 
    218                          wee(n+u_ldown,5,2)*le(n+t_ldown+u_right)*ue(n+t_ldown+u_right,l)     & 
    219                                       ) 
    220          v_e(n+u_ldown,l,:)= ue(n+u_ldown,l)*normal(n+u_ldown,:) + up_e*tangent(n+u_ldown,:) 
    221  
    222           cc(n+u_right,l,:) = xyz_e(n+u_right,:) - v_e(n+u_right,l,:)*0.5*bigt    !! redge is mid point of edge i  
    223           cc(n+u_lup,l,:)   = xyz_e(n+u_lup,:)   - v_e(n+u_lup,l,:)*0.5*bigt  
    224           cc(n+u_ldown,l,:) = xyz_e(n+u_ldown,:) - v_e(n+u_ldown,l,:)*0.5*bigt  
    225        ENDDO 
    226      ENDDO   
    227     END DO 
    228 !123     continue            
    229 !========================================================================================================== 
    230     DO l = 1,llm 
    231         DO j=jj_begin-1,jj_end+1 
    232          DO i=ii_begin-1,ii_end+1 
    233             n=(j-1)*iim+i 
    234         if (ne(n,right)*ue(n+u_right,l)>0) then  
    235           ed = cc(n+u_right,l,:) - xyz_i(n,:) 
    236          qe(n+u_right,l)=qi(n,l)+sum2(gradq3d(n,l,:),ed)  
    237         else 
    238            ed = cc(n+u_right,l,:) - xyz_i(n+t_right,:) 
    239          qe(n+u_right,l)=qi(n+t_right,l)+sum2(gradq3d(n+t_right,l,:),ed)  
    240         endif 
    241         if (ne(n,lup)*ue(n+u_lup,l)>0) then  
    242            ed = cc(n+u_lup,l,:) - xyz_i(n,:) 
    243          qe(n+u_lup,l)=qi(n,l)+sum2(gradq3d(n,l,:),ed) 
    244         else 
    245            ed = cc(n+u_lup,l,:) - xyz_i(n+t_lup,:) 
    246          qe(n+u_lup,l)=qi(n+t_lup,l)+sum2(gradq3d(n+t_lup,l,:),ed)  
    247         endif 
    248         if (ne(n,ldown)*ue(n+u_ldown,l)>0) then  
    249           ed = cc(n+u_ldown,l,:) - xyz_i(n,:) 
    250          qe(n+u_ldown,l)=qi(n,l)+ sum2(gradq3d(n,l,:),ed)  
    251         else 
    252          ed = cc(n+u_ldown,l,:) - xyz_i(n+t_ldown,:) 
    253         qe(n+u_ldown,l)=qi(n+t_ldown,l)+sum2(gradq3d(n+t_ldown,l,:),ed)  
    254         endif 
    255         end do 
    256     end do 
    257    END DO 
    258  
    259  
    260          DO j=jj_begin-1,jj_end+1 
    261            DO i=ii_begin-1,ii_end+1 
    262             n=(j-1)*iim+i  
    263             flxx(n+u_right,:) = he(n+u_right,:)*qe(n+u_right,:)*ne(n,right) 
    264             flxx(n+u_lup,:) = he(n+u_lup,:)*qe(n+u_lup,:)*ne(n,lup)  
    265             flxx(n+u_ldown,:) = he(n+u_ldown,:)*qe(n+u_ldown,:)*ne(n,ldown) 
    266            ENDDO 
    267          ENDDO 
    268                                  
    269      DO j=jj_begin,jj_end 
    270        DO i=ii_begin,ii_end 
    271         n=(j-1)*iim+i   
    272  
    273          dhi(n,:)= -(1/Ai(n))*(he(n+u_right,:)*ne(n,right) + he(n+u_lup,:)*ne(n,lup) & 
    274                        + he(n+u_ldown,:)*ne(n,ldown) + he(n+u_rup,:)*ne(n,rup) & 
    275                        + he(n+u_left,:)*ne(n,left) +  he(n+u_rdown,:)*ne(n,rdown) )  
    276  
    277          dqi(n,:)= -(1/Ai(n))*(flxx(n+u_right,:)+flxx(n+u_lup,:)       &  
    278                          +flxx(n+u_ldown,:) - flxx(n+u_rup,:)    & 
    279                          - flxx(n+u_left,:) - flxx(n+u_rdown,:) )    
    280          him_old(:) = him(n,:)  
    281          him(n,:) = him(n,:) + dhi(n,:)  
    282          qi(n,:) = (qi(n,:)*him_old(:) + dqi(n,:))/him(n,:)   
    283                         
    284        END DO 
    285      END DO 
    286  
    287         CONTAINS 
    288 !==================================================================================== 
    289           REAL FUNCTION sum2(a1,a2) 
    290         IMPLICIT NONE 
    291         REAL,INTENT(IN):: a1(3), a2(3) 
    292                 sum2 = a1(1)*a2(1)+a1(2)*a2(2)+a1(3)*a2(3) 
    293            END FUNCTION sum2 
    294   
    295                 END SUBROUTINE advect2  
    296 !========================================================================== 
     241          n=(j-1)*iim+i   
     242 
     243          dhi(n,:)= -(1/Ai(n))*(he(n+u_right,:)*ne(n,right) + he(n+u_lup,:)*ne(n,lup) & 
     244               + he(n+u_ldown,:)*ne(n,ldown) + he(n+u_rup,:)*ne(n,rup) & 
     245               + he(n+u_left,:)*ne(n,left) +  he(n+u_rdown,:)*ne(n,rdown) )  
     246 
     247          dqi(n,:)= -(1/Ai(n))*(flxx(n+u_right,:)+flxx(n+u_lup,:)       &  
     248               +flxx(n+u_ldown,:) - flxx(n+u_rup,:)    & 
     249               - flxx(n+u_left,:) - flxx(n+u_rdown,:) )    
     250          him_old(:) = him(n,:)  
     251          him(n,:) = him(n,:) + dhi(n,:)  
     252          qi(n,:) = (qi(n,:)*him_old(:) + dqi(n,:))/him(n,:)   
     253 
     254       END DO 
     255    END DO 
     256 
     257  CONTAINS 
     258    !==================================================================================== 
     259    REAL FUNCTION sum2(a1,a2) 
     260      IMPLICIT NONE 
     261      REAL,INTENT(IN):: a1(3), a2(3) 
     262      sum2 = a1(1)*a2(1)+a1(2)*a2(2)+a1(3)*a2(3) 
     263    END FUNCTION sum2 
     264 
     265  END SUBROUTINE compute_advect_horiz 
     266  !========================================================================== 
    297267  SUBROUTINE gradq(n0,n1,n2,n3,q,dq,det) 
    298   USE geometry 
    299   USE metric 
    300   USE dimensions 
    301   IMPLICIT NONE 
     268    USE geometry 
     269    USE metric 
     270    USE dimensions 
     271    IMPLICIT NONE 
    302272    INTEGER, INTENT(IN) :: n0,n1,n2,n3 
    303273    REAL,INTENT(IN)     :: q(iim*jjm) 
    … …  
    313283    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)  
    314284    A(3,1)=xyz_v(n3,1);              A(3,2)= xyz_v(n3,2);             A(3,3)= xyz_v(n3,3) 
    315     
    316    
     285 
     286 
    317287    dq(1) = q(n1)-q(n0) 
    318288    dq(2) = q(n2)-q(n0) 
    319289    dq(3) = 0.0 
    320 !    CALL DGESV(3,1,A,3,IPIV,dq(:),3,info) 
    321         CALL determinant(A(:,1),A(:,2),A(:,3),det) 
    322         CALL determinant(dq,A(:,2),A(:,3),detx) 
    323         CALL determinant(A(:,1),dq,A(:,3),dety) 
    324         CALL determinant(A(:,1),A(:,2),dq,detz) 
    325         dq(1) = detx 
    326         dq(2) = dety 
    327         dq(3) = detz  
     290    !    CALL DGESV(3,1,A,3,IPIV,dq(:),3,info) 
     291    CALL determinant(A(:,1),A(:,2),A(:,3),det) 
     292    CALL determinant(dq,A(:,2),A(:,3),detx) 
     293    CALL determinant(A(:,1),dq,A(:,3),dety) 
     294    CALL determinant(A(:,1),A(:,2),dq,detz) 
     295    dq(1) = detx 
     296    dq(2) = dety 
     297    dq(3) = detz  
    328298  END SUBROUTINE gradq 
    329 !========================================================================== 
    330           SUBROUTINE determinant(a1,a2,a3,det) 
    331         IMPLICIT NONE 
    332         REAL, DIMENSION(3) :: a1, a2,a3 
    333         REAL ::  det,x1,x2,x3 
    334         x1 =  a1(1) * (a2(2) * a3(3) - a2(3) * a3(2)) 
    335         x2 =  a1(2) * (a2(1) * a3(3) - a2(3) * a3(1)) 
    336         x3 =  a1(3) * (a2(1) * a3(2) - a2(2) * a3(1)) 
    337         det =  x1 - x2 + x3 
    338         END SUBROUTINE     
    339         END MODULE   
     299  !========================================================================== 
     300  SUBROUTINE determinant(a1,a2,a3,det) 
     301    IMPLICIT NONE 
     302    REAL, DIMENSION(3) :: a1, a2,a3 
     303    REAL ::  det,x1,x2,x3 
     304    x1 =  a1(1) * (a2(2) * a3(3) - a2(3) * a3(2)) 
     305    x2 =  a1(2) * (a2(1) * a3(3) - a2(3) * a3(1)) 
     306    x3 =  a1(3) * (a2(1) * a3(2) - a2(2) * a3(1)) 
     307    det =  x1 - x2 + x3 
     308  END SUBROUTINE determinant 
     309 
     310END MODULE advect_mod 
Note: See TracChangeset for help on using the changeset viewer.