Changeset 136


Ignore:
Timestamp:
02/15/13 13:42:05 (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

    r69 r136  
    115115 
    116116  !===================================================================================================    
    117   SUBROUTINE compute_advect_horiz(normal,tangent,qi,gradq3d,him,ue,he,bigt) 
     117  SUBROUTINE compute_advect_horiz(update_mass, normal,tangent,qi,gradq3d,him,ue,he,bigt) 
    118118    USE domain_mod 
    119119    USE dimensions 
    … …  
    121121    USE metric 
    122122    IMPLICIT NONE 
     123    LOGICAL, INTENT(IN)      :: update_mass 
    123124    REAL(rstd),INTENT(IN)    :: normal(3*iim*jjm,3) 
    124125    REAL(rstd),INTENT(IN)    :: tangent(3*iim*jjm,3) 
    … …  
    131132 
    132133    REAL(rstd) :: dqi(iim*jjm,llm),dhi(iim*jjm,llm)  
    133     REAL(rstd) :: cc(3*iim*jjm,llm,3)  
     134    REAL(rstd) :: cc(3*iim*jjm,llm,3)   ! barycenter of quadrilateral, where q is evaluated (1-point quadrature) 
    134135    REAL(rstd) :: v_e(3*iim*jjm,llm,3) 
    135136    REAL(rstd) :: up_e 
    … …  
    142143 
    143144 
    144     !go to  123 
     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 
    145147    DO l = 1,llm 
    146148       DO j=jj_begin-1,jj_end+1 
    … …  
    196198       ENDDO 
    197199    END DO 
    198     !123         continue            
    199     !========================================================================================================== 
     200    ! end of tracer-independant computations 
     201 
     202    ! evaluate traver field at point cc usin piecewise linear reconstruction 
     203    ! q(cc)= q0 + gradq.(cc-xyz_i)  with xi centroid of hexagon 
    200204    DO l = 1,llm 
    201205       DO j=jj_begin-1,jj_end+1 
    … …  
    227231    END DO 
    228232 
    229  
     233    ! evaluate q flux as mass flux * qe 
     234    ! TODO : loop over k ! 
     235    ! TODO : merge with previous loop and eliminate unnecessary temporary qe 
    230236    DO j=jj_begin-1,jj_end+1 
    231237       DO i=ii_begin-1,ii_end+1 
    … …  
    237243    ENDDO 
    238244 
     245    ! update q and, if update_mass, mass 
     246    ! TODO : loop over k ! 
    239247    DO j=jj_begin,jj_end 
    240248       DO i=ii_begin,ii_end 

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

    r133 r136  
    22  USE icosa 
    33  PRIVATE 
    4   INTEGER,PARAMETER::iapp_tracvl= 1 
    54 
    65  TYPE(t_field),POINTER :: f_normal(:) 
    … …  
    87  TYPE(t_field),POINTER :: f_gradq3d(:) 
    98 
    10   PUBLIC init_advect_tracer, advect_tracer_rhodz, advect_tracer 
     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(:) 
     18 
     19  REAL(rstd), PARAMETER :: pente_max=2.0 ! for vlz 
     20 
     21  PUBLIC init_advect_tracer, advect_tracer 
    1122 
    1223CONTAINS 
    … …  
    2233    CALL allocate_field(f_tangent,field_u,type_real,3) 
    2334    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) 
    2445 
    2546    DO ind=1,ndomain 
    … …  
    3354  END SUBROUTINE init_advect_tracer 
    3455 
    35   SUBROUTINE advect_tracer(f_ps,f_u,f_q) 
     56  SUBROUTINE advect_tracer(f_hfluxt, f_wfluxt,f_u, f_q,f_rhodz) 
    3657    USE icosa 
    37     USE advect_mod 
    38     USE disvert_mod 
    39     USE mpipara 
    40     IMPLICIT NONE 
    41     TYPE(t_field),POINTER :: f_ps(:) 
    42     TYPE(t_field),POINTER :: f_u(:) 
    43     TYPE(t_field),POINTER :: f_q(:) 
    44     REAL(rstd),POINTER :: q(:,:,:) 
    45     REAL(rstd),POINTER :: u(:,:)  
    46     REAL(rstd),POINTER :: ps(:)  
    47     REAL(rstd),POINTER :: massflx(:,:) 
    48     REAL(rstd),POINTER :: rhodz(:,:) 
    49     TYPE(t_field),POINTER,SAVE :: f_massflxc(:) 
    50     TYPE(t_field),POINTER,SAVE :: f_massflx(:) 
    51     TYPE(t_field),POINTER,SAVE :: f_uc(:) 
    52     TYPE(t_field),POINTER,SAVE :: f_rhodzm1(:) 
    53     TYPE(t_field),POINTER,SAVE :: f_rhodz(:) 
    54     REAL(rstd),POINTER,SAVE :: massflxc(:,:) 
    55     REAL(rstd),POINTER,SAVE :: uc(:,:) 
    56     REAL(rstd),POINTER,SAVE :: rhodzm1(:,:) 
    57     REAL(rstd):: bigt 
    58     INTEGER :: ind,it,i,j,l,ij 
    59     INTEGER,SAVE :: iadvtr=0 
    60     LOGICAL,SAVE:: first=.TRUE.  
    61 !------------------------------------------------------sarvesh 
    62        CALL transfert_request(f_ps,req_i1) 
    63        CALL transfert_request(f_u,req_e1) 
    64        CALL transfert_request(f_u,req_e1) 
    65        CALL transfert_request(f_q,req_i1) 
    66          
    67     IF ( first ) THEN  
    68        CALL allocate_field(f_rhodz,field_t,type_real,llm)  
    69        CALL allocate_field(f_rhodzm1,field_t,type_real,llm)  
    70        CALL allocate_field(f_massflxc,field_u,type_real,llm)  
    71        CALL allocate_field(f_massflx,field_u,type_real,llm)  
    72        CALL allocate_field(f_uc,field_u,type_real,llm)  
    73        first = .FALSE. 
    74     END IF 
    75  
    76     DO ind=1,ndomain 
    77        CALL swap_dimensions(ind) 
    78        CALL swap_geometry(ind) 
    79        rhodz=f_rhodz(ind) 
    80        massflx=f_massflx(ind) 
    81        ps=f_ps(ind) 
    82        u=f_u(ind) 
    83  
    84        DO l = 1, llm 
    85           DO j=jj_begin-1,jj_end+1 
    86              DO i=ii_begin-1,ii_end+1 
    87                 ij=(j-1)*iim+i 
    88                 rhodz(ij,l) = (ap(l) - ap(l+1) + (bp(l)-bp(l+1))*ps(ij))/g  
    89              ENDDO 
    90           ENDDO 
    91        ENDDO 
    92  
    93        DO l = 1, llm 
    94           DO j=jj_begin-1,jj_end+1 
    95              DO i=ii_begin-1,ii_end+1 
    96                 ij=(j-1)*iim+i 
    97                 massflx(ij+u_right,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_right,l))*u(ij+u_right,l)*le(ij+u_right) 
    98                 massflx(ij+u_lup,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_lup,l))*u(ij+u_lup,l)*le(ij+u_lup) 
    99                 massflx(ij+u_ldown,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_ldown,l))*u(ij+u_ldown,l)*le(ij+u_ldown) 
    100              ENDDO 
    101           ENDDO 
    102        ENDDO 
    103     ENDDO 
    104  
    105     IF ( iadvtr == 0 ) THEN  
    106        DO ind=1,ndomain          
    107           CALL swap_dimensions(ind) 
    108           CALL swap_geometry(ind) 
    109           rhodz=f_rhodz(ind)   
    110           rhodzm1 = f_rhodzm1(ind)  
    111           massflxc = f_massflxc(ind)  ! accumulated mass fluxes 
    112           uc = f_uc(ind)              ! time-integrated normal velocity to compute back-trajectory (Miura) 
    113           rhodzm1 = rhodz 
    114           massflxc = 0.0  
    115           uc = 0.0  
    116        END DO 
    117        CALL transfert_request(f_rhodzm1,req_i1)   !----> 
    118        CALL transfert_request(f_massflxc,req_e1) !----> 
    119        CALL transfert_request(f_massflxc,req_e1) !------> 
    120        CALL transfert_request(f_uc,req_e1) !----> 
    121        CALL transfert_request(f_uc,req_e1)  
    122     END IF 
    123  
    124     iadvtr = iadvtr + 1  
    125  
    126     DO ind=1,ndomain 
    127        CALL swap_dimensions(ind) 
    128        CALL swap_geometry(ind) 
    129        massflx  = f_massflx(ind) 
    130        massflxc = f_massflxc(ind)  
    131        uc = f_uc(ind)  
    132        u  = f_u(ind)  
    133        massflxc = massflxc + massflx  
    134        uc = uc + u  
    135     END DO 
    136  
    137     IF ( iadvtr == iapp_tracvl ) THEN  
    138        IF (is_mpi_root) PRINT *, 'Advection scheme' 
    139        bigt = dt*iapp_tracvl 
    140        DO ind=1,ndomain 
    141           CALL swap_dimensions(ind) 
    142           CALL swap_geometry(ind) 
    143           uc = f_uc(ind)  
    144           uc = uc/real(iapp_tracvl)  
    145           massflxc = f_massflxc(ind) 
    146           massflxc = massflxc*dt 
    147           ! now massflx is time-integrated 
    148        END DO 
    149  
    150        CALL vlsplt(f_q,f_rhodzm1,f_massflxc,2.0,f_uc,bigt)  
    151        iadvtr = 0 
    152     END IF 
    153   END SUBROUTINE advect_tracer 
    154  
    155   SUBROUTINE vlsplt(f_q,f_rhodz,f_massflx,pente_max,f_u,bigt) 
    15658    USE field_mod 
    15759    USE domain_mod 
    … …  
    16163    USE metric 
    16264    USE advect_mod 
     65    USE advect_mod 
     66    USE disvert_mod 
     67    USE mpipara 
    16368    IMPLICIT NONE 
    16469 
    165     TYPE(t_field),POINTER :: f_q(:) 
    166     TYPE(t_field),POINTER :: f_u(:) 
    167     TYPE(t_field),POINTER :: f_rhodz(:)  
    168     TYPE(t_field),POINTER :: f_massflx(:) 
    169  
    170     TYPE(t_field),POINTER,SAVE :: f_wg(:) 
    171     TYPE(t_field),POINTER,SAVE :: f_zm(:) 
    172     TYPE(t_field),POINTER,SAVE :: f_zq(:) 
    173  
    174     REAL(rstd)::bigt 
    175     REAL(rstd),POINTER :: q(:,:,:) 
    176     REAL(rstd),POINTER :: u(:,:) 
    177     REAL(rstd),POINTER :: rhodz(:,:)  
    178     REAL(rstd),POINTER :: massflx(:,:)  
    179     REAL(rstd),POINTER :: wg(:,:) 
    180     REAL(rstd),POINTER :: zq(:,:,:)  
    181     REAL(rstd),POINTER :: zm(:,:)  
    182     REAL(rstd),POINTER :: tangent(:,:) 
    183     REAL(rstd),POINTER :: normal(:,:) 
    184     REAL(rstd),POINTER :: gradq3d(:,:,:) 
    185  
    186     REAL(rstd)::pente_max 
    187     LOGICAL,SAVE::first = .true.  
    188     INTEGER :: i,ij,l,j,ind,k 
    189     REAL(rstd) :: zzpbar, zzw  
    190     REAL::qvmax,qvmin  
    191  
    192     IF ( first ) THEN  
    193        CALL allocate_field(f_wg,field_t,type_real,llm) 
    194        CALL allocate_field(f_zm,field_t,type_real,llm) 
    195        CALL allocate_field(f_zq,field_t,type_real,llm,nqtot)  
    196        first = .FALSE. 
    197     END IF 
     70    TYPE(t_field),POINTER :: f_hfluxt(:)   ! time-integrated horizontal mass flux 
     71    TYPE(t_field),POINTER :: f_wfluxt(:)   ! time-integrated vertical mass flux 
     72    TYPE(t_field),POINTER :: f_u(:)        ! velocity (for back-trajectories) 
     73    TYPE(t_field),POINTER :: f_q(:)        ! tracer 
     74    TYPE(t_field),POINTER :: f_rhodz(:)    ! mass field at beginning of macro time step 
     75 
     76    REAL(rstd),POINTER :: q(:,:,:), normal(:,:,:), tangent(:,:,:), gradq3d(:,:,:) 
     77    REAL(rstd),POINTER :: hfluxt(:,:), wfluxt(:,:) 
     78    REAL(rstd),POINTER :: rhodz(:,:), u(:,:)  
     79    INTEGER :: ind 
     80 
     81    CALL transfert_request(f_u,req_e1) 
     82    CALL transfert_request(f_u,req_e1) 
     83    CALL transfert_request(f_q,req_i1) 
     84         
     85    IF (is_mpi_root) PRINT *, 'Advection scheme' 
    19886 
    19987    DO ind=1,ndomain 
    20088       CALL swap_dimensions(ind) 
    20189       CALL swap_geometry(ind) 
    202        q=f_q(ind)  
    203        rhodz=f_rhodz(ind)  
    204        zq=f_zq(ind) 
    205        zm=f_zm(ind)  
    206        zm = rhodz ; zq = q   
    207        wg = f_wg(ind) 
    208        wg = 0.0 
    209        massflx=f_massflx(ind)  
    210        CALL advtrac(massflx,wg) ! compute vertical mass fluxes 
    211     END DO 
    212  
    213     DO ind=1,ndomain 
    214        CALL swap_dimensions(ind) 
    215        CALL swap_geometry(ind) 
    216        zq=f_zq(ind) 
    217        zm=f_zm(ind)  
    218        wg=f_wg(ind) 
    219        wg=wg*0.5 
    220        DO k = 1, nqtot 
    221           CALL vlz(zq(:,:,k),2.,zm,wg) 
    222        END DO 
    223     END DO 
    224  
    225     DO ind=1,ndomain 
    226        CALL swap_dimensions(ind) 
    227        CALL swap_geometry(ind) 
    228        zq=f_zq(ind) 
    229        zq = f_zq(ind) 
    230        zm = f_zm(ind) 
    231        massflx =f_massflx(ind) 
    232        u = f_u(ind) 
     90       q       = f_q(ind) 
     91       rhodz   = f_rhodz(ind) 
     92       hfluxt  = f_hfluxt(ind)  
     93       wfluxt  = f_wfluxt(ind)  
     94       normal  = f_normal(ind) 
    23395       tangent = f_tangent(ind) 
    234        normal = f_normal(ind) 
    23596       gradq3d = f_gradq3d(ind) 
    236  
    237        DO k = 1,nqtot 
    238           CALL compute_gradq3d(zq(:,:,k),gradq3d) 
    239           CALL compute_advect_horiz(tangent,normal,zq(:,:,k),gradq3d,zm,u,massflx,bigt)  
    240        END DO 
    241     END DO 
    242  
    243     DO ind=1,ndomain 
    244        CALL swap_dimensions(ind) 
    245        CALL swap_geometry(ind) 
    246        q = f_q(ind) 
    247        zq = f_zq(ind)  
    248        zm = f_zm(ind)  
    249        wg = f_wg(ind) 
    250        DO k = 1,nqtot 
    251           CALL vlz(zq(:,:,k),2.,zm,wg) 
    252        END DO 
    253        q = zq  
    254     END DO 
    255  
    256   END SUBROUTINE vlsplt 
    257  
    258   !==============================================================================   
    259   SUBROUTINE advtrac(massflx,wgg) 
     97       CALL compute_adv_tracer(tangent,normal,hfluxt,wfluxt,u, q,rhodz,gradq3d) 
     98    END DO 
     99 
     100  END SUBROUTINE advect_tracer 
     101 
     102  SUBROUTINE compute_adv_tracer(tangent,normal,hfluxt,wfluxt,u, q,rhodz,gradq3d) 
     103    USE icosa 
     104    USE field_mod 
    260105    USE domain_mod 
    261106    USE dimensions 
    … …  
    263108    USE geometry 
    264109    USE metric 
     110    USE advect_mod 
     111    USE advect_mod 
    265112    USE disvert_mod 
     113    USE mpipara 
     114    REAL(rstd), INTENT(IN) :: tangent(:,:,:), normal(:,:,:) 
     115    REAL(rstd), INTENT(IN) :: hfluxt(:,:), wfluxt(:,:), u(:,:) 
     116    REAL(rstd), INTENT(INOUT) :: rhodz(:,:) 
     117    REAL(rstd), INTENT(INOUT) :: q(:,:,:) 
     118    REAL(rstd), INTENT(OUT)   :: gradq3d(:,:,:) 
     119    INTEGER :: k 
     120    ! for mass/tracer consistency each transport step also updates the mass field rhodz 
     121    ! mass is updated after all tracers have been updated, when k==nqtot 
     122 
     123    ! 1/2 vertical transport 
     124    DO k = 1, nqtot 
     125       CALL vlz(k==nqtot, 0.5,wfluxt,rhodz,q(:,:,k)) 
     126    END DO 
     127 
     128    ! horizontal transport 
     129    DO k = 1,nqtot 
     130       CALL compute_gradq3d(q(:,:,k),gradq3d) 
     131       CALL compute_advect_horiz(k==nqtot, tangent,normal,q(:,:,k),gradq3d,rhodz,u,hfluxt,dt*itau_adv)  
     132    END DO 
     133 
     134    ! 1/2 vertical transport 
     135    DO k = 1,nqtot 
     136       CALL vlz(k==nqtot, 0.5,wfluxt,rhodz, q(:,:,k)) 
     137    END DO 
     138  END SUBROUTINE compute_adv_tracer 
     139   
     140  SUBROUTINE vlz(update_mass, fac,wfluxt,mass, q) 
     141    ! 
     142    !     Auteurs:   P.Le Van, F.Hourdin, F.Forget, T. Dubos 
     143    ! 
     144    !    ******************************************************************** 
     145    !     Update tracers using vertical mass flux only 
     146    !     Van Leer scheme with minmod limiter 
     147    !     wfluxt >0 for upward transport 
     148    !    ******************************************************************** 
     149    USE icosa 
    266150    IMPLICIT NONE 
    267     REAL(rstd),INTENT(IN) :: massflx(iim*3*jjm,llm) 
    268     REAL(rstd),INTENT(OUT) :: wgg(iim*jjm,llm) 
    269  
    270     INTEGER :: i,j,ij,l 
    271     REAL(rstd) :: convm(iim*jjm,llm)  
    272  
    273     DO l = 1, llm 
    274        DO j=jj_begin,jj_end 
    275           DO i=ii_begin,ii_end 
    276              ij=(j-1)*iim+i 
    277              ! divergence of horizontal flux 
    278              convm(ij,l)= 1/(Ai(ij))*(ne(ij,right)*massflx(ij+u_right,l)   +  & 
    279                   ne(ij,rup)*massflx(ij+u_rup,l)       +  & 
    280                   ne(ij,lup)*massflx(ij+u_lup,l)       +  & 
    281                   ne(ij,left)*massflx(ij+u_left,l)     +  & 
    282                   ne(ij,ldown)*massflx(ij+u_ldown,l)   +  & 
    283                   ne(ij,rdown)*massflx(ij+u_rdown,l)) 
     151    LOGICAL, INTENT(IN)       :: update_mass 
     152    REAL(rstd), INTENT(IN)    :: fac, wfluxt(iim*jjm,llm+1) ! vertical mass flux 
     153    REAL(rstd), INTENT(INOUT) :: mass(iim*jjm,llm) 
     154    REAL(rstd), INTENT(INOUT) :: q(iim*jjm,llm) 
     155 
     156    REAL(rstd) :: dq(iim*jjm,llm), & ! increase of q 
     157         dzqw(iim*jjm,llm),        & ! vertical finite difference of q  
     158         adzqw(iim*jjm,llm),       & ! abs(dzqw) 
     159         dzq(iim*jjm,llm),         & ! limited slope of q 
     160         wq(iim*jjm,llm+1)           ! time-integrated flux of q 
     161    REAL(rstd) :: dzqmax, newmass, sigw, qq, w 
     162    INTEGER :: i,ij,l,j 
     163 
     164    ! finite difference of q 
     165    DO l=2,llm 
     166       DO j=jj_begin,jj_end 
     167          DO i=ii_begin,ii_end 
     168             ij=(j-1)*iim+i 
     169             dzqw(ij,l)=q(ij,l)-q(ij,l-1) 
     170             adzqw(ij,l)=abs(dzqw(ij,l)) 
    284171          ENDDO 
    285172       ENDDO 
    286173    ENDDO 
    287174 
    288     ! accumulate divergence from top of model 
    289     DO  l = llm-1, 1, -1 
    290        DO j=jj_begin,jj_end 
    291           DO i=ii_begin,ii_end 
    292              ij=(j-1)*iim+i 
    293              convm(ij,l) = convm(ij,l) + convm(ij,l+1) 
     175    ! minmod-limited slope of q 
     176    ! dzq = slope*dz, i.e. the reconstructed q varies by dzq inside level l 
     177    DO l=2,llm-1 
     178       DO j=jj_begin,jj_end 
     179          DO i=ii_begin,ii_end 
     180             ij=(j-1)*iim+i 
     181             IF(dzqw(ij,l)*dzqw(ij,l+1).gt.0.) THEN 
     182                dzq(ij,l) = 0.5*( dzqw(ij,l)+dzqw(ij,l+1) ) 
     183                dzqmax    = pente_max * min( adzqw(ij,l),adzqw(ij,l+1) ) 
     184                dzq(ij,l) = sign( min(abs(dzq(ij,l)),dzqmax) , dzq(ij,l) )  ! NB : sign(a,b)=a*sign(b) 
     185             ELSE 
     186                dzq(ij,l)=0. 
     187             ENDIF 
    294188          ENDDO 
    295189       ENDDO 
    296190    ENDDO 
    297191 
    298 !!! Compute vertical velocity 
    299     DO l = 1,llm-1 
    300        DO j=jj_begin,jj_end 
    301           DO i=ii_begin,ii_end 
    302              ij=(j-1)*iim+i 
    303              wgg( ij, l+1 ) = (convm( ij, l+1 ) - bp(l+1) * convm( ij, 1 ))  
    304           ENDDO 
    305        ENDDO 
    306     ENDDO 
    307  
     192    ! 0 slope in top and bottom layers 
    308193    DO j=jj_begin,jj_end 
    309194       DO i=ii_begin,ii_end 
    310195          ij=(j-1)*iim+i 
    311           wgg(ij,1)  = 0. 
     196          dzq(ij,1)=0. 
     197          dzq(ij,llm)=0. 
    312198       ENDDO 
    313199    ENDDO 
    314   END SUBROUTINE advtrac 
    315  
    316   SUBROUTINE vlz(q,pente_max,masse,wgg) 
    317     !c 
    318     !c     Auteurs:   P.Le Van, F.Hourdin, F.Forget  
    319     !c 
    320     !c    ******************************************************************** 
    321     !c     Shema  d'advection " pseudo amont " . 
    322     !c    ******************************************************************** 
    323     USE icosa 
    324     IMPLICIT NONE 
    325     !c 
    326     !c   Arguments: 
    327     !c   ---------- 
    328     REAL masse(iim*jjm,llm),pente_max 
    329     REAL q(iim*jjm,llm) 
    330     REAL wgg(iim*jjm,llm),w(iim*jjm,llm+1)  
    331     REAL dq(iim*jjm,llm) 
    332     INTEGER i,ij,l,j,ii 
    333     !c 
    334     REAL wq(iim*jjm,llm+1),newmasse 
    335  
    336     REAL dzq(iim*jjm,llm),dzqw(iim*jjm,llm),adzqw(iim*jjm,llm),dzqmax 
    337     REAL sigw 
    338  
    339     REAL      SSUM 
    340  
    341  
    342     w(:,1:llm) = -wgg(:,:)  ! w>0 for downward transport 
    343     w(:,llm+1) = 0.0  
    344  
    345     !c    On oriente tout dans le sens de la pression c'est a dire dans le 
    346     !c    sens de W 
    347  
    348     DO l=2,llm 
    349        DO j=jj_begin,jj_end 
    350           DO i=ii_begin,ii_end 
    351              ij=(j-1)*iim+i 
    352              dzqw(ij,l)=q(ij,l-1)-q(ij,l) 
    353              adzqw(ij,l)=abs(dzqw(ij,l)) 
    354           ENDDO 
    355        ENDDO 
    356     ENDDO 
    357  
    358     DO l=2,llm-1 
    359        DO j=jj_begin,jj_end 
    360           DO i=ii_begin,ii_end 
    361              ij=(j-1)*iim+i 
    362              IF(dzqw(ij,l)*dzqw(ij,l+1).gt.0.) THEN 
    363                 dzq(ij,l)=0.5*(dzqw(ij,l)+dzqw(ij,l+1)) 
     200 
     201    ! sigw = fraction of mass that leaves level l/l+1 
     202    ! then amount of q leaving level l/l+1 = wq = w * qq 
     203    DO l = 1,llm-1 
     204       DO j=jj_begin,jj_end 
     205          DO i=ii_begin,ii_end 
     206             ij=(j-1)*iim+i 
     207             IF(wfluxt(ij,l+1).gt.0.) THEN 
     208                ! downward transport, sigw<0 
     209                ! qq = q if sigw=-1 , qq = q-dzq/2 if sigw=0                 
     210                w          = fac*wfluxt(ij,l+1) 
     211                sigw       = w/mass(ij,l+1) 
     212                qq         = q(ij,l+1) - 0.5*(1.+sigw)*dzq(ij,l+1) 
     213                wq(ij,l+1) = w*qq 
    364214             ELSE 
    365                 dzq(ij,l)=0. 
     215                ! upward transport, sigw>0 
     216                ! qq = q if sigw=1 , qq = q+dzq/2 if sigw=0 
     217                w          = fac*wfluxt(ij,l) 
     218                sigw       = w/mass(ij,l) 
     219                qq         = q(ij,l)+0.5*(1.-sigw)*dzq(ij,l) 
     220                wq(ij,l+1) = w*qq 
    366221             ENDIF 
    367              dzqmax=pente_max*min(adzqw(ij,l),adzqw(ij,l+1)) 
    368              dzq(ij,l)=sign(min(abs(dzq(ij,l)),dzqmax),dzq(ij,l)) 
    369           ENDDO 
    370        ENDDO 
    371     ENDDO 
    372  
    373     DO l=2,llm-1 
    374        DO j=jj_begin,jj_end 
    375           DO i=ii_begin,ii_end 
    376              ij=(j-1)*iim+i 
    377              dzq(ij,1)=0. 
    378              dzq(ij,llm)=0. 
    379           ENDDO 
    380        ENDDO 
    381     ENDDO 
    382     !c --------------------------------------------------------------- 
    383     !c   .... calcul des termes d'advection verticale  ....... 
    384     !c --------------------------------------------------------------- 
    385  
    386     !c calcul de  - d( q   * w )/ d(sigma)    qu'on ajoute a  dq pour calculer dq 
    387  
    388     DO l = 1,llm-1 
    389        DO j=jj_begin,jj_end 
    390           DO i=ii_begin,ii_end 
    391              ij=(j-1)*iim+i 
    392              IF(w(ij,l+1).gt.0.) THEN 
    393                 ! upwind only if downward transport 
    394                 sigw=w(ij,l+1)/masse(ij,l+1) 
    395                 wq(ij,l+1)=w(ij,l+1)*(q(ij,l+1)+0.5*(1.-sigw)*dzq(ij,l+1)) 
    396              ELSE 
    397                 ! upwind only if upward transport 
    398                 sigw=w(ij,l+1)/masse(ij,l) 
    399                 wq(ij,l+1)=w(ij,l+1)*(q(ij,l)-0.5*(1.+sigw)*dzq(ij,l)) 
    400              ENDIF 
    401           ENDDO 
    402        ENDDO 
    403     END DO 
    404  
     222          ENDDO 
     223       ENDDO 
     224    END DO 
     225    ! wq = 0 at top and bottom 
    405226    DO j=jj_begin,jj_end 
    406227       DO i=ii_begin,ii_end 
    … …  
    411232    END DO 
    412233 
     234    ! update q, mass is updated only after all q's have been updated 
    413235    DO l=1,llm 
    414236       DO j=jj_begin,jj_end 
    415237          DO i=ii_begin,ii_end 
    416238             ij=(j-1)*iim+i 
    417              ! masse -= dw/dz but w>0 <=> downward 
    418              newmasse=masse(ij,l)+(w(ij,l+1)-w(ij,l))  
    419 !             dq(ij,l) = (wq(ij,l+1)-wq(ij,l)) !================>>>> 
    420              q(ij,l)=(q(ij,l)*masse(ij,l)+wq(ij,l+1)-wq(ij,l))/newmasse     
    421 !             masse(ij,l)=newmasse 
    422           ENDDO 
    423        ENDDO 
    424     END DO 
    425     RETURN 
     239             newmass = mass(ij,l) + fac*(wfluxt(ij,l)-wfluxt(ij,l+1)) 
     240             q(ij,l) = ( q(ij,l)*mass(ij,l) + wq(ij,l)-wq(ij,l+1) ) / newmass 
     241             IF(update_mass) mass(ij,l)=newmass 
     242          ENDDO 
     243       ENDDO 
     244    END DO 
     245 
    426246  END SUBROUTINE vlz 
    427247 
Note: See TracChangeset for help on using the changeset viewer.