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

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

File:
1 edited

Legend:

Unmodified
Added
Removed

  • codes/icosagcm/trunk/src/advect_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 
Note: See TracChangeset for help on using the changeset viewer.