Index: /trunk/LMDZ.PLUTO/deftank/xios/context_pcm_physics.xml
===================================================================
--- /trunk/LMDZ.PLUTO/deftank/xios/context_pcm_physics.xml	(revision 3508)
+++ /trunk/LMDZ.PLUTO/deftank/xios/context_pcm_physics.xml	(revision 3508)
@@ -0,0 +1,150 @@
+<context id="LMDZ">
+
+    <!-- Definition of a calendar is done in code via wxios_set_cal -->
+
+    <!-- Definition of domains and groups of domains -->
+    <domain_definition>
+        <domain_group id="dom_glo" data_dim="1">
+          <domain id="dom_glo" />
+        </domain_group>
+
+    </domain_definition>
+
+    <domain_definition>
+        <domain id="dom_64_48" type="rectilinear" ni_glo="64" nj_glo="48"   >
+          <generate_rectilinear_domain/>
+          <interpolate_domain order="1"/>
+        </domain>
+
+        <domain id="dom_96_95" type="rectilinear"  ni_glo="96" nj_glo="95"  >
+          <generate_rectilinear_domain/>
+          <interpolate_domain order="1"/>
+        </domain>
+
+        <domain id="dom_144_142" type="rectilinear" ni_glo="144" nj_glo="142"   >
+          <generate_rectilinear_domain/>
+          <interpolate_domain order="1"/>
+        </domain>
+
+        <domain id="dom_512_360" type="rectilinear" ni_glo="512" nj_glo="360"   >
+          <generate_rectilinear_domain/>
+          <interpolate_domain order="1"/>
+        </domain>
+
+        <domain id="dom_720_360" type="rectilinear" ni_glo="720" nj_glo="360" >
+          <generate_rectilinear_domain/>
+          <interpolate_domain order="1"/>
+        </domain>
+
+        <domain id="dom_128_96" type="rectilinear" ni_glo="128" nj_glo="96" >
+          <generate_rectilinear_domain/>
+          <interpolate_domain order="1"/>
+        </domain>
+
+        <domain id="dom_256_192" type="rectilinear" ni_glo="256" nj_glo="192" >
+          <generate_rectilinear_domain/>
+          <interpolate_domain order="1"/>
+        </domain>
+
+        <domain id="dom_out" domain_ref="dom_128_96"/>
+    </domain_definition>
+    <!-- <domain_definition src="./domain_def.xml"/> -->
+
+    <!-- Definition of vertical axes
+                           NB: these must be defined and set in code via
+         xios_output_mod -->
+    <axis_definition>
+        <axis id="presnivs"
+              standard_name="Pseudo-pressure of model vertical levels"
+              unit="Pa">
+        </axis>
+        <axis id="altitude"
+              standard_name="Pseudo-altitude of model vertical levels"
+              unit="km">
+        </axis>
+        <axis id="soil_layers"
+              standard_name="soil layer depth at intermediate level"
+              unit="m">
+        </axis>
+        <axis id="interlayer"
+              standard_name="Levels between layers (altitude+1)"
+              unit="km">
+        </axis>
+        <axis id="controle_axe"
+              standard_name="Controle"
+              n_glo="100">
+        </axis>
+        <!--Added for diagspecIR-->
+        <axis id="IR_Wavenumber"
+               standard_name="Band mid frequency in the infrared"
+               unit="cm-1">
+        </axis>
+        <!--Added for diagspecVI-->
+        <axis id="VI_Wavenumber"
+              standard_name="Band mid frequency in the visible"
+              unit=
+            "cm-1">
+        </axis>
+
+    </axis_definition>
+
+    <scalar_definition>
+           <scalar id="scalar" />
+    </scalar_definition>
+
+    <grid_definition>
+        <!-- toggle axis id below to change output vertical axis -->
+        <grid id="grid_3D">
+            <domain domain_ref="dom_glo" />
+            <axis axis_ref="altitude" />
+        </grid>
+        <grid id="grid_2D">
+            <domain domain_ref="dom_glo" />
+        </grid>
+        <!--Added for diagspecIR-->
+        <grid id='grid_IR_3D'>
+            <domain domain_ref="dom_glo" />
+            <axis axis_ref='IR_Wavenumber'/>
+        </grid>
+        <!-- Added for diagspecVI-->
+        <grid id="grid_VI_3D">
+          <domain domain_ref="dom_glo" />
+          <axis axis_ref="VI_Wavenumber"/>
+        </grid>
+        <!-- grid for interpolated lon-lat outputs -->
+        <!-- to output scalars -->
+        <grid id="scalarpoint" >
+            <scalar scalar_ref="scalar" />
+        </grid>
+        <!-- output grids -->
+        <grid id="grid_3D_out">
+            <domain domain_ref="dom_out" />
+            <axis axis_ref="altitude" />
+        </grid>
+        <grid id="grid_2D_out">
+            <domain domain_ref="dom_out" />
+        </grid>
+        <!--added for diagspecIR-->
+        <grid id="grid_IR_3D_out">
+            <domain domain_ref="dom_out"/>
+            <axis axis_ref="IR_Wavenumber"/>
+        </grid>
+        <grid id="grid_IR">
+            <axis axis_ref="IR_Wavenumber"/>
+        </grid>
+        <!--Added for diagspecVI-->
+        <grid id="grid_VI_3D_out">
+            <domain domain_ref="dom_out"/>
+            <axis axis_ref="VI_Wavenumber"/>
+        </axis>
+        <grid id="grid_VI">
+          <axis axis_ref="VI_Wavenumber"/>
+        </grid>
+    </grid_definition>
+
+    <!-- Definition of model variables -->
+    <field_definition src="./field_def_physics.xml" />
+
+    <!-- Definition of output files -->
+    <file_definition src="./file_def_physics.xml" />
+</context>
Index: /trunk/LMDZ.PLUTO/deftank/xios/field_def_physics.xml
===================================================================
--- /trunk/LMDZ.PLUTO/deftank/xios/field_def_physics.xml	(revision 3508)
+++ /trunk/LMDZ.PLUTO/deftank/xios/field_def_physics.xml	(revision 3508)
@@ -0,0 +1,362 @@
+<!-- ============================================================== -->
+<!-- field_def_physics.xml                                          -->
+<!-- Definition of all existing variables that can be outputed
+     from Generic PCM physics                                       -->
+<!-- ============================================================== -->
+
+    <field_definition prec="4"
+                      freq_op="1ts"
+                      enabled="true">
+
+        <!-- 0D variables -->
+        <field_group id="fields_0D" grid_ref="scalarpoint">
+            <field id="Ls"
+                   long_name="Solar longitude"
+                   unit="degrees" />
+            <field id="RA"
+                   long_name="Right Ascension"
+                   unit="degrees" />
+            <field id="Declin"
+                   long_name="Solar declination"
+                   unit="degrees" />
+            <field id="dist_star"
+                   long_name="Distance to star"
+                   unit="AU" />
+        </field_group>
+
+        <field_group id="fields_controle" grid_ref="controle_grid">
+            <field id="controle"
+                   long_name="controle" />
+        </field_group>
+
+        <field_group id="fields_altitude" grid_ref="altitude_grid">
+            <field id="aps"
+                   long_name="aps" />
+            <field id="bps"
+                   long_name="bps" />
+        </field_group>
+
+        <field_group id="fields_interlayer" grid_ref="interlayer_grid">
+            <field id="ap"
+                   long_name="ap" />
+            <field id="bp"
+                   long_name="bp" />
+        </field_group>
+
+        <!-- 1D variables (spectral)-->
+        <field_group id="IR_fields" grid_ref="grid_IR">
+              <field id="IR_Bandwidth"
+                     long_name="Bandwidth in the infrared"
+                     unit="cm-1"/>
+        </field_group>
+        <field_group id="VI_field" grid_ref="grid_VI">
+              <field id="VI_Bandwidth"
+                     long_name="Bandwith in the visible"
+                     unit="cm-1"/>
+       </field_group>
+        <!-- 2D variables -->
+        <field_group id="fields_2D" grid_ref="grid_2D">
+            <field id="area"
+                   long_name="Mesh area"
+                   unit="m2" />
+            <field id="phisfi"
+                   long_name="Surface geopotential"
+                   unit="m2.s-2" />
+            <field id="emis"
+                   long_name="Surface emissivity"
+                   unit="" />
+            <field id="ALB"
+                   long_name="Surface Albedo"
+                    unit="" />
+
+            <field id="ps"
+                   long_name="Surface Pressure"
+                   unit="Pa" />
+            <field id="tsurf"
+                   long_name="Surface Temperature"
+                   unit="K" />
+
+            <!-- Tracers -->
+            <field id="n2_surf"
+                   long_name="N2 ice thickness"
+                   unit="kg/m2" />
+            <field id="ch4_gas_surf"
+                   long_name="CH4 vapor at surface"
+                   unit="kg/m2" />
+            <field id="ch4_ice_surf"
+                   long_name="CH4 ice at surface"
+                   unit="kg/m2" />
+            <field id="co_gas_surf"
+                   long_name="co vapor at surface"
+                   unit="kg/m2" />
+            <field id="co_ice_surf"
+                   long_name="co ice at surface"
+                   unit="kg/m2" />
+            <field id="haze_surf"
+                   long_name="haze at surface"
+                   unit="kg/m2" />
+            <field id="prec_haze_surf"
+                   long_name="precursor haze at surface"
+                   unit="kg/m2" />
+
+            <field id="n2_col"
+                   long_name="N2 vapor column"
+                   unit="kg/m2" />
+            <field id="ch4_gas_col"
+                   long_name="CH4 vapor column"
+                   unit="kg/m2" />
+            <field id="ch4_ice_col"
+                   long_name="CH4 ice column"
+                   unit="kg/m2" />
+            <field id="co_gas_col"
+                   long_name="co vapor column"
+                   unit="kg/m2" />
+            <field id="co_ice_col"
+                   long_name="co ice column"
+                   unit="kg/m2" />
+            <field id="haze_col"
+                   long_name="haze column"
+                   unit="kg/m2" />
+            <field id="prec_haze_col"
+                   long_name="precursor haze column"
+                   unit="kg/m2" />
+
+            <field id="n2_iceflux"
+                   long_name="N2 flux"
+                   unit="kg/m2/s" />
+            <field id="ch4_iceflux"
+                   long_name="CH4 ice flux"
+                   unit="kg/m2/s" />
+            <field id="co_iceflux"
+                   long_name="co ice flux"
+                   unit="kg/m2/s" />
+
+            <field id="vmr_ch4"
+                   long_name="CH4 volume mixing ratio"
+                   unit="" />
+            <field id="vmr_co"
+                   long_name="CO volume mixing ratio"
+                   unit="" />
+
+            <field id="tau_col"
+                   long_name="Total aerosol optical depth"
+                   unit="opacity" />
+
+            <field id="zdqsdif_ch4_ice"
+                   long_name="zdqsdif ch4"
+                   unit="" />
+            <field id="zdqssed_ch4"
+                   long_name="zdqssed ch4"
+                   unit="" />
+            <field id="zdqssed_haze"
+                   long_name="zdqssed haze"
+                   unit="" />
+
+            <!-- Radiation -->
+            <field id="ASR"
+                   long_name="Absorbed stellar radiation"
+                    unit="W/m2" />
+            <field id="ISR"
+                   long_name="incoming solar radiation"
+                   unit="W/m2" />
+            <field id="OLR"
+                   long_name="outgoing longwave radiation"
+                   unit="W/m2" />
+            <field id="GND"
+                   long_name="Heat flux from ground"
+                   unit="W/m2" />
+            <field id="DYN"
+                   long_name="Dynamical heat input"
+                   unit="W/m2" />
+            <field id="fluxsurf_lw"
+                   long_name="Longwave radiation at the surface"
+                   unit="W.m-2" />
+            <field id="sensibFlux"
+                   long_name="Sensible heat flux"
+                   unit="W/m2" />
+            <field id="lscaledE"
+                   long_name="heat from largescale"
+                   unit="W/m2" />
+            <field id="madjdE"
+                   long_name="heat from moistadj"
+                   unit="W/m2" />
+            <field id="genericconddE"
+                   long_name="heat from generic condensation"
+                   unit="W/m2" />
+
+        </field_group>
+
+        <!-- 3D variables -->
+        <field_group id="fields_3D" grid_ref="grid_3D">
+
+            <field id="zzlay"
+                   long_name="Midlayer altitude"
+                   unit="m" />
+            <field id="zzlev"
+                   long_name="Interlayer altitude"
+                   unit="m" />
+            <field id="pphi"
+                   long_name="Geopotential"
+                   unit="m2.s-2" />
+
+            <field id="temperature"
+                   long_name="Atmospheric temperature"
+                   unit="K" />
+            <field id="teta"
+                   long_name="Potential temperature"
+                   unit="K" />
+            <field id="u"
+                   long_name="Eastward Zonal Wind"
+                   unit="m/s" />
+            <field id="v"
+                   long_name="Northward Meridional Wind"
+                   unit="m/s" />
+            <field id="w"
+                   long_name="Downward Vertical Wind"
+                   unit="m/s" />
+            <field id="omega"
+                   long_name="Downward omega velocity"
+                   unit="Pa/s" />
+            <field id="p"
+                   long_name="Atmospheric pressure"
+                   unit="Pa" />
+            <field id="detr"
+                   long_name="Detrainment"
+                   unit="kg/m2/s" />
+            <field id="entr"
+                   long_name="Entrainment"
+                   unit="kg/m2/s" />
+            <field id="w_plm"
+                   long_name="Plume vertical velocity"
+                   unit="m/s" />
+            <field id="q2"
+                   long_name="turbulent kinetic energy"
+                   unit="J.kg^-1" />
+            <field id="dEzradsw"
+                   long_name="radiative heating"
+                   unit="w.m^-2" />
+            <field id="dEzradlw"
+                   long_name="radiative heating"
+                   unit="w.m^-2" />
+            <field id="qsatatm"
+                   long_name="atm qsat"
+                   unit=" " />
+            <field id="dt_generic_condensation"
+                   long_name="heating from generic condensation"
+                   unit="K s-1" />
+            <field id="zdtsw"
+                   long_name="SW heating"
+                   unit="T s-1" />
+            <field id="zdtlw"
+                   long_name="LW heating"
+                   unit="T s-1" />
+            <field id="dtrad"
+                   long_name="radiative heating"
+                   unit="K s-1" />
+            <field id="zdtdyn"
+                   long_name="Dyn. heating"
+                   unit="T s-1" />
+            <field id="dtmoistadj"
+                   long_name="moist adj heating"
+                   unit="K s-1" />
+            <field id="dtdryadj"
+                   long_name="dry adj heating"
+                   unit="K s-1" />
+
+            <!--Tracers-->
+            <field id="n2"
+                   long_name="N2 mass mixing ratio"
+                   unit="kg/kg" />
+            <field id="ch4_gas"
+                   long_name="CH4 vapor mass mixing ratio"
+                   unit="kg/kg" />
+            <field id="ch4_ice"
+                   long_name="CH4 ice mass mixing ratio"
+                   unit="kg/kg" />
+            <field id="co_gas"
+                   long_name="co vapor mass mixing ratio"
+                   unit="kg/kg" />
+            <field id="co_ice"
+                   long_name="co ice mass mixing ratio"
+                   unit="kg/kg" />
+            <field id="haze"
+                   long_name="haze mass mixing ratio"
+                   unit="kg/kg" />
+            <field id="prec_haze"
+                   long_name="precursor haze mass mixing ratio"
+                   unit="kg/kg" />
+            <field id="zdtc"
+                   long_name="tendancy T cond N2"
+                   unit="K" />
+
+            <field id="zrho_ch4"
+                   long_name="zrho_ch4"
+                   unit="kg/m3" />
+            <field id="zrho_co"
+                   long_name="zrho_co"
+                   unit="kg/m3" />
+            <field id="zrho_haze"
+                   long_name="zrho_haze"
+                   unit="kg/m3" />
+            <field id="zdqcn2_ch4"
+                   long_name="zdq condn2 ch4"
+                   unit="" />
+            <field id="zdqdif_ch4"
+                   long_name="zdqdif_ch4"
+                   unit="" />
+            <field id="zdqadj_ch4"
+                   long_name="zdqadj_ch4"
+                   unit="" />
+            <field id="zdqsed_ch4"
+                   long_name="zdqsed_ch4"
+                   unit="" />
+            <field id="zdtch4cloud"
+                   long_name="ch4 cloud"
+                   unit="T s-1" />
+            <field id="zdqch4cloud"
+                   long_name="ch4 cloud"
+                   unit="T s-1" />
+            <field id="zdqrho_photprec"
+                   long_name="zdqrho_photprec"
+                   unit="kg/m3/s" />
+            <field id="zdqphot_prec"
+                   long_name="zdqphot_prec"
+                   unit="" />
+            <field id="zdqphot_ch4"
+                   long_name="zdqphot_ch4"
+                   unit="" />
+
+            <field id="zdqhaze_ch4"
+                   long_name="zdqhaze_ch4"
+                   unit="" />
+            <field id="zdqhaze_prec"
+                   long_name="zdqhaze_prec"
+                   unit="" />
+            <field id="zdqhaze_haze"
+                   long_name="zdqhaze_haze"
+                   unit="" />
+            <field id="zdqconv_prec"
+                   long_name="zdqconv_prec"
+                   unit="" />
+
+        </field_group>
+
+        <!--Spectral IR 3D variables-->
+        <field_group id="spectral_IR_3D" grid_ref="grid_IR_3D">
+              <field id="OLR3D"
+                     long_name="OLR(lon,lat,band)"
+                     unit="W/m^2/cm^-1" />
+        </field_group>
+
+        <!-- SPectral VI 3D variables-->
+        <field_group id="spectral_VI_3D" grid_ref="grid_VI_3D">
+              <field id="OSR3D"
+                     long_name="OSR(lon,lat,band"
+                     unit="W/m^2/cm^-1"/>
+              <field id="GSR3D"
+                     long_name="GSR(lon,lat,band)"
+                     unit="w/m^2/cm^-1"/>
+        </field_group>
+    </field_definition>
+
+
Index: /trunk/LMDZ.PLUTO/deftank/xios/file_def_physics.xml
===================================================================
--- /trunk/LMDZ.PLUTO/deftank/xios/file_def_physics.xml	(revision 3508)
+++ /trunk/LMDZ.PLUTO/deftank/xios/file_def_physics.xml	(revision 3508)
@@ -0,0 +1,123 @@
+<!-- =================================================================== -->
+<!-- file_def_physics.xml                                                -->
+<!-- Definition of output files for Generic PCM physics                  -->
+<!-- =================================================================== -->
+
+    <file_definition>
+
+        <!-- Instantaneous outputs; Here output_freq in physics time steps -->
+        <!--
+           1 ts == 1 physical ts
+           if physical ts is half a day,
+           10ts   below means every 5 days
+           2000ts below means every 1000 days
+          -->
+        <file id="histins"
+              name="Xhistins"
+              output_freq="1ts"
+              sync_freq="1ts"
+              type="one_file"
+              time_units="days"
+              enabled="true">
+
+
+            <!-- VARS 0D -->
+            <field_group operation="instant"
+                         grid_ref="scalarpoint"
+                         freq_op="1ts">
+                <field field_ref="Ls" />
+            </field_group>
+
+            <!--<field_group domain_ref="dom_glo"> -->
+            <!--<field_group grid_ref="grid_2D_out">-->
+            <!-- toggle domain_ref to output on native or lon-lat domain -->
+
+              <!-- VARS 2D -->
+              <field_group operation="instant" prec="4"
+                           grid_ref="grid_2D_out"
+                           freq_op="1ts">
+                  <field field_ref="area" operation="once" />
+                  <field field_ref="ps" enabled="true" />
+                  <field field_ref="ISR" enabled="true" />
+                  <field field_ref="OLR" enabled="true" />
+                  <field field_ref="ASR" enabled="true" />
+              </field_group>
+
+              <!-- VARS 3D -->
+              <field_group operation="instant"
+                           grid_ref="grid_3D_out"
+                           freq_op="1ts">
+                  <field field_ref="u" />
+                  <field field_ref="v" />
+                  <field field_ref="omega"/>
+                  <field field_ref="temperature" />
+                  <field field_ref="n2" />
+                  <field field_ref="ch4_gas" />
+                  <field field_ref="ch4_ice" />
+                  <field field_ref="co_gas" />
+                  <field field_ref="co_ice" />
+                  <field field_ref="haze" />
+
+                  <!--field field_ref="h2o_vap" />
+                  <field field_ref="h2o_ice" />
+                  <field field_ref="w_plm" />
+                  <field field_ref="detr" />
+                  <field field_ref="entr" /-->
+              </field_group>
+
+            <!--</field_group>-->
+        </file>
+       <!--XspecIR-->
+      <file id="specIR"
+            name="XspecIR"
+            output_freq="400ts"
+            type="one_file"
+            time_units="days"
+            enabled="false">
+          <!--Var 1D ?-->
+          <field_group operation="instant"
+                        grid_ref="grid_IR"
+                        freq_op="1ts">
+              <field field_ref="IR_Bandwidth" operation="once"/>
+          </field_group>
+          <!--Var 2D-->
+          <field_group operation="instant" prec="4"
+                        grid_ref="grid_2D_out"
+                        freq_op="1ts">
+              <field field_ref="area" operation="once" />
+          </field_group>
+          <!--Var 3D-->
+          <field_group operation="instant"
+                        grid_ref="grid_IR_3D_out"
+                        freq_op="1ts">
+              <field field_ref="OLR3D" />
+          </field_group>
+      </file>
+       <!--XspecVI-->
+      <file id="specVI"
+            name="XspecVI"
+            output_freq="400ts"
+            type="one_file"
+            time_units="days"
+            enabled="true">
+          <!--Var 1D -->
+          <field_group operation="instant"
+                        grid_ref="grid_VI"
+                        freq_op="1ts">
+              <field field_ref="VI_Bandwidth" operation="once"/>
+          </field_group>
+          <!--Var 2D-->
+          <field_group operation="instant" prec="4"
+                        grid_ref="grid_2D_out"
+                        freq_op="1ts">
+              <field field_ref="area" operation="once" />
+          </field_group>
+          <!--Var 3D-->
+          <field_group operation="instant"
+                        grid_ref="grid_VI_3D_out"
+                        freq_op="1ts">
+              <field field_ref="OSR3D" />
+              <field field_ref="GSR3D" />
+          </field_group>
+      </file>
+    </file_definition>
Index: /trunk/LMDZ.PLUTO/deftank/xios/iodef.xml
===================================================================
--- /trunk/LMDZ.PLUTO/deftank/xios/iodef.xml	(revision 3508)
+++ /trunk/LMDZ.PLUTO/deftank/xios/iodef.xml	(revision 3508)
@@ -0,0 +1,27 @@
+<?xml version="1.0"?>
+<simulation>
+        <context  id="xios">
+            <variable_definition>
+                <!-- set print_file to true to have XIOS write its
+                     outputs in separate files -->
+                <variable id="print_file" type="bool"> true </variable>
+
+                <variable id="buffer_size_factor" type="double">4</variable>
+                <variable_group id="buffer">
+                buffer_size = 85000000
+                buffer_server_factor_size = 2
+                </variable_group>
+
+                <variable_group id="parameters" >
+                    <!-- set using_server to true to run with XIOS
+                         in detached server/client mode -->
+                    <variable id="using_server" type="bool">false</variable>
+                    <!-- info_level is the verbosity level of XIOS
+                         0 for no messages; 100 for lots -->
+                    <variable id="info_level" type="int">10</variable>
+                </variable_group>
+            </variable_definition>
+        </context>
+
+        <context id="LMDZ" src="./context_pcm_physics.xml"/>
+</simulation>
Index: /trunk/LMDZ.PLUTO/libf/phypluto/physiq_mod.F90
===================================================================
--- /trunk/LMDZ.PLUTO/libf/phypluto/physiq_mod.F90	(revision 3507)
+++ /trunk/LMDZ.PLUTO/libf/phypluto/physiq_mod.F90	(revision 3508)
@@ -100,4 +100,5 @@
       use wxios, only: wxios_context_init, xios_context_finalize
 #endif
+      use write_output_mod, only: write_output
 
       implicit none
@@ -1038,5 +1039,4 @@
             inertia_min=minval(inertiedat(ig,:))
             inertia_max=maxval(inertiedat(ig,:))
-            print*, "inertia min max" , inertia_min, inertia_max
             ! diurnal and annual skin depth
             diurnal_skin=(inertia_min/volcapa)*sqrt(daysec/pi)
@@ -2334,6 +2334,6 @@
             call wstats(ngrid,trim(noms(iq))//'_col',trim(noms(iq))//'_col',  &
                         'kg m^-2',2,qcol(1,iq) )
-            call wstats(ngrid,trim(noms(iq))//'_reff',trim(noms(iq))//'_reff',&
-                        'm',3,reffrad(1,1,iq))
+            ! call wstats(ngrid,trim(noms(iq))//'_reff',trim(noms(iq))//'_reff',&
+            !             'm',3,reffrad(1,1,iq)) ! bug (2): Subscript #3 of the array REFFRAD has value 2 which is greater than the upper bound of 1
          end do
 
@@ -2357,80 +2357,63 @@
 !------------------------------------------------------------------------------
 
-      call writediagfi(ngrid,"Ls","solar longitude","deg",0,zls*180./pi)
-      call writediagfi(ngrid,"Lss","sub solar longitude","deg",0,zlss*180./pi)
-      call writediagfi(ngrid,"RA","right ascension","deg",0,right_ascen*180./pi)
-      call writediagfi(ngrid,"Declin","solar declination","deg",0,declin*180./pi)
-      call writediagfi(ngrid,"tsurf","Surface temperature","K",2,tsurf)
-      call writediagfi(ngrid,"ps","Surface pressure","Pa",2,ps)
-      call writediagfi(ngrid,"emis","Emissivity","",2,emis)
-
-      !! Pluto outputs
-      call writediagfi(ngrid,"dist_star","dist_star","AU",0,dist_star)
+      call write_output("Ls","solar longitude","deg",zls*180./pi)
+      ! call write_output("Lss","sub solar longitude","deg",zlss*180./pi)
+      call write_output("RA","right ascension","deg",right_ascen*180./pi)
+      call write_output("Declin","solar declination","deg",declin*180./pi)
+      call write_output("dist_star","dist_star","AU",dist_star)
+
+      call write_output("tsurf","Surface temperature","K",tsurf)
+      call write_output("ps","Surface pressure","Pa",ps)
+      call write_output("emis","Emissivity","",emis)
 
       if (fast) then
-         call writediagfi(ngrid,"globave","surf press","Pa",0,globave)
-         call writediagfi(ngrid,"fluxrad","fluxrad","W m-2",2,fluxrad)
-         call writediagfi(ngrid,"fluxgrd","fluxgrd","W m-2",2,fluxgrd)
-         call writediagfi(ngrid,"capcal","capcal","W.s m-2 K-1",2,capcal)
-         ! call writediagfi(ngrid,"dplanck","dplanck","W.s m-2 K-1",2,dplanck)
-         call writediagfi(ngrid,"tsoil","tsoil","K",3,tsoil)
+         call write_output("globave","surf press","Pa",globave)
+         call write_output("fluxrad","fluxrad","W m-2",fluxrad)
+         call write_output("fluxgrd","fluxgrd","W m-2",fluxgrd)
+         call write_output("capcal","capcal","W.s m-2 K-1",capcal)
+         ! call write_output("dplanck","dplanck","W.s m-2 K-1",dplanck)
+         call write_output("tsoil","tsoil","K",tsoil)
       else
          if (check_physics_outputs) then
             ! Check the validity of updated fields at the end of the physics step
-            call check_physics_fields("HERE physiq:", zt, zu, zv, pplev, zq)
-         endif
-
-         call writediagfi(ngrid,"temp","temperature","K",3,zt)
-         call writediagfi(ngrid,"teta","potential temperature","K",3,zh)
-         call writediagfi(ngrid,"u","Zonal wind","m.s-1",3,zu)
-         call writediagfi(ngrid,"v","Meridional wind","m.s-1",3,zv)
-         call writediagfi(ngrid,"w","Vertical wind","m.s-1",3,pw)
-         call writediagfi(ngrid,"p","Pressure","Pa",3,pplay)
+            call check_physics_fields("physiq:", zt, zu, zv, pplev, zq)
+         endif
+
+         call write_output("zzlay","Midlayer altitude", "m",zzlay(:,:))
+         call write_output("zzlev","Interlayer altitude", "m",zzlev(:,1:nlayer))
+         !call write_output('pphi','Geopotential',' ',pphi)
+
+         call write_output("temperature","temperature","K",zt)
+         call write_output("teta","potential temperature","K",zh)
+         call write_output("u","Zonal wind","m.s-1",zu)
+         call write_output("v","Meridional wind","m.s-1",zv)
+         call write_output("w","Vertical wind","m.s-1",pw)
+         call write_output("p","Pressure","Pa",pplay)
       endif
-
-!     Subsurface temperatures
-!        call writediagsoil(ngrid,"tsurf","Surface temperature","K",2,tsurf)
-!        call writediagsoil(ngrid,"temp","temperature","K",3,tsoil)
 
       ! Total energy balance diagnostics
       if(callrad)then
-
-         call writediagfi(ngrid,"ALB","Surface albedo"," ",2,albedo_equivalent)
-         call writediagfi(ngrid,"ALB_1st","First Band Surface albedo"," ",2,albedo(:,1))
-         call writediagfi(ngrid,"ISR","incoming stellar rad.","W m-2",2,fluxtop_dn)
-         call writediagfi(ngrid,"ASR","absorbed stellar rad.","W m-2",2,fluxabs_sw)
-         call writediagfi(ngrid,"OLR","outgoing longwave rad.","W m-2",2,fluxtop_lw)
-         call writediagfi(ngrid,"shad","rings"," ", 2, fract)
-
-!           call writediagfi(ngrid,"fluxsurfsw","sw surface flux.","W m-2",2,fluxsurf_sw)
-!           call writediagfi(ngrid,"fluxsurflw","lw back radiation.","W m-2",2,fluxsurf_lw)
-
-         call writediagfi(ngrid,"GND","heat flux from ground","W m-2",2,fluxgrd)
-         ! endif
-
-         call writediagfi(ngrid,"DYN","dynamical heat input","W m-2",2,fluxdyn)
-
+         call write_output("ALB","Surface albedo"," ",albedo_equivalent)
+         call write_output("ASR","absorbed stellar rad.","W m-2",fluxabs_sw)
+         call write_output("ISR","incoming stellar rad.","W m-2",fluxtop_dn)
+         call write_output("OLR","outgoing longwave rad.","W m-2",fluxtop_lw)
+         call write_output("GND","heat flux from ground","W m-2",fluxgrd)
+         call write_output("DYN","dynamical heat input","W m-2",fluxdyn)
       endif ! end of 'callrad'
 
       if(enertest) then
          if (calldifv) then
-
-            call writediagfi(ngrid,"q2","turbulent kinetic energy","J.kg^-1",3,q2)
-            call writediagfi(ngrid,"sensibFlux","sensible heat flux","w.m^-2",2,sensibFlux)
-
-!             call writediagfi(ngrid,"dEzdiff","turbulent diffusion heating (-sensible flux)","w.m^-2",3,dEzdiff)
-!             call writediagfi(ngrid,"dEdiff","integrated turbulent diffusion heating (-sensible flux)","w.m^-2",2,dEdiff)
-!             call writediagfi(ngrid,"dEdiffs","In TurbDiff (correc rad+latent heat) surf nrj change","w.m^-2",2,dEdiffs)
-
+            call write_output("q2","turbulent kinetic energy","J.kg^-1",q2)
+            call write_output("sensibFlux","sensible heat flux","w.m^-2",sensibFlux)
          endif
 
          if (corrk) then
-            call writediagfi(ngrid,"dEzradsw","radiative heating","w.m^-2",3,dEzradsw)
-            call writediagfi(ngrid,"dEzradlw","radiative heating","w.m^-2",3,dEzradlw)
+            call write_output("dEzradsw","radiative heating","w.m^-2",dEzradsw)
+            call write_output("dEzradlw","radiative heating","w.m^-2",dEzradlw)
          endif
 
          if (generic_condensation) then
-            call writediagfi(ngrid,"genericconddE","heat from generic condensation","W m-2",2,genericconddE)
-            call writediagfi(ngrid,"dt_generic_condensation","heating from generic condensation","K s-1",3,dt_generic_condensation)
+            call write_output("genericconddE","heat from generic condensation","W m-2",genericconddE)
+            call write_output("dt_generic_condensation","heating from generic condensation","K s-1",dt_generic_condensation)
          endif
 
@@ -2442,81 +2425,73 @@
             do nw=1,L_NSPECTV
                write(str2,'(i2.2)') nw
-               call writediagfi(ngrid,'dtauv'//str2,'Layer optical thickness attenuation in VI band '//str2,'',1,int_dtauv(:,nlayer:1:-1,nw))
+               call write_output('dtauv'//str2,'Layer optical thickness attenuation in VI band '//str2,'',int_dtauv(:,nlayer:1:-1,nw))
             enddo
             do nw=1,L_NSPECTI
                write(str2,'(i2.2)') nw
-               call writediagfi(ngrid,'dtaui'//str2,'Layer optical thickness attenuation in IR band '//str2,'',1,int_dtaui(:,nlayer:1:-1,nw))
+               call write_output('dtaui'//str2,'Layer optical thickness attenuation in IR band '//str2,'',int_dtaui(:,nlayer:1:-1,nw))
             enddo
          endif
 
          ! Temporary inclusions for heating diagnostics.
-         call writediagfi(ngrid,"zdtsw","SW heating","T s-1",3,zdtsw)
-         call writediagfi(ngrid,"zdtlw","LW heating","T s-1",3,zdtlw)
-         call writediagfi(ngrid,"dtrad","radiative heating","K s-1",3,dtrad)
-         call writediagfi(ngrid,"zdtdyn","Dyn. heating","T s-1",3,zdtdyn)
+         call write_output("zdtsw","SW heating","T s-1",zdtsw)
+         call write_output("zdtlw","LW heating","T s-1",zdtlw)
+         call write_output("dtrad","radiative heating","K s-1",dtrad)
+         call write_output("zdtdyn","Dyn. heating","T s-1",zdtdyn)
 
          ! For Debugging.
-         !call writediagfi(ngrid,'rnat','Terrain type',' ',2,real(rnat))
-         !call writediagfi(ngrid,'pphi','Geopotential',' ',3,pphi)
+         !call write_output('rnat','Terrain type',' ',real(rnat))
 
       ! Output tracers.
       if (tracer) then
-         ! call writediagfi(ngrid,"zdtc","tendancy T cond N2","K",3,zdtc)
+         ! call write_output("zdtc","tendancy T cond N2","K",zdtc)
 
          do iq=1,nq
-            call writediagfi(ngrid,noms(iq),noms(iq),'kg/kg',3,zq(1,1,iq))
-            ! call writediagfi(ngrid,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf',  &
-         !                    'kg m^-2',2,qsurf_hist(1,iq) )
-           call writediagfi(ngrid,trim(noms(iq))//'_col',trim(noms(iq))//'_col',    &
-                           'kg m^-2',2,qcol(1,iq) )
-            call writediagfi(ngrid,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf',  &
-                         'kg m^-2',2,qsurf(1,iq) )
-
-            if(generic_condensation)then
-               call writediagfi(ngrid,"rneb_generic","GCS cloud fraction (generic condensation)"," ",3,rneb_generic)
-               call writediagfi(ngrid,"CLF","GCS cloud fraction"," ",3,cloudfrac)
-               call writediagfi(ngrid,"RH_generic","GCS relative humidity"," ",3,RH_generic)
-            endif
-           ! call writediagfi(ngrid,"tau_col","Total aerosol optical depth","[]",2,tau_col)
+            call write_output(noms(iq),noms(iq),'kg/kg',zq(:,:,iq))
+            ! call write_output(trim(noms(iq))//'_surf',trim(noms(iq))//'_surf',  &
+         !                    'kg m^-2',qsurf_hist(1,iq) )
+           call write_output(trim(noms(iq))//'_col',trim(noms(iq))//'_col',    &
+                           'kg m^-2',qcol(:,iq) )
+            call write_output(trim(noms(iq))//'_surf',trim(noms(iq))//'_surf',  &
+                         'kg m^-2',qsurf(:,iq) )
          enddo ! end of 'nq' loop
 
          !Pluto specific
-         call writediagfi(ngrid,'n2_iceflux','n2_iceflux',"kg m^-2 s^-1",2,flusurf(1,igcm_n2) )
+         call write_output('n2_iceflux','n2_iceflux',"kg m^-2 s^-1",flusurf(1,igcm_n2) )
          if (haze_radproffix)then
-            call writediagfi(ngrid,'haze_reff','haze_reff','m',3,reffrad(1,1,1))
+            call write_output('haze_reff','haze_reff','m',reffrad(1,1,1))
          end if
          if (methane) then
-            ! call writediagfi(ngrid,"rice_ch4","ch4 ice mass mean radius","m",3,rice_ch4)
-            ! call writediagfi(ngrid,"zq1temp_ch4"," "," ",2,zq1temp_ch4)
-            ! call writediagfi(ngrid,"qsat_ch4"," "," ",2,qsat_ch4)
-            ! call writediagfi(ngrid,"qsat_ch4_l1"," "," ",2,qsat_ch4_l1)
-
-            call writediagfi(ngrid,'ch4_iceflux','ch4_iceflux',&
-                              "kg m^-2 s^-1",2,flusurf(1,igcm_ch4_ice) )
-            call writediagfi(ngrid,"vmr_ch4","vmr_ch4","%",2,vmr_ch4)
+            ! call write_output("rice_ch4","ch4 ice mass mean radius","m",rice_ch4)
+            ! call write_output("zq1temp_ch4"," "," ",zq1temp_ch4)
+            ! call write_output("qsat_ch4"," "," ",qsat_ch4)
+            ! call write_output("qsat_ch4_l1"," "," ",qsat_ch4_l1)
+
+            call write_output('ch4_iceflux','ch4_iceflux',&
+                              "kg m^-2 s^-1",flusurf(1,igcm_ch4_ice) )
+            call write_output("vmr_ch4","vmr_ch4","%",vmr_ch4)
             if (.not.fast) then
-               call writediagfi(ngrid,"zrho_ch4","zrho_ch4","kg.m-3",3,zrho_ch4(:,:))
+               call write_output("zrho_ch4","zrho_ch4","kg.m-3",zrho_ch4(:,:))
             endif
 
             ! Tendancies
-            call writediagfi(ngrid,"zdqcn2_ch4","zdq condn2 ch4","",&
-                           3,zdqc(:,:,igcm_ch4_gas))
-            call writediagfi(ngrid,"zdqdif_ch4","zdqdif ch4","",&
-                           3,zdqdif(:,:,igcm_ch4_gas))
-            call writediagfi(ngrid,"zdqsdif_ch4_ice","zdqsdif ch4","",&
-                           2,zdqsdif(:,igcm_ch4_ice))
-            call writediagfi(ngrid,"zdqadj_ch4","zdqadj ch4","",&
-                           3,zdqadj(:,:,igcm_ch4_gas))
+            call write_output("zdqcn2_ch4","zdq condn2 ch4","",&
+                           zdqc(:,:,igcm_ch4_gas))
+            call write_output("zdqdif_ch4","zdqdif ch4","",&
+                           zdqdif(:,:,igcm_ch4_gas))
+            call write_output("zdqsdif_ch4_ice","zdqsdif ch4","",&
+                           zdqsdif(:,igcm_ch4_ice))
+            call write_output("zdqadj_ch4","zdqadj ch4","",&
+                           zdqadj(:,:,igcm_ch4_gas))
             if (sedimentation) then
-               call writediagfi(ngrid,"zdqsed_ch4","zdqsed ch4","",&
-                              3,zdqsed(:,:,igcm_ch4_gas))
-               call writediagfi(ngrid,"zdqssed_ch4","zdqssed ch4","",&
-                              2,zdqssed(:,igcm_ch4_gas))
+               call write_output("zdqsed_ch4","zdqsed ch4","",&
+                              zdqsed(:,:,igcm_ch4_gas))
+               call write_output("zdqssed_ch4","zdqssed ch4","",&
+                              zdqssed(:,igcm_ch4_gas))
             endif
             if (metcloud.and.(.not.fast)) then
-               call writediagfi(ngrid,"zdtch4cloud","ch4 cloud","T s-1",&
-                           3,zdtch4cloud)
-               call writediagfi(ngrid,"zdqch4cloud","ch4 cloud","T s-1",&
-                           3,zdqch4cloud(1,1,igcm_ch4_gas))
+               call write_output("zdtch4cloud","ch4 cloud","T s-1",&
+                           zdtch4cloud)
+               call write_output("zdqch4cloud","ch4 cloud","T s-1",&
+                           zdqch4cloud(1,1,igcm_ch4_gas))
             endif
 
@@ -2524,42 +2499,42 @@
 
          if (carbox) then
-            ! call writediagfi(ngrid,"zdtcocloud","tendancy T cocloud","K",3,zdtcocloud)
-            call writediagfi(ngrid,'co_iceflux','co_iceflux',&
-                               "kg m^-2 s^-1",2,flusurf(1,igcm_co_ice) )
-            call writediagfi(ngrid,"vmr_co","vmr_co","%",2,vmr_co)
+            ! call write_output("zdtcocloud","tendancy T cocloud","K",zdtcocloud)
+            call write_output('co_iceflux','co_iceflux',&
+                               "kg m^-2 s^-1",flusurf(1,igcm_co_ice) )
+            call write_output("vmr_co","vmr_co","%",vmr_co)
             if (.not.fast) THEN
-               call writediagfi(ngrid,"zrho_co","zrho_co","kg.m-3",3,zrho_co(:,:))
+               call write_output("zrho_co","zrho_co","kg.m-3",zrho_co(:,:))
             endif
          endif
 
          if (haze) then
-   !         call writediagfi(ngrid,"zrho_haze","zrho_haze","kg.m-3",3,zrho_haze(:,:))
-            call writediagfi(ngrid,"zdqrho_photprec","zdqrho_photprec",&
-                        "kg.m-3.s-1",3,zdqrho_photprec(:,:))
-            call writediagfi(ngrid,"zdqphot_prec","zdqphot_prec","",&
-                                                3,zdqphot_prec(:,:))
-            call writediagfi(ngrid,"zdqhaze_ch4","zdqhaze_ch4","",&
-                     3,zdqhaze(:,:,igcm_ch4_gas))
-            call writediagfi(ngrid,"zdqhaze_prec","zdqhaze_prec","",&
-                     3,zdqhaze(:,:,igcm_prec_haze))
+   !         call write_output("zrho_haze","zrho_haze","kg.m-3",zrho_haze(:,:))
+            call write_output("zdqrho_photprec","zdqrho_photprec",&
+                        "kg.m-3.s-1",zdqrho_photprec(:,:))
+            call write_output("zdqphot_prec","zdqphot_prec","",&
+                                                zdqphot_prec(:,:))
+            call write_output("zdqhaze_ch4","zdqhaze_ch4","",&
+                     zdqhaze(:,:,igcm_ch4_gas))
+            call write_output("zdqhaze_prec","zdqhaze_prec","",&
+                     zdqhaze(:,:,igcm_prec_haze))
             if (igcm_haze.ne.0) then
-               call writediagfi(ngrid,"zdqhaze_haze","zdqhaze_haze","",&
-                        3,zdqhaze(:,:,igcm_haze))
+               call write_output("zdqhaze_haze","zdqhaze_haze","",&
+                        zdqhaze(:,:,igcm_haze))
                if (sedimentation) then
-                  call writediagfi(ngrid,"zdqssed_haze","zdqssed haze",&
-                     "kg/m2/s",2,zdqssed(:,igcm_haze))
+                  call write_output("zdqssed_haze","zdqssed haze",&
+                     "kg/m2/s",zdqssed(:,igcm_haze))
                endif
             endif
-            call writediagfi(ngrid,"zdqphot_ch4","zdqphot_ch4","",&
-                                                3,zdqphot_ch4(:,:))
-            call writediagfi(ngrid,"zdqconv_prec","zdqconv_prec","",&
-                                                3,zdqconv_prec(:,:))
-   !         call writediagfi(ngrid,"zdqhaze_col","zdqhaze col","kg/m2/s",
-   !     &                   2,zdqhaze_col(:))
+            call write_output("zdqphot_ch4","zdqphot_ch4","",&
+                                                zdqphot_ch4(:,:))
+            call write_output("zdqconv_prec","zdqconv_prec","",&
+                                                zdqconv_prec(:,:))
+   !         call write_output("zdqhaze_col","zdqhaze col","kg/m2/s",
+   !     &                   zdqhaze_col(:))
          endif
 
          if (aerohaze) then
-            call writediagfi(ngrid,"tau_col",&
-               "Total aerosol optical depth","opacity",2,tau_col)
+            call write_output("tau_col",&
+               "Total aerosol optical depth","opacity",tau_col)
          endif
 
