Index: codes/icosagcm/trunk/src/caldyn_gcm.f90
===================================================================
--- codes/icosagcm/trunk/src/caldyn_gcm.f90	(revision 386)
+++ codes/icosagcm/trunk/src/caldyn_gcm.f90	(revision 387)
@@ -133,5 +133,5 @@
     
     REAL(rstd),POINTER :: ps(:), dps(:)
-    REAL(rstd),POINTER :: mass(:,:), theta_rhodz(:,:), dtheta_rhodz(:,:)
+    REAL(rstd),POINTER :: mass(:,:), theta_rhodz(:,:,:), dtheta_rhodz(:,:,:)
     REAL(rstd),POINTER :: u(:,:), du(:,:), hflux(:,:), wflux(:,:)
     REAL(rstd),POINTER :: qu(:,:)
@@ -194,5 +194,5 @@
           qu=f_qu(ind)
           qv=f_qv(ind)
-          CALL compute_pvort(ps,u,theta_rhodz, mass,theta,qu,qv) ! COM00 COM01 COM02
+          CALL compute_pvort(ps,u,theta_rhodz(:,:,1), mass,theta,qu,qv) ! COM00 COM01 COM02
        ENDDO
 !       CALL checksum(f_mass)
@@ -208,5 +208,4 @@
           ps=f_ps(ind)
           u=f_u(ind)
-          theta_rhodz=f_theta_rhodz(ind)
           mass=f_mass(ind)
           theta = f_theta(ind)
@@ -219,10 +218,10 @@
           dtheta_rhodz=f_dtheta_rhodz(ind)
           du=f_du(ind)
-          CALL compute_caldyn_horiz(u,mass,qu,theta,pk,geopot, hflux,convm,dtheta_rhodz,du)
+          CALL compute_caldyn_horiz(u,mass,qu,theta,pk,geopot, hflux,convm,dtheta_rhodz(:,:,1),du)
           IF(caldyn_eta==eta_mass) THEN
              wflux=f_wflux(ind)
              wwuu=f_wwuu(ind)
              dps=f_dps(ind)
-             CALL compute_caldyn_vert(u,theta,mass,convm, wflux,wwuu, dps, dtheta_rhodz, du)
+             CALL compute_caldyn_vert(u,theta,mass,convm, wflux,wwuu, dps, dtheta_rhodz(:,:,1), du)
           END IF
        ENDDO       
@@ -245,5 +244,5 @@
           qu=f_qu(ind)
           qv=f_qv(ind)
-          CALL compute_pvort(ps,u,theta_rhodz, mass,theta,qu,qv)
+          CALL compute_pvort(ps,u,theta_rhodz(:,:,1), mass,theta,qu,qv)
           pk = f_pk(ind)
           geopot = f_geopot(ind)  
@@ -253,5 +252,5 @@
           dtheta_rhodz=f_dtheta_rhodz(ind)
           du=f_du(ind)
-          CALL compute_caldyn_horiz(u,mass,qu,theta,pk,geopot, hflux,convm,dtheta_rhodz,du)
+          CALL compute_caldyn_horiz(u,mass,qu,theta,pk,geopot, hflux,convm,dtheta_rhodz(:,:,1),du)
           IF(caldyn_eta==eta_mass) THEN
              wflux=f_wflux(ind)
Index: codes/icosagcm/trunk/src/caldyn_hevi.f90
===================================================================
--- codes/icosagcm/trunk/src/caldyn_hevi.f90	(revision 386)
+++ codes/icosagcm/trunk/src/caldyn_hevi.f90	(revision 387)
@@ -49,5 +49,5 @@
     
     REAL(rstd),POINTER :: ps(:), dps(:)
-    REAL(rstd),POINTER :: mass(:,:), theta_rhodz(:,:), dtheta_rhodz(:,:)
+    REAL(rstd),POINTER :: mass(:,:), theta_rhodz(:,:,:), dtheta_rhodz(:,:,:)
     REAL(rstd),POINTER :: du(:,:), dW(:,:), dPhi(:,:), hflux(:,:), wflux(:,:)
     REAL(rstd),POINTER :: u(:,:), w(:,:), qu(:,:), qv(:,:)
@@ -107,5 +107,5 @@
        mass=f_mass(ind)
        theta = f_theta(ind)
-       CALL compute_theta(ps,theta_rhodz, mass,theta)
+       CALL compute_theta(ps,theta_rhodz(:,:,1), mass,theta)
        pk = f_pk(ind)
        geopot = f_geopot(ind)
@@ -163,5 +163,5 @@
           CALL compute_caldyn_slow_NH(u,mass,geopot,W, hflux,du,dPhi,dW)
        END IF
-       CALL compute_caldyn_Coriolis(hflux,theta,qu, convm,dtheta_rhodz,du)
+       CALL compute_caldyn_Coriolis(hflux,theta,qu, convm,dtheta_rhodz(:,:,1),du)
        IF(caldyn_eta==eta_mass) THEN
           wflux=f_wflux(ind)
Index: codes/icosagcm/trunk/src/check_conserve.f90
===================================================================
--- codes/icosagcm/trunk/src/check_conserve.f90	(revision 386)
+++ codes/icosagcm/trunk/src/check_conserve.f90	(revision 387)
@@ -256,5 +256,5 @@
     REAL(rstd), POINTER :: ue(:,:)
     REAL(rstd), POINTER :: p(:,:) 
-    REAL(rstd), POINTER :: theta_rhodz(:,:) 
+    REAL(rstd), POINTER :: theta_rhodz(:,:,:) 
     REAL(rstd), POINTER :: pk(:,:) 
     REAL(rstd), POINTER :: phis(:) 
@@ -275,5 +275,5 @@
       pk=f_pk(ind) 
       phis=f_phis(ind) 
-      CALL compute_energy(ind,ue,p,rhodz,theta_rhodz,pk,phis, &
+      CALL compute_energy(ind,ue,p,rhodz,theta_rhodz(:,:,1),pk,phis, &
            e, s, AAM_mass, AAM_vel, rmsv)
     END DO 
Index: codes/icosagcm/trunk/src/dissip_gcm.f90
===================================================================
--- codes/icosagcm/trunk/src/dissip_gcm.f90	(revision 386)
+++ codes/icosagcm/trunk/src/dissip_gcm.f90	(revision 387)
@@ -241,6 +241,6 @@
     IF (is_master) PRINT *, 'mean T-cell edge size (km)', 1.45*radius/iim_glo/1000., &
                               'effective T-cell half-edge size (km)', dumax**(-.5/nitergdiv)/1000
-    IF (is_master)  PRINT *, 'Max. time step assuming c=340 m/s and Courant number=4 :', &
-                               4./340.*dumax**(-.5/nitergdiv)
+    IF (is_master)  PRINT *, 'Max. time step assuming c=340 m/s and Courant number=3 (ARK2.3) :', &
+                               3./340.*dumax**(-.5/nitergdiv)
 
     cgraddiv=dumax**(-1./nitergdiv)
@@ -469,5 +469,5 @@
     REAL(rstd),POINTER         :: due_diss1(:,:)
     REAL(rstd),POINTER         :: due_diss2(:,:)
-    REAL(rstd),POINTER         :: dtheta_rhodz(:,:)
+    REAL(rstd),POINTER         :: dtheta_rhodz(:,:,:)
     REAL(rstd),POINTER         :: dtheta_rhodz_diss(:,:)
 
@@ -508,5 +508,5 @@
             due(ij+u_ldown,l) = -0.5*( due_diss1(ij+u_ldown,l)/tau_graddiv(l) + due_diss2(ij+u_ldown,l)/tau_gradrot(l))*itau_dissip
 
-            dtheta_rhodz(ij,l) = -0.5*dtheta_rhodz_diss(ij,l)/tau_divgrad(l)*itau_dissip
+            dtheta_rhodz(ij,l,1) = -0.5*dtheta_rhodz_diss(ij,l)/tau_divgrad(l)*itau_dissip
         ENDDO
       ENDDO
@@ -715,5 +715,5 @@
     
     REAL(rstd),POINTER :: mass(:,:)
-    REAL(rstd),POINTER :: theta_rhodz(:,:)
+    REAL(rstd),POINTER :: theta_rhodz(:,:,:)
     REAL(rstd),POINTER :: dtheta_rhodz(:,:)
 
@@ -733,5 +733,5 @@
 !$SIMD
         DO ij=ij_begin,ij_end
-            dtheta_rhodz(ij,l) = theta_rhodz(ij,l) / mass(ij,l)
+            dtheta_rhodz(ij,l) = theta_rhodz(ij,l,1) / mass(ij,l)
         ENDDO
       ENDDO
Index: codes/icosagcm/trunk/src/etat0.f90
===================================================================
--- codes/icosagcm/trunk/src/etat0.f90	(revision 386)
+++ codes/icosagcm/trunk/src/etat0.f90	(revision 387)
@@ -137,5 +137,5 @@
     REAL(rstd),POINTER :: mass(:,:)
     REAL(rstd),POINTER :: phis(:)
-    REAL(rstd),POINTER :: theta_rhodz(:,:)
+    REAL(rstd),POINTER :: theta_rhodz(:,:,:)
     REAL(rstd),POINTER :: temp(:,:)
     REAL(rstd),POINTER :: u(:,:)
@@ -162,7 +162,7 @@
 
       IF( TRIM(etat0_type)=='williamson91.6' ) THEN
-       CALL compute_etat0_collocated(ps,mass, phis, theta_rhodz, u, geopot, W, q)
+         CALL compute_etat0_collocated(ps,mass, phis, theta_rhodz(:,:,1), u, geopot, W, q)
       ELSE
-       CALL compute_etat0_collocated(ps,mass, phis, temp, u, geopot, W, q)
+         CALL compute_etat0_collocated(ps,mass, phis, temp, u, geopot, W, q)
       ENDIF
     ENDDO
Index: codes/icosagcm/trunk/src/euler_scheme.f90
===================================================================
--- codes/icosagcm/trunk/src/euler_scheme.f90	(revision 386)
+++ codes/icosagcm/trunk/src/euler_scheme.f90	(revision 387)
@@ -33,5 +33,5 @@
     REAL(rstd),POINTER :: u(:,:) , du(:,:)
     REAL(rstd),POINTER :: mass(:,:), dmass(:,:)
-    REAL(rstd),POINTER :: theta_rhodz(:,:), dtheta_rhodz(:,:)
+    REAL(rstd),POINTER :: theta_rhodz(:,:,:), dtheta_rhodz(:,:,:)
     REAL(rstd),POINTER :: hflux(:,:),wflux(:,:),hfluxt(:,:),wfluxt(:,:)
     INTEGER :: ind
@@ -77,5 +77,5 @@
              u(ij+u_lup,l)=u(ij+u_lup,l)+dt*du(ij+u_lup,l)
              u(ij+u_ldown,l)=u(ij+u_ldown,l)+dt*du(ij+u_ldown,l)
-             theta_rhodz(ij,l)=theta_rhodz(ij,l)+dt*dtheta_rhodz(ij,l)
+             theta_rhodz(ij,l,1)=theta_rhodz(ij,l,1)+dt*dtheta_rhodz(ij,l,1)
           ENDDO
        ENDDO
Index: codes/icosagcm/trunk/src/grid_param.f90
===================================================================
--- codes/icosagcm/trunk/src/grid_param.f90	(revision 386)
+++ codes/icosagcm/trunk/src/grid_param.f90	(revision 387)
@@ -4,5 +4,6 @@
   INTEGER,PARAMETER  :: nb_face=10
   INTEGER  :: llm=19
-  INTEGER  :: nqtot
+  INTEGER  :: nqtot ! number of tracers handled by advection scheme
+  INTEGER  :: nqdyn ! number of dynamical tracers : 1 if dry, more if moist
 
 CONTAINS
Index: codes/icosagcm/trunk/src/hevi_scheme.f90
===================================================================
--- codes/icosagcm/trunk/src/hevi_scheme.f90	(revision 386)
+++ codes/icosagcm/trunk/src/hevi_scheme.f90	(revision 387)
@@ -80,14 +80,14 @@
              CALL update_2D(bjl(l,j), f_ps, f_dps_slow(:,l))
           ELSE
-             CALL update(bjl(l,j), f_mass, f_dmass_slow(:,l))
+             CALL update_3D(bjl(l,j), f_mass, f_dmass_slow(:,l))
           END IF
-          CALL update(bjl(l,j), f_theta_rhodz, f_dtheta_rhodz_slow(:,l))
-          CALL update(bjl(l,j), f_u, f_du_slow(:,l))
-          CALL update(cjl(l,j), f_u, f_du_fast(:,l))
+          CALL update_4D(bjl(l,j), f_theta_rhodz, f_dtheta_rhodz_slow(:,l))
+          CALL update_3D(bjl(l,j), f_u, f_du_slow(:,l))
+          CALL update_3D(cjl(l,j), f_u, f_du_fast(:,l))
           IF(.NOT. hydrostatic) THEN
-             CALL update(bjl(l,j), f_W, f_dW_slow(:,l))
-             CALL update(cjl(l,j), f_W, f_dW_fast(:,l))
-             CALL update(bjl(l,j), f_geopot, f_dPhi_slow(:,l))
-             CALL update(cjl(l,j), f_geopot, f_dPhi_fast(:,l))
+             CALL update_3D(bjl(l,j), f_W, f_dW_slow(:,l))
+             CALL update_3D(cjl(l,j), f_W, f_dW_fast(:,l))
+             CALL update_3D(bjl(l,j), f_geopot, f_dPhi_slow(:,l))
+             CALL update_3D(cjl(l,j), f_geopot, f_dPhi_fast(:,l))
           END IF
        END DO
@@ -96,5 +96,24 @@
   END SUBROUTINE HEVI_scheme
   
-  SUBROUTINE update(w, f_y, f_dy)
+  SUBROUTINE update_4D(w, f_y, f_dy)
+    USE dimensions
+    USE grid_param, ONLY : nqdyn
+    REAL(rstd) :: w
+    TYPE(t_field) :: f_y(:), f_dy(:)
+    REAL(rstd), POINTER :: y(:,:,:), dy(:,:,:)
+    INTEGER :: ind, iq
+    IF(w /= 0.) THEN
+       DO ind=1,ndomain
+          IF (.NOT. assigned_domain(ind)) CYCLE
+          CALL swap_dimensions(ind)
+          dy=f_dy(ind); y=f_y(ind)
+          DO iq=1,nqdyn
+             CALL compute_update_3D(w,y(:,:,iq),dy(:,:,iq))
+          END DO
+       ENDDO
+    END IF
+  END SUBROUTINE update_4D
+    
+  SUBROUTINE update_3D(w, f_y, f_dy)
     USE dimensions
     REAL(rstd) :: w
@@ -107,10 +126,10 @@
           CALL swap_dimensions(ind)
           dy=f_dy(ind); y=f_y(ind)
-          CALL compute_update(w,y,dy)
+          CALL compute_update_3D(w,y,dy)
        ENDDO
     END IF
-  END SUBROUTINE update
+  END SUBROUTINE update_3D
     
-  SUBROUTINE compute_update(w, y, dy)
+  SUBROUTINE compute_update_3D(w, y, dy)
     USE omp_para
     USE disvert_mod
@@ -123,5 +142,5 @@
        y(:,l)=y(:,l)+w*dy(:,l)
     ENDDO
-  END SUBROUTINE compute_update
+  END SUBROUTINE compute_update_3D
   
   SUBROUTINE update_2D(w, f_y, f_dy)
Index: codes/icosagcm/trunk/src/observable.f90
===================================================================
--- codes/icosagcm/trunk/src/observable.f90	(revision 386)
+++ codes/icosagcm/trunk/src/observable.f90	(revision 387)
@@ -214,5 +214,5 @@
     TYPE(t_field), POINTER :: f_ps(:), f_phis(:), f_theta_rhodz(:), &  ! IN
          f_pks(:), f_p(:), f_theta(:), f_pk(:), f_phi(:)               ! OUT
-    REAL(rstd),POINTER :: pk(:,:), p(:,:), theta(:,:), theta_rhodz(:,:), &
+    REAL(rstd),POINTER :: pk(:,:), p(:,:), theta(:,:), theta_rhodz(:,:,:), &
          phi(:,:), phis(:), ps(:), pks(:)
     INTEGER :: ind
@@ -233,5 +233,5 @@
        theta_rhodz = f_theta_rhodz(ind)
        theta = f_theta(ind)
-       CALL compute_theta_rhodz2theta(ps, theta_rhodz,theta,0)
+       CALL compute_theta_rhodz2theta(ps, theta_rhodz(:,:,1),theta,0)
        phis = f_phis(ind)
        phi = f_phi(ind)
Index: codes/icosagcm/trunk/src/physics.f90
===================================================================
--- codes/icosagcm/trunk/src/physics.f90	(revision 386)
+++ codes/icosagcm/trunk/src/physics.f90	(revision 387)
@@ -94,9 +94,4 @@
     TYPE(t_field),POINTER :: f_wflux(:)
     TYPE(t_field),POINTER :: f_q(:)
-    REAL(rstd),POINTER :: phis(:)
-    REAL(rstd),POINTER :: ps(:)
-    REAL(rstd),POINTER :: theta_rhodz(:,:)
-    REAL(rstd),POINTER :: ue(:,:)
-    REAL(rstd),POINTER :: q(:,:,:)
 
     LOGICAL:: firstcall,lastcall
@@ -153,5 +148,4 @@
     REAL(rstd),POINTER :: ps(:)
     REAL(rstd),POINTER :: temp(:,:)
-    REAL(rstd),POINTER :: theta_rhodz(:,:)
     REAL(rstd),POINTER :: ue(:,:)
     REAL(rstd),POINTER :: dulon(:,:)
Index: codes/icosagcm/trunk/src/theta_rhodz.f90
===================================================================
--- codes/icosagcm/trunk/src/theta_rhodz.f90	(revision 386)
+++ codes/icosagcm/trunk/src/theta_rhodz.f90	(revision 387)
@@ -58,5 +58,5 @@
   
     REAL(rstd), POINTER :: ps(:)
-    REAL(rstd), POINTER :: theta_rhodz(:,:)
+    REAL(rstd), POINTER :: theta_rhodz(:,:,:)
     REAL(rstd), POINTER :: temp(:,:)
     REAL(rstd), POINTER :: p(:,:)
@@ -81,5 +81,5 @@
       CALL compute_exner(ps,p,pks,pk,0)
 !$OMP BARRIER
-      CALL compute_theta_rhodz2temperature(p, pk, theta_rhodz,temp,0)
+      CALL compute_theta_rhodz2temperature(p, pk, theta_rhodz(:,:,1),temp,0)
     ENDDO
 !$OMP BARRIER
@@ -97,5 +97,5 @@
   
     REAL(rstd), POINTER :: ps(:)
-    REAL(rstd), POINTER :: theta_rhodz(:,:)
+    REAL(rstd), POINTER :: theta_rhodz(:,:,:)
     REAL(rstd), POINTER :: temp(:,:)
     REAL(rstd), POINTER :: p(:,:)
@@ -120,5 +120,5 @@
       CALL compute_exner(ps,p,pks,pk,0)
 !$OMP BARRIER
-      CALL compute_temperature2theta_rhodz(p, pk, temp, theta_rhodz, 0)
+      CALL compute_temperature2theta_rhodz(p, pk, temp, theta_rhodz(:,:,1), 0)
     ENDDO
 !$OMP BARRIER
Index: codes/icosagcm/trunk/src/timeloop_gcm.f90
===================================================================
--- codes/icosagcm/trunk/src/timeloop_gcm.f90	(revision 386)
+++ codes/icosagcm/trunk/src/timeloop_gcm.f90	(revision 387)
@@ -41,4 +41,6 @@
        STOP
     END IF
+
+    nqdyn = 1 ! one dynamical tracer = theta for the moment
 
     def='ARK2.3'
@@ -87,5 +89,5 @@
     CALL allocate_field(f_mass,field_t,type_real,llm,name='mass')
     CALL allocate_field(f_rhodz,field_t,type_real,llm,name='rhodz')
-    CALL allocate_field(f_theta_rhodz,field_t,type_real,llm,name='theta_rhodz')
+    CALL allocate_field(f_theta_rhodz,field_t,type_real,llm,nqdyn,name='theta_rhodz')
     CALL allocate_field(f_u,field_u,type_real,llm,name='u')
     CALL allocate_field(f_geopot,field_t,type_real,llm+1,name='geopot')
@@ -103,10 +105,10 @@
        CALL allocate_field(f_dps,field_t,type_real,name='dps')
        CALL allocate_field(f_dmass,field_t,type_real,llm, name='dmass')
-       CALL allocate_field(f_dtheta_rhodz,field_t,type_real,llm,name='dtheta_rhodz')
+       CALL allocate_field(f_dtheta_rhodz,field_t,type_real,llm,nqdyn,name='dtheta_rhodz')
        CALL allocate_field(f_du,field_u,type_real,llm,name='du')
        ! Model state at previous time step (RK/MLF)
        CALL allocate_field(f_psm1,field_t,type_real,name='psm1')
        CALL allocate_field(f_massm1,field_t,type_real,llm, name='massm1')
-       CALL allocate_field(f_theta_rhodzm1,field_t,type_real,llm,name='theta_rhodzm1')
+       CALL allocate_field(f_theta_rhodzm1,field_t,type_real,llm,nqdyn,name='theta_rhodzm1')
        CALL allocate_field(f_um1,field_u,type_real,llm,name='um1')
     CASE(hevi)
@@ -114,5 +116,5 @@
        CALL allocate_fields(nb_stage,f_dps_slow, field_t,type_real,name='dps_slow')
        CALL allocate_fields(nb_stage,f_dmass_slow, field_t,type_real,llm, name='dmass_slow')
-       CALL allocate_fields(nb_stage,f_dtheta_rhodz_slow, field_t,type_real,llm,name='dtheta_rhodz_fast')
+       CALL allocate_fields(nb_stage,f_dtheta_rhodz_slow, field_t,type_real,llm,nqdyn,name='dtheta_rhodz_fast')
        CALL allocate_fields(nb_stage,f_du_slow, field_u,type_real,llm,name='du_slow')
        CALL allocate_fields(nb_stage,f_du_fast, field_u,type_real,llm,name='du_fast')
@@ -131,5 +133,5 @@
        CALL allocate_field(f_psm2,field_t,type_real)
        CALL allocate_field(f_massm2,field_t,type_real,llm)
-       CALL allocate_field(f_theta_rhodzm2,field_t,type_real,llm)
+       CALL allocate_field(f_theta_rhodzm2,field_t,type_real,llm,nqdyn)
        CALL allocate_field(f_um2,field_u,type_real,llm)
     END SELECT
@@ -208,5 +210,5 @@
     fluxt_zero=.TRUE.
 
-    IF(positive_theta) CALL copy_theta_to_q1(f_theta_rhodz,f_rhodz,f_q)
+    IF(positive_theta) CALL copy_theta_to_q(f_theta_rhodz,f_rhodz,f_q)
 
     !$OMP MASTER
@@ -303,5 +305,5 @@
              END DO
           ENDIF
-          IF(positive_theta) CALL copy_q1_to_theta(f_theta_rhodz,f_rhodz,f_q)
+          IF(positive_theta) CALL copy_q_to_theta(f_theta_rhodz,f_rhodz,f_q)
        END IF
        
@@ -360,8 +362,8 @@
   END SUBROUTINE print_iteration
 
-  SUBROUTINE copy_theta_to_q1(f_theta_rhodz,f_rhodz,f_q)
+  SUBROUTINE copy_theta_to_q(f_theta_rhodz,f_rhodz,f_q)
     TYPE(t_field),POINTER :: f_theta_rhodz(:),f_rhodz(:), f_q(:)
-    REAL(rstd), POINTER :: theta_rhodz(:,:), rhodz(:,:), q(:,:,:)
-    INTEGER :: ind
+    REAL(rstd), POINTER :: theta_rhodz(:,:,:), rhodz(:,:), q(:,:,:)
+    INTEGER :: ind, iq
     DO ind=1, ndomain
        IF (.NOT. assigned_domain(ind)) CYCLE
@@ -371,12 +373,14 @@
        rhodz=f_rhodz(ind)
        q=f_q(ind)
-       q(:,:,1)  = theta_rhodz(:,:)/rhodz(:,:)
+       DO iq=1, nqdyn
+          q(:,:,iq)  = theta_rhodz(:,:,iq)/rhodz(:,:)
+       END DO
     END DO
-  END SUBROUTINE copy_theta_to_q1
-
-  SUBROUTINE copy_q1_to_theta(f_theta_rhodz,f_rhodz,f_q)
+  END SUBROUTINE copy_theta_to_q
+
+  SUBROUTINE copy_q_to_theta(f_theta_rhodz,f_rhodz,f_q)
     TYPE(t_field),POINTER :: f_theta_rhodz(:),f_rhodz(:), f_q(:)
-    REAL(rstd), POINTER :: theta_rhodz(:,:), rhodz(:,:), q(:,:,:)
-    INTEGER :: ind
+    REAL(rstd), POINTER :: theta_rhodz(:,:,:), rhodz(:,:), q(:,:,:)
+    INTEGER :: ind, iq
     DO ind=1, ndomain
        IF (.NOT. assigned_domain(ind)) CYCLE
@@ -386,7 +390,9 @@
        rhodz=f_rhodz(ind)
        q=f_q(ind)
-       theta_rhodz(:,:) = rhodz(:,:)*q(:,:,1)
+       DO iq=1,nqdyn
+          theta_rhodz(:,:,iq) = rhodz(:,:)*q(:,:,iq)
+       END DO
     END DO
-  END SUBROUTINE copy_q1_to_theta
+  END SUBROUTINE copy_q_to_theta
 
 END MODULE timeloop_gcm_mod
