Index: trunk/LMDZ.MARS/libf/phymars/physiq_mod.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/physiq_mod.F	(revision 2930)
+++ trunk/LMDZ.MARS/libf/phymars/physiq_mod.F	(revision 2931)
@@ -2512,26 +2512,19 @@
 
 #ifndef MESOSCALE
-            IF (ngrid .eq. 1) THEN
-               dimout=0
-            ELSE
-               dimout=2
-            ENDIF
             DO n=1,n_out
                write(zstring, '(F8.6)') z_out(n)
                call WRITEDIAGFI(ngrid,'T_out_'//trim(zstring),
-     &   'potential temperature at z_out','K',dimout,T_out(:,n))
+     &   'potential temperature at z_out','K',2,T_out(:,n))
                call WRITEDIAGFI(ngrid,'u_out_'//trim(zstring),
-     &   'horizontal velocity norm at z_out','m/s',dimout,u_out(:,n))
+     &   'horizontal velocity norm at z_out','m/s',2,u_out(:,n))
             ENDDO
             call WRITEDIAGFI(ngrid,'u_star',
-     &   'friction velocity','m/s',dimout,ustar)
+     &   'friction velocity','m/s',2,ustar)
             call WRITEDIAGFI(ngrid,'teta_star',
-     &   'friction potential temperature','K',dimout,tstar)
-!            call WRITEDIAGFI(ngrid,'L',
-!     &   'Monin Obukhov length','m',dimout,L_mo)
+     &   'friction potential temperature','K',2,tstar)
             call WRITEDIAGFI(ngrid,'vvv',
-     &   'Vertical velocity variance at zout','m',dimout,vvv)
+     &   'Vertical velocity variance at zout','m',2,vvv)
             call WRITEDIAGFI(ngrid,'vhf',
-     &   'Vertical heat flux at zout','m',dimout,vhf)
+     &   'Vertical heat flux at zout','m',2,vhf)
 #else
          T_out1(:)=T_out(:,1)
@@ -2610,5 +2603,5 @@
 #endif
 
-         IF (ngrid.NE.1) then
+c         IF (ngrid.NE.1) then
 
 c        -------------------------------------------------------------------
@@ -2640,8 +2633,10 @@
            if(doubleq) then
               do ig=1,ngrid 
+         IF (sedimentation) THEN 
                 dqdustsurf(ig) = 
      &                zdqssed(ig,igcm_dust_mass)*tauscaling(ig)
                 dndustsurf(ig) = 
      &                zdqssed(ig,igcm_dust_number)*tauscaling(ig)
+         ENDIF
                 ndust(ig,:) =
      &                zq(ig,:,igcm_dust_number)*tauscaling(ig)
@@ -2651,8 +2646,10 @@
               if (scavenging) then
                 do ig=1,ngrid 
+         IF (sedimentation) THEN 
                   dqdustsurf(ig) = dqdustsurf(ig) + 
      &                     zdqssed(ig,igcm_ccn_mass)*tauscaling(ig)
                   dndustsurf(ig) = dndustsurf(ig) + 
      &                     zdqssed(ig,igcm_ccn_number)*tauscaling(ig)
+         ENDIF
                   nccn(ig,:) =
      &                     zq(ig,:,igcm_ccn_number)*tauscaling(ig)
@@ -2853,10 +2850,4 @@
           call wstats(ngrid,"emis","Surface emissivity","w.m-1",2,
      &                emis(:,iflat))
-c          call wstats(ngrid,"ssurf","Surface stress","N.m-2",
-c    &                2,zstress)
-c          call wstats(ngrid,"sw_htrt","sw heat.rate",
-c    &                 "W.m-2",3,zdtsw)
-c          call wstats(ngrid,"lw_htrt","lw heat.rate",
-c    &                 "W.m-2",3,zdtlw)
           call wstats(ngrid,"fluxsurf_dir_dn_sw",
      &                "Direct incoming SW flux at surface",
@@ -2950,14 +2941,4 @@
              
              if (doubleq) then
-c            call wstats(ngrid,'qsurf','qsurf',
-c     &                       'kg.m-2',2,qsurf(1,igcm_dust_mass))
-c            call wstats(ngrid,'Nsurf','N particles',
-c     &                       'N.m-2',2,qsurf(1,igcm_dust_number))
-c            call wstats(ngrid,'dqsdev','ddevil lift',
-c    &                       'kg.m-2.s-1',2,zdqsdev(1,1))
-c            call wstats(ngrid,'dqssed','sedimentation',
-c     &                       'kg.m-2.s-1',2,zdqssed(1,1))
-c            call wstats(ngrid,'dqsdif','diffusion',
-c     &                       'kg.m-2.s-1',2,zdqsdif(1,1))
                call wstats(ngrid,'dqsdust',
      &                        'deposited surface dust mass',
@@ -3159,6 +3140,4 @@
      &    "Surface emissivity","w.m-1",2, emis(:,islope))
          ENDDO
-c        call WRITEDIAGFI(ngrid,"pplay","Pressure","Pa",3,zplay)
-c        call WRITEDIAGFI(ngrid,"pplev","Pressure","Pa",3,zplev)
          call WRITEDIAGFI(ngrid,"zzlev","Interlayer altitude",
      &                    "m",3,zzlev(:,1:nlayer))
@@ -3223,9 +3202,6 @@
          call WRITEDIAGFI(ngrid,"w","Vertical wind","m.s-1",3,pw)
          call WRITEDIAGFI(ngrid,"rho","density","kg.m-3",3,rho)
-c        call WRITEDIAGFI(ngrid,"q2","q2","kg.m-3",3,q2)
-c        call WRITEDIAGFI(ngrid,'Teta','T potentielle','K',3,zh)
          call WRITEDIAGFI(ngrid,"pressure","Pressure","Pa",3,zplay)
-c        call WRITEDIAGFI(ngrid,"ssurf","Surface stress","N.m-2",2,
-c    &                  zstress)
+         call WRITEDIAGFI(ngrid,"pplev","Pressure","Pa",3,zplev)
         call WRITEDIAGFI(ngrid,'sw_htrt','sw heat. rate',
      &                   'K/s',3,zdtsw)
@@ -3501,19 +3477,4 @@
         end if !(water)
 
-        if (water.and..not.photochem) then
-          iq = nq
-c          write(str2(1:2),'(i2.2)') iq
-c          call WRITEDIAGFI(ngrid,'dqs'//str2,'dqscloud',
-c    &                      'kg.m-2',2,zdqscloud(1,iq))
-c          call WRITEDIAGFI(ngrid,'dqch'//str2,'var chim',
-c    &                      'kg/kg',3,zdqchim(1,1,iq))
-c          call WRITEDIAGFI(ngrid,'dqd'//str2,'var dif',
-c    &                      'kg/kg',3,zdqdif(1,1,iq))
-c          call WRITEDIAGFI(ngrid,'dqa'//str2,'var adj',
-c    &                      'kg/kg',3,zdqadj(1,1,iq))
-c          call WRITEDIAGFI(ngrid,'dqc'//str2,'var c',
-c    &                      'kg/kg',3,zdqc(1,1,iq))
-        end if  !(water.and..not.photochem)
-
 c        ----------------------------------------------------------
 c        Outputs of the dust cycle
@@ -3540,18 +3501,4 @@
 #ifndef MESOINI
            if (doubleq) then
-c            call WRITEDIAGFI(ngrid,'qsurf','qsurf',
-c     &                       'kg.m-2',2,qsurf(1,igcm_dust_mass))
-c            call WRITEDIAGFI(ngrid,'Nsurf','N particles',
-c     &                       'N.m-2',2,qsurf(1,igcm_dust_number))
-c            call WRITEDIAGFI(ngrid,'dqsdev','ddevil lift',
-c    &                       'kg.m-2.s-1',2,zdqsdev(1,1))
-c            call WRITEDIAGFI(ngrid,'dqssed','sedimentation',
-c     &                       'kg.m-2.s-1',2,zdqssed(1,1))
-c            call WRITEDIAGFI(ngrid,'dqsdif','diffusion',
-c     &                       'kg.m-2.s-1',2,zdqsdif(1,1))
-c             call WRITEDIAGFI(ngrid,'sedice','sedimented ice',
-c     &                       'kg.m-2.s-1',2,zdqssed(:,igcm_h2o_ice))
-c             call WRITEDIAGFI(ngrid,'subice','sublimated ice',
-c     &                       'kg.m-2.s-1',2,zdqsdif(:,igcm_h2o_ice))
              call WRITEDIAGFI(ngrid,'dqsdust',
      &                        'deposited surface dust mass',
@@ -3648,4 +3595,7 @@
      &          ,'sedimentation q','kg.m-2.s-1',3,
      &          zdqsed(:,:,igcm_stormdust_mass))
+             call WRITEDIAGFI(ngrid,'zdqsed_dustN'
+     &          ,'sedimentation N','Nbr.m-2.s-1',3,
+     &          zdqsed(:,:,igcm_dust_number))
              call WRITEDIAGFI(ngrid,'rdust','rdust',
      &                        'm',3,rdust)
@@ -3714,9 +3664,4 @@
          ENDDO
            endif ! (scavenging)
-
-c          if (submicron) then
-c            call WRITEDIAGFI(ngrid,'dustsubm','subm mass mr',
-c    &                        'kg/kg',3,pq(1,1,igcm_dust_submicron))
-c          endif ! (submicron)
 
 #else
@@ -3789,20 +3734,4 @@
      $           0,[PhiEscD])
 
-!            call wstats(ngrid,"PhiH","H escape flux","s-1",
-!     $           0,[PhiEscH])
-!            call wstats(ngrid,"PhiH2","H2 escape flux","s-1",
-!     $           0,[PhiEscH2])
-!            call wstats(ngrid,"PhiD","D escape flux","s-1",
-!     $           0,[PhiEscD])
-            
-!            call wstats(ngrid,"q15um","15 um cooling","K/s",
-!     $           3,zdtnlte)
-!            call wstats(ngrid,"quv","UV heating","K/s",
-!     $           3,zdteuv)
-!            call wstats(ngrid,"cond","Thermal conduction","K/s",
-!     $           3,zdtconduc)
-!            call wstats(ngrid,"qnir","NIR heating","K/s",
-!     $           3,zdtnirco2)
-
          endif  !(callthermos)
 
@@ -3822,19 +3751,7 @@
 c        ----------------------------------------------------------
       if (calltherm) then
-!        call WRITEDIAGFI(ngrid,'dtke',
-!     &                   'tendance tke thermiques','m**2/s**2',
-!     &                   3,dtke_th)
-!        call WRITEDIAGFI(ngrid,'d_u_ajs',
-!     &                   'tendance u thermiques','m/s',
-!     &                   3,pdu_th*ptimestep)
-!        call WRITEDIAGFI(ngrid,'d_v_ajs',
-!     &                   'tendance v thermiques','m/s',
-!     &                   3,pdv_th*ptimestep)
-!          if (nq .eq. 2) then
-!            call WRITEDIAGFI(ngrid,'deltaq_th',
-!     &                       'delta q thermiques','kg/kg',
-!     &                       3,ptimestep*pdq_th(:,:,2))
-!          end if
-
+        call WRITEDIAGFI(ngrid,'lmax_th',
+     &              'hauteur du thermique','point',
+     &                         2,lmax_th_out)
         call WRITEDIAGFI(ngrid,'zmax_th',
      &                   'hauteur du thermique','m',
@@ -3869,389 +3786,35 @@
      &                     3,tsoil(:,:,islope))
          ENDDO
-         ! Write surface temperature
-!        call writediagsoil(ngrid,"tsurf","Surface temperature","K",
-!    &                     2,tsurf)
-
-c        ==========================================================
-c        END OF WRITEDIAGFI
-c        ==========================================================
-#endif
-! of ifdef MESOSCALE
-
-      ELSE     ! if(ngrid.eq.1)
-
-#ifndef MESOSCALE
-         write(*,
-     &    '("Ls =",f11.6," tau_pref_scenario(",f4.0," Pa) =",f9.6)')
-     &    zls*180./pi,odpref,tau_pref_scenario
-c      ----------------------------------------------------------------------
-c      Output in grads file "g1d" (ONLY when using testphys1d)
-c      (output at every X physical timestep)
-c      ----------------------------------------------------------------------
-c
-c        CALL writeg1d(ngrid,1,fluxsurf_lw,'Fs_ir','W.m-2')
-c        CALL writeg1d(ngrid,1,tsurf,'tsurf','K')
-c        CALL writeg1d(ngrid,1,ps,'ps','Pa')
-c        CALL writeg1d(ngrid,nlayer,zt,'T','K')
-c        CALL writeg1d(ngrid,nlayer,pu,'u','m.s-1')
-c        CALL writeg1d(ngrid,nlayer,pv,'v','m.s-1')
-c        CALL writeg1d(ngrid,nlayer,pw,'w','m.s-1')
-        call writediagsoil(ngrid,"soiltemp","Soil temperature","K",
-     &                     3,tsoil(:,:,iflat))
-         do islope=1,nslope
-          write(str2(1:2),'(i2.2)') islope
-        call writediagsoil(ngrid,"soiltemp_slope"//str2,
-     &                     "Soil temperature","K",
-     &                     3,tsoil(:,:,islope))
-         ENDDO
-
-! THERMALS STUFF 1D
-         if(calltherm) then
-
-        call WRITEDIAGFI(ngrid,'lmax_th',
-     &              'hauteur du thermique','point',
-     &                         0,lmax_th_out)
-        call WRITEDIAGFI(ngrid,'zmax_th',
-     &              'hauteur du thermique','m',
-     &                         0,zmax_th)
-        call WRITEDIAGFI(ngrid,'hfmax_th',
-     &              'maximum TH heat flux','K.m/s',
-     &                         0,hfmax_th)
-        call WRITEDIAGFI(ngrid,'wstar',
-     &              'maximum TH vertical velocity','m/s',
-     &                         0,wstar)
-
-         end if ! of if (calltherm)
-
-         call WRITEDIAGFI(ngrid,'w','vertical velocity'
-     &                              ,'m/s',1,pw)
-         call WRITEDIAGFI(ngrid,"q2","q2","kg.m-3",1,q2)
-         call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",0,
-     &                  tsurf(:,iflat))
-         do islope=1,nslope
-          write(str2(1:2),'(i2.2)') islope
-         call WRITEDIAGFI(ngrid,"tsurf_slope"//str2,
-     &             "Surface temperature","K",0,
-     &                  tsurf(:,islope))
-         ENDDO
-         call WRITEDIAGFI(ngrid,"u","u wind","m/s",1,zu)
-         call WRITEDIAGFI(ngrid,"v","v wind","m/s",1,zv)
-
-         call WRITEDIAGFI(ngrid,"pplay","Pressure","Pa",1,zplay)
-         call WRITEDIAGFI(ngrid,"pplev","Pressure","Pa",1,zplev)
-         call WRITEDIAGFI(ngrid,"rho","rho","kg.m-3",1,rho)
+
+!PREVIOUSLY IN 1D ONLY
          call WRITEDIAGFI(ngrid,"dtrad","rad. heat. rate",
-     &              "K.s-1",1,dtrad)
-         call WRITEDIAGFI(ngrid,'sw_htrt','sw heat. rate',
-     &                   'w.m-2',1,zdtsw)
-         call WRITEDIAGFI(ngrid,'lw_htrt','lw heat. rate',
-     &                   'w.m-2',1,zdtlw)
-         call WRITEDIAGFI(ngrid,"co2ice","co2 ice thickness"
-     &       ,"kg.m-2",0,qsurf(:,igcm_co2,iflat))
-         do islope=1,nslope
-          write(str2(1:2),'(i2.2)') islope
-         call WRITEDIAGFI(ngrid,"co2ice_slope"//str2,
-     &                  "co2 ice thickness"
-     &       ,"kg.m-2",0,qsurf(:,igcm_co2,islope))
-         ENDDO
-
-         if (igcm_co2.ne.0) then
-          call co2sat(ngrid*nlayer,zt,zqsatco2)
-           do ig=1,ngrid 
-            do l=1,nlayer
-               satuco2(ig,l) = zq(ig,l,igcm_co2)* 
-     &              (mmean(ig,l)/44.01)*zplay(ig,l)/zqsatco2(ig,l)
-                  
-c               write(*,*) "In PHYSIQMOD, pt,zt,time ",pt(ig,l)
-c     &              ,zt(ig,l),ptime
-            enddo
-           enddo
-          endif
-
-         call WRITEDIAGFI(ngrid,'ps','Surface pressure','Pa',0,ps)
-         call WRITEDIAGFI(ngrid,'temp','Temperature ',
-     &                       'K JA',1,zt)
-c         call WRITEDIAGFI(ngrid,'temp2','Temperature ',
-c     &        'K JA2',1,pt)
-
-c           CALL writeg1d(ngrid,1,tau,'tau','SI')
-            do iq=1,nq
-c              CALL writeg1d(ngrid,nlayer,zq(1,1,iq),noms(iq),'kg/kg') 
-               call WRITEDIAGFI(ngrid,trim(noms(iq)),
-     &              trim(noms(iq)),'kg/kg',1,zq(1,1,iq))
-            end do
-           if (doubleq) then
-             call WRITEDIAGFI(ngrid,'rdust','rdust',
-     &                        'm',1,rdust)
-           endif ! doubleq 1D
+     &              "K.s-1",3,dtrad)
+
            if (rdstorm) then
              call writediagfi(1,'aerosol_dust','opacity of env. dust',''
-     &                        ,1,aerosol(:,:,iaer_dust_doubleq))
+     &                        ,3,aerosol(:,:,iaer_dust_doubleq))
              call writediagfi(1,'aerosol_stormdust',
      &                         'opacity of storm dust',''
-     &                        ,1,aerosol(:,:,iaer_stormdust_doubleq))
+     &                        ,3,aerosol(:,:,iaer_stormdust_doubleq))
              call WRITEDIAGFI(ngrid,'dqsdifdustq','diffusion',
-     &          'kg.m-2.s-1',0,zdqsdif(1,igcm_dust_mass,iflat))
+     &          'kg.m-2.s-1',2,zdqsdif(:,igcm_dust_mass,iflat))
          do islope=1,nslope
           write(str2(1:2),'(i2.2)') islope
              call WRITEDIAGFI(ngrid,'dqsdifdustq_slope'//str2,
      &             'diffusion',
-     &          'kg.m-2.s-1',0,zdqsdif(1,igcm_dust_mass,islope))
+     &          'kg.m-2.s-1',2,zdqsdif(:,igcm_dust_mass,islope))
          ENDDO
              call WRITEDIAGFI(ngrid,'dqsdifrdsq','diffusion',
-     &          'kg.m-2.s-1',0,zdqsdif(1,igcm_stormdust_mass,iflat))
+     &          'kg.m-2.s-1',2,zdqsdif(:,igcm_stormdust_mass,iflat))
          do islope=1,nslope
           write(str2(1:2),'(i2.2)') islope
              call WRITEDIAGFI(ngrid,'dqsdifrdsq_slope'//str2,
      &           'diffusion',
-     &          'kg.m-2.s-1',0,zdqsdif(1,igcm_stormdust_mass,islope))
+     &          'kg.m-2.s-1',2,zdqsdif(:,igcm_stormdust_mass,islope))
          ENDDO
-             call WRITEDIAGFI(ngrid,'mstormdtot',
-     &                        'total mass of stormdust only',
-     &                        'kg.m-2',0,mstormdtot)
-             call WRITEDIAGFI(ngrid,'mdusttot',
-     &                        'total mass of dust only',
-     &                        'kg.m-2',0,mdusttot)
-             call WRITEDIAGFI(ngrid,'tau_pref_scenario',
-     &                        'Prescribed dust ref opt depth at 610 Pa',
-     &                        'NU',0,tau_pref_scenario)
-             call WRITEDIAGFI(ngrid,'tau_pref_gcm',
-     &                        'Dust ref opt depth at 610 Pa in the GCM',
-     &                        'NU',0,tau_pref_gcm)
-             call WRITEDIAGFI(ngrid,'rdsdqsdust',
-     &                        'deposited surface stormdust mass',
-     &                        'kg.m-2.s-1',0,rdsdqdustsurf)
-             call WRITEDIAGFI(ngrid,'rdsdustq','storm Dust mass mr',
-     &                        'kg/kg',1,rdsqdust)
-             call WRITEDIAGFI(ngrid,"stormfract", 
-     &                         "fraction of the mesh,with stormdust",
-     &                         "none",0,totstormfract)
-             call WRITEDIAGFI(ngrid,'rdsqsurf',
-     &                 'stormdust at surface',
-     &                 'kg.m-2',0,qsurf(:,igcm_stormdust_mass,iflat))
-         do islope=1,nslope
-          write(str2(1:2),'(i2.2)') islope
-             call WRITEDIAGFI(ngrid,'rdsqsurf_slope'//str2,
-     &                 'stormdust at surface',
-     &                 'kg.m-2',0,qsurf(:,igcm_stormdust_mass,islope))
-         ENDDO
-             call WRITEDIAGFI(ngrid,'qsurf',
-     &                  'dust mass at surface',
-     &                  'kg.m-2',0,qsurf(:,igcm_dust_mass,iflat))
-         do islope=1,nslope
-          write(str2(1:2),'(i2.2)') islope
-             call WRITEDIAGFI(ngrid,'qsurf_slope'//str2,
-     &                  'dust mass at surface',
-     &                  'kg.m-2',0,qsurf(:,igcm_dust_mass,islope))
-         ENDDO
-             call WRITEDIAGFI(ngrid,'wspeed','vertical speed stormdust',
-     &                        'm/s',1,wspeed)
-             call WRITEDIAGFI(ngrid,'totaldustq','total dust mass',
-     &                        'kg/kg',1,qdusttotal)
-             call WRITEDIAGFI(ngrid,'dsords',
-     &                       'density scaled opacity of stormdust',
-     &                       'm2.kg-1',1,dsords)
-             call WRITEDIAGFI(ngrid,'zdqsed_dustq'
-     &          ,'sedimentation q','kg.m-2.s-1',1,
-     &          zdqsed(1,:,igcm_dust_mass))
-             call WRITEDIAGFI(ngrid,'zdqsed_rdsq'
-     &          ,'sedimentation q','kg.m-2.s-1',1,
-     &          zdqsed(1,:,igcm_stormdust_mass))
-           endif !(rdstorm 1D)
-
-           if (water.AND.tifeedback) then
-             call WRITEDIAGFI(ngrid,"soiltemp",
-     &                              "Soil temperature","K",
-     &                              1,tsoil)
-             call WRITEDIAGFI(ngrid,'soilti',
-     &                       'Soil Thermal Inertia',
-     &                       'J.s-1/2.m-2.K-1',1,inertiesoil)
-           endif
-         
-cccccccccccccccccc scavenging & water outputs 1D TN ccccccccccccccc
-ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-         IF (water) THEN
-
-           if (.not.activice) then
-
-             tauTES=0
-             do l=1,nlayer
-               Qabsice = min(
-     &             max(0.4e6*rice(1,l)*(1.+nuice_ref)-0.05 ,0.),1.2
-     &                      )
-               opTES(1,l)= 0.75 * Qabsice *
-     &             zq(1,l,igcm_h2o_ice) *
-     &             (zplev(1,l) - zplev(1,l+1)) / g
-     &             / (rho_ice * rice(1,l) * (1.+nuice_ref))
-               tauTES=tauTES+ opTES(1,l)
-             enddo
-             CALL WRITEDIAGFI(ngrid,'tauTESap',
-     &                         'tau abs 825 cm-1',
-     &                         '',0,tauTES)
-           else
-
-             CALL WRITEDIAGFI(ngrid,'tauTES',
-     &                         'tau abs 825 cm-1',
-     &                         '',0,taucloudtes)
-           endif
-     
-           mtot = 0
-           icetot = 0
-           h2otot = qsurf(1,igcm_h2o_ice,iflat)
-           if (hdo) THEN
-           mtotD = 0
-           icetotD = 0
-           hdotot = qsurf(1,igcm_hdo_ice,iflat)
-           ENDIF !hdo
-
-           do l=1,nlayer
-             mtot = mtot +  zq(1,l,igcm_h2o_vap) 
-     &                 * (zplev(1,l) - zplev(1,l+1)) / g
-             icetot = icetot +  zq(1,l,igcm_h2o_ice) 
-     &                 * (zplev(1,l) - zplev(1,l+1)) / g
-               if (hdo) THEN
-                 mtotD = mtotD +  zq(1,l,igcm_hdo_vap)
-     &                 * (zplev(1,l) - zplev(1,l+1)) / g
-                 icetotD = icetotD +  zq(1,l,igcm_hdo_ice)
-     &                 * (zplev(1,l) - zplev(1,l+1)) / g
-               ENDIF !hdo
-           end do
-           h2otot = h2otot+mtot+icetot
-               IF (hdo) then
-                   hdotot = hdotot+mtotD+icetotD
-               ENDIF ! hdo
-
-
-             CALL WRITEDIAGFI(ngrid,'h2otot',
-     &                         'h2otot',
-     &                         'kg/m2',0,h2otot)
-             CALL WRITEDIAGFI(ngrid,'mtot',
-     &                         'mtot',
-     &                         'kg/m2',0,mtot)
-             CALL WRITEDIAGFI(ngrid,'icetot',
-     &                         'icetot',
-     &                         'kg/m2',0,icetot)
-             CALL WRITEDIAGFI(ngrid,'albedo',
-     &                         'albedo',
-     &                         '',2,albedo(1,1,iflat))
-         do islope=1,nslope
-          write(str2(1:2),'(i2.2)') islope
-             CALL WRITEDIAGFI(ngrid,'albedo_slope'//str2,
-     &                         'albedo',
-     &                         '',2,albedo(1,1,islope))
-         ENDDO
-             IF (hdo) THEN
-             CALL WRITEDIAGFI(ngrid,'mtotD',
-     &                         'mtotD',
-     &                         'kg/m2',0,mtotD)
-             CALL WRITEDIAGFI(ngrid,'icetotD',
-     &                         'icetotD',
-     &                         'kg/m2',0,icetotD)
-             CALL WRITEDIAGFI(ngrid,'hdotot',
-     &                         'hdotot',
-     &                         'kg/m2',0,hdotot)
-
-C           Calculation of the D/H ratio
-            do l=1,nlayer
-                if (zq(1,l,igcm_h2o_vap).gt.qperemin) then
-                    DoH_vap(1,l) = 0.5*( zq(1,l,igcm_hdo_vap)/
-     &                zq(1,l,igcm_h2o_vap) )/155.76e-6
-                else
-                    DoH_vap(1,l) = 0.
-                endif
-            enddo
-
-            do l=1,nlayer
-                if (zq(1,l,igcm_h2o_ice).gt.qperemin) then
-                    DoH_ice(1,l) = 0.5*( zq(1,l,igcm_hdo_ice)/
-     &                  zq(1,l,igcm_h2o_ice) )/155.76e-6
-                else
-                    DoH_ice(1,l) = 0.
-                endif
-            enddo
-
-            CALL WRITEDIAGFI(ngrid,'DoH_vap',
-     &                       'D/H ratio in vapor',
-     &                       ' ',1,DoH_vap)
-            CALL WRITEDIAGFI(ngrid,'DoH_ice',
-     &                       'D/H ratio in ice',
-     &                       '',1,DoH_ice)
-
-             ENDIF !Hdo
-
-
-           if (scavenging) then
-
-             rave = 0
-             do l=1,nlayer
-cccc Column integrated effective ice radius
-cccc is weighted by total ice surface area (BETTER) 
-             rave = rave + tauscaling(1) *
-     &              zq(1,l,igcm_ccn_number) *
-     &              (zplev(1,l) - zplev(1,l+1)) / g * 
-     &              rice(1,l) * rice(1,l)*  (1.+nuice_ref)
-             enddo
-             rave=icetot*0.75/max(rave*pi*rho_ice,1.e-30) ! surface weight
-
-              Nccntot= 0
-              call watersat(ngrid*nlayer,zt,zplay,zqsat)
-               do l=1,nlayer
-                   Nccntot = Nccntot + 
-     &              zq(1,l,igcm_ccn_number)*tauscaling(1)
-     &              *(zplev(1,l) - zplev(1,l+1)) / g
-                   satu(1,l) = zq(1,l,igcm_h2o_vap)/zqsat(1,l)
-                   satu(1,l) = (max(satu(1,l),float(1))-1)
-!     &                      * zq(1,l,igcm_h2o_vap) *
-!     &                      (zplev(1,l) - zplev(1,l+1)) / g
-               enddo
-               call WRITEDIAGFI(ngrid,"satu","vap in satu","kg/kg",1,
-     &                    satu)
-               CALL WRITEDIAGFI(ngrid,'Nccntot',
-     &                         'Nccntot',
-     &                         'nbr/m2',0,Nccntot)
-
-             call WRITEDIAGFI(ngrid,'zdqsed_dustq'
-     & ,'sedimentation q','kg.m-2.s-1',1,zdqsed(1,:,igcm_dust_mass))
-             call WRITEDIAGFI(ngrid,'zdqsed_dustN'
-     &,'sedimentation N','Nbr.m-2.s-1',1,
-     &                                 zdqsed(1,:,igcm_dust_number))
-
-           else ! of if (scavenging)
-
-cccc Column integrated effective ice radius 
-cccc is weighted by total ice mass         (LESS GOOD)
-             rave = 0
-             do l=1,nlayer
-               rave = rave + zq(1,l,igcm_h2o_ice)
-     &              * (zplev(1,l) - zplev(1,l+1)) / g
-     &              *  rice(1,l) * (1.+nuice_ref)
-             enddo
-             rave=max(rave/max(icetot,1.e-30),1.e-30) ! mass weight
-           endif ! of if (scavenging)
-
-
-           CALL WRITEDIAGFI(ngrid,'reffice',
-     &                      'reffice',
-     &                      'm',0,rave)
-
-          !Alternative A. Pottier weighting
-           rave2 = 0.
-           totrave2 = 0.
-           do l=1,nlayer
-              rave2 =rave2+ zq(1,l,igcm_h2o_ice)*rice(1,l)
-              totrave2 = totrave2 + zq(1,l,igcm_h2o_ice)
-           end do
-           rave2=max(rave2/max(totrave2,1.e-30),1.e-30)
-          CALL WRITEDIAGFI(ngrid,'rmoym',
-     &                      'reffice',
-     &                      'm',0,rave2)
-
-           do iq=1,nq
-             call WRITEDIAGFI(ngrid,trim(noms(iq))//'_s',
-     &       trim(noms(iq))//'_s','kg/kg',0,qsurf(1,iq,iflat))
-           end do
-     
-        call WRITEDIAGFI(ngrid,"watercap","Water ice thickness"
-     &                                  ,"kg.m-2",0,watercap)
+          endif !(rdstorm)
+
+         if(water) then
+           if (.not.scavenging) then
         call WRITEDIAGFI(ngrid,'zdqcloud_ice','cloud ice',
      &            'kg.m-2.s-1',1,zdqcloud(1,:,igcm_h2o_ice))
@@ -4266,38 +3829,26 @@
         call WRITEDIAGFI(ngrid,'zdqcloud_vapD','cloud vap hdo',
      &            'kg.m-2.s-1',1,zdqcloud(1,:,igcm_hdo_vap))
-
         ENDIF ! hdo
-
-        call WRITEDIAGFI(ngrid,"rice","ice radius","m",1,
-     &                    rice)
-              
-              if (CLFvarying) then
-                 call WRITEDIAGFI(ngrid,'totcloudfrac',
-     &                'Total cloud fraction',
-     &                       ' ',0,totcloudfrac)
-              endif !clfvarying
-
+           endif !not.scavenging
         ENDIF ! of IF (water)
-         
-ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-
-         zlocal(1)=-log(zplay(1,1)/zplev(1,1))* Rnew(1,1)*zt(1,1)/g
-
-         do l=2,nlayer-1
-            tmean=zt(1,l)
-            if(zt(1,l).ne.zt(1,l-1))
-     &        tmean=(zt(1,l)-zt(1,l-1))/log(zt(1,l)/zt(1,l-1))
-              zlocal(l)= zlocal(l-1)
-     &        -log(zplay(1,l)/zplay(1,l-1))*rnew(1,l)*tmean/g
-         enddo
-         zlocal(nlayer)= zlocal(nlayer-1)-
-     &                   log(zplay(1,nlayer)/zplay(1,nlayer-1))*
-     &                   rnew(1,nlayer)*tmean/g
-
-
+!PREVIOUSLY IN 1D ONLY
+
+c        ==========================================================
+c        END OF WRITEDIAGFI
+c        ==========================================================
 #endif
-
-      END IF       ! if(ngrid.ne.1)
+! of ifdef MESOSCALE
+
+c      ELSE     ! if(ngrid.eq.1)
+
+c#ifndef MESOSCALE
+c         write(*,
+c     &    '("Ls =",f11.6," tau_pref_scenario(",f4.0," Pa) =",f9.6)')
+c     &    zls*180./pi,odpref,tau_pref_scenario
+c#endif
+
+c      END IF       ! if(ngrid.ne.1)
+
+
 ! test for co2 conservation with co2 microphysics
       if (igcm_co2_ice.ne.0) then
Index: trunk/LMDZ.MARS/libf/phymars/topmons_mod.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/topmons_mod.F90	(revision 2930)
+++ trunk/LMDZ.MARS/libf/phymars/topmons_mod.F90	(revision 2931)
@@ -642,48 +642,10 @@
 ! WRITEDIAGFI
 ! **********************************************************************
-        IF (ngrid.eq.1) THEN
-               CALL WRITEDIAGFI(ngrid,'wup_top', &
-                'wup_top','',1,wup(:,:))
-               CALL WRITEDIAGFI(ngrid,'wrad_top', &
-                'wrad_top','',1,wrad(:,:))
-               CALL WRITEDIAGFI(ngrid,'wfin_top', &
-                'wfin_top','',1,wfin(:,:))
-               CALL WRITEDIAGFI(ngrid,'wlwmax_top', &
-                'wlwmax_top','',0,wup(:,lwmax(1)))
-               CALL WRITEDIAGFI(ngrid,'w0_top', &
-                'w0_top','',0,wup(:,lsummit(1)+1))
-               CALL WRITEDIAGFI(ngrid,'alpha_hmons', &
-                'alpha_hmons','',0,alpha_hmons)
-               CALL WRITEDIAGFI(ngrid,'zt_top', &
-                'zt_top','',1,t_top(:,:))
-               CALL WRITEDIAGFI(ngrid,'dt_top', &
-                'dt_top','',0,dt_top(:))
-               CALL WRITEDIAGFI(ngrid,'envtemp', &
-                'envtemp','',1,zt(:,:))
-               CALL WRITEDIAGFI(ngrid,'t_env', &
-                't_env','',1,t_env(:,:))
-               CALL WRITEDIAGFI(ngrid,'newzt_top', &
-                'newzt_top','',1,newzt(:,:))
-               CALL WRITEDIAGFI(ngrid,'deltahr_top', &
-                'deltahr_top','',1,deltahr(:,:))
-               CALL WRITEDIAGFI(ngrid,'rholev', &
-                'rholev','',1,rho(:,:))
-               CALL WRITEDIAGFI(ngrid,'rhobarz', &
-                'rhobarz','',1,rhobarz(:,:))
-               CALL WRITEDIAGFI(ngrid,'zlsummit', &
-                'zlsummit','',0,pzlay(:,lsummit(1)))
-               CALL WRITEDIAGFI(ngrid,'zlwmax', &
-                'zlwmax','',0,pzlay(:,lwmax(1)))
-               CALL WRITEDIAGFI(ngrid,'pzlev_top', &
-                'pzlev_top','',1,pzlev(:,:))
-               CALL WRITEDIAGFI(ngrid,'coefdetrain', &
-                'coefdetrain','',1,coefdetrain(:,:))
-               CALL WRITEDIAGFI(ngrid,'zdtvert', &
-                'zdtvert','',1,zdtvert(:,:))
-               CALL WRITEDIAGFI(ngrid,'entr', &
-                'entr','',0,entr(:))
-        ELSE
                CALL WRITEDIAGFI(ngrid,'wup_top', &
                 'wup_top','',3,wup(:,:))
+               CALL WRITEDIAGFI(ngrid,'wlwmax_top', &
+                'wlwmax_top','',2,wup(:,lwmax(1)))
+               CALL WRITEDIAGFI(ngrid,'w0_top', &
+                'w0_top','',2,wup(:,lsummit(1)+1))
                CALL WRITEDIAGFI(ngrid,'wup2_top', &
                 'wup2_top','',3,wup2(:,:))
@@ -700,4 +662,6 @@
                CALL WRITEDIAGFI(ngrid,'zt_top', &
                 'zt_top','',3,t_top(:,:))
+               CALL WRITEDIAGFI(ngrid,'dt_top', &
+                'dt_top','',2,dt_top(:))
                CALL WRITEDIAGFI(ngrid,'ztemp', &
                 'envtemp','',3,zt(:,:))
@@ -719,6 +683,12 @@
                               'detrain_rate', &
                               '',3,detrain_rate(:,:))
-
-        ENDIF
+               CALL WRITEDIAGFI(ngrid,'rholev', &
+                'rholev','',3,rho(:,:))
+               CALL WRITEDIAGFI(ngrid,'zlsummit', &
+                'zlsummit','',2,pzlay(:,lsummit(1)))
+               CALL WRITEDIAGFI(ngrid,'pzlev_top', &
+                'pzlev_top','',3,pzlev(:,:))
+               CALL WRITEDIAGFI(ngrid,'entr', &
+                'entr','',2,entr(:))
         
         END SUBROUTINE topmons
