Changeset 138 for codes/icosagcm/trunk


Ignore:
Timestamp:
02/16/13 17:03:57 (11 years ago)
Author:
dubos
Message:

Transport now working again - tested with dcmip4.1.0

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

Legend:

Unmodified
Added
Removed

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

    r137 r138  
    33  IMPLICIT NONE  
    44 
     5  PRIVATE 
     6   
     7  PUBLIC :: init_advect, compute_backward_traj, compute_gradq3d, compute_advect_horiz 
     8 
    59CONTAINS 
    610 
    … …  
    812 
    913  SUBROUTINE init_advect(normal,tangent) 
    10     USE domain_mod 
    11     USE dimensions 
    12     USE geometry 
    13     USE metric 
    14     USE vector 
    1514    IMPLICIT NONE 
    1615    REAL(rstd),INTENT(OUT) :: normal(3*iim*jjm,3) 
    … …  
    4847 
    4948  SUBROUTINE compute_gradq3d(qi,gradq3d) 
    50     USE domain_mod 
    51     USE dimensions 
    52     USE geometry 
    53     USE metric 
    5449    IMPLICIT NONE 
    5550    REAL(rstd),INTENT(IN)  :: qi(iim*jjm,llm) 
    … …  
    6257    REAL(rstd) :: ar 
    6358    INTEGER :: i,j,n,k,ind,l 
     59 
     60    ! TODO : precompute ar, drop arr as output argument of gradq ? 
     61 
    6462    !========================================================================================== GRADIENT  
     63    ! Compute gradient at triangles solving a linear system 
     64    ! arr = area of triangle joining centroids of hexagons 
    6565    Do l = 1,llm  
    6666       DO j=jj_begin-1,jj_end+1 
    … …  
    114114  END SUBROUTINE compute_gradq3d 
    115115 
    116   ! Backward trajectories, used in Miura approach 
     116  ! Backward trajectories, for use with Miura approach 
    117117  SUBROUTINE compute_backward_traj(normal,tangent,ue,tau, cc) 
    118     USE domain_mod 
    119     USE dimensions 
    120     USE geometry 
    121     USE metric 
    122118    IMPLICIT NONE 
    123119    REAL(rstd),INTENT(IN)    :: normal(3*iim*jjm,3) 
    … …  
    127123    REAL(rstd),INTENT(IN)    :: tau 
    128124 
    129     REAL(rstd) :: v_e(3), up_e 
    130     REAL(rstd) :: qe 
    131     REAL(rstd) :: ed(3)     
     125    REAL(rstd) :: v_e(3), up_e, qe, ed(3)     
    132126    INTEGER :: i,j,n,l 
    133127 
     128    ! TODO : compute normal displacement ue*tau as hfluxt / mass(upwind) then reconstruct tangential displacement 
     129     
    134130    ! reconstruct tangential wind then 3D wind at edge then cc = edge midpoint - u*tau 
    135131    DO l = 1,llm 
    … …  
    189185  ! Horizontal transport (S. Dubey, T. Dubos) 
    190186  ! Slope-limited van Leer approach with hexagons 
    191   SUBROUTINE compute_advect_horiz(update_mass, normal,tangent,qi,gradq3d,him,ue,hfluxt,bigt) 
    192     USE domain_mod 
    193     USE dimensions 
    194     USE geometry 
    195     USE metric 
    196     IMPLICIT NONE 
    197     LOGICAL, INTENT(IN)      :: update_mass 
    198     REAL(rstd),INTENT(IN)    :: normal(3*iim*jjm,3) 
    199     REAL(rstd),INTENT(IN)    :: tangent(3*iim*jjm,3) 
    200     REAL(rstd),INTENT(INOUT) :: qi(iim*jjm,llm) 
    201     REAL(rstd),INTENT(IN)    :: gradq3d(iim*jjm,llm,3)  
    202     REAL(rstd),INTENT(INOUT) :: him(iim*jjm,llm) 
    203     REAL(rstd),INTENT(IN)    :: ue(iim*3*jjm,llm) 
    204     REAL(rstd),INTENT(IN)    :: hfluxt(3*iim*jjm,llm) ! mass flux 
    205     REAL(rstd),INTENT(IN)    :: bigt 
    206  
    207     REAL(rstd) :: dqi(iim*jjm,llm),dhi(iim*jjm,llm)  
    208     REAL(rstd) :: cc(3*iim*jjm,llm,3)   ! barycenter of quadrilateral, where q is evaluated (1-point quadrature) 
    209     REAL(rstd) :: v_e(3*iim*jjm,llm,3) 
    210     REAL(rstd) :: up_e 
    211  
    212 !    REAL(rstd) :: qe(3*iim*jjm,llm) 
    213     REAL(rstd) :: qe 
    214     REAL(rstd) :: ed(3)     
    215     REAL(rstd) :: flxx(3*iim*jjm,llm) 
     187  SUBROUTINE compute_advect_horiz(update_mass,hfluxt,cc,gradq3d, mass,qi)  
     188    IMPLICIT NONE 
     189    LOGICAL, INTENT(IN)       :: update_mass 
     190    REAL(rstd), INTENT(IN)    :: gradq3d(iim*jjm,llm,3)  
     191    REAL(rstd), INTENT(IN)    :: hfluxt(3*iim*jjm,llm) ! mass flux 
     192    REAL(rstd), INTENT(IN)    :: cc(3*iim*jjm,llm,3)   ! barycenter of quadrilateral, where q is evaluated (1-point quadrature) 
     193    REAL(rstd), INTENT(INOUT) :: mass(iim*jjm,llm) 
     194    REAL(rstd), INTENT(INOUT) :: qi(iim*jjm,llm) 
     195 
     196    REAL(rstd) :: dq,dmass,qe,ed(3), newmass 
     197    REAL(rstd) :: qflux(3*iim*jjm,llm) 
    216198    INTEGER :: i,j,n,k,l 
    217     REAL(rstd):: him_old(llm)  
    218  
    219  
    220     ! reconstruct tangential wind then 3D wind at edge then cc = edge midpoint - u*dt/2 
    221     ! TODO : this does not depend on q and should be done only once, not for each tracer 
     199 
     200    ! evaluate tracer field at point cc using piecewise linear reconstruction 
     201    ! q(cc)= q0 + gradq.(cc-xyz_i)  with xi centroid of hexagon 
     202    ! ne*hfluxt>0 iff outgoing 
    222203    DO l = 1,llm 
    223204       DO j=jj_begin-1,jj_end+1 
    224205          DO i=ii_begin-1,ii_end+1 
    225206             n=(j-1)*iim+i 
    226  
    227              up_e =1/de(n+u_right)*(                                                   & 
    228                   wee(n+u_right,1,1)*le(n+u_rup)*ue(n+u_rup,l)+                        & 
    229                   wee(n+u_right,2,1)*le(n+u_lup)*ue(n+u_lup,l)+                        & 
    230                   wee(n+u_right,3,1)*le(n+u_left)*ue(n+u_left,l)+                      & 
    231                   wee(n+u_right,4,1)*le(n+u_ldown)*ue(n+u_ldown,l)+                    & 
    232                   wee(n+u_right,5,1)*le(n+u_rdown)*ue(n+u_rdown,l)+                    &  
    233                   wee(n+u_right,1,2)*le(n+t_right+u_ldown)*ue(n+t_right+u_ldown,l)+    & 
    234                   wee(n+u_right,2,2)*le(n+t_right+u_rdown)*ue(n+t_right+u_rdown,l)+    & 
    235                   wee(n+u_right,3,2)*le(n+t_right+u_right)*ue(n+t_right+u_right,l)+    & 
    236                   wee(n+u_right,4,2)*le(n+t_right+u_rup)*ue(n+t_right+u_rup,l)+        & 
    237                   wee(n+u_right,5,2)*le(n+t_right+u_lup)*ue(n+t_right+u_lup,l)         & 
    238                   ) 
    239              v_e(n+u_right,l,:)= ue(n+u_right,l)*normal(n+u_right,:) + up_e*tangent(n+u_right,:) 
    240  
    241              up_e=1/de(n+u_lup)*(                                                        & 
    242                   wee(n+u_lup,1,1)*le(n+u_left)*ue(n+u_left,l)+                        & 
    243                   wee(n+u_lup,2,1)*le(n+u_ldown)*ue(n+u_ldown,l)+                      & 
    244                   wee(n+u_lup,3,1)*le(n+u_rdown)*ue(n+u_rdown,l)+                      & 
    245                   wee(n+u_lup,4,1)*le(n+u_right)*ue(n+u_right,l)+                      & 
    246                   wee(n+u_lup,5,1)*le(n+u_rup)*ue(n+u_rup,l)+                          &  
    247                   wee(n+u_lup,1,2)*le(n+t_lup+u_right)*ue(n+t_lup+u_right,l)+          & 
    248                   wee(n+u_lup,2,2)*le(n+t_lup+u_rup)*ue(n+t_lup+u_rup,l)+              & 
    249                   wee(n+u_lup,3,2)*le(n+t_lup+u_lup)*ue(n+t_lup+u_lup,l)+              & 
    250                   wee(n+u_lup,4,2)*le(n+t_lup+u_left)*ue(n+t_lup+u_left,l)+            & 
    251                   wee(n+u_lup,5,2)*le(n+t_lup+u_ldown)*ue(n+t_lup+u_ldown,l)           & 
    252                   ) 
    253              v_e(n+u_lup,l,:)= ue(n+u_lup,l)*normal(n+u_lup,:) + up_e*tangent(n+u_lup,:) 
    254  
    255              up_e=1/de(n+u_ldown)*(                                                    & 
    256                   wee(n+u_ldown,1,1)*le(n+u_rdown)*ue(n+u_rdown,l)+                    & 
    257                   wee(n+u_ldown,2,1)*le(n+u_right)*ue(n+u_right,l)+                    & 
    258                   wee(n+u_ldown,3,1)*le(n+u_rup)*ue(n+u_rup,l)+                        & 
    259                   wee(n+u_ldown,4,1)*le(n+u_lup)*ue(n+u_lup,l)+                        & 
    260                   wee(n+u_ldown,5,1)*le(n+u_left)*ue(n+u_left,l)+                      &  
    261                   wee(n+u_ldown,1,2)*le(n+t_ldown+u_lup)*ue(n+t_ldown+u_lup,l)+        & 
    262                   wee(n+u_ldown,2,2)*le(n+t_ldown+u_left)*ue(n+t_ldown+u_left,l)+      & 
    263                   wee(n+u_ldown,3,2)*le(n+t_ldown+u_ldown)*ue(n+t_ldown+u_ldown,l)+    & 
    264                   wee(n+u_ldown,4,2)*le(n+t_ldown+u_rdown)*ue(n+t_ldown+u_rdown,l)+    & 
    265                   wee(n+u_ldown,5,2)*le(n+t_ldown+u_right)*ue(n+t_ldown+u_right,l)     & 
    266                   ) 
    267              v_e(n+u_ldown,l,:)= ue(n+u_ldown,l)*normal(n+u_ldown,:) + up_e*tangent(n+u_ldown,:) 
    268  
    269              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  
    270              cc(n+u_lup,l,:)   = xyz_e(n+u_lup,:)   - v_e(n+u_lup,l,:)*0.5*bigt  
    271              cc(n+u_ldown,l,:) = xyz_e(n+u_ldown,:) - v_e(n+u_ldown,l,:)*0.5*bigt  
    272           ENDDO 
    273        ENDDO 
    274     END DO 
    275     ! end of tracer-independant computations 
    276  
    277     ! evaluate traver field at point cc usin piecewise linear reconstruction 
    278     ! q(cc)= q0 + gradq.(cc-xyz_i)  with xi centroid of hexagon 
    279     DO l = 1,llm 
    280        DO j=jj_begin-1,jj_end+1 
    281           DO i=ii_begin-1,ii_end+1 
    282              n=(j-1)*iim+i 
    283              if (ne(n,right)*ue(n+u_right,l)>0) then  
     207             if (ne(n,right)*hfluxt(n+u_right,l)>0) then  
    284208                ed = cc(n+u_right,l,:) - xyz_i(n,:) 
    285                 qe=qi(n,l)+sum2(gradq3d(n,l,:),ed)  
     209                qe = qi(n,l)+sum2(gradq3d(n,l,:),ed)  
    286210             else 
    287211                ed = cc(n+u_right,l,:) - xyz_i(n+t_right,:) 
    288                 qe=qi(n+t_right,l)+sum2(gradq3d(n+t_right,l,:),ed)  
     212                qe = qi(n+t_right,l)+sum2(gradq3d(n+t_right,l,:),ed)  
    289213             endif 
    290              flxx(n+u_right,l) = hfluxt(n+u_right,:)*qe*ne(n,right) 
     214             qflux(n+u_right,l) = hfluxt(n+u_right,l)*qe 
    291215              
    292              if (ne(n,lup)*ue(n+u_lup,l)>0) then  
     216             if (ne(n,lup)*hfluxt(n+u_lup,l)>0) then  
    293217                ed = cc(n+u_lup,l,:) - xyz_i(n,:) 
    294                 qe(n+u_lup,l)=qi(n,l)+sum2(gradq3d(n,l,:),ed) 
     218                qe = qi(n,l)+sum2(gradq3d(n,l,:),ed) 
    295219             else 
    296220                ed = cc(n+u_lup,l,:) - xyz_i(n+t_lup,:) 
    297                 qe(n+u_lup,l)=qi(n+t_lup,l)+sum2(gradq3d(n+t_lup,l,:),ed)  
     221                qe = qi(n+t_lup,l)+sum2(gradq3d(n+t_lup,l,:),ed)  
    298222             endif 
    299              flxx(n+u_lup,:) = hfluxt(n+u_lup,:)*qe(n+u_lup,:)*ne(n,lup)  
     223             qflux(n+u_lup,l) = hfluxt(n+u_lup,l)*qe  
    300224              
    301              if (ne(n,ldown)*ue(n+u_ldown,l)>0) then  
     225             if (ne(n,ldown)*hfluxt(n+u_ldown,l)>0) then  
    302226                ed = cc(n+u_ldown,l,:) - xyz_i(n,:) 
    303                 qe(n+u_ldown,l)=qi(n,l)+ sum2(gradq3d(n,l,:),ed)  
     227                qe = qi(n,l)+sum2(gradq3d(n,l,:),ed)  
    304228             else 
    305229                ed = cc(n+u_ldown,l,:) - xyz_i(n+t_ldown,:) 
    306                 qe(n+u_ldown,l)=qi(n+t_ldown,l)+sum2(gradq3d(n+t_ldown,l,:),ed)  
     230                qe = qi(n+t_ldown,l)+sum2(gradq3d(n+t_ldown,l,:),ed)  
    307231             endif 
    308              flxx(n+u_ldown,:) = hfluxt(n+u_ldown,:)*qe(n+u_ldown,:)*ne(n,ldown) 
     232             qflux(n+u_ldown,l) = hfluxt(n+u_ldown,l)*qe 
    309233          end do 
    310234       end do 
    311235    END DO 
    312236 
    313     ! evaluate q flux as mass flux * qe 
    314     ! TODO : loop over k ! 
    315     ! TODO : merge with previous loop and eliminate unnecessary temporary qe 
    316     DO j=jj_begin-1,jj_end+1 
    317        DO i=ii_begin-1,ii_end+1 
    318           n=(j-1)*iim+i  
    319        ENDDO 
    320     ENDDO 
    321  
    322     ! update q and, if update_mass, mass 
    323     ! TODO : loop over k ! 
    324     DO j=jj_begin,jj_end 
    325        DO i=ii_begin,ii_end 
    326           n=(j-1)*iim+i   
    327  
    328           dhi(n,:)= -(1/Ai(n))*(he(n+u_right,:)*ne(n,right) + he(n+u_lup,:)*ne(n,lup) & 
    329                + he(n+u_ldown,:)*ne(n,ldown) + he(n+u_rup,:)*ne(n,rup) & 
    330                + he(n+u_left,:)*ne(n,left) +  he(n+u_rdown,:)*ne(n,rdown) )  
    331  
    332           dqi(n,:)= -(1/Ai(n))*(flxx(n+u_right,:)+flxx(n+u_lup,:)       &  
    333                +flxx(n+u_ldown,:) - flxx(n+u_rup,:)    & 
    334                - flxx(n+u_left,:) - flxx(n+u_rdown,:) )    
    335           him_old(:) = him(n,:)  
    336           him(n,:) = him(n,:) + dhi(n,:)  
    337           qi(n,:) = (qi(n,:)*him_old(:) + dqi(n,:))/him(n,:)   
    338  
     237    ! update q and, if update_mass, update mass 
     238    DO l = 1,llm 
     239       DO j=jj_begin,jj_end 
     240          DO i=ii_begin,ii_end 
     241             n=(j-1)*iim+i   
     242             ! sign convention as Ringler et al. (2010) eq. 21 
     243             dmass =   hfluxt(n+u_right,l)*ne(n,right)   & 
     244                     + hfluxt(n+u_lup,l)  *ne(n,lup)     & 
     245                     + hfluxt(n+u_ldown,l)*ne(n,ldown)   & 
     246                     + hfluxt(n+u_rup,l)  *ne(n,rup)     & 
     247                     + hfluxt(n+u_left,l) *ne(n,left)    & 
     248                     + hfluxt(n+u_rdown,l)*ne(n,rdown) 
     249 
     250             dq    =   qflux(n+u_right,l) *ne(n,right)   & 
     251                     + qflux(n+u_lup,l)   *ne(n,lup)     & 
     252                     + qflux(n+u_ldown,l) *ne(n,ldown)   & 
     253                     + qflux(n+u_rup,l)   *ne(n,rup)     & 
     254                     + qflux(n+u_left,l)  *ne(n,left)    & 
     255                     + qflux(n+u_rdown,l) *ne(n,rdown) 
     256 
     257             newmass = mass(n,l) - dmass/Ai(n) 
     258             qi(n,l) = (qi(n,l)*mass(n,l) - dq/Ai(n) ) / newmass 
     259             IF(update_mass) mass(n,l) = newmass 
     260          END DO 
    339261       END DO 
    340262    END DO 
    341263 
    342264  CONTAINS 
     265 
    343266    !==================================================================================== 
    344     REAL FUNCTION sum2(a1,a2) 
     267    PURE REAL(rstd) FUNCTION sum2(a1,a2) 
    345268      IMPLICIT NONE 
    346       REAL,INTENT(IN):: a1(3), a2(3) 
     269      REAL(rstd),INTENT(IN):: a1(3), a2(3) 
    347270      sum2 = a1(1)*a2(1)+a1(2)*a2(2)+a1(3)*a2(3) 
     271      !      sum2 = 0. ! Godunov scheme 
    348272    END FUNCTION sum2 
    349273 
    350274  END SUBROUTINE compute_advect_horiz 
     275 
    351276  !========================================================================== 
    352   SUBROUTINE gradq(n0,n1,n2,n3,q,dq,det) 
    353     USE geometry 
    354     USE metric 
    355     USE dimensions 
     277  PURE SUBROUTINE gradq(n0,n1,n2,n3,q,dq,det) 
    356278    IMPLICIT NONE 
    357279    INTEGER, INTENT(IN) :: n0,n1,n2,n3 
    358     REAL,INTENT(IN)     :: q(iim*jjm) 
    359     REAL,INTENT(OUT)    :: dq(3) 
    360     REAL(rstd)  ::det,detx,dety,detz 
     280    REAL(rstd), INTENT(IN)     :: q(iim*jjm) 
     281    REAL(rstd), INTENT(OUT)    :: dq(3), det 
     282     
     283    REAL(rstd)  ::detx,dety,detz 
    361284    INTEGER    :: info 
    362285    INTEGER    :: IPIV(3) 
    … …  
    365288    REAL(rstd) :: B(3) 
    366289 
    367     A(1,1)=xyz_i(n1,1) -xyz_i(n0,1); A(1,2)=xyz_i(n1,2)- xyz_i(n0,2); A(1,3)=xyz_i(n1,3)  - xyz_i(n0,3)  
    368     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)  
    369     A(3,1)=xyz_v(n3,1);              A(3,2)= xyz_v(n3,2);             A(3,3)= xyz_v(n3,3) 
    370  
     290    ! TODO : replace A by A1,A2,A3 
     291    A(1,1)=xyz_i(n1,1)-xyz_i(n0,1);  A(1,2)=xyz_i(n1,2)-xyz_i(n0,2); A(1,3)=xyz_i(n1,3)-xyz_i(n0,3)  
     292    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)  
     293    A(3,1)=xyz_v(n3,1);              A(3,2)= xyz_v(n3,2);            A(3,3)=xyz_v(n3,3) 
    371294 
    372295    dq(1) = q(n1)-q(n0) 
    … …  
    382305    dq(3) = detz  
    383306  END SUBROUTINE gradq 
     307 
    384308  !========================================================================== 
    385   SUBROUTINE determinant(a1,a2,a3,det) 
    386     IMPLICIT NONE 
    387     REAL, DIMENSION(3) :: a1, a2,a3 
    388     REAL ::  det,x1,x2,x3 
     309  PURE SUBROUTINE determinant(a1,a2,a3,det) 
     310    IMPLICIT NONE 
     311    REAL(rstd), DIMENSION(3), INTENT(IN) :: a1,a2,a3 
     312    REAL(rstd), INTENT(OUT) :: det 
     313    REAL(rstd) ::  x1,x2,x3 
    389314    x1 =  a1(1) * (a2(2) * a3(3) - a2(3) * a3(2)) 
    390315    x2 =  a1(2) * (a2(1) * a3(3) - a2(3) * a3(1)) 

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

    r137 r138  
    11MODULE advect_tracer_mod 
    22  USE icosa 
     3  IMPLICIT NONE 
    34  PRIVATE 
    45 
    … …  
    89  TYPE(t_field),POINTER :: f_cc(:)  ! starting point of backward-trajectory (Miura approach) 
    910 
    10 !  TYPE(t_field),POINTER :: f_rhodzm1(:) 
    11 !  TYPE(t_field),POINTER :: f_rhodz(:) 
    12  
    1311  REAL(rstd), PARAMETER :: pente_max=2.0 ! for vlz 
    1412 
    … …  
    1917  SUBROUTINE init_advect_tracer 
    2018    USE advect_mod 
    21     IMPLICIT NONE 
    2219    REAL(rstd),POINTER :: tangent(:,:) 
    2320    REAL(rstd),POINTER :: normal(:,:) 
    2421    INTEGER :: ind 
    2522 
    26     CALL allocate_field(f_normal,field_u,type_real,3) 
    27     CALL allocate_field(f_tangent,field_u,type_real,3) 
    28     CALL allocate_field(f_gradq3d,field_t,type_real,llm,3) 
    29     CALL allocate_field(f_cc,field_u,type_real,llm,3) 
    30  
    31 !    CALL allocate_field(f_rhodz,field_t,type_real,llm)  
    32 !    CALL allocate_field(f_rhodzm1,field_t,type_real,llm)  
     23    CALL allocate_field(f_normal,field_u,type_real,3, name='normal') 
     24    CALL allocate_field(f_tangent,field_u,type_real,3, name='tangent') 
     25    CALL allocate_field(f_gradq3d,field_t,type_real,llm,3, name='gradq3d') 
     26    CALL allocate_field(f_cc,field_u,type_real,llm,3, name='cc') 
    3327 
    3428    DO ind=1,ndomain 
    … …  
    4337 
    4438  SUBROUTINE advect_tracer(f_hfluxt, f_wfluxt,f_u, f_q,f_rhodz) 
    45     USE icosa 
    46     USE field_mod 
    47     USE domain_mod 
    48     USE dimensions 
    49     USE grid_param 
    50     USE geometry 
    51     USE metric 
    5239    USE advect_mod 
    53     USE advect_mod 
    54     USE disvert_mod 
    5540    USE mpipara 
    5641    IMPLICIT NONE 
    … …  
    6247    TYPE(t_field),POINTER :: f_rhodz(:)    ! mass field at beginning of macro time step 
    6348 
    64     REAL(rstd),POINTER :: q(:,:,:), normal(:,:,:), tangent(:,:,:), gradq3d(:,:,:) 
     49    REAL(rstd),POINTER :: q(:,:,:), normal(:,:), tangent(:,:), gradq3d(:,:,:), cc(:,:,:) 
    6550    REAL(rstd),POINTER :: hfluxt(:,:), wfluxt(:,:) 
    6651    REAL(rstd),POINTER :: rhodz(:,:), u(:,:)  
    67     INTEGER :: ind 
     52    INTEGER :: ind,k 
    6853 
    6954    CALL transfert_request(f_u,req_e1) 
    70     CALL transfert_request(f_u,req_e1) 
     55!    CALL transfert_request(f_hfluxt,req_e1)  ! BUG : This (unnecessary) transfer makes the computation go wrong 
     56    CALL transfert_request(f_wfluxt,req_i1) 
    7157    CALL transfert_request(f_q,req_i1) 
     58    CALL transfert_request(f_rhodz,req_i1) 
    7259         
    7360    IF (is_mpi_root) PRINT *, 'Advection scheme' 
    7461 
     62!    DO ind=1,ndomain 
     63!       CALL swap_dimensions(ind) 
     64!       CALL swap_geometry(ind) 
     65!       normal  = f_normal(ind) 
     66!       tangent = f_tangent(ind) 
     67!       cc      = f_cc(ind) 
     68!       u       = f_u(ind) 
     69!       q       = f_q(ind) 
     70!       rhodz   = f_rhodz(ind) 
     71!       hfluxt  = f_hfluxt(ind)  
     72!       wfluxt  = f_wfluxt(ind)  
     73!       gradq3d = f_gradq3d(ind) 
     74! 
     75!       ! 1/2 vertical transport 
     76!       DO k = 1, nqtot 
     77!          CALL vlz(k==nqtot,0.5, wfluxt,rhodz,q(:,:,k)) 
     78!       END DO 
     79! 
     80!       ! horizontal transport 
     81!       CALL compute_backward_traj(tangent,normal,u,0.5*dt*itau_adv, cc)  
     82!       DO k = 1,nqtot 
     83!          CALL compute_gradq3d(q(:,:,k),gradq3d) 
     84!          CALL compute_advect_horiz(k==nqtot,hfluxt,cc,gradq3d, rhodz,q(:,:,k)) 
     85!       END DO 
     86! 
     87!       ! 1/2 vertical transport 
     88!       DO k = 1,nqtot 
     89!          CALL vlz(k==nqtot, 0.5,wfluxt,rhodz, q(:,:,k)) 
     90!       END DO 
     91!    END DO 
     92 
     93    ! 1/2 vertical transport + back-trajectories 
     94    DO ind=1,ndomain 
     95       CALL swap_dimensions(ind) 
     96       CALL swap_geometry(ind) 
     97       normal  = f_normal(ind) 
     98       tangent = f_tangent(ind) 
     99       cc      = f_cc(ind) 
     100       u       = f_u(ind) 
     101       q       = f_q(ind) 
     102       rhodz   = f_rhodz(ind) 
     103       wfluxt  = f_wfluxt(ind)  
     104       DO k = 1, nqtot 
     105          CALL vlz(k==nqtot,0.5, wfluxt,rhodz,q(:,:,k)) 
     106       END DO 
     107       CALL compute_backward_traj(tangent,normal,u,0.5*dt*itau_adv, cc)  
     108    END DO 
     109 
     110    CALL transfert_request(f_q,req_i1)      ! necessary ? 
     111    CALL transfert_request(f_rhodz,req_i1)  ! necessary ? 
     112 
     113    ! horizontal transport - split in two to place transfer of gradq3d 
     114 
     115    DO k = 1, nqtot 
     116       DO ind=1,ndomain 
     117          CALL swap_dimensions(ind) 
     118          CALL swap_geometry(ind) 
     119          q       = f_q(ind) 
     120          gradq3d = f_gradq3d(ind) 
     121          CALL compute_gradq3d(q(:,:,k),gradq3d) 
     122       END DO 
     123 
     124       CALL transfert_request(f_gradq3d,req_i1) 
     125 
     126       DO ind=1,ndomain 
     127          CALL swap_dimensions(ind) 
     128          CALL swap_geometry(ind) 
     129          cc      = f_cc(ind) 
     130          q       = f_q(ind) 
     131          rhodz   = f_rhodz(ind) 
     132          hfluxt  = f_hfluxt(ind)  
     133          gradq3d = f_gradq3d(ind) 
     134          CALL compute_advect_horiz(k==nqtot,hfluxt,cc,gradq3d, rhodz,q(:,:,k)) 
     135       END DO 
     136    END DO  
     137 
     138    CALL transfert_request(f_q,req_i1)      ! necessary ? 
     139    CALL transfert_request(f_rhodz,req_i1)  ! necessary ? 
     140 
     141    ! 1/2 vertical transport 
    75142    DO ind=1,ndomain 
    76143       CALL swap_dimensions(ind) 
    … …  
    78145       q       = f_q(ind) 
    79146       rhodz   = f_rhodz(ind) 
    80        hfluxt  = f_hfluxt(ind)  
    81147       wfluxt  = f_wfluxt(ind)  
    82        normal  = f_normal(ind) 
    83        tangent = f_tangent(ind) 
    84        gradq3d = f_gradq3d(ind) 
    85        CALL compute_adv_tracer(tangent,normal,hfluxt,wfluxt,u, q,rhodz,gradq3d) 
     148       DO k = 1,nqtot 
     149          CALL vlz(k==nqtot, 0.5,wfluxt,rhodz, q(:,:,k)) 
     150       END DO 
    86151    END DO 
    87152 
    88153  END SUBROUTINE advect_tracer 
    89154 
    90   SUBROUTINE compute_adv_tracer(tangent,normal,hfluxt,wfluxt,u, q,rhodz,gradq3d) 
    91     USE icosa 
    92     USE field_mod 
    93     USE domain_mod 
    94     USE dimensions 
    95     USE grid_param 
    96     USE geometry 
    97     USE metric 
    98     USE advect_mod 
    99     USE advect_mod 
    100     USE disvert_mod 
    101     USE mpipara 
    102     REAL(rstd), INTENT(IN) :: tangent(:,:,:), normal(:,:,:) 
    103     REAL(rstd), INTENT(IN) :: hfluxt(:,:), wfluxt(:,:), u(:,:) 
    104     REAL(rstd), INTENT(INOUT) :: rhodz(:,:) 
    105     REAL(rstd), INTENT(INOUT) :: q(:,:,:) 
    106     REAL(rstd), INTENT(OUT)   :: gradq3d(:,:,:) 
    107     INTEGER :: k 
    108     ! for mass/tracer consistency each transport step also updates the mass field rhodz 
    109     ! mass is updated after all tracers have been updated, when k==nqtot 
    110  
    111     ! 1/2 vertical transport 
    112     DO k = 1, nqtot 
    113        CALL vlz(k==nqtot, 0.5,wfluxt,rhodz,q(:,:,k)) 
    114     END DO 
    115  
    116     ! horizontal transport 
    117     DO k = 1,nqtot 
    118        CALL compute_gradq3d(q(:,:,k),gradq3d) 
    119        CALL compute_advect_horiz(k==nqtot, tangent,normal,q(:,:,k),gradq3d,rhodz,u,hfluxt,dt*itau_adv)  
    120     END DO 
    121  
    122     ! 1/2 vertical transport 
    123     DO k = 1,nqtot 
    124        CALL vlz(k==nqtot, 0.5,wfluxt,rhodz, q(:,:,k)) 
    125     END DO 
    126   END SUBROUTINE compute_adv_tracer 
    127    
    128155  SUBROUTINE vlz(update_mass, fac,wfluxt,mass, q) 
    129156    ! 
    … …  
    135162    !     wfluxt >0 for upward transport 
    136163    !    ******************************************************************** 
    137     USE icosa 
    138164    IMPLICIT NONE 
    139165    LOGICAL, INTENT(IN)       :: update_mass 
    … …  
    152178    ! finite difference of q 
    153179    DO l=2,llm 
    154        DO j=jj_begin,jj_end 
    155           DO i=ii_begin,ii_end 
     180       DO j=jj_begin-1,jj_end+1 
     181          DO i=ii_begin-1,ii_end+1 
    156182             ij=(j-1)*iim+i 
    157183             dzqw(ij,l)=q(ij,l)-q(ij,l-1) 
    … …  
    164190    ! dzq = slope*dz, i.e. the reconstructed q varies by dzq inside level l 
    165191    DO l=2,llm-1 
    166        DO j=jj_begin,jj_end 
    167           DO i=ii_begin,ii_end 
     192       DO j=jj_begin-1,jj_end+1 
     193          DO i=ii_begin-1,ii_end+1 
    168194             ij=(j-1)*iim+i 
    169195             IF(dzqw(ij,l)*dzqw(ij,l+1).gt.0.) THEN 
    … …  
    179205 
    180206    ! 0 slope in top and bottom layers 
    181     DO j=jj_begin,jj_end 
    182        DO i=ii_begin,ii_end 
     207    DO j=jj_begin-1,jj_end+1 
     208       DO i=ii_begin-1,ii_end+1 
    183209          ij=(j-1)*iim+i 
    184210          dzq(ij,1)=0. 
    … …  
    190216    ! then amount of q leaving level l/l+1 = wq = w * qq 
    191217    DO l = 1,llm-1 
    192        DO j=jj_begin,jj_end 
    193           DO i=ii_begin,ii_end 
    194              ij=(j-1)*iim+i 
    195              IF(wfluxt(ij,l+1).gt.0.) THEN 
    196                 ! downward transport, sigw<0 
    197                 ! qq = q if sigw=-1 , qq = q-dzq/2 if sigw=0                 
    198                 w          = fac*wfluxt(ij,l+1) 
     218       DO j=jj_begin-1,jj_end+1 
     219          DO i=ii_begin-1,ii_end+1 
     220             ij=(j-1)*iim+i 
     221             w = fac*wfluxt(ij,l+1) 
     222             IF(w>0.) THEN  ! upward transport, upwind side is at level l  
     223                sigw       = w/mass(ij,l) 
     224                qq         = q(ij,l)+0.5*(1.-sigw)*dzq(ij,l) ! qq = q if sigw=1 , qq = q+dzq/2 if sigw=0 
     225             ELSE           ! downward transport, upwind side is at level l+1 
    199226                sigw       = w/mass(ij,l+1) 
    200                 qq         = q(ij,l+1) - 0.5*(1.+sigw)*dzq(ij,l+1) 
    201                 wq(ij,l+1) = w*qq 
    202              ELSE 
    203                 ! upward transport, sigw>0 
    204                 ! qq = q if sigw=1 , qq = q+dzq/2 if sigw=0 
    205                 w          = fac*wfluxt(ij,l) 
    206                 sigw       = w/mass(ij,l) 
    207                 qq         = q(ij,l)+0.5*(1.-sigw)*dzq(ij,l) 
    208                 wq(ij,l+1) = w*qq 
     227                qq         = q(ij,l+1)-0.5*(1.+sigw)*dzq(ij,l+1) ! qq = q if sigw=-1 , qq = q-dzq/2 if sigw=0                 
    209228             ENDIF 
     229             wq(ij,l+1) = w*qq 
    210230          ENDDO 
    211231       ENDDO 
    212232    END DO 
    213233    ! wq = 0 at top and bottom 
    214     DO j=jj_begin,jj_end 
    215        DO i=ii_begin,ii_end 
     234    DO j=jj_begin-1,jj_end+1 
     235       DO i=ii_begin-1,ii_end+1 
    216236          ij=(j-1)*iim+i 
    217237          wq(ij,llm+1)=0. 
    … …  
    222242    ! update q, mass is updated only after all q's have been updated 
    223243    DO l=1,llm 
    224        DO j=jj_begin,jj_end 
    225           DO i=ii_begin,ii_end 
     244       DO j=jj_begin-1,jj_end+1 
     245          DO i=ii_begin-1,ii_end+1 
    226246             ij=(j-1)*iim+i 
    227247             newmass = mass(ij,l) + fac*(wfluxt(ij,l)-wfluxt(ij,l+1)) 

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

    r134 r138  
    297297    REAL(rstd),INTENT(IN)  :: qu(iim*3*jjm,llm) 
    298298 
    299     REAL(rstd),INTENT(OUT) :: du(iim*3*jjm,llm), hflux(iim*3*jjm,llm) 
     299    REAL(rstd),INTENT(OUT) :: du(iim*3*jjm,llm), hflux(iim*3*jjm,llm) ! hflux in kg/s 
    300300    REAL(rstd),INTENT(OUT) :: dtheta_rhodz(iim*jjm,llm) 
    301301    REAL(rstd),INTENT(OUT) :: dps(iim*jjm) 
    302     REAL(rstd),INTENT(OUT) :: wflux(iim*jjm,llm+1) ! vertical mass flux 
     302    REAL(rstd),INTENT(OUT) :: wflux(iim*jjm,llm+1) ! vertical mass flux (kg/m2/s) 
    303303     
    304304    INTEGER :: i,j,ij,l 
    … …  
    312312    REAL(rstd),ALLOCATABLE,SAVE :: phi(:,:)    ! geopotential 
    313313    REAL(rstd),ALLOCATABLE,SAVE :: Ftheta(:,:) ! theta flux 
    314     REAL(rstd),ALLOCATABLE,SAVE :: convm(:,:)  ! mass flux convergence 
     314    REAL(rstd),ALLOCATABLE,SAVE :: divm(:,:)  ! mass flux divergence 
    315315    REAL(rstd),ALLOCATABLE,SAVE :: berni(:,:)  ! Bernouilli function 
    316316     
    … …  
    328328    ALLOCATE(phi(iim*jjm,llm))      ! geopotential 
    329329    ALLOCATE(Ftheta(3*iim*jjm,llm)) ! theta flux    
    330     ALLOCATE(convm(iim*jjm,llm))    ! mass flux convergence 
     330    ALLOCATE(divm(iim*jjm,llm))    ! mass flux convergence 
    331331    ALLOCATE(berni(iim*jjm,llm))    ! bernouilli term 
    332332 
    … …  
    363363      DO i=ii_begin,ii_end 
    364364        ij=(j-1)*iim+i   
    365         ! convm = +div(mass flux), sign convention as in Ringler et al. 2012, eq. 21 
    366         convm(ij,l)= 1./Ai(ij)*(ne(ij,right)*hflux(ij+u_right,l)   +  & 
     365        ! divm = +div(mass flux), sign convention as in Ringler et al. 2012, eq. 21 
     366        divm(ij,l)= 1./Ai(ij)*(ne(ij,right)*hflux(ij+u_right,l)   +  & 
    367367                                ne(ij,rup)*hflux(ij+u_rup,l)       +  &   
    368368                                ne(ij,lup)*hflux(ij+u_lup,l)       +  &   
    … …  
    383383  ENDDO 
    384384    
    385 !!! cumulate mass flux convergence from top to bottom 
     385!!! cumulate mass flux divergence from top to bottom 
    386386  DO  l = llm-1, 1, -1 
    387387!$OMP DO 
    … …  
    389389      DO i=ii_begin,ii_end 
    390390        ij=(j-1)*iim+i 
    391         convm(ij,l) = convm(ij,l) + convm(ij,l+1) 
     391        divm(ij,l) = divm(ij,l) + divm(ij,l+1) 
    392392      ENDDO 
    393393    ENDDO 
    … …  
    402402        ! w = int(z,ztop,div(flux)dz) + B(eta)dps/dt 
    403403        ! => w>0 for upward transport 
    404         wflux( ij, l+1 ) = convm( ij, l+1 ) - bp(l+1) * convm( ij, 1 ) 
     404        wflux( ij, l+1 ) = divm( ij, l+1 ) - bp(l+1) * divm( ij, 1 ) 
    405405      ENDDO 
    406406    ENDDO 
    407407  ENDDO 
    408408 
    409 ! compute dps, vertical mass flux at the surface = 0 
     409! compute dps, set vertical mass flux at the surface to 0 
    410410!$OMP DO 
    411411  DO j=jj_begin,jj_end 
    … …  
    414414      wflux(ij,1)  = 0. 
    415415      ! dps/dt = -int(div flux)dz 
    416       dps(ij)=-convm(ij,1) * g  
     416      dps(ij)=-divm(ij,1) * g  
    417417    ENDDO 
    418418  ENDDO 
    … …  
    644644    DEALLOCATE(phi)     ! geopotential 
    645645    DEALLOCATE(Ftheta)  ! theta flux    
    646     DEALLOCATE(convm)   ! mass flux convergence 
     646    DEALLOCATE(divm)    ! mass flux divergence 
    647647    DEALLOCATE(berni)   ! bernouilli term 
    648648!!$OMP END MASTER 

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

    r80 r138  
    174174     
    175175      q(:,:,1)=theta(:,:) 
     176      IF(nqtot>2) q(:,:,3)=1. 
    176177 
    177178      DO l=1,llm 
    … …  
    200201                                                            -cos(latc)*sin(lat)*cos(lon-lonc))/(sqrt(1-K**2))) 
    201202             q(ij,l,2)=ABS(g/preff*(-1./radius*duodeta*dthetaodlat-(2*sin(lat)*omega+vort)*dthetaodeta)) 
     203             IF(nqtot>3) q(ij,l,4)=cos(lon)*cos(lat) 
    202204           ENDDO 
    203205         ENDDO 

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

    r29 r138  
    1111     
    1212  TYPE t_field 
     13    CHARACTER(30)      :: name 
    1314    REAL(rstd),POINTER :: rval2d(:) 
    1415    REAL(rstd),POINTER :: rval3d(:,:) 
    … …  
    4546CONTAINS 
    4647 
    47   SUBROUTINE allocate_field(field,field_type,data_type,dim1,dim2) 
     48  SUBROUTINE allocate_field(field,field_type,data_type,dim1,dim2,name) 
    4849  USE domain_mod 
    4950  IMPLICIT NONE 
    … …  
    5253    INTEGER,INTENT(IN) :: data_type 
    5354    INTEGER,OPTIONAL :: dim1,dim2 
     55    CHARACTER(*), OPTIONAL :: name 
    5456    INTEGER :: ind 
    5557    INTEGER :: ii_size,jj_size 
    5658 
    57     ALLOCATE(field(ndomain))     
     59    ALLOCATE(field(ndomain)) 
    5860 
    5961    DO ind=1,ndomain 
    60    
     62      IF(PRESENT(name)) THEN 
     63         field(ind)%name = name 
     64      ELSE 
     65         field(ind)%name = '(unkown)' 
     66      END IF 
     67 
    6168      IF (PRESENT(dim2)) THEN 
    6269        field(ind)%ndim=4  
    … …  
    224231    TYPE(t_field),INTENT(IN) :: field 
    225232     
    226     IF (field%ndim/=2 .OR. field%data_type/=type_real) STOP 'get_val_r2d : bad pointer assignation'         
     233    IF (field%ndim/=2 .OR. field%data_type/=type_real) THEN 
     234       PRINT *, 'get_val_r2d : bad pointer assignation with ' // TRIM(field%name)  
     235       STOP 
     236    END IF 
    227237    field_pt=>field%rval2d 
    228238  END SUBROUTINE  getval_r2d 
    … …  
    233243    TYPE(t_field),INTENT(IN) :: field 
    234244     
    235     IF (field%ndim/=3 .OR. field%data_type/=type_real) STOP 'get_val_r3d : bad pointer assignation'         
     245    IF (field%ndim/=3 .OR. field%data_type/=type_real) THEN 
     246       PRINT *, 'get_val_r3d : bad pointer assignation with ' // TRIM(field%name)  
     247       STOP 
     248    END IF 
    236249    field_pt=>field%rval3d 
    237250  END SUBROUTINE  getval_r3d 
    … …  
    242255    TYPE(t_field),INTENT(IN) :: field 
    243256     
    244     IF (field%ndim/=4 .OR. field%data_type/=type_real) STOP 'get_val_r4d : bad pointer assignation'         
     257    IF (field%ndim/=4 .OR. field%data_type/=type_real) THEN 
     258       PRINT *, 'get_val_r4d : bad pointer assignation with ' // TRIM(field%name) 
     259       STOP 
     260    END IF 
    245261    field_pt=>field%rval4d 
    246   END SUBROUTINE  getval_r4d 
    247    
     262  END SUBROUTINE  getval_r4d   
    248263 
    249264   

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

    r135 r138  
    4646  INTEGER :: it,i,j,n, nb_stage, stage, matsuno_period, scheme 
    4747  CHARACTER(len=255) :: scheme_name 
    48   LOGICAL :: fluxt_zero ! set to .TRUE. to start accumulating fluxes in time 
     48  LOGICAL :: fluxt_zero(ndomain) ! set to .TRUE. to start accumulating fluxes in time 
    4949!  INTEGER :: itaumax 
    5050!  REAL(rstd) ::write_period 
    … …  
    141141     CALL swap_geometry(ind) 
    142142     rhodz=f_rhodz(ind); ps=f_ps(ind) 
    143      CALL compute_rhodz(ps,rhodz) ! save rhodz for transport scheme before dynamics update ps 
     143     CALL compute_rhodz(.TRUE., ps,rhodz) ! save rhodz for transport scheme before dynamics update ps 
    144144  END DO 
    145   fluxt_zero=.FALSE. 
    146  
     145  fluxt_zero=.TRUE. 
     146 
     147  ! check that rhodz is consistent with ps 
     148  CALL transfert_request(f_rhodz,req_i1) 
     149  CALL transfert_request(f_ps,req_i1) 
     150  DO ind=1,ndomain 
     151     CALL swap_dimensions(ind) 
     152     CALL swap_geometry(ind) 
     153     rhodz=f_rhodz(ind); ps=f_ps(ind) 
     154     CALL compute_rhodz(.FALSE., ps, rhodz)    
     155  END DO 
     156   
    147157  DO it=0,itaumax 
    148158 
    … …  
    182192 
    183193    IF(MOD(it+1,itau_adv)==0) THEN 
     194       CALL transfert_request(f_wfluxt,req_i1) ! FIXME 
     195!       CALL transfert_request(f_hfluxt,req_e1) ! FIXME 
     196 
    184197       CALL advect_tracer(f_hfluxt,f_wfluxt,f_u, f_q,f_rhodz)  ! update q and rhodz after RK step 
    185198       fluxt_zero=.TRUE. 
     199 
     200       ! FIXME : check that rhodz is consistent with ps 
     201       CALL transfert_request(f_rhodz,req_i1) 
     202       CALL transfert_request(f_ps,req_i1) 
     203       CALL transfert_request(f_dps,req_i1) ! FIXME 
     204       CALL transfert_request(f_wflux,req_i1) ! FIXME 
     205       DO ind=1,ndomain 
     206          CALL swap_dimensions(ind) 
     207          CALL swap_geometry(ind) 
     208          rhodz=f_rhodz(ind); ps=f_ps(ind); dps=f_dps(ind);  
     209          wflux=f_wflux(ind); wfluxt=f_wfluxt(ind) 
     210          CALL compute_rhodz(.FALSE., ps, rhodz)    
     211       END DO 
     212 
    186213    END IF 
    187214 
    … …  
    195222    LOGICAL :: with_dps 
    196223    INTEGER :: ind 
    197        
    198224    DO ind=1,ndomain 
     225       CALL swap_dimensions(ind) 
     226       CALL swap_geometry(ind) 
    199227       IF(with_dps) THEN 
    200228          ps=f_ps(ind) ; dps=f_dps(ind) ;  
    201229          ps(:)=ps(:)+dt*dps(:) 
    202230          hflux=f_hflux(ind);     hfluxt=f_hfluxt(ind) 
    203           wflux=f_hflux(ind);     wfluxt=f_wfluxt(ind) 
    204           CALL accumulate_fluxes(hflux,wflux,hfluxt,wfluxt,dt,fluxt_zero) 
     231          wflux=f_wflux(ind);     wfluxt=f_wfluxt(ind) 
     232          CALL accumulate_fluxes(hflux,wflux,hfluxt,wfluxt,dt,fluxt_zero(ind)) 
    205233       END IF 
    206234       u=f_u(ind) ; theta_rhodz=f_theta_rhodz(ind) 
    … …  
    221249       
    222250      DO ind=1,ndomain 
     251         CALL swap_dimensions(ind) 
     252         CALL swap_geometry(ind) 
    223253         ps=f_ps(ind)   ; u=f_u(ind)   ; theta_rhodz=f_theta_rhodz(ind) 
    224254         psm1=f_psm1(ind) ; um1=f_um1(ind) ; theta_rhodzm1=f_theta_rhodzm1(ind) 
    … …  
    234264         IF(stage==nb_stage) THEN ! accumulate mass fluxes at last stage 
    235265            hflux=f_hflux(ind);     hfluxt=f_hfluxt(ind) 
    236             wflux=f_hflux(ind);     wfluxt=f_wfluxt(ind) 
    237             CALL accumulate_fluxes(hflux,wflux,hfluxt,wfluxt,dt,fluxt_zero) 
     266            wflux=f_wflux(ind);     wfluxt=f_wfluxt(ind) 
     267            CALL accumulate_fluxes(hflux,wflux, hfluxt,wfluxt, dt,fluxt_zero(ind)) 
    238268         END IF 
    239269      END DO 
    … …  
    245275 
    246276      DO ind=1,ndomain 
     277        CALL swap_dimensions(ind) 
     278        CALL swap_geometry(ind) 
    247279        ps=f_ps(ind)   ; u=f_u(ind)   ; theta_rhodz=f_theta_rhodz(ind) 
    248280        psm1=f_psm1(ind) ; um1=f_um1(ind) ; theta_rhodzm1=f_theta_rhodzm1(ind) 
    … …  
    277309    tau = dt/nb_stage 
    278310      DO ind=1,ndomain 
     311        CALL swap_dimensions(ind) 
     312        CALL swap_geometry(ind) 
     313 
    279314        ps=f_ps(ind)   ; u=f_u(ind)   ; theta_rhodz=f_theta_rhodz(ind) 
    280315        psm1=f_psm1(ind) ; um1=f_um1(ind) ; theta_rhodzm1=f_theta_rhodzm1(ind) 
    … …  
    320355 
    321356      DO ind=1,ndomain 
     357        CALL swap_dimensions(ind) 
     358        CALL swap_geometry(ind) 
    322359        ps=f_ps(ind)   ; u=f_u(ind)   ; theta_rhodz=f_theta_rhodz(ind) 
    323360        dps=f_dps(ind) ; du=f_du(ind) ; dtheta_rhodz=f_dtheta_rhodz(ind) 
    … …  
    343380  END SUBROUTINE timeloop     
    344381 
    345   SUBROUTINE compute_rhodz(ps, rhodz) 
     382  SUBROUTINE compute_rhodz(comp, ps, rhodz) 
    346383    USE icosa 
    347384    USE disvert_mod 
     385    LOGICAL, INTENT(IN) :: comp ! .TRUE. to compute, .FALSE. to check 
    348386    REAL(rstd), INTENT(IN) :: ps(iim*jjm) 
    349387    REAL(rstd), INTENT(OUT) :: rhodz(iim*jjm,llm) 
    350     INTEGER :: l,i,j,ij 
     388    REAL(rstd) :: m, err 
     389    INTEGER :: l,i,j,ij,dd 
     390    err=0. 
     391    IF(comp) THEN  
     392       dd=1 
     393    ELSE 
     394!       dd=-1 
     395       dd=0 
     396    END IF 
     397 
    351398    DO l = 1, llm 
    352        DO j=jj_begin-1,jj_end+1 
    353           DO i=ii_begin-1,ii_end+1 
     399       DO j=jj_begin-dd,jj_end+dd 
     400          DO i=ii_begin-dd,ii_end+dd 
    354401             ij=(j-1)*iim+i 
    355              rhodz(ij,l) = (ap(l) - ap(l+1) + (bp(l)-bp(l+1))*ps(ij))/g  
     402             m = ( ap(l) - ap(l+1) + (bp(l)-bp(l+1))*ps(ij) )/g  
     403             IF(comp) THEN 
     404                rhodz(ij,l) = m 
     405             ELSE 
     406                err = MAX(err,abs(m-rhodz(ij,l))) 
     407             END IF 
    356408          ENDDO 
    357409       ENDDO 
    358410    ENDDO 
     411 
     412    IF(.NOT. comp) THEN 
     413       IF(err>1e-10) THEN 
     414          PRINT *, 'Discrepancy between ps and rhodz detected', err 
     415          STOP 
     416       ELSE 
     417!          PRINT *, 'No discrepancy between ps and rhodz detected' 
     418       END IF 
     419    END IF 
     420 
    359421  END SUBROUTINE compute_rhodz 
    360422 
    361   SUBROUTINE accumulate_fluxes(hflux,wflux,hfluxt,wfluxt,tau,fluxt_zero) 
     423  SUBROUTINE accumulate_fluxes(hflux,wflux, hfluxt,wfluxt, tau,fluxt_zero) 
    362424    USE icosa 
    363425    REAL(rstd), INTENT(IN)    :: hflux(3*iim*jjm,llm), wflux(iim*jjm,llm+1) 
    … …  
    366428    LOGICAL, INTENT(INOUT) :: fluxt_zero 
    367429    IF(fluxt_zero) THEN 
     430!       PRINT *, 'Accumulating fluxes (first)' 
    368431       fluxt_zero=.FALSE. 
    369432       hfluxt = tau*hflux 
    370433       wfluxt = tau*wflux 
    371434    ELSE 
     435!       PRINT *, 'Accumulating fluxes (next)' 
    372436       hfluxt = hfluxt + tau*hflux 
    373437       wfluxt = wfluxt + tau*wflux 
Note: See TracChangeset for help on using the changeset viewer.