Index: codes/icosagcm/trunk/src/observable.f90
===================================================================
--- codes/icosagcm/trunk/src/observable.f90	(revision 404)
+++ codes/icosagcm/trunk/src/observable.f90	(revision 413)
@@ -16,5 +16,5 @@
 
   PUBLIC init_observable, write_output_fields_basic, f_theta
-  LOGICAL,SAVE :: first_output=.TRUE.
+
 !$OMP THREADPRIVATE(first_output)
   
@@ -36,6 +36,8 @@
     CALL allocate_field(f_pmid,  field_t,type_real,llm,  name='pmid')       ! mid layer pressure
   END SUBROUTINE init_observable
-  
-  SUBROUTINE write_output_fields_basic(f_phis, f_ps, f_mass, f_geopot, f_theta_rhodz, f_u, f_W, f_q)
+
+  SUBROUTINE write_output_fields_basic(init, f_phis, f_ps, f_mass, f_geopot, f_theta_rhodz, f_u, f_W, f_q)
+    USE xios
+    USE disvert_mod
     USE wind_mod
     USE output_field_mod
@@ -43,99 +45,115 @@
     USE time_mod
     USE xios
-    USE disvert_mod
     USE earth_const
     USE pression_mod
     USE vertical_interp_mod
     USE theta2theta_rhodz_mod
-    USE wind_mod
     USE omega_mod
-    
-    TYPE(t_field),POINTER :: f_phis(:), f_ps(:), f_mass(:), f_geopot(:), f_theta_rhodz(:), f_u(:), f_W(:), f_q(:)
-!    IF (is_master) PRINT *,'CALL write_output_fields_basic'
+    LOGICAL, INTENT(IN) :: init
+    INTEGER :: l
     REAL :: scalar(1)
     REAL :: mid_ap(llm)
     REAL :: mid_bp(llm)
-    INTEGER :: l
-
-    IF (first_output) THEN
-      scalar(1)=dt
-      CALL xios_send_field("timestep", scalar)
-      scalar(1)=preff
-      CALL xios_send_field("preff", scalar)
-      CALL xios_send_field("ap",ap)
-      CALL xios_send_field("bp",bp)
-      DO l=1,llm
-        mid_ap(l)=(ap(l)+ap(l+1))/2
-        mid_bp(l)=(bp(l)+bp(l+1))/2
-      ENDDO
-      CALL xios_send_field("mid_ap",mid_ap)
-      CALL xios_send_field("mid_bp",mid_bp)
-
-      CALL output_field("phis",f_phis)
-        
-      first_output=.FALSE.
-    ENDIF
+
+    TYPE(t_field),POINTER :: f_phis(:), f_ps(:), f_mass(:), f_geopot(:), f_theta_rhodz(:), f_u(:), f_W(:), f_q(:)
+!    IF (is_master) PRINT *,'CALL write_output_fields_basic'
+
+    IF(init) THEN
+       scalar(1)=dt
+       CALL xios_send_field("timestep", scalar)
+       scalar(1)=preff
+       CALL xios_send_field("preff", scalar)
+       CALL xios_send_field("ap",ap)
+       CALL xios_send_field("bp",bp)
+       DO l=1,llm
+          mid_ap(l)=(ap(l)+ap(l+1))/2
+          mid_bp(l)=(bp(l)+bp(l+1))/2
+       ENDDO
+       CALL xios_send_field("mid_ap",mid_ap)
+       CALL xios_send_field("mid_bp",mid_bp)
+
+       CALL output_field("phis",f_phis)
+       CALL output_field("Ai",geom%Ai)       
+    END IF
+
+    CALL divide_by_mass(1, f_mass, f_theta_rhodz, f_buf_i)
+    IF(init) THEN
+       CALL output_field("theta_init",f_buf_i)
+    ELSE
+       CALL output_field("theta",f_buf_i)
+    END IF
+
+    IF(nqdyn>1) THEN
+       CALL divide_by_mass(2, f_mass, f_theta_rhodz, f_buf_i)
+       IF(init) THEN
+          CALL output_field("dyn_q_init",f_buf_i)
+       ELSE
+          CALL output_field("dyn_q",f_buf_i)
+       END IF
+    END IF
+
+    CALL theta_rhodz2temperature(f_ps,f_theta_rhodz,f_buf_i) ; 
+    CALL Tv2T(f_buf_i,f_q,f_buf1_i)
+    IF(init) THEN
+       CALL output_field("temp_init",f_buf_i)
+    ELSE
+       CALL output_field("temp",f_buf_i)
+       CALL vertical_interp(f_ps,f_buf_i,f_buf_s,85000.)
+       CALL output_field("t850",f_buf_s)
+       CALL vertical_interp(f_ps,f_buf_i,f_buf_s,50000.)
+       CALL output_field("t500",f_buf_s)
+       CALL vertical_interp(f_ps,f_buf_i,f_buf_s,preff)
+       CALL output_field("SST",f_buf_s)       
+    END IF
     
     CALL progonostic_vel_to_horiz(f_geopot, f_ps, f_mass, f_u, f_W, f_buf_uh, f_buf_i)
     CALL transfert_request(f_buf_uh,req_e1_vect) 
-    CALL output_field("uz",f_buf_i)
-    CALL vertical_interp(f_ps,f_buf_i,f_buf_s,85000.)
-    CALL output_field("w850",f_buf_s)
-    CALL vertical_interp(f_ps,f_buf_i,f_buf_s,50000.)
-    CALL output_field("w500",f_buf_s)
-    
-    
     CALL un2ulonlat(f_buf_uh, f_buf_ulon, f_buf_ulat)
-    CALL output_field("ulon",f_buf_ulon)
-    CALL output_field("ulat",f_buf_ulat)
-    CALL output_field("ps",f_ps)
-    CALL output_field("Ai",geom%Ai)
-
-    !       CALL output_field("dps",f_dps)
-    CALL output_field("mass",f_mass)
-    CALL output_field("geopot",f_geopot)
-    !       CALL output_field("dmass",f_dmass)
-    !       CALL output_field("vort",f_qv)
-    
-    
-    CALL theta_rhodz2temperature(f_ps,f_theta_rhodz,f_buf_i) ; 
-    CALL Tv2T(f_buf_i,f_q,f_buf1_i) 
-    CALL output_field("temp",f_buf_i)
-    CALL vertical_interp(f_ps,f_buf_i,f_buf_s,85000.)
-    CALL output_field("t850",f_buf_s)
-    CALL vertical_interp(f_ps,f_buf_i,f_buf_s,50000.)
-    CALL output_field("t500",f_buf_s)
-    CALL vertical_interp(f_ps,f_buf_i,f_buf_s,preff)
-    CALL output_field("SST",f_buf_s)
-
-
-    CALL extract_slice(f_theta_rhodz, f_buf_i,1)
-    CALL output_field("theta",f_buf_i)
-           
-    !       CALL output_field("exner",f_pk)
-    !       CALL output_field("pv",f_qv)
-    CALL output_field("q",f_q)
     CALL pression_mid(f_ps, f_pmid)
-    CALL output_field("p",f_pmid)
-
-    CALL vertical_interp(f_ps,f_buf_ulon,f_buf_s,85000.)
-    CALL output_field("u850",f_buf_s)
-    CALL vertical_interp(f_ps,f_buf_ulon,f_buf_s,50000.)
-    CALL output_field("u500",f_buf_s)
-
-    CALL vertical_interp(f_ps,f_buf_ulat,f_buf_s,85000.)
-    CALL output_field("v850",f_buf_s)
-    CALL vertical_interp(f_ps,f_buf_ulat,f_buf_s,50000.)
-    CALL output_field("v500",f_buf_s)
-    
-    CALL w_omega(f_ps, f_u, f_buf_i)
-    CALL output_field("omega",f_buf_i)
-    CALL vertical_interp(f_ps,f_buf_i,f_buf_s,85000.)
-    CALL output_field("omega850",f_buf_s)
-    CALL vertical_interp(f_ps,f_buf_i,f_buf_s,50000.)
-    CALL output_field("omega500",f_buf_s)
-
-
-    
+    IF(init) THEN
+       CALL output_field("uz_init",f_buf_i)
+       CALL output_field("ulon_init",f_buf_ulon)
+       CALL output_field("ulat_init",f_buf_ulat)
+       CALL output_field("p_init",f_pmid)
+       CALL output_field("ps_init",f_ps)
+       CALL output_field("mass_init",f_mass)
+       CALL output_field("geopot_init",f_geopot)
+       CALL output_field("q_init",f_q)
+    ELSE
+       CALL output_field("uz",f_buf_i)
+       CALL output_field("ulon",f_buf_ulon)
+       CALL output_field("ulat",f_buf_ulat)
+       CALL output_field("p",f_pmid)
+       CALL output_field("ps",f_ps)
+       CALL output_field("mass",f_mass)
+       CALL output_field("geopot",f_geopot)
+       CALL output_field("q",f_q)
+
+       !       CALL output_field("exner",f_pk)
+       !       CALL output_field("pv",f_qv)
+       
+       CALL vertical_interp(f_ps,f_buf_ulon,f_buf_s,85000.)
+       CALL output_field("u850",f_buf_s)
+       CALL vertical_interp(f_ps,f_buf_ulon,f_buf_s,50000.)
+       CALL output_field("u500",f_buf_s)
+       
+       CALL vertical_interp(f_ps,f_buf_ulat,f_buf_s,85000.)
+       CALL output_field("v850",f_buf_s)
+       CALL vertical_interp(f_ps,f_buf_ulat,f_buf_s,50000.)
+       CALL output_field("v500",f_buf_s)
+
+       CALL vertical_interp(f_ps,f_buf_i,f_buf_s,85000.)
+       CALL output_field("w850",f_buf_s)
+       CALL vertical_interp(f_ps,f_buf_i,f_buf_s,50000.)
+       CALL output_field("w500",f_buf_s)    
+
+       CALL w_omega(f_ps, f_u, f_buf_i)
+       CALL output_field("omega",f_buf_i)
+       CALL vertical_interp(f_ps,f_buf_i,f_buf_s,85000.)
+       CALL output_field("omega850",f_buf_s)
+       CALL vertical_interp(f_ps,f_buf_i,f_buf_s,50000.)
+       CALL output_field("omega500",f_buf_s)
+    END IF
+
   END SUBROUTINE write_output_fields_basic
   
@@ -344,4 +362,20 @@
     
   END SUBROUTINE Tv2T
-  
+
+  SUBROUTINE divide_by_mass(iq, f_mass, f_theta_rhodz, f_theta)
+    INTEGER, INTENT(IN) :: iq
+    TYPE(t_field), POINTER :: f_mass(:), f_theta_rhodz(:), f_theta(:)
+    REAL(rstd), POINTER :: mass(:,:), theta_rhodz(:,:,:), theta(:,:)
+    INTEGER :: ind
+    DO ind=1,ndomain
+       IF (.NOT. assigned_domain(ind)) CYCLE
+       CALL swap_dimensions(ind)
+       CALL swap_geometry(ind)
+       mass=f_mass(ind)
+       theta_rhodz=f_theta_rhodz(ind)
+       theta=f_theta(ind)
+       theta(:,:) = theta_rhodz(:,:,iq) / mass(:,:)
+    END DO
+  END SUBROUTINE divide_by_mass
+    
 END MODULE observable_mod
