Index: codes/icosagcm/trunk/src/caldyn_gcm.f90
===================================================================
--- codes/icosagcm/trunk/src/caldyn_gcm.f90	(revision 353)
+++ codes/icosagcm/trunk/src/caldyn_gcm.f90	(revision 354)
@@ -10,9 +10,5 @@
   !$OMP THREADPRIVATE(out_u, p, qu)
 
-  TYPE(t_field),POINTER :: f_buf_i(:), f_buf_ulon(:), f_buf_ulat(:), f_buf_u3d(:)
-  TYPE(t_field),POINTER :: f_buf_v(:), f_buf_s(:), f_buf_p(:)
-
 ! temporary shared variable for caldyn
-  TYPE(t_field),POINTER :: f_theta(:)
   TYPE(t_field),POINTER :: f_pk(:)
   TYPE(t_field),POINTER :: f_wwuu(:)
@@ -24,6 +20,5 @@
   TYPE(t_message) :: req_ps, req_mass, req_theta_rhodz, req_u, req_qu
 
-  PUBLIC init_caldyn, caldyn_BC, caldyn, write_output_fields, &
-       req_ps, req_mass
+  PUBLIC init_caldyn, caldyn_BC, caldyn, req_ps, req_mass
 
 CONTAINS
@@ -31,5 +26,6 @@
   SUBROUTINE init_caldyn
     USE icosa
-    USE exner_mod
+    USE observable_mod
+!    USE exner_mod
     USE mpipara
     USE omp_para
@@ -72,14 +68,4 @@
     CALL allocate_field(f_qu,field_u,type_real,llm) 
     CALL allocate_field(f_qv,field_z,type_real,llm) 
-  
-    CALL allocate_field(f_buf_i,   field_t,type_real,llm,name="buffer_i")
-    CALL allocate_field(f_buf_p,   field_t,type_real,llm+1) 
-    CALL allocate_field(f_buf_u3d, field_t,type_real,3,llm)  ! 3D vel at cell centers
-    CALL allocate_field(f_buf_ulon,field_t,type_real,llm)
-    CALL allocate_field(f_buf_ulat,field_t,type_real,llm)
-    CALL allocate_field(f_buf_v,   field_z,type_real,llm)
-    CALL allocate_field(f_buf_s,   field_t,type_real)
-
-    CALL allocate_field(f_theta, field_t,type_real,llm,  name='theta')   ! potential temperature
     CALL allocate_field(f_pk,    field_t,type_real,llm,  name='pk')
     CALL allocate_field(f_wwuu,  field_u,type_real,llm+1,name='wwuu')
@@ -134,4 +120,5 @@
        f_geopot, f_hflux, f_wflux, f_dps, f_dmass, f_dtheta_rhodz, f_du)
     USE icosa
+    USE observable_mod
     USE disvert_mod, ONLY : caldyn_eta, eta_mass
     USE vorticity_mod
@@ -294,25 +281,4 @@
 
 !$OMP BARRIER
-    IF (write_out) THEN
-
-       IF (is_master) PRINT *,'CALL write_output_fields'
-
-! ---> for openMP test to fix later
-!       CALL write_output_fields(f_ps, f_phis, f_dps, f_u, f_theta_rhodz, f_q, &
-!            f_buf_i, f_buf_v, f_buf_u3d, f_buf_ulon, f_buf_ulat, f_buf_s, f_buf_p)
-       CALL un2ulonlat(f_u, 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("dps",f_dps)
-!       CALL output_field("mass",f_mass)
-!       CALL output_field("dmass",f_dmass)
-       CALL output_field("vort",f_qv)
-       CALL output_field("theta",f_theta)
-!       CALL output_field("exner",f_pk)
-!       CALL output_field("pv",f_qv)
- 
-    END IF
-    
     !    CALL check_mass_conservation(f_ps,f_dps)
     CALL trace_end("caldyn")
@@ -946,135 +912,3 @@
   END SUBROUTINE check_mass_conservation  
   
-  SUBROUTINE write_output_fields(f_ps, f_phis, f_dps, f_u, f_theta_rhodz, f_q, &
-       f_buf_i, f_buf_v, f_buf_i3, f_buf1_i, f_buf2_i, f_buf_s, f_buf_p)
-    USE icosa
-    USE vorticity_mod
-    USE theta2theta_rhodz_mod
-    USE pression_mod
-    USE omega_mod
-    USE write_field_mod
-    USE vertical_interp_mod
-    USE wind_mod
-    TYPE(t_field),POINTER :: f_ps(:), f_phis(:), f_u(:), f_theta_rhodz(:), f_q(:), f_dps(:), &
-         f_buf_i(:), f_buf_v(:), f_buf_i3(:), f_buf1_i(:), f_buf2_i(:), f_buf_s(:), f_buf_p(:)
-    
-    REAL(rstd) :: out_pression_level
-    CHARACTER(LEN=255) :: str_pression
-    CHARACTER(LEN=255) :: physics_type
-    
-    out_pression_level=0.
-    CALL getin("out_pression_level",out_pression_level) 
-    WRITE(str_pression,*) INT(out_pression_level/100)
-    str_pression=ADJUSTL(str_pression)
-    
-    CALL writefield("ps",f_ps)
-    CALL writefield("dps",f_dps)
-    CALL writefield("phis",f_phis)
-    CALL vorticity(f_u,f_buf_v)
-    CALL writefield("vort",f_buf_v)
-
-    CALL w_omega(f_ps, f_u, f_buf_i)
-    CALL writefield('omega', f_buf_i)
-    IF (out_pression_level<=preff .AND. out_pression_level > 0) THEN
-      CALL vertical_interp(f_ps,f_buf_i,f_buf_s,out_pression_level)
-      CALL writefield("omega"//TRIM(str_pression),f_buf_s)
-    ENDIF
-    
-    ! Temperature
-!    CALL theta_rhodz2temperature(f_ps,f_theta_rhodz,f_buf_i) ; ! FIXME
-     
-    CALL getin('physics',physics_type)
-    IF (TRIM(physics_type)=='dcmip') THEN
-      CALL Tv2T(f_buf_i,f_q,f_buf1_i) 
-      CALL writefield("T",f_buf1_i)
-      IF (out_pression_level<=preff .AND. out_pression_level > 0) THEN
-        CALL vertical_interp(f_ps,f_buf1_i,f_buf_s,out_pression_level)
-        CALL writefield("T"//TRIM(str_pression),f_buf_s)
-      ENDIF
-    ELSE
-      CALL writefield("T",f_buf_i)
-      IF (out_pression_level<=preff .AND. out_pression_level > 0) THEN
-        CALL vertical_interp(f_ps,f_buf_i,f_buf_s,out_pression_level)
-        CALL writefield("T"//TRIM(str_pression),f_buf_s)
-      ENDIF
-    ENDIF
-   
-    ! velocity components
-    CALL un2ulonlat(f_u, f_buf1_i, f_buf2_i)
-    CALL writefield("ulon",f_buf1_i)
-    CALL writefield("ulat",f_buf2_i)
-
-    IF (out_pression_level<=preff .AND. out_pression_level > 0) THEN
-      CALL vertical_interp(f_ps,f_buf1_i,f_buf_s,out_pression_level)
-      CALL writefield("ulon"//TRIM(str_pression),f_buf_s)
-      CALL vertical_interp(f_ps,f_buf2_i,f_buf_s,out_pression_level)
-      CALL writefield("ulat"//TRIM(str_pression),f_buf_s)
-    ENDIF
-    
-    ! geopotential ! FIXME
-    CALL thetarhodz2geopot(f_ps,f_phis,f_theta_rhodz, f_buf_s,f_buf_p,f_buf1_i,f_buf2_i,f_buf_i)
-    CALL writefield("p",f_buf_p)
-!    CALL writefield("phi",f_geopot)   ! geopotential
-    CALL writefield("theta",f_buf1_i) ! potential temperature
-    CALL writefield("pk",f_buf2_i)    ! Exner pressure
-  
-  END SUBROUTINE write_output_fields
-  
-  SUBROUTINE thetarhodz2geopot(f_ps,f_phis,f_theta_rhodz, f_pks,f_p,f_theta,f_pk,f_phi) 
-    USE field_mod
-    USE pression_mod
-    USE exner_mod
-    USE geopotential_mod
-    USE theta2theta_rhodz_mod
-    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(:,:), &
-         phi(:,:), phis(:), ps(:), pks(:)
-    INTEGER :: ind
-
-    DO ind=1,ndomain
-       IF (.NOT. assigned_domain(ind)) CYCLE
-       CALL swap_dimensions(ind)
-       CALL swap_geometry(ind)
-       ps = f_ps(ind)
-       p  = f_p(ind)
-!$OMP BARRIER
-       CALL compute_pression(ps,p,0)
-       pk = f_pk(ind)
-       pks = f_pks(ind)
-!$OMP BARRIER
-       CALL compute_exner(ps,p,pks,pk,0)
-!$OMP BARRIER
-       theta_rhodz = f_theta_rhodz(ind)
-       theta = f_theta(ind)
-       CALL compute_theta_rhodz2theta(ps, theta_rhodz,theta,0)
-       phis = f_phis(ind)
-       phi = f_phi(ind)
-       CALL compute_geopotential(phis,pks,pk,theta,phi,0)
-    END DO
-
-  END SUBROUTINE thetarhodz2geopot
-  
-  SUBROUTINE Tv2T(f_Tv, f_q, f_T)
-  USE icosa
-  IMPLICIT NONE
-    TYPE(t_field), POINTER :: f_TV(:)
-    TYPE(t_field), POINTER :: f_q(:)
-    TYPE(t_field), POINTER :: f_T(:)
-    
-    REAL(rstd),POINTER :: Tv(:,:), q(:,:,:), T(:,:)
-    INTEGER :: ind
-    
-    DO ind=1,ndomain
-       IF (.NOT. assigned_domain(ind)) CYCLE
-       CALL swap_dimensions(ind)
-       CALL swap_geometry(ind)
-       Tv=f_Tv(ind)
-       q=f_q(ind)
-       T=f_T(ind)
-       T=Tv/(1+0.608*q(:,:,1))
-    END DO
-    
-  END SUBROUTINE Tv2T
-  
 END MODULE caldyn_gcm_mod
Index: codes/icosagcm/trunk/src/geopotential_mod.f90
===================================================================
--- codes/icosagcm/trunk/src/geopotential_mod.f90	(revision 353)
+++ codes/icosagcm/trunk/src/geopotential_mod.f90	(revision 354)
@@ -1,9 +1,11 @@
 MODULE geopotential_mod
+  IMPLICIT NONE
+  PRIVATE
 
+  PUBLIC :: compute_geopotential
 CONTAINS
 
-  SUBROUTINE geopotential(f_phis,f_pks,f_pk,f_theta,f_phi)
+  SUBROUTINE geopotential(f_phis,f_pks,f_pk,f_theta,f_phi) ! ORPHAN
   USE icosa
-  IMPLICIT NONE
     TYPE(t_field), POINTER :: f_phis(:)
     TYPE(t_field), POINTER :: f_pks(:)
Index: codes/icosagcm/trunk/src/observable.f90
===================================================================
--- codes/icosagcm/trunk/src/observable.f90	(revision 354)
+++ codes/icosagcm/trunk/src/observable.f90	(revision 354)
@@ -0,0 +1,179 @@
+MODULE observable_mod
+  USE icosa
+  IMPLICIT NONE
+  PRIVATE
+
+  TYPE(t_field),POINTER, SAVE :: f_buf_i(:), f_buf_ulon(:), f_buf_ulat(:), f_buf_u3d(:)
+  TYPE(t_field),POINTER, SAVE :: f_buf1_i(:), f_buf2_i(:)
+  TYPE(t_field),POINTER, SAVE :: f_buf_v(:), f_buf_s(:), f_buf_p(:)
+
+! temporary shared variable for caldyn
+  TYPE(t_field),POINTER, SAVE :: f_theta(:)
+
+  PUBLIC init_observable, write_output_fields_basic, f_theta
+
+CONTAINS
+  
+  SUBROUTINE init_observable
+    CALL allocate_field(f_buf_i,   field_t,type_real,llm,name="buffer_i")
+    CALL allocate_field(f_buf_p,   field_t,type_real,llm+1) 
+    CALL allocate_field(f_buf_u3d, field_t,type_real,3,llm)  ! 3D vel at cell centers
+    CALL allocate_field(f_buf_ulon,field_t,type_real,llm, name="buf_ulon")
+    CALL allocate_field(f_buf_ulat,field_t,type_real,llm, name="buf_ulat")
+    CALL allocate_field(f_buf_v,   field_z,type_real,llm)
+    CALL allocate_field(f_buf_s,   field_t,type_real)
+
+    CALL allocate_field(f_theta, field_t,type_real,llm,  name='theta')   ! potential temperature
+  END SUBROUTINE init_observable
+  
+  SUBROUTINE write_output_fields_basic(f_ps, f_u, f_q)
+    USE wind_mod
+    USE output_field_mod
+    USE omp_para
+    TYPE(t_field),POINTER :: f_ps(:), f_u(:), f_q(:)
+    IF (is_master) PRINT *,'CALL write_output_fields_basic'
+    CALL un2ulonlat(f_u, 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("dps",f_dps)
+    !CALL output_field("mass",f_mass)
+    !       CALL output_field("dmass",f_dmass)
+    !       CALL output_field("vort",f_qv)
+    CALL output_field("theta",f_theta)
+    !       CALL output_field("exner",f_pk)
+    !       CALL output_field("pv",f_qv)
+    CALL output_field("q",f_q)
+  END SUBROUTINE write_output_fields_basic
+  
+  SUBROUTINE write_output_fields(f_ps, f_phis, f_dps, f_u, f_theta_rhodz, f_q, &
+       f_buf_i, f_buf_v, f_buf_i3, f_buf1_i, f_buf2_i, f_buf_s, f_buf_p)
+    USE icosa
+    USE vorticity_mod
+    USE theta2theta_rhodz_mod
+    USE pression_mod
+    USE omega_mod
+    USE write_field_mod
+    USE vertical_interp_mod
+    USE wind_mod
+    TYPE(t_field),POINTER :: f_ps(:), f_phis(:), f_u(:), f_theta_rhodz(:), f_q(:), f_dps(:), &
+         f_buf_i(:), f_buf_v(:), f_buf_i3(:), f_buf1_i(:), f_buf2_i(:), f_buf_s(:), f_buf_p(:)
+    
+    REAL(rstd) :: out_pression_level
+    CHARACTER(LEN=255) :: str_pression
+    CHARACTER(LEN=255) :: physics_type
+    
+    out_pression_level=0.
+    CALL getin("out_pression_level",out_pression_level) 
+    WRITE(str_pression,*) INT(out_pression_level/100)
+    str_pression=ADJUSTL(str_pression)
+    
+    CALL writefield("ps",f_ps)
+    CALL writefield("dps",f_dps)
+    CALL writefield("phis",f_phis)
+    CALL vorticity(f_u,f_buf_v)
+    CALL writefield("vort",f_buf_v)
+
+    CALL w_omega(f_ps, f_u, f_buf_i)
+    CALL writefield('omega', f_buf_i)
+    IF (out_pression_level<=preff .AND. out_pression_level > 0) THEN
+      CALL vertical_interp(f_ps,f_buf_i,f_buf_s,out_pression_level)
+      CALL writefield("omega"//TRIM(str_pression),f_buf_s)
+    ENDIF
+    
+    ! Temperature
+!    CALL theta_rhodz2temperature(f_ps,f_theta_rhodz,f_buf_i) ; ! FIXME
+     
+    CALL getin('physics',physics_type)
+    IF (TRIM(physics_type)=='dcmip') THEN
+      CALL Tv2T(f_buf_i,f_q,f_buf1_i) 
+      CALL writefield("T",f_buf1_i)
+      IF (out_pression_level<=preff .AND. out_pression_level > 0) THEN
+        CALL vertical_interp(f_ps,f_buf1_i,f_buf_s,out_pression_level)
+        CALL writefield("T"//TRIM(str_pression),f_buf_s)
+      ENDIF
+    ELSE
+      CALL writefield("T",f_buf_i)
+      IF (out_pression_level<=preff .AND. out_pression_level > 0) THEN
+        CALL vertical_interp(f_ps,f_buf_i,f_buf_s,out_pression_level)
+        CALL writefield("T"//TRIM(str_pression),f_buf_s)
+      ENDIF
+    ENDIF
+   
+    ! velocity components
+    CALL un2ulonlat(f_u, f_buf1_i, f_buf2_i)
+    CALL writefield("ulon",f_buf1_i)
+    CALL writefield("ulat",f_buf2_i)
+
+    IF (out_pression_level<=preff .AND. out_pression_level > 0) THEN
+      CALL vertical_interp(f_ps,f_buf1_i,f_buf_s,out_pression_level)
+      CALL writefield("ulon"//TRIM(str_pression),f_buf_s)
+      CALL vertical_interp(f_ps,f_buf2_i,f_buf_s,out_pression_level)
+      CALL writefield("ulat"//TRIM(str_pression),f_buf_s)
+    ENDIF
+    
+    ! geopotential ! FIXME
+    CALL thetarhodz2geopot(f_ps,f_phis,f_theta_rhodz, f_buf_s,f_buf_p,f_buf1_i,f_buf2_i,f_buf_i)
+    CALL writefield("p",f_buf_p)
+!    CALL writefield("phi",f_geopot)   ! geopotential
+    CALL writefield("theta",f_buf1_i) ! potential temperature
+    CALL writefield("pk",f_buf2_i)    ! Exner pressure
+  
+  END SUBROUTINE write_output_fields
+  
+  SUBROUTINE thetarhodz2geopot(f_ps,f_phis,f_theta_rhodz, f_pks,f_p,f_theta,f_pk,f_phi) 
+    USE field_mod
+    USE pression_mod
+    USE exner_mod
+    USE geopotential_mod
+    USE theta2theta_rhodz_mod
+    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(:,:), &
+         phi(:,:), phis(:), ps(:), pks(:)
+    INTEGER :: ind
+    
+    DO ind=1,ndomain
+       IF (.NOT. assigned_domain(ind)) CYCLE
+       CALL swap_dimensions(ind)
+       CALL swap_geometry(ind)
+       ps = f_ps(ind)
+       p  = f_p(ind)
+!$OMP BARRIER
+       CALL compute_pression(ps,p,0)
+       pk = f_pk(ind)
+       pks = f_pks(ind)
+!$OMP BARRIER
+       CALL compute_exner(ps,p,pks,pk,0)
+!$OMP BARRIER
+       theta_rhodz = f_theta_rhodz(ind)
+       theta = f_theta(ind)
+       CALL compute_theta_rhodz2theta(ps, theta_rhodz,theta,0)
+       phis = f_phis(ind)
+       phi = f_phi(ind)
+       CALL compute_geopotential(phis,pks,pk,theta,phi,0)
+    END DO
+
+  END SUBROUTINE thetarhodz2geopot
+  
+  SUBROUTINE Tv2T(f_Tv, f_q, f_T)
+    TYPE(t_field), POINTER :: f_TV(:)
+    TYPE(t_field), POINTER :: f_q(:)
+    TYPE(t_field), POINTER :: f_T(:)
+    
+    REAL(rstd),POINTER :: Tv(:,:), q(:,:,:), T(:,:)
+    INTEGER :: ind
+    
+    DO ind=1,ndomain
+       IF (.NOT. assigned_domain(ind)) CYCLE
+       CALL swap_dimensions(ind)
+       CALL swap_geometry(ind)
+       Tv=f_Tv(ind)
+       q=f_q(ind)
+       T=f_T(ind)
+       T=Tv/(1+0.608*q(:,:,1))
+    END DO
+    
+  END SUBROUTINE Tv2T
+  
+END MODULE observable_mod
Index: codes/icosagcm/trunk/src/timeloop_gcm.f90
===================================================================
--- codes/icosagcm/trunk/src/timeloop_gcm.f90	(revision 353)
+++ codes/icosagcm/trunk/src/timeloop_gcm.f90	(revision 354)
@@ -30,4 +30,5 @@
   USE icosa
   USE dissip_gcm_mod
+  USE observable_mod
   USE caldyn_mod
   USE etat0_mod
@@ -131,4 +132,5 @@
     CALL init_dissip
     CALL init_sponge
+    CALL init_observable
     CALL init_caldyn
     CALL init_guided
@@ -157,4 +159,5 @@
   USE caldyn_mod
   USE caldyn_gcm_mod, ONLY : req_ps, req_mass
+  USE observable_mod
   USE etat0_mod
   USE guided_mod
@@ -242,5 +245,6 @@
     IF (mod(it,itau_out)==0 ) THEN
       CALL update_time_counter(dt*it)
-      CALL output_field("q",f_q)
+      CALL write_output_fields_basic(f_ps, f_u, f_q)
+!      CALL output_field("q",f_q)
     ENDIF
     
