Changeset 137


Ignore:
Timestamp:
02/15/13 13:59:06 (11 years ago)
Author:
dubos
Message:

Heavily modified advect_tracer.f90

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

Legend:

Unmodified
Added
Removed

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

    r136 r137  
    114114  END SUBROUTINE compute_gradq3d 
    115115 
    116   !===================================================================================================    
    117   SUBROUTINE compute_advect_horiz(update_mass, normal,tangent,qi,gradq3d,him,ue,he,bigt) 
     116  ! Backward trajectories, used in Miura approach 
     117  SUBROUTINE compute_backward_traj(normal,tangent,ue,tau, cc) 
    118118    USE domain_mod 
    119119    USE dimensions 
    … …  
    121121    USE metric 
    122122    IMPLICIT NONE 
    123     LOGICAL, INTENT(IN)      :: update_mass 
    124123    REAL(rstd),INTENT(IN)    :: normal(3*iim*jjm,3) 
    125124    REAL(rstd),INTENT(IN)    :: tangent(3*iim*jjm,3) 
    126     REAL(rstd),INTENT(INOUT) :: qi(iim*jjm,llm) 
    127     REAL(rstd),INTENT(IN)    :: gradq3d(iim*jjm,llm,3)  
    128     REAL(rstd),INTENT(INOUT) :: him(iim*jjm,llm) 
    129125    REAL(rstd),INTENT(IN)    :: ue(iim*3*jjm,llm) 
    130     REAL(rstd),INTENT(IN)    :: he(3*iim*jjm,llm) ! mass flux 
    131     REAL(rstd),INTENT(IN)    :: bigt 
    132  
    133     REAL(rstd) :: dqi(iim*jjm,llm),dhi(iim*jjm,llm)  
    134     REAL(rstd) :: cc(3*iim*jjm,llm,3)   ! barycenter of quadrilateral, where q is evaluated (1-point quadrature) 
    135     REAL(rstd) :: v_e(3*iim*jjm,llm,3) 
    136     REAL(rstd) :: up_e 
    137  
    138     REAL(rstd) :: qe(3*iim*jjm,llm) 
     126    REAL(rstd),INTENT(OUT)   :: cc(3*iim*jjm,llm,3) ! start of backward trajectory 
     127    REAL(rstd),INTENT(IN)    :: tau 
     128 
     129    REAL(rstd) :: v_e(3), up_e 
     130    REAL(rstd) :: qe 
    139131    REAL(rstd) :: ed(3)     
    140     REAL(rstd) :: flxx(3*iim*jjm,llm) 
    141     INTEGER :: i,j,n,k,l 
    142     REAL(rstd):: him_old(llm)  
    143  
    144  
    145     ! reconstruct tangential wind then 3D wind at edge then cc = edge midpoint - u*dt/2 
    146     ! TODO : this does not depend on q and should be done only once, not for each tracer 
     132    INTEGER :: i,j,n,l 
     133 
     134    ! reconstruct tangential wind then 3D wind at edge then cc = edge midpoint - u*tau 
    147135    DO l = 1,llm 
    148136       DO j=jj_begin-1,jj_end+1 
    … …  
    162150                  wee(n+u_right,5,2)*le(n+t_right+u_lup)*ue(n+t_right+u_lup,l)         & 
    163151                  ) 
     152             v_e = ue(n+u_right,l)*normal(n+u_right,:) + up_e*tangent(n+u_right,:) 
     153             cc(n+u_right,l,:) = xyz_e(n+u_right,:) - v_e*tau 
     154 
     155             up_e=1/de(n+u_lup)*(                                                      & 
     156                  wee(n+u_lup,1,1)*le(n+u_left)*ue(n+u_left,l)+                        & 
     157                  wee(n+u_lup,2,1)*le(n+u_ldown)*ue(n+u_ldown,l)+                      & 
     158                  wee(n+u_lup,3,1)*le(n+u_rdown)*ue(n+u_rdown,l)+                      & 
     159                  wee(n+u_lup,4,1)*le(n+u_right)*ue(n+u_right,l)+                      & 
     160                  wee(n+u_lup,5,1)*le(n+u_rup)*ue(n+u_rup,l)+                          &  
     161                  wee(n+u_lup,1,2)*le(n+t_lup+u_right)*ue(n+t_lup+u_right,l)+          & 
     162                  wee(n+u_lup,2,2)*le(n+t_lup+u_rup)*ue(n+t_lup+u_rup,l)+              & 
     163                  wee(n+u_lup,3,2)*le(n+t_lup+u_lup)*ue(n+t_lup+u_lup,l)+              & 
     164                  wee(n+u_lup,4,2)*le(n+t_lup+u_left)*ue(n+t_lup+u_left,l)+            & 
     165                  wee(n+u_lup,5,2)*le(n+t_lup+u_ldown)*ue(n+t_lup+u_ldown,l)           & 
     166                  ) 
     167             v_e = ue(n+u_lup,l)*normal(n+u_lup,:) + up_e*tangent(n+u_lup,:) 
     168             cc(n+u_lup,l,:) = xyz_e(n+u_lup,:) - v_e*tau 
     169 
     170             up_e=1/de(n+u_ldown)*(                                                    & 
     171                  wee(n+u_ldown,1,1)*le(n+u_rdown)*ue(n+u_rdown,l)+                    & 
     172                  wee(n+u_ldown,2,1)*le(n+u_right)*ue(n+u_right,l)+                    & 
     173                  wee(n+u_ldown,3,1)*le(n+u_rup)*ue(n+u_rup,l)+                        & 
     174                  wee(n+u_ldown,4,1)*le(n+u_lup)*ue(n+u_lup,l)+                        & 
     175                  wee(n+u_ldown,5,1)*le(n+u_left)*ue(n+u_left,l)+                      &  
     176                  wee(n+u_ldown,1,2)*le(n+t_ldown+u_lup)*ue(n+t_ldown+u_lup,l)+        & 
     177                  wee(n+u_ldown,2,2)*le(n+t_ldown+u_left)*ue(n+t_ldown+u_left,l)+      & 
     178                  wee(n+u_ldown,3,2)*le(n+t_ldown+u_ldown)*ue(n+t_ldown+u_ldown,l)+    & 
     179                  wee(n+u_ldown,4,2)*le(n+t_ldown+u_rdown)*ue(n+t_ldown+u_rdown,l)+    & 
     180                  wee(n+u_ldown,5,2)*le(n+t_ldown+u_right)*ue(n+t_ldown+u_right,l)     & 
     181                  ) 
     182             v_e = ue(n+u_ldown,l)*normal(n+u_ldown,:) + up_e*tangent(n+u_ldown,:) 
     183             cc(n+u_ldown,l,:) = xyz_e(n+u_ldown,:) - v_e*tau 
     184          ENDDO 
     185       ENDDO 
     186    END DO 
     187  END SUBROUTINE compute_backward_traj 
     188 
     189  ! Horizontal transport (S. Dubey, T. Dubos) 
     190  ! 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) 
     216    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 
     222    DO l = 1,llm 
     223       DO j=jj_begin-1,jj_end+1 
     224          DO i=ii_begin-1,ii_end+1 
     225             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                  ) 
    164239             v_e(n+u_right,l,:)= ue(n+u_right,l)*normal(n+u_right,:) + up_e*tangent(n+u_right,:) 
    165240 
    … …  
    208283             if (ne(n,right)*ue(n+u_right,l)>0) then  
    209284                ed = cc(n+u_right,l,:) - xyz_i(n,:) 
    210                 qe(n+u_right,l)=qi(n,l)+sum2(gradq3d(n,l,:),ed)  
     285                qe=qi(n,l)+sum2(gradq3d(n,l,:),ed)  
    211286             else 
    212287                ed = cc(n+u_right,l,:) - xyz_i(n+t_right,:) 
    213                 qe(n+u_right,l)=qi(n+t_right,l)+sum2(gradq3d(n+t_right,l,:),ed)  
     288                qe=qi(n+t_right,l)+sum2(gradq3d(n+t_right,l,:),ed)  
    214289             endif 
     290             flxx(n+u_right,l) = hfluxt(n+u_right,:)*qe*ne(n,right) 
     291              
    215292             if (ne(n,lup)*ue(n+u_lup,l)>0) then  
    216293                ed = cc(n+u_lup,l,:) - xyz_i(n,:) 
    … …  
    220297                qe(n+u_lup,l)=qi(n+t_lup,l)+sum2(gradq3d(n+t_lup,l,:),ed)  
    221298             endif 
     299             flxx(n+u_lup,:) = hfluxt(n+u_lup,:)*qe(n+u_lup,:)*ne(n,lup)  
     300              
    222301             if (ne(n,ldown)*ue(n+u_ldown,l)>0) then  
    223302                ed = cc(n+u_ldown,l,:) - xyz_i(n,:) 
    … …  
    227306                qe(n+u_ldown,l)=qi(n+t_ldown,l)+sum2(gradq3d(n+t_ldown,l,:),ed)  
    228307             endif 
     308             flxx(n+u_ldown,:) = hfluxt(n+u_ldown,:)*qe(n+u_ldown,:)*ne(n,ldown) 
    229309          end do 
    230310       end do 
    … …  
    237317       DO i=ii_begin-1,ii_end+1 
    238318          n=(j-1)*iim+i  
    239           flxx(n+u_right,:) = he(n+u_right,:)*qe(n+u_right,:)*ne(n,right) 
    240           flxx(n+u_lup,:) = he(n+u_lup,:)*qe(n+u_lup,:)*ne(n,lup)  
    241           flxx(n+u_ldown,:) = he(n+u_ldown,:)*qe(n+u_ldown,:)*ne(n,ldown) 
    242319       ENDDO 
    243320    ENDDO 

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

    r136 r137  
    66  TYPE(t_field),POINTER :: f_tangent(:) 
    77  TYPE(t_field),POINTER :: f_gradq3d(:) 
    8  
    9   TYPE(t_field),POINTER :: f_wg(:) 
    10   TYPE(t_field),POINTER :: f_zm(:) 
    11   TYPE(t_field),POINTER :: f_zq(:) 
    12  
    13   TYPE(t_field),POINTER :: f_massflxc(:) 
    14   TYPE(t_field),POINTER :: f_massflx(:) 
    15   TYPE(t_field),POINTER :: f_uc(:) 
    16   TYPE(t_field),POINTER :: f_rhodzm1(:) 
    17   TYPE(t_field),POINTER :: f_rhodz(:) 
     8  TYPE(t_field),POINTER :: f_cc(:)  ! starting point of backward-trajectory (Miura approach) 
     9 
     10!  TYPE(t_field),POINTER :: f_rhodzm1(:) 
     11!  TYPE(t_field),POINTER :: f_rhodz(:) 
    1812 
    1913  REAL(rstd), PARAMETER :: pente_max=2.0 ! for vlz 
    … …  
    3327    CALL allocate_field(f_tangent,field_u,type_real,3) 
    3428    CALL allocate_field(f_gradq3d,field_t,type_real,llm,3) 
    35  
    36     CALL allocate_field(f_wg,field_t,type_real,llm) 
    37     CALL allocate_field(f_zm,field_t,type_real,llm) 
    38     CALL allocate_field(f_zq,field_t,type_real,llm,nqtot)  
    39  
    40     CALL allocate_field(f_rhodz,field_t,type_real,llm)  
    41     CALL allocate_field(f_rhodzm1,field_t,type_real,llm)  
    42     CALL allocate_field(f_massflxc,field_u,type_real,llm)  
    43     CALL allocate_field(f_massflx,field_u,type_real,llm)  
    44     CALL allocate_field(f_uc,field_u,type_real,llm) 
     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)  
    4533 
    4634    DO ind=1,ndomain 
Note: See TracChangeset for help on using the changeset viewer.