Changeset 22


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

Location:
codes/icosagcm/trunk/src
Files:
4 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 

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

    r19 r22  
    44  INTEGER,PARAMETER::iapp_tracvl= 3 
    55  REAL(rstd),SAVE :: dt 
    6    
     6 
     7  TYPE(t_field),POINTER :: f_normal(:) 
     8  TYPE(t_field),POINTER :: f_tangent(:) 
     9  TYPE(t_field),POINTER :: f_gradq3d(:) 
     10 
    711  PUBLIC init_advect_tracer, advect_tracer 
    812 
    913CONTAINS 
    10    
     14 
    1115  SUBROUTINE init_advect_tracer(dt_in) 
    12   USE advect_mod 
    13   IMPLICIT NONE 
     16    USE advect_mod 
     17    IMPLICIT NONE 
    1418    REAL(rstd),INTENT(IN) :: dt_in 
    15      
     19    REAL(rstd),POINTER :: tangent(:,:) 
     20    REAL(rstd),POINTER :: normal(:,:) 
     21 
    1622    dt=dt_in 
    17      
    18     CALL init_advect 
    19      
     23    CALL allocate_field(f_normal,field_u,type_real,3) 
     24    CALL allocate_field(f_tangent,field_u,type_real,3) 
     25    CALL allocate_field(f_gradq3d,field_t,type_real,llm,3) 
     26 
     27    DO ind=1,ndomain 
     28       CALL swap_dimensions(ind) 
     29       CALL swap_geometry(ind) 
     30       normal=f_normal(ind) 
     31       tangent=f_tangent(ind) 
     32       CALL init_advect(normal,tangent) 
     33    END DO 
     34 
    2035  END SUBROUTINE init_advect_tracer 
    21    
    22    
    23   SUBROUTINE advect_tracer(f_ps,f_u, f_q) 
    24   USE icosa 
    25   USE advect_mod 
    26   USE disvert_mod 
    27   IMPLICIT NONE 
    28   TYPE(t_field),POINTER :: f_ps(:) 
    29   TYPE(t_field),POINTER :: f_u(:) 
    30   TYPE(t_field),POINTER :: f_q(:) 
    31   REAL(rstd),POINTER :: q(:,:,:) 
    32   REAL(rstd),POINTER :: u(:,:)  
    33   REAL(rstd),POINTER :: ps(:)  
    34   REAL(rstd),POINTER :: massflx(:,:) 
    35   REAL(rstd),POINTER :: rhodz(:,:) 
    36   TYPE(t_field),POINTER,SAVE :: f_massflxc(:) 
    37   TYPE(t_field),POINTER,SAVE :: f_massflx(:) 
    38   TYPE(t_field),POINTER,SAVE :: f_uc(:) 
    39   TYPE(t_field),POINTER,SAVE :: f_rhodzm1(:) 
    40   TYPE(t_field),POINTER,SAVE :: f_rhodz(:) 
    41   REAL(rstd),POINTER,SAVE :: massflxc(:,:) 
    42   REAL(rstd),POINTER,SAVE :: uc(:,:) 
    43   REAL(rstd),POINTER,SAVE :: rhodzm1(:,:) 
    44   REAL(rstd):: bigt  
    45   INTEGER :: ind,it,iapptrac,i,j,l,ij 
    46   INTEGER,SAVE :: iadvtr=0 
    47   LOGICAL,SAVE:: first=.TRUE.  
    48 !------------------------------------------------------sarvesh 
     36 
     37  SUBROUTINE advect_tracer(f_ps,f_u,f_q) 
     38    USE icosa 
     39    USE advect_mod 
     40    USE disvert_mod 
     41    IMPLICIT NONE 
     42    TYPE(t_field),POINTER :: f_ps(:) 
     43    TYPE(t_field),POINTER :: f_u(:) 
     44    TYPE(t_field),POINTER :: f_q(:) 
     45    REAL(rstd),POINTER :: q(:,:,:) 
     46    REAL(rstd),POINTER :: u(:,:)  
     47    REAL(rstd),POINTER :: ps(:)  
     48    REAL(rstd),POINTER :: massflx(:,:) 
     49    REAL(rstd),POINTER :: rhodz(:,:) 
     50    TYPE(t_field),POINTER,SAVE :: f_massflxc(:) 
     51    TYPE(t_field),POINTER,SAVE :: f_massflx(:) 
     52    TYPE(t_field),POINTER,SAVE :: f_uc(:) 
     53    TYPE(t_field),POINTER,SAVE :: f_rhodzm1(:) 
     54    TYPE(t_field),POINTER,SAVE :: f_rhodz(:) 
     55    REAL(rstd),POINTER,SAVE :: massflxc(:,:) 
     56    REAL(rstd),POINTER,SAVE :: uc(:,:) 
     57    REAL(rstd),POINTER,SAVE :: rhodzm1(:,:) 
     58    REAL(rstd):: bigt 
     59    INTEGER :: ind,it,i,j,l,ij 
     60    INTEGER,SAVE :: iadvtr=0 
     61    LOGICAL,SAVE:: first=.TRUE.  
     62    !------------------------------------------------------sarvesh 
    4963    IF ( first ) THEN  
    50       CALL allocate_field(f_rhodz,field_t,type_real,llm)  
    51       CALL allocate_field(f_rhodzm1,field_t,type_real,llm)  
    52       CALL allocate_field(f_massflxc,field_u,type_real,llm)  
    53       CALL allocate_field(f_massflx,field_u,type_real,llm)  
    54       CALL allocate_field(f_uc,field_u,type_real,llm)  
    55       first = .FALSE. 
    56     END IF  
    57  
    58     DO ind=1,ndomain 
    59       CALL swap_dimensions(ind) 
    60       CALL swap_geometry(ind) 
    61       rhodz=f_rhodz(ind) 
    62       massflx=f_massflx(ind) 
    63       ps=f_ps(ind) 
    64       u=f_u(ind) 
    65        
    66       DO l = 1, llm 
    67         DO j=jj_begin-1,jj_end+1 
    68           DO i=ii_begin-1,ii_end+1 
    69             ij=(j-1)*iim+i 
    70             rhodz(ij,l) = (ap(l) - ap(l+1) + (bp(l)-bp(l+1))*ps(ij))/g  
    71           ENDDO 
    72         ENDDO 
    73       ENDDO 
    74  
    75       DO l = 1, llm 
    76         DO j=jj_begin-1,jj_end+1 
    77           DO i=ii_begin-1,ii_end+1 
    78             ij=(j-1)*iim+i 
    79             massflx(ij+u_right,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_right,l))*u(ij+u_right,l)*le(ij+u_right) 
    80             massflx(ij+u_lup,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_lup,l))*u(ij+u_lup,l)*le(ij+u_lup) 
    81             massflx(ij+u_ldown,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_ldown,l))*u(ij+u_ldown,l)*le(ij+u_ldown) 
    82           ENDDO 
    83         ENDDO 
    84       ENDDO 
    85      
    86     ENDDO 
    87      
    88  
    89      
    90         IF ( iadvtr == 0 ) THEN  
    91          DO ind=1,ndomain          
    92          CALL swap_dimensions(ind) 
    93          CALL swap_geometry(ind) 
    94            rhodz=f_rhodz(ind)   
    95            rhodzm1 = f_rhodzm1(ind)  
    96            massflxc = f_massflxc(ind)  
    97            rhodzm1 = rhodz 
    98            massflxc = 0.0  
    99            uc = f_uc(ind) 
    100            uc = 0.0  
    101          END DO 
    102          CALL transfert_request(f_rhodzm1,req_i1)   !----> 
    103          CALL transfert_request(f_massflxc,req_e1) !----> 
    104          CALL transfert_request(f_massflxc,req_e1) !------> 
    105          CALL transfert_request(f_uc,req_e1) !----> 
    106          CALL transfert_request(f_uc,req_e1)  
    107         END IF 
    108  
    109         iadvtr = iadvtr + 1  
    110         iapptrac = iadvtr   
    111  
    112       DO ind=1,ndomain 
    113        CALL swap_dimensions(ind) 
    114        CALL swap_geometry(ind) 
    115             massflx=f_massflx(ind) 
    116             rhodzm1 = f_rhodzm1(ind)  
    117             massflxc = f_massflxc(ind)  
    118             massflxc = massflxc + massflx  
    119              uc = f_uc(ind)  
    120              u  = f_u(ind)  
    121              uc = uc + u  
     64       CALL allocate_field(f_rhodz,field_t,type_real,llm)  
     65       CALL allocate_field(f_rhodzm1,field_t,type_real,llm)  
     66       CALL allocate_field(f_massflxc,field_u,type_real,llm)  
     67       CALL allocate_field(f_massflx,field_u,type_real,llm)  
     68       CALL allocate_field(f_uc,field_u,type_real,llm)  
     69       first = .FALSE. 
     70    END IF 
     71 
     72    DO ind=1,ndomain 
     73       CALL swap_dimensions(ind) 
     74       CALL swap_geometry(ind) 
     75       rhodz=f_rhodz(ind) 
     76       massflx=f_massflx(ind) 
     77       ps=f_ps(ind) 
     78       u=f_u(ind) 
     79 
     80       DO l = 1, llm 
     81          DO j=jj_begin-1,jj_end+1 
     82             DO i=ii_begin-1,ii_end+1 
     83                ij=(j-1)*iim+i 
     84                rhodz(ij,l) = (ap(l) - ap(l+1) + (bp(l)-bp(l+1))*ps(ij))/g  
     85             ENDDO 
     86          ENDDO 
     87       ENDDO 
     88 
     89       DO l = 1, llm 
     90          DO j=jj_begin-1,jj_end+1 
     91             DO i=ii_begin-1,ii_end+1 
     92                ij=(j-1)*iim+i 
     93                massflx(ij+u_right,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_right,l))*u(ij+u_right,l)*le(ij+u_right) 
     94                massflx(ij+u_lup,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_lup,l))*u(ij+u_lup,l)*le(ij+u_lup) 
     95                massflx(ij+u_ldown,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_ldown,l))*u(ij+u_ldown,l)*le(ij+u_ldown) 
     96             ENDDO 
     97          ENDDO 
     98       ENDDO 
     99    ENDDO 
     100 
     101    IF ( iadvtr == 0 ) THEN  
     102       DO ind=1,ndomain          
     103          CALL swap_dimensions(ind) 
     104          CALL swap_geometry(ind) 
     105          rhodz=f_rhodz(ind)   
     106          rhodzm1 = f_rhodzm1(ind)  
     107          massflxc = f_massflxc(ind)  ! accumulated mass fluxes 
     108          uc = f_uc(ind)              ! time-integrated normal velocity to compute back-trajectory (Miura) 
     109          rhodzm1 = rhodz 
     110          massflxc = 0.0  
     111          uc = 0.0  
    122112       END DO 
    123  
    124      IF ( iadvtr == iapp_tracvl ) THEN  
    125          bigt = dt*iapp_tracvl 
    126          DO ind=1,ndomain 
    127            CALL swap_dimensions(ind) 
    128            CALL swap_geometry(ind) 
    129               uc = f_uc(ind)  
    130               uc = uc/real(iapp_tracvl)  
    131          END DO 
     113       CALL transfert_request(f_rhodzm1,req_i1)   !----> 
     114       CALL transfert_request(f_massflxc,req_e1) !----> 
     115       CALL transfert_request(f_massflxc,req_e1) !------> 
     116       CALL transfert_request(f_uc,req_e1) !----> 
     117       CALL transfert_request(f_uc,req_e1)  
     118    END IF 
     119 
     120    iadvtr = iadvtr + 1  
     121 
     122    DO ind=1,ndomain 
     123       CALL swap_dimensions(ind) 
     124       CALL swap_geometry(ind) 
     125       massflx  = f_massflx(ind) 
     126       massflxc = f_massflxc(ind)  
     127       uc = f_uc(ind)  
     128       u  = f_u(ind)  
     129       massflxc = massflxc + massflx  
     130       uc = uc + u  
     131    END DO 
     132 
     133    IF ( iadvtr == iapp_tracvl ) THEN  
     134       bigt = dt*iapp_tracvl 
     135       DO ind=1,ndomain 
     136          CALL swap_dimensions(ind) 
     137          CALL swap_geometry(ind) 
     138          uc = f_uc(ind)  
     139          uc = uc/real(iapp_tracvl)  
     140          massflxc = f_massflx(ind) 
     141          massflxc = massflxc*dt 
     142          ! now massflx is time-integrated 
     143       END DO 
    132144 
    133145       CALL vlsplt(f_q,f_rhodzm1,f_massflxc,2.0,f_uc,bigt)  
    134             iadvtr = 0 
    135      END IF 
     146       iadvtr = 0 
     147    END IF 
    136148  END SUBROUTINE advect_tracer 
    137    
    138 !==============================================================================   
    139   SUBROUTINE advtrac(massflx,wgg) 
    140   USE domain_mod 
    141   USE dimensions 
    142   USE grid_param 
    143   USE geometry 
    144   USE metric 
    145   USE disvert_mod 
    146   IMPLICIT NONE 
    147     REAL(rstd),INTENT(IN) :: massflx(iim*3*jjm,llm) 
    148     REAL(rstd),INTENT(OUT) :: wgg(iim*jjm,llm) 
    149       
    150     INTEGER :: i,j,ij,l 
    151     REAL(rstd) :: convm(iim*jjm,llm)  
    152     
    153      DO l = 1, llm 
    154       DO j=jj_begin,jj_end 
    155        DO i=ii_begin,ii_end 
    156         ij=(j-1)*iim+i 
    157         convm(ij,l)= 1/(Ai(ij))*(ne(ij,right)*massflx(ij+u_right,l)   +  & 
    158                                 ne(ij,rup)*massflx(ij+u_rup,l)       +  & 
    159                                 ne(ij,lup)*massflx(ij+u_lup,l)       +  & 
    160                                 ne(ij,left)*massflx(ij+u_left,l)     +  & 
    161                                 ne(ij,ldown)*massflx(ij+u_ldown,l)   +  & 
    162                                 ne(ij,rdown)*massflx(ij+u_rdown,l)) 
    163         ENDDO 
    164       ENDDO 
    165      ENDDO 
    166  
    167     DO  l = llm-1, 1, -1 
    168      DO j=jj_begin,jj_end 
    169       DO i=ii_begin,ii_end 
    170         ij=(j-1)*iim+i 
    171         convm(ij,l) = convm(ij,l) + convm(ij,l+1) 
    172       ENDDO 
    173      ENDDO 
    174     ENDDO 
    175  
    176     !!! Compute vertical velocity 
    177   DO l = 1,llm-1 
    178     DO j=jj_begin,jj_end 
    179       DO i=ii_begin,ii_end 
    180         ij=(j-1)*iim+i 
    181         wgg( ij, l+1 ) = (convm( ij, l+1 ) - bp(l+1) * convm( ij, 1 ))  
    182       ENDDO 
    183     ENDDO 
    184   ENDDO 
    185  
    186     DO j=jj_begin,jj_end 
    187      DO i=ii_begin,ii_end 
    188       ij=(j-1)*iim+i 
    189       wgg(ij,1)  = 0. 
    190      ENDDO 
    191    ENDDO 
    192  END SUBROUTINE advtrac  
    193  
    194         SUBROUTINE vlsplt(f_q,f_rhodz,f_massflx,pente_max,f_u,bigt)  
    195   USE field_mod 
    196   USE domain_mod 
    197   USE dimensions 
    198   USE grid_param 
    199   USE geometry 
    200   USE metric 
    201   USE advect_mod 
    202   IMPLICIT NONE 
    203  
    204   TYPE(t_field),POINTER :: f_q(:) 
    205   TYPE(t_field),POINTER :: f_u(:) 
    206   TYPE(t_field),POINTER :: f_rhodz(:)  
    207   TYPE(t_field),POINTER :: f_massflx(:) 
    208   TYPE(t_field),POINTER,SAVE :: f_wg(:) 
    209   TYPE(t_field),POINTER,SAVE :: f_zm(:) 
    210   TYPE(t_field),POINTER,SAVE :: f_zq(:) 
    211  
    212   REAL(rstd)::bigt,dt 
    213   REAL(rstd),POINTER :: q(:,:,:) 
    214   REAL(rstd),POINTER :: u(:,:) 
    215   REAL(rstd),POINTER :: rhodz(:,:)  
    216   REAL(rstd),POINTER :: massflx(:,:)  
    217   REAL(rstd),POINTER,SAVE :: wg(:,:) 
    218   REAL(rstd),POINTER,SAVE::zq(:,:,:)  
    219   REAL(rstd),POINTER,SAVE::zm(:,:)  
    220  
    221   REAL(rstd)::pente_max 
    222   LOGICAL,SAVE::first = .true.  
    223   INTEGER :: i,ij,l,j,ind,k 
    224   REAL(rstd) :: zzpbar, zzw  
    225   REAL::qvmax,qvmin  
    226    
    227      IF ( first ) THEN  
     149 
     150  SUBROUTINE vlsplt(f_q,f_rhodz,f_massflx,pente_max,f_u,bigt) 
     151    USE field_mod 
     152    USE domain_mod 
     153    USE dimensions 
     154    USE grid_param 
     155    USE geometry 
     156    USE metric 
     157    USE advect_mod 
     158    IMPLICIT NONE 
     159 
     160    TYPE(t_field),POINTER :: f_q(:) 
     161    TYPE(t_field),POINTER :: f_u(:) 
     162    TYPE(t_field),POINTER :: f_rhodz(:)  
     163    TYPE(t_field),POINTER :: f_massflx(:) 
     164 
     165    TYPE(t_field),POINTER,SAVE :: f_wg(:) 
     166    TYPE(t_field),POINTER,SAVE :: f_zm(:) 
     167    TYPE(t_field),POINTER,SAVE :: f_zq(:) 
     168 
     169    REAL(rstd)::bigt,dt 
     170    REAL(rstd),POINTER :: q(:,:,:) 
     171    REAL(rstd),POINTER :: u(:,:) 
     172    REAL(rstd),POINTER :: rhodz(:,:)  
     173    REAL(rstd),POINTER :: massflx(:,:)  
     174    REAL(rstd),POINTER :: wg(:,:) 
     175    REAL(rstd),POINTER :: zq(:,:,:)  
     176    REAL(rstd),POINTER :: zm(:,:)  
     177    REAL(rstd),POINTER :: tangent(:,:) 
     178    REAL(rstd),POINTER :: normal(:,:) 
     179    REAL(rstd),POINTER :: gradq3d(:,:,:) 
     180 
     181    REAL(rstd)::pente_max 
     182    LOGICAL,SAVE::first = .true.  
     183    INTEGER :: i,ij,l,j,ind,k 
     184    REAL(rstd) :: zzpbar, zzw  
     185    REAL::qvmax,qvmin  
     186 
     187    IF ( first ) THEN  
    228188       CALL allocate_field(f_wg,field_t,type_real,llm) 
    229189       CALL allocate_field(f_zm,field_t,type_real,llm) 
    230190       CALL allocate_field(f_zq,field_t,type_real,llm,nqtot)  
    231191       first = .FALSE. 
    232       END IF  
    233  
    234      DO ind=1,ndomain 
    235        CALL swap_dimensions(ind) 
    236        CALL swap_geometry(ind) 
    237          q=f_q(ind)  
    238          rhodz=f_rhodz(ind)  
    239          zq=f_zq(ind) 
    240         zm=f_zm(ind)  
    241          zm = rhodz ; zq = q   
    242          wg = f_wg(ind) 
    243         wg = 0.0 
    244          massflx=f_massflx(ind)  
    245        CALL advtrac(massflx,wg)  
    246      END DO 
    247  
    248 !   CALL transfert_request(f_wg,req_i1)  
    249  
    250       DO ind=1,ndomain 
    251         CALL swap_dimensions(ind) 
    252         CALL swap_geometry(ind) 
    253          zq=f_zq(ind) 
    254         zm=f_zm(ind)  
    255          wg=f_wg(ind) 
    256         wg=wg*0.5 
     192    END IF 
     193 
     194    DO ind=1,ndomain 
     195       CALL swap_dimensions(ind) 
     196       CALL swap_geometry(ind) 
     197       q=f_q(ind)  
     198       rhodz=f_rhodz(ind)  
     199       zq=f_zq(ind) 
     200      zm=f_zm(ind)  
     201       zm = rhodz ; zq = q   
     202       wg = f_wg(ind) 
     203      wg = 0.0 
     204       massflx=f_massflx(ind)  
     205       CALL advtrac(massflx,wg) ! compute vertical mass fluxes 
     206    END DO 
     207 
     208    !   CALL transfert_request(f_wg,req_i1)  
     209 
     210    DO ind=1,ndomain 
     211       CALL swap_dimensions(ind) 
     212       CALL swap_geometry(ind) 
     213       zq=f_zq(ind) 
     214      zm=f_zm(ind)  
     215       wg=f_wg(ind) 
     216      wg=wg*0.5 
    257217       DO k = 1, nqtot 
    258         CALL vlz(zq(:,:,k),2.,zm,wg) 
     218          CALL vlz(zq(:,:,k),2.,zm,wg) 
    259219       END DO 
    260       END DO 
    261  
    262         DO ind=1,ndomain 
    263           CALL swap_dimensions(ind) 
    264           CALL swap_geometry(ind) 
    265           CALL swap_advect(ind) 
    266           zq=f_zq(ind) 
    267           zq = f_zq(ind) 
    268           zm = f_zm(ind) 
    269           massflx =f_massflx(ind) 
    270           u = f_u(ind)  
    271          DO k = 1,nqtot 
    272            CALL advect1(zq(:,:,k)) 
    273            CALL advect2(zq(:,:,k),zm,u,massflx,bigt)  
    274          END DO 
    275         END DO 
    276  
    277         DO ind=1,ndomain 
    278            CALL swap_dimensions(ind) 
    279            CALL swap_geometry(ind) 
    280             q = f_q(ind) 
    281            zq = f_zq(ind)  
    282            zm = f_zm(ind)  
    283            wg = f_wg(ind) 
    284           DO k = 1,nqtot 
    285             CALL vlz(zq(:,:,k),2.,zm,wg) 
    286           END DO  
    287             q = zq  
    288         END DO 
    289  
    290   END SUBROUTINE vlsplt  
    291  
    292          
    293      SUBROUTINE vlz(q,pente_max,masse,wgg) 
    294 !c 
    295 !c     Auteurs:   P.Le Van, F.Hourdin, F.Forget  
    296 !c 
    297 !c    ******************************************************************** 
    298 !c     Shema  d'advection " pseudo amont " . 
    299 !c    ******************************************************************** 
    300       USE icosa 
    301       IMPLICIT NONE 
    302 !c 
    303 !c   Arguments: 
    304 !c   ---------- 
    305       REAL masse(iim*jjm,llm),pente_max 
    306       REAL q(iim*jjm,llm) 
    307       REAL wgg(iim*jjm,llm),w(iim*jjm,llm+1)  
    308       REAL dq(iim*jjm,llm) 
    309       INTEGER i,ij,l,j,ii 
    310 !c 
    311       REAL wq(iim*jjm,llm+1),newmasse 
    312  
    313       REAL dzq(iim*jjm,llm),dzqw(iim*jjm,llm),adzqw(iim*jjm,llm),dzqmax 
    314       REAL sigw 
    315  
    316       REAL      SSUM 
    317  
    318  
    319          w(:,1:llm) = wgg(:,:) 
    320          w(:,llm+1) = 0.0  
    321  
    322 !c    On oriente tout dans le sens de la pression c'est a dire dans le 
    323 !c    sens de W 
    324  
    325        DO l=2,llm 
    326         DO j=jj_begin,jj_end 
    327          DO i=ii_begin,ii_end 
    328             ij=(j-1)*iim+i 
    329             dzqw(ij,l)=q(ij,l-1)-q(ij,l) 
    330             adzqw(ij,l)=abs(dzqw(ij,l)) 
    331          ENDDO 
    332         ENDDO 
    333        ENDDO 
    334  
    335        DO l=2,llm-1 
    336         DO j=jj_begin,jj_end 
    337          DO i=ii_begin,ii_end 
    338             ij=(j-1)*iim+i 
    339            IF(dzqw(ij,l)*dzqw(ij,l+1).gt.0.) THEN 
     220    END DO 
     221 
     222    DO ind=1,ndomain 
     223       CALL swap_dimensions(ind) 
     224       CALL swap_geometry(ind) 
     225       zq=f_zq(ind) 
     226       zq = f_zq(ind) 
     227       zm = f_zm(ind) 
     228       massflx =f_massflx(ind) 
     229       u = f_u(ind) 
     230       tangent = f_tangent(ind) 
     231       normal = f_normal(ind) 
     232       gradq3d = f_gradq3d(ind) 
     233 
     234       DO k = 1,nqtot 
     235          CALL compute_gradq3d(zq(:,:,k),gradq3d) 
     236          CALL compute_advect_horiz(zq(:,:,k),zm,u,massflx,bigt)  
     237       END DO 
     238    END DO 
     239 
     240    DO ind=1,ndomain 
     241       CALL swap_dimensions(ind) 
     242       CALL swap_geometry(ind) 
     243       q = f_q(ind) 
     244       zq = f_zq(ind)  
     245       zm = f_zm(ind)  
     246       wg = f_wg(ind) 
     247       DO k = 1,nqtot 
     248          CALL vlz(zq(:,:,k),2.,zm,wg) 
     249       END DO 
     250       q = zq  
     251    END DO 
     252 
     253  END SUBROUTINE vlsplt 
     254 
     255  !==============================================================================   
     256  SUBROUTINE advtrac(massflx,wgg) 
     257    USE domain_mod 
     258    USE dimensions 
     259    USE grid_param 
     260    USE geometry 
     261    USE metric 
     262    USE disvert_mod 
     263    IMPLICIT NONE 
     264    REAL(rstd),INTENT(IN) :: massflx(iim*3*jjm,llm) 
     265    REAL(rstd),INTENT(OUT) :: wgg(iim*jjm,llm) 
     266 
     267    INTEGER :: i,j,ij,l 
     268    REAL(rstd) :: convm(iim*jjm,llm)  
     269 
     270    DO l = 1, llm 
     271       DO j=jj_begin,jj_end 
     272          DO i=ii_begin,ii_end 
     273             ij=(j-1)*iim+i 
     274             ! divergence of horizontal flux 
     275             convm(ij,l)= 1/(Ai(ij))*(ne(ij,right)*massflx(ij+u_right,l)   +  & 
     276                  ne(ij,rup)*massflx(ij+u_rup,l)       +  & 
     277                  ne(ij,lup)*massflx(ij+u_lup,l)       +  & 
     278                  ne(ij,left)*massflx(ij+u_left,l)     +  & 
     279                  ne(ij,ldown)*massflx(ij+u_ldown,l)   +  & 
     280                  ne(ij,rdown)*massflx(ij+u_rdown,l)) 
     281          ENDDO 
     282       ENDDO 
     283    ENDDO 
     284 
     285    ! accumulate divergence from top of model 
     286    DO  l = llm-1, 1, -1 
     287       DO j=jj_begin,jj_end 
     288          DO i=ii_begin,ii_end 
     289             ij=(j-1)*iim+i 
     290             convm(ij,l) = convm(ij,l) + convm(ij,l+1) 
     291          ENDDO 
     292       ENDDO 
     293    ENDDO 
     294 
     295!!! Compute vertical velocity 
     296    DO l = 1,llm-1 
     297       DO j=jj_begin,jj_end 
     298          DO i=ii_begin,ii_end 
     299             ij=(j-1)*iim+i 
     300             wgg( ij, l+1 ) = (convm( ij, l+1 ) - bp(l+1) * convm( ij, 1 ))  
     301          ENDDO 
     302       ENDDO 
     303    ENDDO 
     304 
     305    DO j=jj_begin,jj_end 
     306       DO i=ii_begin,ii_end 
     307          ij=(j-1)*iim+i 
     308          wgg(ij,1)  = 0. 
     309       ENDDO 
     310    ENDDO 
     311  END SUBROUTINE advtrac 
     312 
     313  SUBROUTINE vlz(q,pente_max,masse,wgg) 
     314    !c 
     315    !c     Auteurs:   P.Le Van, F.Hourdin, F.Forget  
     316    !c 
     317    !c    ******************************************************************** 
     318    !c     Shema  d'advection " pseudo amont " . 
     319    !c    ******************************************************************** 
     320    USE icosa 
     321    IMPLICIT NONE 
     322    !c 
     323    !c   Arguments: 
     324    !c   ---------- 
     325    REAL masse(iim*jjm,llm),pente_max 
     326    REAL q(iim*jjm,llm) 
     327    REAL wgg(iim*jjm,llm),w(iim*jjm,llm+1)  
     328    REAL dq(iim*jjm,llm) 
     329    INTEGER i,ij,l,j,ii 
     330    !c 
     331    REAL wq(iim*jjm,llm+1),newmasse 
     332 
     333    REAL dzq(iim*jjm,llm),dzqw(iim*jjm,llm),adzqw(iim*jjm,llm),dzqmax 
     334    REAL sigw 
     335 
     336    REAL      SSUM 
     337 
     338 
     339    w(:,1:llm) = -wgg(:,:)  ! w>0 for downward transport 
     340    w(:,llm+1) = 0.0  
     341 
     342    !c    On oriente tout dans le sens de la pression c'est a dire dans le 
     343    !c    sens de W 
     344 
     345    DO l=2,llm 
     346       DO j=jj_begin,jj_end 
     347          DO i=ii_begin,ii_end 
     348             ij=(j-1)*iim+i 
     349             dzqw(ij,l)=q(ij,l-1)-q(ij,l) 
     350             adzqw(ij,l)=abs(dzqw(ij,l)) 
     351          ENDDO 
     352       ENDDO 
     353    ENDDO 
     354 
     355    DO l=2,llm-1 
     356       DO j=jj_begin,jj_end 
     357          DO i=ii_begin,ii_end 
     358             ij=(j-1)*iim+i 
     359             IF(dzqw(ij,l)*dzqw(ij,l+1).gt.0.) THEN 
    340360                dzq(ij,l)=0.5*(dzqw(ij,l)+dzqw(ij,l+1)) 
    341            ELSE 
     361             ELSE 
    342362                dzq(ij,l)=0. 
    343            ENDIF 
    344             dzqmax=pente_max*min(adzqw(ij,l),adzqw(ij,l+1)) 
    345             dzq(ij,l)=sign(min(abs(dzq(ij,l)),dzqmax),dzq(ij,l)) 
    346          ENDDO 
    347         ENDDO 
    348         ENDDO 
    349  
    350        DO l=2,llm-1 
    351         DO j=jj_begin,jj_end 
    352          DO i=ii_begin,ii_end 
    353             ij=(j-1)*iim+i 
    354             dzq(ij,1)=0. 
    355            dzq(ij,llm)=0. 
    356          ENDDO 
    357         ENDDO 
    358       ENDDO 
    359 !c --------------------------------------------------------------- 
    360 !c   .... calcul des termes d'advection verticale  ....... 
    361 !c --------------------------------------------------------------- 
    362  
    363 !c calcul de  - d( q   * w )/ d(sigma)    qu'on ajoute a  dq pour calculer dq 
    364  
    365         DO l = 1,llm-1 
    366            DO j=jj_begin,jj_end 
    367            DO i=ii_begin,ii_end 
     363             ENDIF 
     364             dzqmax=pente_max*min(adzqw(ij,l),adzqw(ij,l+1)) 
     365             dzq(ij,l)=sign(min(abs(dzq(ij,l)),dzqmax),dzq(ij,l)) 
     366          ENDDO 
     367       ENDDO 
     368    ENDDO 
     369 
     370    DO l=2,llm-1 
     371       DO j=jj_begin,jj_end 
     372          DO i=ii_begin,ii_end 
     373             ij=(j-1)*iim+i 
     374             dzq(ij,1)=0. 
     375             dzq(ij,llm)=0. 
     376          ENDDO 
     377      ENDDO 
     378    ENDDO 
     379    !c --------------------------------------------------------------- 
     380    !c   .... calcul des termes d'advection verticale  ....... 
     381    !c --------------------------------------------------------------- 
     382 
     383    !c calcul de  - d( q   * w )/ d(sigma)    qu'on ajoute a  dq pour calculer dq 
     384 
     385    DO l = 1,llm-1 
     386       DO j=jj_begin,jj_end 
     387          DO i=ii_begin,ii_end 
    368388             ij=(j-1)*iim+i 
    369389             IF(w(ij,l+1).gt.0.) THEN 
    370              sigw=w(ij,l+1)/masse(ij,l+1) 
    371              wq(ij,l+1)=w(ij,l+1)*(q(ij,l+1)+0.5*(1.-sigw)*dzq(ij,l+1)) 
     390                ! upwind only if downward transport 
     391                sigw=w(ij,l+1)/masse(ij,l+1) 
     392                wq(ij,l+1)=w(ij,l+1)*(q(ij,l+1)+0.5*(1.-sigw)*dzq(ij,l+1)) 
    372393             ELSE 
    373              sigw=w(ij,l+1)/masse(ij,l) 
    374              wq(ij,l+1)=w(ij,l+1)*(q(ij,l)-0.5*(1.+sigw)*dzq(ij,l)) 
    375             ENDIF 
    376            ENDDO 
    377          ENDDO 
    378          END DO 
    379  
    380           DO j=jj_begin,jj_end 
    381            DO i=ii_begin,ii_end 
    382              ij=(j-1)*iim+i 
    383              wq(ij,llm+1)=0. 
    384              wq(ij,1)=0. 
    385            ENDDO 
    386           END DO 
    387  
    388       DO l=1,llm 
    389        DO j=jj_begin,jj_end 
    390          DO i=ii_begin,ii_end 
    391              ij=(j-1)*iim+i 
    392             newmasse=masse(ij,l)+(w(ij,l+1)-w(ij,l)) 
    393             dq(ij,l) = (wq(ij,l+1)-wq(ij,l)) !================>>>> 
    394             q(ij,l)=(q(ij,l)*masse(ij,l)+wq(ij,l+1)-wq(ij,l))/newmasse     
    395                           
    396             masse(ij,l)=newmasse 
    397            ENDDO 
    398         ENDDO 
    399       END DO     
    400       RETURN 
    401       END SUBROUTINE vlz 
     394                ! upwind only if upward transport 
     395                sigw=w(ij,l+1)/masse(ij,l) 
     396                wq(ij,l+1)=w(ij,l+1)*(q(ij,l)-0.5*(1.+sigw)*dzq(ij,l)) 
     397             ENDIF 
     398          ENDDO 
     399       ENDDO 
     400    END DO 
     401 
     402    DO j=jj_begin,jj_end 
     403       DO i=ii_begin,ii_end 
     404          ij=(j-1)*iim+i 
     405          wq(ij,llm+1)=0. 
     406          wq(ij,1)=0. 
     407       ENDDO 
     408    END DO 
     409 
     410    DO l=1,llm 
     411       DO j=jj_begin,jj_end 
     412          DO i=ii_begin,ii_end 
     413             ij=(j-1)*iim+i 
     414             ! masse -= dw/dz but w>0 <=> downward 
     415             newmasse=masse(ij,l)+(w(ij,l+1)-w(ij,l))  
     416!             dq(ij,l) = (wq(ij,l+1)-wq(ij,l)) !================>>>> 
     417             q(ij,l)=(q(ij,l)*masse(ij,l)+wq(ij,l+1)-wq(ij,l))/newmasse     
     418!             masse(ij,l)=newmasse 
     419          ENDDO 
     420       ENDDO 
     421    END DO 
     422    RETURN 
     423  END SUBROUTINE vlz 
    402424 
    403425END MODULE advect_tracer_mod 

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

    r21 r22  
    438438        ij=(j-1)*iim+i   
    439439! signe ? attention d (rho theta dz) 
     440        ! dtheta_rhodz =  -div(flux.theta) 
    440441        dtheta_rhodz(ij,l)=-1./Ai(ij)*(ne(ij,right)*Ftheta(ij+u_right,l)    +  & 
    441442                                 ne(ij,rup)*Ftheta(ij+u_rup,l)        +  &   
    … …  
    458459      DO i=ii_begin,ii_end 
    459460        ij=(j-1)*iim+i 
    460 !signe ? 
     461        ! convm = +div(mass flux), sign convention as in Ringler et al. 2012, eq. 21 
    461462        convm(ij,l)= 1./Ai(ij)*(ne(ij,right)*Fe(ij+u_right,l)   +  & 
    462463                                ne(ij,rup)*Fe(ij+u_rup,l)       +  &   
    … …  
    498499    DO i=ii_begin,ii_end 
    499500      ij=(j-1)*iim+i 
     501      ! dps/dt = -int(div flux)dz 
    500502      dps(ij)=-convm(ij,1) * g  
    501503    ENDDO 
    … …  
    510512      DO i=ii_begin,ii_end 
    511513        ij=(j-1)*iim+i 
     514        ! w = int(z,ztop,div(flux)dz) + B(eta)dps/dt 
     515        ! => w>0 for upward transport 
    512516        w( ij, l+1 ) = convm( ij, l+1 ) - bp(l+1) * convm( ij, 1 ) 
    513517      ENDDO 
    … …  
    682686    DO j=jj_begin,jj_end 
    683687      DO i=ii_begin,ii_end 
     688         ! ww>0 <=> upward transport 
    684689        ij=(j-1)*iim+i 
    685690        ww            =   0.5 * w(ij,l+1) * (theta(ij,l) + theta(ij,l+1) ) 
    686         dtheta_rhodz(ij, l ) = dtheta_rhodz(ij, l )  -  ww  
     691        dtheta_rhodz(ij, l ) = dtheta_rhodz(ij, l )  -  ww 
    687692        dtheta_rhodz(ij,l+1) = dtheta_rhodz(ij,l+1)  +  ww 
    688693      ENDDO 

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

    r19 r22  
    145145 
    146146 
    147     DO it=1,500 
     147    DO it=1,20 
    148148 
    149149      CALL transfert_request(f_u,req_dissip) 
    … …  
    212212 
    213213 
    214     DO it=1,500 
     214    DO it=1,20 
    215215  
    216216      CALL transfert_request(f_u,req_dissip) 
    … …  
    279279    ENDDO 
    280280 
    281     DO it=1,500 
     281    DO it=1,20 
    282282 
    283283      CALL transfert_request(f_theta,req_i1) 
Note: See TracChangeset for help on using the changeset viewer.