Changeset 599


Ignore:
Timestamp:
10/19/17 17:04:26 (7 years ago)
Author:
dubos
Message:

trunk : backported commits r582-r598 (transport diagnostics)

Location:
codes/icosagcm/trunk/src
Files:
5 added
1 deleted
8 edited
2 moved

Legend:

Unmodified
Added
Removed

  • codes/icosagcm/trunk/src/diagnostics/check_conserve.f90

    r548 r599  
    299299    REAL(rstd),INTENT(INOUT) :: e, s, AAM_mass, AAM_vel, rmsv 
    300300     
    301     REAL(rstd) :: ucenter(iim*jjm,3,llm) 
     301    REAL(rstd) :: ucenter(iim*jjm,llm,3) 
    302302    REAL(rstd) :: ulon(iim*jjm,llm) 
    303303    REAL(rstd) :: ulat(iim*jjm,llm) 

  • codes/icosagcm/trunk/src/diagnostics/observable.f90

    r581 r599  
    66  TYPE(t_field),POINTER, SAVE :: f_buf_i(:), & 
    77       f_buf_Fel(:), f_buf_uh(:), & ! horizontal velocity, different from prognostic velocity if NH 
    8        f_buf_ulon(:), f_buf_ulat(:), & 
    9        f_buf_u3d(:) ! unused, remove ? 
     8       f_buf_ulon(:), f_buf_ulat(:) 
    109  TYPE(t_field),POINTER, SAVE :: f_buf1_i(:), f_buf2_i(:) 
    1110  TYPE(t_field),POINTER, SAVE :: f_buf_v(:), f_buf_s(:), f_buf_p(:) 
    … …  
    2423    CALL allocate_field(f_buf2_i,   field_t,type_real,llm,name="buffer2_i") 
    2524    CALL allocate_field(f_buf_p,   field_t,type_real,llm+1)  
    26     CALL allocate_field(f_buf_u3d, field_t,type_real,3,llm)  ! 3D vel at cell centers 
    2725    CALL allocate_field(f_buf_ulon,field_t,type_real,llm, name="buf_ulon") 
    2826    CALL allocate_field(f_buf_ulat,field_t,type_real,llm, name="buf_ulat") 
    … …  
    4947    USE theta2theta_rhodz_mod 
    5048    USE omega_mod 
     49    USE diagflux_mod 
    5150    LOGICAL, INTENT(IN) :: init 
    5251    INTEGER :: l 
    … …  
    160159    END IF 
    161160 
     161    IF(.NOT. init) THEN 
     162       IF(diagflux_on) THEN 
     163          CALL flux_centered_lonlat(1./(itau_out*dt) , f_massfluxt, f_buf_ulon, f_buf_ulat) 
     164          CALL output_field("mass_t", f_masst) 
     165          CALL output_field("massflux_lon",f_buf_ulon) 
     166          CALL output_field("massflux_lat",f_buf_ulat) 
     167 
     168          CALL transfert_request(f_epotfluxt,req_e1_vect)  
     169          CALL flux_centered_lonlat(1./(itau_out*dt) , f_epotfluxt, f_buf_ulon, f_buf_ulat) 
     170          CALL output_field("epot_t", f_epot) 
     171          CALL output_field("epotflux_lon",f_buf_ulon) 
     172          CALL output_field("epotflux_lat",f_buf_ulat) 
     173 
     174          CALL transfert_request(f_ekinfluxt,req_e1_vect)  
     175          CALL flux_centered_lonlat(1./(itau_out*dt) , f_ekinfluxt, f_buf_ulon, f_buf_ulat) 
     176          CALL output_field("ekin_t", f_ekin) 
     177          CALL output_field("ekinflux_lon",f_buf_ulon) 
     178          CALL output_field("ekinflux_lat",f_buf_ulat) 
     179 
     180          CALL transfert_request(f_enthalpyfluxt,req_e1_vect)  
     181          CALL flux_centered_lonlat(1./(itau_out*dt) , f_enthalpyfluxt, f_buf_ulon, f_buf_ulat) 
     182          CALL output_field("enthalpy_t", f_enthalpy) 
     183          CALL output_field("enthalpyflux_lon",f_buf_ulon) 
     184          CALL output_field("enthalpyflux_lat",f_buf_ulat) 
     185 
     186          CALL qflux_centered_lonlat(1./(itau_out*dt) , f_qfluxt, f_qfluxt_lon, f_qfluxt_lat) 
     187          CALL output_field("qmass_t", f_qmasst) 
     188          CALL output_field("qflux_lon",f_qfluxt_lon) 
     189          CALL output_field("qflux_lat",f_qfluxt_lat) 
     190          CALL zero_qfluxt  ! restart accumulating fluxes 
     191       END IF 
     192    END IF 
    162193  END SUBROUTINE write_output_fields_basic 
    163194   

  • codes/icosagcm/trunk/src/diagnostics/wind.F90

    r581 r599  
    11MODULE wind_mod 
    2  
    3 CONTAINS 
    4  
    5   SUBROUTINE un2ulonlat(f_u, f_ulon, f_ulat) 
     2  USE omp_para 
    63  USE icosa 
    74  IMPLICIT NONE 
     5  PRIVATE 
     6 
     7  PUBLIC :: compute_wind_centered, compute_flux_centered, & 
     8       compute_wind_centered_lonlat_compound, compute_wind2d_perp_from_lonlat_compound, & 
     9       compute_wind_centered_from_lonlat_compound, compute_wind_perp_from_lonlat_compound, & 
     10       un2ulonlat, ulonlat2un 
     11        
     12CONTAINS 
     13 
     14  SUBROUTINE un2ulonlat(f_u, f_ulon, f_ulat) 
    815    TYPE(t_field), POINTER :: f_u(:) ! IN  : normal velocity components on edges 
    9     TYPE(t_field), POINTER :: f_ulon(:), f_ulat(:) ! OUT : velocity reconstructed at hexagons 
    10      
     16    TYPE(t_field), POINTER :: f_ulon(:), f_ulat(:) ! OUT : velocity reconstructed at hexagons     
    1117    REAL(rstd),POINTER :: u(:,:),  ulon(:,:), ulat(:,:) 
    1218    INTEGER :: ind 
    … …  
    2531 
    2632  SUBROUTINE ulonlat2un(f_ulon, f_ulat,f_u) 
    27   USE icosa 
    28   IMPLICIT NONE 
    2933    TYPE(t_field), POINTER :: f_ulon(:), f_ulat(:) ! IN : velocity reconstructed at hexagons 
    3034    TYPE(t_field), POINTER :: f_u(:) ! OUT  : normal velocity components on edges 
    … …  
    4650  
    4751  SUBROUTINE compute_wind_centered(ue,ucenter) 
    48   USE icosa 
    49   USE omp_para 
    50   IMPLICIT NONE 
    5152  REAL(rstd) :: ue(3*iim*jjm,llm) 
    52   REAL(rstd) :: ucenter(iim*jjm,3,llm) 
    53   INTEGER :: i,j,ij,l     
    54   
     53  REAL(rstd) :: ucenter(iim*jjm,llm,3) 
     54  INTEGER :: ij,l 
     55  REAL(rstd), PARAMETER :: scale=1. 
     56  REAL(rstd) :: fac, ue_le, cx,cy,cz, ux,uy,uz 
     57#include "../kernels/wind_centered.k90" 
     58 END SUBROUTINE compute_wind_centered 
     59  
     60  SUBROUTINE compute_flux_centered(scale,ue,ucenter) 
     61  REAL(rstd), INTENT(IN) :: scale 
     62  REAL(rstd) :: ue(3*iim*jjm,llm) 
     63  REAL(rstd) :: ucenter(iim*jjm,llm,3) 
     64  INTEGER :: ij,l 
     65  REAL(rstd) :: fac, ue_le, cx,cy,cz, ux,uy,uz 
     66#include "../kernels/flux_centered.k90" 
     67  END SUBROUTINE compute_flux_centered 
     68  
     69   
     70 SUBROUTINE compute_wind_on_edge(ue,uedge) 
     71  REAL(rstd) :: ue(3*iim*jjm,llm) 
     72  REAL(rstd) :: uedge(3*iim*jjm,llm,3) 
     73 
     74  REAL(rstd) :: ut(3*iim*jjm,llm) 
     75  INTEGER :: i,j,ij,l      
     76     
     77    CALL compute_tangential_compound(ue,ut) 
     78   
    5579    DO l=ll_begin,ll_end 
    5680      DO j=jj_begin,jj_end 
    5781        DO i=ii_begin,ii_end 
    5882          ij=(j-1)*iim+i 
    59           ucenter(ij,:,l)=1/Ai(ij)*                                                                                                & 
    60                         ( ne(ij,right)*ue(ij+u_right,l)*le(ij+u_right)*((xyz_v(ij+z_rdown,:)+xyz_v(ij+z_rup,:))/2-centroid(ij,:))  & 
    61                          + ne(ij,rup)*ue(ij+u_rup,l)*le(ij+u_rup)*((xyz_v(ij+z_rup,:)+xyz_v(ij+z_up,:))/2-centroid(ij,:))          & 
    62                          + ne(ij,lup)*ue(ij+u_lup,l)*le(ij+u_lup)*((xyz_v(ij+z_up,:)+xyz_v(ij+z_lup,:))/2-centroid(ij,:))          & 
    63                          + ne(ij,left)*ue(ij+u_left,l)*le(ij+u_left)*((xyz_v(ij+z_lup,:)+xyz_v(ij+z_ldown,:))/2-centroid(ij,:))    & 
    64                          + ne(ij,ldown)*ue(ij+u_ldown,l)*le(ij+u_ldown)*((xyz_v(ij+z_ldown,:)+xyz_v(ij+z_down,:))/2-centroid(ij,:))& 
    65                          + ne(ij,rdown)*ue(ij+u_rdown,l)*le(ij+u_rdown)*((xyz_v(ij+z_down,:)+xyz_v(ij+z_rdown,:))/2-centroid(ij,:))) 
    66         ENDDO 
    67       ENDDO 
    68     ENDDO 
    69   
    70  END SUBROUTINE compute_wind_centered 
    71   
    72    
    73  SUBROUTINE compute_wind_on_edge(ue,uedge) 
    74   USE icosa 
    75   USE omp_para 
    76      
    77   IMPLICIT NONE 
    78   REAL(rstd) :: ue(3*iim*jjm,llm) 
    79   REAL(rstd) :: uedge(3*iim*jjm,3,llm) 
    80  
    81   REAL(rstd) :: ut(3*iim*jjm,llm) 
    82   INTEGER :: i,j,ij,l      
    83      
    84     CALL compute_tangential_compound(ue,ut) 
    85    
    86     DO l=ll_begin,ll_end 
    87       DO j=jj_begin,jj_end 
    88         DO i=ii_begin,ii_end 
    89           ij=(j-1)*iim+i 
    90           uedge(ij+u_right,:,l)=ue(ij+u_right,l)*ep_e(ij+u_right,:)*ne(ij,right) + ut(ij+u_right,l)*et_e(ij+u_right,:)*ne(ij,right)  
    91           uedge(ij+u_lup,:,l)=ue(ij+u_lup,l)*ep_e(ij+u_lup,:)*ne(ij,lup) + ut(ij+u_lup,l)*et_e(ij+u_lup,:)*ne(ij,lup) 
    92           uedge(ij+u_ldown,:,l)=ue(ij+u_ldown,l)*ep_e(ij+u_ldown,:)*ne(ij,ldown) + ut(ij+u_ldown,l)*et_e(ij+u_ldown,:)*ne(ij,ldown) 
     83          uedge(ij+u_right,l,:)=ue(ij+u_right,l)*ep_e(ij+u_right,:)*ne(ij,right) + ut(ij+u_right,l)*et_e(ij+u_right,:)*ne(ij,right)  
     84          uedge(ij+u_lup,l,:)=ue(ij+u_lup,l)*ep_e(ij+u_lup,:)*ne(ij,lup) + ut(ij+u_lup,l)*et_e(ij+u_lup,:)*ne(ij,lup) 
     85          uedge(ij+u_ldown,l,:)=ue(ij+u_ldown,l)*ep_e(ij+u_ldown,:)*ne(ij,ldown) + ut(ij+u_ldown,l)*et_e(ij+u_ldown,:)*ne(ij,ldown) 
    9386        ENDDO 
    9487      ENDDO 
    … …  
    10093  
    10194 SUBROUTINE compute_tangential_compound(ue,ut) 
    102   USE icosa   
    103   USE omp_para 
    104   IMPLICIT NONE 
    10595  REAL(rstd) :: ue(3*iim*jjm,llm) 
    10696  REAL(rstd) :: ut(3*iim*jjm,llm) 
    … …  
    155145 END SUBROUTINE compute_tangential_compound 
    156146  
    157  SUBROUTINE compute_wind_lonlat_compound(u, ulon, ulat) 
    158   USE icosa   
    159   USE omp_para 
    160      
    161   IMPLICIT NONE 
    162   REAL(rstd) :: u(3*iim*jjm,3,llm) 
    163   REAL(rstd) :: ulon(3*iim*jjm,3,llm) 
    164   REAL(rstd) :: ulat(3*iim*jjm,3,llm) 
    165  
    166   INTEGER :: i,j,ij,l      
    167      
     147! SUBROUTINE compute_wind_lonlat_compound(u, ulon, ulat) 
     148!  REAL(rstd) :: u(3*iim*jjm,3,llm) 
     149!  REAL(rstd) :: ulon(3*iim*jjm,3,llm) 
     150!  REAL(rstd) :: ulat(3*iim*jjm,3,llm) 
     151! 
     152!  INTEGER :: i,j,ij,l      
     153!     
     154!   
     155!    DO l=ll_begin,ll_end 
     156!      DO j=jj_begin-1,jj_end+1 
     157!        DO i=ii_begin-1,ii_end+1 
     158!          ij=(j-1)*iim+i 
     159!          ulon(ij+u_right,:,l)=sum(u(ij+u_right,:,l)*elon_e(ij+u_right,:))*elon_e(ij+u_right,:)  
     160!          ulon(ij+u_lup,:,l)=sum(u(ij+u_lup,:,l)*elon_e(ij+u_lup,:))*elon_e(ij+u_lup,:) 
     161!          ulon(ij+u_ldown,:,l)=sum(u(ij+u_ldown,:,l)*elon_e(ij+u_ldown,:))*elon_e(ij+u_ldown,:) 
     162!           
     163!          ulat(ij+u_right,:,l)=sum(u(ij+u_right,:,l)*elat_e(ij+u_right,:))*elat_e(ij+u_right,:)  
     164!          ulat(ij+u_lup,:,l)=sum(u(ij+u_lup,:,l)*elat_e(ij+u_lup,:))*elat_e(ij+u_lup,:)  
     165!          ulat(ij+u_ldown,:,l)=sum(u(ij+u_ldown,:,l)*elat_e(ij+u_ldown,:))*elat_e(ij+u_ldown,:)  
     166!           
     167!        ENDDO 
     168!      ENDDO 
     169!    ENDDO 
     170!  
     171! END SUBROUTINE compute_wind_lonlat_compound 
     172  
     173  SUBROUTINE compute_wind_from_lonlat_compound(ulon, ulat, u) 
     174  REAL(rstd) :: u(3*iim*jjm,llm,3) 
     175  REAL(rstd) :: ulon(3*iim*jjm,llm) 
     176  REAL(rstd) :: ulat(3*iim*jjm,llm) 
     177 
     178  INTEGER :: i,j,ij,l      
    168179   
    169180    DO l=ll_begin,ll_end 
    … …  
    171182        DO i=ii_begin-1,ii_end+1 
    172183          ij=(j-1)*iim+i 
    173           ulon(ij+u_right,:,l)=sum(u(ij+u_right,:,l)*elon_e(ij+u_right,:))*elon_e(ij+u_right,:)  
    174           ulon(ij+u_lup,:,l)=sum(u(ij+u_lup,:,l)*elon_e(ij+u_lup,:))*elon_e(ij+u_lup,:) 
    175           ulon(ij+u_ldown,:,l)=sum(u(ij+u_ldown,:,l)*elon_e(ij+u_ldown,:))*elon_e(ij+u_ldown,:) 
    176            
    177           ulat(ij+u_right,:,l)=sum(u(ij+u_right,:,l)*elat_e(ij+u_right,:))*elat_e(ij+u_right,:)  
    178           ulat(ij+u_lup,:,l)=sum(u(ij+u_lup,:,l)*elat_e(ij+u_lup,:))*elat_e(ij+u_lup,:)  
    179           ulat(ij+u_ldown,:,l)=sum(u(ij+u_ldown,:,l)*elat_e(ij+u_ldown,:))*elat_e(ij+u_ldown,:)  
    180            
    181         ENDDO 
    182       ENDDO 
    183     ENDDO 
    184   
    185  END SUBROUTINE compute_wind_lonlat_compound 
    186   
    187   SUBROUTINE compute_wind_from_lonlat_compound(ulon, ulat, u) 
    188   USE icosa   
    189   USE omp_para 
    190      
    191   IMPLICIT NONE 
    192   REAL(rstd) :: u(3*iim*jjm,3,llm) 
    193   REAL(rstd) :: ulon(3*iim*jjm,llm) 
    194   REAL(rstd) :: ulat(3*iim*jjm,llm) 
    195  
    196   INTEGER :: i,j,ij,l      
    197    
    198     DO l=ll_begin,ll_end 
    199       DO j=jj_begin-1,jj_end+1 
    200         DO i=ii_begin-1,ii_end+1 
    201           ij=(j-1)*iim+i 
    202           u(ij+u_right,:,l)=ulon(ij+u_right,l)*elon_e(ij+u_right,:)+ ulat(ij+u_right,l)*elat_e(ij+u_right,:) 
    203           u(ij+u_lup,:,l)=ulon(ij+u_lup,l)*elon_e(ij+u_lup,:) + ulat(ij+u_lup,l)*elat_e(ij+u_lup,:) 
    204           u(ij+u_ldown,:,l)=ulon(ij+u_ldown,l)*elon_e(ij+u_ldown,:) + ulat(ij+u_ldown,l)*elat_e(ij+u_ldown,:) 
     184          u(ij+u_right,l,:)=ulon(ij+u_right,l)*elon_e(ij+u_right,:)+ ulat(ij+u_right,l)*elat_e(ij+u_right,:) 
     185          u(ij+u_lup,l,:)=ulon(ij+u_lup,l)*elon_e(ij+u_lup,:) + ulat(ij+u_lup,l)*elat_e(ij+u_lup,:) 
     186          u(ij+u_ldown,l,:)=ulon(ij+u_ldown,l)*elon_e(ij+u_ldown,:) + ulat(ij+u_ldown,l)*elat_e(ij+u_ldown,:) 
    205187        ENDDO 
    206188      ENDDO 
    … …  
    210192  
    211193  SUBROUTINE compute_wind_centered_from_lonlat_compound(ulon, ulat, u) 
    212   USE icosa   
    213   USE omp_para 
    214      
    215   IMPLICIT NONE 
    216   REAL(rstd) :: u(iim*jjm,3,llm) 
     194  REAL(rstd) :: u(iim*jjm,llm,3) 
    217195  REAL(rstd) :: ulon(iim*jjm,llm) 
    218196  REAL(rstd) :: ulat(iim*jjm,llm) 
    219  
    220197  INTEGER :: i,j,ij,l      
    221198  DO l=ll_begin,ll_end 
    … …  
    223200        DO i=ii_begin-1,ii_end+1 
    224201          ij=(j-1)*iim+i 
    225           u(ij,:,l)=ulon(ij,l)*elon_i(ij,:)+ ulat(ij,l)*elat_i(ij,:) 
    226         ENDDO 
    227       ENDDO 
    228     ENDDO 
    229   
     202          u(ij,l,:)=ulon(ij,l)*elon_i(ij,:)+ ulat(ij,l)*elat_i(ij,:) 
     203        ENDDO 
     204      ENDDO 
     205    ENDDO  
    230206  END SUBROUTINE compute_wind_centered_from_lonlat_compound 
    231207  
    232208  SUBROUTINE compute_wind2D_from_lonlat_compound(ulon, ulat, u) 
    233   USE icosa   
    234    
    235   IMPLICIT NONE 
    236209  REAL(rstd) :: u(3*iim*jjm,3) 
    237210  REAL(rstd) :: ulon(3*iim*jjm) 
    … …  
    252225  
    253226  SUBROUTINE compute_wind_perp_from_lonlat_compound(ulon, ulat, up) 
    254   USE icosa   
    255   USE omp_para 
    256      
    257   IMPLICIT NONE 
    258227  REAL(rstd) :: up(3*iim*jjm,llm) 
    259228  REAL(rstd) :: ulon(3*iim*jjm,llm) 
    260229  REAL(rstd) :: ulat(3*iim*jjm,llm) 
    261   REAL(rstd) :: u(3*iim*jjm,3,llm) 
     230  REAL(rstd) :: u(3*iim*jjm,llm,3) 
    262231 
    263232  INTEGER :: i,j,ij,l      
    … …  
    269238        DO i=ii_begin-1,ii_end+1 
    270239          ij=(j-1)*iim+i 
    271           up(ij+u_right,l)=sum(u(ij+u_right,:,l)*ep_e(ij+u_right,:)) 
    272           up(ij+u_lup,l)=sum(u(ij+u_lup,:,l)*ep_e(ij+u_lup,:)) 
    273           up(ij+u_ldown,l)=sum(u(ij+u_ldown,:,l)*ep_e(ij+u_ldown,:)) 
     240          up(ij+u_right,l)=sum(u(ij+u_right,l,:)*ep_e(ij+u_right,:)) 
     241          up(ij+u_lup,l)=sum(u(ij+u_lup,l,:)*ep_e(ij+u_lup,:)) 
     242          up(ij+u_ldown,l)=sum(u(ij+u_ldown,l,:)*ep_e(ij+u_ldown,:)) 
    274243        ENDDO 
    275244      ENDDO 
    … …  
    279248    
    280249  SUBROUTINE compute_wind2D_perp_from_lonlat_compound(ulon, ulat, up) 
    281   USE icosa   
    282      
    283   IMPLICIT NONE 
    284250  REAL(rstd) :: up(3*iim*jjm) 
    285251  REAL(rstd) :: ulon(3*iim*jjm) 
    … …  
    302268    
    303269 SUBROUTINE compute_wind_centered_lonlat_compound(uc, ulon, ulat) 
    304   USE icosa   
    305   USE omp_para 
    306      
    307   IMPLICIT NONE 
    308   REAL(rstd) :: uc(iim*jjm,3,llm) 
     270  REAL(rstd) :: uc(iim*jjm,llm,3) 
    309271  REAL(rstd) :: ulon(iim*jjm,llm) 
    310272  REAL(rstd) :: ulat(iim*jjm,llm) 
    311273 
    312274  INTEGER :: i,j,ij,l      
    313      
    314275   
    315276    DO l=ll_begin,ll_end 
    … …  
    317278        DO i=ii_begin,ii_end 
    318279          ij=(j-1)*iim+i 
    319           ulon(ij,l)=sum(uc(ij,:,l)*elon_i(ij,:)) 
    320           ulat(ij,l)=sum(uc(ij,:,l)*elat_i(ij,:))  
     280          ulon(ij,l)=sum(uc(ij,l,:)*elon_i(ij,:)) 
     281          ulat(ij,l)=sum(uc(ij,l,:)*elat_i(ij,:))  
    321282        ENDDO 
    322283      ENDDO 
    … …  
    326287 
    327288 SUBROUTINE compute_wind_centered_from_wind_lonlat_centered(ulon, ulat,uc) 
    328   USE icosa   
    329   USE omp_para 
    330      
    331   IMPLICIT NONE 
    332289  REAL(rstd) :: ulon(iim*jjm,llm) 
    333290  REAL(rstd) :: ulat(iim*jjm,llm) 
    334   REAL(rstd) :: uc(iim*jjm,3,llm) 
     291  REAL(rstd) :: uc(iim*jjm,llm,3) 
    335292 
    336293  INTEGER :: i,j,ij,l      
    … …  
    341298        DO i=ii_begin,ii_end 
    342299          ij=(j-1)*iim+i 
    343           uc(ij,:,l)=ulon(ij,l)*elon_i(ij,:)+ulat(ij,l)*elat_i(ij,:) 
     300          uc(ij,l,:)=ulon(ij,l)*elon_i(ij,:)+ulat(ij,l)*elat_i(ij,:) 
    344301        ENDDO 
    345302      ENDDO 
    … …  
    348305 END SUBROUTINE compute_wind_centered_from_wind_lonlat_centered 
    349306 
    350  
    351  
    352307 SUBROUTINE compute_wind_perp_from_wind_centered(uc,un) 
    353   USE icosa   
    354   USE omp_para 
    355      
     308 
    356309  IMPLICIT NONE 
    357   REAL(rstd),INTENT(IN)   :: uc(iim*jjm,3,llm) 
     310  REAL(rstd),INTENT(IN)   :: uc(iim*jjm,llm,3) 
    358311  REAL(rstd),INTENT(OUT)  :: un(3*iim*jjm,llm) 
    359312 
    … …  
    365318        DO i=ii_begin,ii_end 
    366319          ij=(j-1)*iim+i 
    367           un(ij+u_right,l) = sum(0.5*(uc(ij,:,l) + uc(ij+t_right,:,l))*ep_e(ij+u_right,:)) 
    368           un(ij+u_lup,l) = sum(0.5*(uc(ij,:,l) + uc(ij+t_lup,:,l))*ep_e(ij+u_lup,:)) 
    369           un(ij+u_ldown,l) = sum(0.5*(uc(ij,:,l) + uc(ij+t_ldown,:,l))*ep_e(ij+u_ldown,:)) 
     320          un(ij+u_right,l) = sum(0.5*(uc(ij,l,:) + uc(ij+t_right,l,:))*ep_e(ij+u_right,:)) 
     321          un(ij+u_lup,l) = sum(0.5*(uc(ij,l,:) + uc(ij+t_lup,l,:))*ep_e(ij+u_lup,:)) 
     322          un(ij+u_ldown,l) = sum(0.5*(uc(ij,l,:) + uc(ij+t_ldown,l,:))*ep_e(ij+u_ldown,:)) 
    370323         ENDDO 
    371324      ENDDO 
    … …  
    376329 
    377330 SUBROUTINE compute_un2ulonlat(un, ulon, ulat) 
    378   USE icosa   
    379      
    380   IMPLICIT NONE 
    381331  REAL(rstd),INTENT(IN)  :: un(3*iim*jjm,llm) 
    382332  REAL(rstd),INTENT(OUT) :: ulon(iim*jjm,llm) 
    383333  REAL(rstd),INTENT(OUT) :: ulat(iim*jjm,llm) 
    384334 
    385   REAL(rstd)             :: uc(iim*jjm,3,llm) 
     335  REAL(rstd)             :: uc(iim*jjm,llm,3) 
    386336     
    387337  CALL compute_wind_centered(un,uc)   
    … …  
    391341 
    392342 SUBROUTINE compute_ulonlat2un(ulon, ulat,un) 
    393   USE icosa   
    394      
    395   IMPLICIT NONE 
    396343  REAL(rstd),INTENT(IN) :: ulon(iim*jjm,llm) 
    397344  REAL(rstd),INTENT(IN) :: ulat(iim*jjm,llm) 
    398345  REAL(rstd),INTENT(OUT)  :: un(3*iim*jjm,llm) 
    399346 
    400   REAL(rstd)             :: uc(iim*jjm,3,llm) 
     347  REAL(rstd)             :: uc(iim*jjm,llm,3) 
    401348     
    402349    CALL compute_wind_centered_from_wind_lonlat_centered(ulon, ulat, uc)   

  • codes/icosagcm/trunk/src/dynamics/caldyn_gcm.F90

    r580 r599  
    55  USE caldyn_kernels_base_mod 
    66  USE caldyn_kernels_mod 
     7  USE omp_para 
     8  USE mpipara 
    79  IMPLICIT NONE 
    810  PRIVATE 
    … …  
    1315   
    1416  SUBROUTINE init_caldyn 
    15     USE icosa 
    16     USE observable_mod 
    17     USE mpipara 
    18     USE omp_para 
    19     IMPLICIT NONE 
    2017    CHARACTER(len=255) :: def 
    2118    INTEGER            :: ind 
    … …  
    124121 
    125122  SUBROUTINE allocate_caldyn 
    126   USE icosa 
    127   IMPLICIT NONE 
    128  
    129123    CALL allocate_field(f_out_u,field_u,type_real,llm)  
    130124    CALL allocate_field(f_qu,field_u,type_real,llm)  
    … …  
    142136 
    143137  SUBROUTINE caldyn_BC(f_phis, f_geopot, f_wflux) 
    144     USE icosa 
    145     USE mpipara 
    146     USE omp_para 
    147138    TYPE(t_field),POINTER :: f_phis(:) 
    148139    TYPE(t_field),POINTER :: f_geopot(:) 
    … …  
    186177  SUBROUTINE caldyn(write_out,f_phis, f_ps, f_mass, f_theta_rhodz, f_u, f_q, & 
    187178       f_geopot, f_hflux, f_wflux, f_dps, f_dmass, f_dtheta_rhodz, f_du) 
    188     USE icosa 
    189179    USE observable_mod 
    190180    USE disvert_mod, ONLY : caldyn_eta, eta_mass 
    191     USE vorticity_mod 
    192     USE kinetic_mod 
    193     USE theta2theta_rhodz_mod 
    194     USE wind_mod 
    195     USE mpipara 
    196181    USE trace 
    197     USE omp_para 
    198     USE output_field_mod 
    199     USE checksum_mod 
    200     IMPLICIT NONE 
    201182    LOGICAL,INTENT(IN)    :: write_out 
    202183    TYPE(t_field),POINTER :: f_phis(:) 
    … …  
    375356 
    376357  SUBROUTINE check_mass_conservation(f_ps,f_dps) 
    377   USE icosa 
    378   USE mpipara 
    379   IMPLICIT NONE 
    380358    TYPE(t_field),POINTER :: f_ps(:) 
    381359    TYPE(t_field),POINTER :: f_dps(:) 

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

    r554 r599  
    11MODULE icosa_init_mod 
    2  
    3  
    42 
    53CONTAINS 
    … …  
    2220  USE restart_mod 
    2321  USE etat0_mod 
     22  USE diagflux_mod 
    2423  IMPLICIT NONE 
    2524   
    … …  
    5554    CALL init_timeloop 
    5655    CALL init_physics 
    57     
     56 
     57    CALL init_diagflux 
    5858    CALL timeloop 
    5959    CALL switch_omp_no_distrib_level 

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

    r548 r599  
    247247 
    248248    REAL(rstd) :: p(iim*jjm,llm+1) 
    249     REAL(rstd) :: uc(iim*jjm,3,llm) 
     249    REAL(rstd) :: uc(iim*jjm,llm,3) 
    250250    REAL(rstd) :: ulon(iim*jjm,llm) 
    251251    REAL(rstd) :: ulat(iim*jjm,llm) 
    … …  
    311311    REAL(rstd) :: dulat(iim*jjm,llm) 
    312312    REAL(rstd) :: ue(3*iim*jjm,llm) 
    313     REAL(rstd) :: duc(iim*jjm,3,llm) 
     313    REAL(rstd) :: duc(iim*jjm,llm,3) 
    314314    REAL(rstd) :: dt2, due 
    315315    INTEGER :: i,j,ij,l 
    … …  
    321321        DO i=ii_begin,ii_end 
    322322          ij=(j-1)*iim+i 
    323           due = sum( (duc(ij,:,l) + duc(ij+t_right,:,l))*ep_e(ij+u_right,:) ) 
     323          due = sum( (duc(ij,l,:) + duc(ij+t_right,l,:))*ep_e(ij+u_right,:) ) 
    324324          ue(ij+u_right,l) = ue(ij+u_right,l) + dt2*due 
    325325 
    326           due = sum( (duc(ij,:,l) + duc(ij+t_lup,:,l))*ep_e(ij+u_lup,:) ) 
     326          due = sum( (duc(ij,l,:) + duc(ij+t_lup,l,:))*ep_e(ij+u_lup,:) ) 
    327327          ue(ij+u_lup,l)=ue(ij+u_lup,l) + dt2*due 
    328328 
    329           due = sum( (duc(ij,:,l) + duc(ij+t_ldown,:,l))*ep_e(ij+u_ldown,:) ) 
     329          due = sum( (duc(ij,l,:) + duc(ij+t_ldown,l,:))*ep_e(ij+u_ldown,:) ) 
    330330          ue(ij+u_ldown,l)=ue(ij+u_ldown,l) + dt2*due 
    331331        ENDDO 

  • codes/icosagcm/trunk/src/physics/physics_lmdz_generic.F90

    r548 r599  
    7878    CALL allocate_field(f_dq,field_t,type_real,llm,nqtot) 
    7979    CALL allocate_field(f_dps,field_t,type_real) 
    80     CALL allocate_field(f_duc,field_t,type_real,3,llm)     
     80    CALL allocate_field(f_duc,field_t,type_real,llm,3) 
    8181!$OMP END PARALLEL     
    8282 
    … …  
    362362          DO i=ii_begin,ii_end 
    363363            ij=(j-1)*iim+i 
    364             duc(ij,:,l)=dulon(ij,l)*elon_i(ij,:)+dulat(ij,l)*elat_i(ij,:) 
     364            duc(ij,l,:)=dulon(ij,l)*elon_i(ij,:)+dulat(ij,l)*elat_i(ij,:) 
    365365          ENDDO 
    366366        ENDDO 
    … …  
    371371          DO i=ii_begin,ii_end 
    372372            ij=(j-1)*iim+i 
    373             u(ij+u_right,l) = u(ij+u_right,l) + dtphy * sum( 0.5*(duc(ij,:,l) + duc(ij+t_right,:,l))*ep_e(ij+u_right,:) ) 
    374             u(ij+u_lup,l) = u(ij+u_lup,l) + dtphy * sum( 0.5*(duc(ij,:,l) + duc(ij+t_lup,:,l))*ep_e(ij+u_lup,:) ) 
    375             u(ij+u_ldown,l) = u(ij+u_ldown,l) + dtphy*sum( 0.5*(duc(ij,:,l) + duc(ij+t_ldown,:,l))*ep_e(ij+u_ldown,:) ) 
     373            u(ij+u_right,l) = u(ij+u_right,l) + dtphy * sum( 0.5*(duc(ij,l,:) + duc(ij+t_right,l,:))*ep_e(ij+u_right,:) ) 
     374            u(ij+u_lup,l) = u(ij+u_lup,l) + dtphy * sum( 0.5*(duc(ij,l,:) + duc(ij+t_lup,l,:))*ep_e(ij+u_lup,:) ) 
     375            u(ij+u_ldown,l) = u(ij+u_ldown,l) + dtphy*sum( 0.5*(duc(ij,l,:) + duc(ij+t_ldown,l,:))*ep_e(ij+u_ldown,:) ) 
    376376          ENDDO 
    377377        ENDDO 

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

    r581 r599  
    168168    USE caldyn_mod 
    169169    USE advect_tracer_mod 
     170    USE diagflux_mod 
    170171    USE physics_mod 
    171172    USE mpipara 
    … …  
    181182    REAL(rstd),POINTER :: rhodz(:,:), mass(:,:), ps(:) 
    182183 
    183     INTEGER :: ind 
    184     INTEGER :: it,i,j,l,n,  stage 
    185     LOGICAL :: fluxt_zero(ndomain) ! set to .TRUE. to start accumulating fluxes in time 
     184    REAL(rstd) :: adv_over_out ! ratio itau_adv/itau_out 
     185    INTEGER :: ind, it,i,j,l,n,  stage 
     186    LOGICAL :: fluxt_zero(ndomain) ! set to .TRUE. to start accumulating mass fluxes in time 
    186187    LOGICAL, PARAMETER :: check_rhodz=.FALSE. 
    187188    INTEGER :: start_clock, stop_clock, rate_clock 
    … …  
    210211 
    211212    IF(positive_theta) CALL copy_theta_to_q(f_theta_rhodz,f_rhodz,f_q) 
     213    IF(diagflux_on) THEN 
     214       adv_over_out = itau_adv*(1./itau_out) 
     215    ELSE 
     216       adv_over_out = 0. 
     217    END IF 
    212218 
    213219    CALL check_conserve(f_ps,f_dps,f_u,f_theta_rhodz,f_phis,itau0)   
    214220 
    215     CALL trace_on 
     221    Call trace_on 
    216222 
    217223    IF (xios_output) THEN ! we must call update_calendar before any XIOS output 
    … …  
    267273          CALL HEVI_scheme(it, fluxt_zero) 
    268274       END SELECT 
    269         
     275 
    270276       IF (MOD(it,itau_dissip)==0) THEN 
    271277           
    … …  
    298304        
    299305       IF(MOD(it,itau_adv)==0) THEN 
    300           CALL advect_tracer(f_hfluxt,f_wfluxt,f_u, f_q,f_rhodz)  ! update q and rhodz after RK step 
    301           fluxt_zero=.TRUE. 
    302           ! FIXME : check that rhodz is consistent with ps 
    303           IF (check_rhodz) THEN 
     306          CALL advect_tracer(f_hfluxt,f_wfluxt,f_u, f_q,f_rhodz, & ! update q and rhodz after RK step 
     307               adv_over_out, f_masst,f_qmasst,f_massfluxt, f_qfluxt)  ! accumulate mass and fluxes if diagflux_on 
     308          fluxt_zero=.TRUE. ! restart accumulation of hfluxt and qfluxt at next call to explicit_scheme / HEVI_scheme 
     309          ! At this point advect_tracer has obtained the halos of u and rhodz, 
     310          ! needed for correct computation of kinetic energy 
     311          IF(diagflux_on) CALL diagflux_energy(adv_over_out, f_phis,f_rhodz,f_theta_rhodz,f_u, f_geopot,f_theta, f_hfluxt) 
     312 
     313          IF (check_rhodz) THEN ! check that rhodz is consistent with ps 
    304314             DO ind=1,ndomain 
    305315                IF (.NOT. assigned_domain(ind)) CYCLE 

  • codes/icosagcm/trunk/src/transport/advect.F90

    r581 r599  
    425425  ! Horizontal transport (S. Dubey, T. Dubos) 
    426426  ! Slope-limited van Leer approach with hexagons 
    427   SUBROUTINE compute_advect_horiz(update_mass,hfluxt,cc,gradq3d, mass,qi)  
     427  SUBROUTINE compute_advect_horiz(update_mass,diagflux_on, hfluxt,cc,gradq3d, mass,qi,qfluxt) 
    428428    USE trace 
    429429    USE omp_para 
    430430    IMPLICIT NONE 
    431     LOGICAL, INTENT(IN)       :: update_mass 
     431    LOGICAL, INTENT(IN)       :: update_mass, diagflux_on 
    432432    REAL(rstd), INTENT(IN)    :: gradq3d(iim*jjm,llm,3)  
    433433    REAL(rstd), INTENT(IN)    :: hfluxt(3*iim*jjm,llm) ! mass flux 
    … …  
    435435    REAL(rstd), INTENT(INOUT) :: mass(iim*jjm,llm) 
    436436    REAL(rstd), INTENT(INOUT) :: qi(iim*jjm,llm) 
     437    REAL(rstd), INTENT(INOUT) :: qfluxt(3*iim*jjm,llm) ! time-integrated tracer flux 
    437438 
    438439    REAL(rstd) :: dq,dmass,qe,ed(3), newmass 
    … …  
    441442 
    442443    CALL trace_start("compute_advect_horiz") 
    443  
    444     ! evaluate tracer field at point cc using piecewise linear reconstruction 
    445     ! q(cc)= q0 + gradq.(cc-xyz_i)  with xi centroid of hexagon 
    446     ! ne*hfluxt>0 iff outgoing 
    447     DO l = ll_begin,ll_end 
    448  
    449 !DIR$ SIMD 
    450       DO ij=ij_begin_ext,ij_end_ext 
    451  
    452              IF (ne(ij,right)*hfluxt(ij+u_right,l)>0) THEN  
    453                 ed = cc(ij+u_right,l,:) - xyz_i(ij,:) 
    454 !                qe = qi(ij,l)+sum2(gradq3d(ij,l,:),ed)  
    455                  qe = qi(ij,l)+gradq3d(ij,l,1)*ed(1)+gradq3d(ij,l,2)*ed(2)+gradq3d(ij,l,3)*ed(3)  
    456              ELSE 
    457                 ed = cc(ij+u_right,l,:) - xyz_i(ij+t_right,:) 
    458 !                qe = qi(ij+t_right,l)+sum2(gradq3d(ij+t_right,l,:),ed)  
    459                 qe = qi(ij+t_right,l) + gradq3d(ij+t_right,l,1)*ed(1)+gradq3d(ij+t_right,l,2)*ed(2)+gradq3d(ij+t_right,l,3)*ed(3)  
    460              ENDIF 
    461              qflux(ij+u_right,l) = hfluxt(ij+u_right,l)*qe 
    462               
    463              IF (ne(ij,lup)*hfluxt(ij+u_lup,l)>0) THEN  
    464                 ed = cc(ij+u_lup,l,:) - xyz_i(ij,:) 
    465 !                qe = qi(ij,l)+sum2(gradq3d(ij,l,:),ed) 
    466                 qe = qi(ij,l)  + gradq3d(ij,l,1)*ed(1)+gradq3d(ij,l,2)*ed(2)+gradq3d(ij,l,3)*ed(3)  
    467              ELSE 
    468                 ed = cc(ij+u_lup,l,:) - xyz_i(ij+t_lup,:) 
    469 !                qe = qi(ij+t_lup,l)+sum2(gradq3d(ij+t_lup,l,:),ed)  
    470                 qe = qi(ij+t_lup,l) + gradq3d(ij+t_lup,l,1)*ed(1)+gradq3d(ij+t_lup,l,2)*ed(2)+gradq3d(ij+t_lup,l,3)*ed(3)  
    471              ENDIF 
    472              qflux(ij+u_lup,l) = hfluxt(ij+u_lup,l)*qe  
    473  
    474              IF (ne(ij,ldown)*hfluxt(ij+u_ldown,l)>0) THEN  
    475                 ed = cc(ij+u_ldown,l,:) - xyz_i(ij,:) 
    476 !                qe = qi(ij,l)+sum2(gradq3d(ij,l,:),ed)  
    477                 qe = qi(ij,l) + gradq3d(ij,l,1)*ed(1)+gradq3d(ij,l,2)*ed(2)+gradq3d(ij,l,3)*ed(3)  
    478              ELSE 
    479                 ed = cc(ij+u_ldown,l,:) - xyz_i(ij+t_ldown,:) 
    480 !                qe = qi(ij+t_ldown,l)+sum2(gradq3d(ij+t_ldown,l,:),ed)  
    481                 qe = qi(ij+t_ldown,l)+gradq3d(ij+t_ldown,l,1)*ed(1)+gradq3d(ij+t_ldown,l,2)*ed(2)+gradq3d(ij+t_ldown,l,3)*ed(3)  
    482              ENDIF 
    483              qflux(ij+u_ldown,l) = hfluxt(ij+u_ldown,l)*qe 
    484        END DO 
    485     END DO 
    486  
    487     ! update q and, if update_mass, update mass 
    488     DO l = ll_begin,ll_end 
    489 !DIR$ SIMD 
    490       DO ij=ij_begin,ij_end 
    491              ! sign convention as Ringler et al. (2010) eq. 21 
    492              dmass =   hfluxt(ij+u_right,l)*ne(ij,right)   & 
    493                      + hfluxt(ij+u_lup,l)  *ne(ij,lup)     & 
    494                      + hfluxt(ij+u_ldown,l)*ne(ij,ldown)   & 
    495                      + hfluxt(ij+u_rup,l)  *ne(ij,rup)     & 
    496                      + hfluxt(ij+u_left,l) *ne(ij,left)    & 
    497                      + hfluxt(ij+u_rdown,l)*ne(ij,rdown) 
    498  
    499              dq    =   qflux(ij+u_right,l) *ne(ij,right)   & 
    500                      + qflux(ij+u_lup,l)   *ne(ij,lup)     & 
    501                      + qflux(ij+u_ldown,l) *ne(ij,ldown)   & 
    502                      + qflux(ij+u_rup,l)   *ne(ij,rup)     & 
    503                      + qflux(ij+u_left,l)  *ne(ij,left)    & 
    504                      + qflux(ij+u_rdown,l) *ne(ij,rdown) 
    505  
    506               
    507              newmass = mass(ij,l) - dmass/Ai(ij) 
    508              qi(ij,l) = (qi(ij,l)*mass(ij,l) - dq/Ai(ij) ) / newmass 
    509              IF(update_mass) mass(ij,l) = newmass 
    510  
    511        END DO 
    512     END DO 
    513  
     444#include "../kernels/advect_horiz.k90" 
    514445    CALL trace_end("compute_advect_horiz") 
    515   CONTAINS 
    516  
    517     !==================================================================================== 
    518     PURE REAL(rstd) FUNCTION sum2(a1,a2) 
    519       IMPLICIT NONE 
    520       REAL(rstd),INTENT(IN):: a1(3), a2(3) 
    521       sum2 = a1(1)*a2(1)+a1(2)*a2(2)+a1(3)*a2(3) 
    522       ! sum2 = 0. ! Godunov scheme 
    523     END FUNCTION sum2 
    524  
    525446  END SUBROUTINE compute_advect_horiz 
    526447 

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

    r548 r599  
    11MODULE advect_tracer_mod 
    22  USE icosa 
     3  USE advect_mod 
    34  IMPLICIT NONE 
    45  PRIVATE 
    … …  
    2526 
    2627  SUBROUTINE init_advect_tracer 
    27     USE advect_mod 
    2828    USE omp_para 
    2929    REAL(rstd),POINTER :: tangent(:,:) 
    … …  
    5454  END SUBROUTINE init_advect_tracer 
    5555 
    56   SUBROUTINE advect_tracer(f_hfluxt, f_wfluxt,f_u, f_q,f_rhodz) 
    57     USE advect_mod 
    58     USE mpipara 
     56  SUBROUTINE advect_tracer(f_hfluxt, f_wfluxt,f_u, f_q,f_rhodz,& 
     57       frac, f_masst,f_qmasst,f_massfluxt,f_qfluxt) 
     58    USE omp_para 
    5959    USE trace 
    6060    USE write_field_mod 
    6161    USE tracer_mod 
    62     IMPLICIT NONE 
    63      
    6462    TYPE(t_field),POINTER :: f_hfluxt(:)   ! time-integrated horizontal mass flux 
    6563    TYPE(t_field),POINTER :: f_wfluxt(:)   ! time-integrated vertical mass flux 
    … …  
    6765    TYPE(t_field),POINTER :: f_q(:)        ! tracer 
    6866    TYPE(t_field),POINTER :: f_rhodz(:)    ! mass field at beginning of macro time step 
     67    REAL(rstd), INTENT(in):: frac          ! ratio itau_adv/itau_out or 0. if not diagflux_on 
     68    TYPE(t_field),POINTER :: f_masst(:)    ! time-integrated mass 
     69    TYPE(t_field),POINTER :: f_qmasst(:)   ! time-integrated tracer mass 
     70    TYPE(t_field),POINTER :: f_massfluxt(:)! time-integrated horizontal mass flux 
     71    TYPE(t_field),POINTER :: f_qfluxt(:)   ! time-integrated horizontal tracer flux 
    6972 
    7073    REAL(rstd),POINTER :: q(:,:,:), normal(:,:), tangent(:,:), sqrt_leng(:), gradq3d(:,:,:), cc(:,:,:) 
    71     REAL(rstd),POINTER :: hfluxt(:,:), wfluxt(:,:) 
     74    REAL(rstd),POINTER :: hfluxt(:,:), wfluxt(:,:), masst(:,:), qmasst(:,:,:), massfluxt(:,:), qfluxt(:,:,:) 
    7275    REAL(rstd),POINTER :: rhodz(:,:), u(:,:)  
    7376! temporary shared variable for vlz 
    … …  
    7780    REAL(rstd),POINTER ::  wq(:,:)           ! time-integrated flux of q 
    7881     
    79      INTEGER :: ind,k, nq_last 
     82    INTEGER :: ind,k, nq_last 
    8083    LOGICAL,SAVE :: first=.TRUE. 
    8184!$OMP THREADPRIVATE(first) 
    … …  
    141144    CALL send_message(f_cc,req_cc) 
    142145 
    143  
    144146    ! horizontal transport - split in two to place transfer of gradq3d 
    145147    DO k = 1, nqtot 
    … …  
    154156          sqrt_leng=f_sqrt_leng(ind) 
    155157          CALL compute_gradq3d(q(:,:,k),sqrt_leng,gradq3d,xyz_i,xyz_v) 
    156  
    157158        END DO 
    158159 
    … …  
    160161        CALL wait_message(req_cc) 
    161162        CALL wait_message(req_gradq3d) 
    162  
    163163 
    164164        DO ind=1,ndomain 
    … …  
    170170          rhodz   = f_rhodz(ind) 
    171171          hfluxt  = f_hfluxt(ind)  
     172          qfluxt  = f_qfluxt(ind)  
    172173          gradq3d = f_gradq3d(ind) 
    173           CALL compute_advect_horiz(k==nq_last,hfluxt,cc,gradq3d, rhodz,q(:,:,k)) 
     174 
     175          IF(frac>0.) THEN ! accumulate mass, mass flux and tracer mass 
     176             qmasst  = f_qmasst(ind)  
     177             qmasst(:,ll_begin:ll_end,k) = qmasst(:,ll_begin:ll_end,k) + & 
     178                  frac*rhodz(:,ll_begin:ll_end)*q(:,ll_begin:ll_end,k) 
     179             IF(k==nq_last) THEN 
     180                masst  = f_masst(ind) 
     181                massfluxt  = f_massfluxt(ind) 
     182                masst(:,ll_begin:ll_end) = masst(:,ll_begin:ll_end)+frac*rhodz(:,ll_begin:ll_end) 
     183                massfluxt(:,ll_begin:ll_end) = massfluxt(:,ll_begin:ll_end)+hfluxt(:,ll_begin:ll_end) 
     184             END IF 
     185          END IF 
     186          CALL compute_advect_horiz(k==nq_last,frac>0., hfluxt,cc,gradq3d, rhodz, q(:,:,k), qfluxt(:,:,k)) 
    174187        END DO 
     188 
    175189      ENDIF 
    176190    END DO  
    … …  
    214228    USE trace 
    215229    USE omp_para 
    216     IMPLICIT NONE 
    217230    LOGICAL, INTENT(IN)       :: update_mass 
    218231    REAL(rstd), INTENT(IN)    :: fac, wfluxt(iim*jjm,llm+1) ! vertical mass flux 
Note: See TracChangeset for help on using the changeset viewer.