Index: trunk/LMDZ.MARS/README
===================================================================
--- trunk/LMDZ.MARS/README	(revision 2906)
+++ trunk/LMDZ.MARS/README	(revision 2907)
@@ -3912,2 +3912,14 @@
 and
 - "at least one opacity is not a missing value but is 0" (in which case we set tau=0)
+
+== 28/02/2023 == RV
+Further implementation of subslope parametrisation.
+Carefull ! This is a devolpment revision and it still need improvements and tests. 
+However, this commit should not change anything for nslope=1. Only nslope=1 is possible for now!
+Sloped variables in the startfi are outputed along the nslope dimension.
+Possibility to output variables on a specific subslope in diagfi using VAR_slopeXX in the diagfi.def,
+with VAR the variable name (ex: tsurf, co2ice...) and XX the slope number (ex: 04, 07...).
+Without any specific mention to slope, variable named VAR in the diagfi.def will correspond to the variable in the flat slope, this can change in the future.
+Other code changes for nslope.gt.1 (sometimes the grid mesh average is used instead of the value of the subslope)
+
+
Index: trunk/LMDZ.MARS/libf/phymars/iostart.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/iostart.F90	(revision 2906)
+++ trunk/LMDZ.MARS/libf/phymars/iostart.F90	(revision 2907)
@@ -692,9 +692,19 @@
             ierr=NF90_REDEF(nid_restart)
 #ifdef NC_DOUBLE
+            if(field_name.eq. "tauscaling" .or. field_name.eq. "totcloudfrac" .or. field_name.eq. "wstar" ) then
             ierr=NF90_DEF_VAR(nid_restart,field_name,NF90_DOUBLE,&
                               (/idim2,idim7/),nvarid)
-#else
+            else
+            ierr=NF90_DEF_VAR(nid_restart,field_name,NF90_DOUBLE,&
+                              (/idim2,idim8,idim7/),nvarid)
+            endif
+#else
+            if(field_name.eq. "tauscaling" .or. field_name.eq. "totcloudfrac" .or. field_name.eq. "wstar" ) then
             ierr=NF90_DEF_VAR(nid_restart,field_name,NF90_FLOAT,&
                               (/idim2,idim7/),nvarid)
+            else
+            ierr=NF90_DEF_VAR(nid_restart,field_name,NF90_DOUBLE,&
+                              (/idim2,idim8,idim7/),nvarid)
+            endif
 #endif
             if (ierr.ne.NF90_NOERR) then
Index: trunk/LMDZ.MARS/libf/phymars/physiq_mod.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/physiq_mod.F	(revision 2906)
+++ trunk/LMDZ.MARS/libf/phymars/physiq_mod.F	(revision 2907)
@@ -815,4 +815,7 @@
       dsotop(:,:)=0.
       dwatercap(:,:)=0
+
+      call compute_meshgridavg(ngrid,nq,albedo,emis,tsurf,qsurf,
+     &       albedo_meshavg,emis_meshavg,tsurf_meshavg,qsurf_meshavg)
       
       ! Dust scenario conversion coefficient from IRabs to VISext
@@ -1029,6 +1032,6 @@
            clearsky=.false. ! part with clouds for both cases CLFvarying true/false
            CALL callradite(icount,ngrid,nlayer,nq,zday,zls,pq,
-     &     albedo(:,:,1),
-     &     emis(:,1),mu0,zplev,zplay,pt,tsurf(:,1),fract,dist_sol,igout,
+     &     albedo_meshavg,emis_meshavg,
+     &     mu0,zplev,zplay,pt,tsurf(:,1),fract,dist_sol,igout,
      &     zdtlw,zdtsw,fluxsurf_lw(:,iflat),fluxsurf_dn_sw(:,:,iflat),
      &     fluxsurf_up_sw,
@@ -1037,5 +1040,5 @@
      &     tau,aerosol,dsodust,tauscaling,dust_rad_adjust,IRtoVIScoef,
      &     taucloudtes,rdust,rice,nuice,riceco2,nuiceco2,
-     &     qsurf(:,igcm_co2,1),rstormdust,rtopdust,totstormfract,
+     &     qsurf_meshavg(:,igcm_co2),rstormdust,rtopdust,totstormfract,
      &     clearatm,dsords,dsotop,nohmons,clearsky,totcloudfrac)
 
@@ -1053,5 +1056,5 @@
                clearsky=.true.  
                CALL callradite(icount,ngrid,nlayer,nq,zday,zls,pq,
-     &              albedo(:,:,1),emis(:,1),mu0,zplev,zplay,pt,
+     &              albedo_meshavg,emis_meshavg,mu0,zplev,zplay,pt,
      &              tsurf(:,1),fract,
      &              dist_sol,igout,zdtlwclf,zdtswclf,
@@ -1061,5 +1064,6 @@
      &              dsodust,tauscaling,dust_rad_adjust,IRtoVIScoef,
      &              taucloudtesclf,rdust,
-     &              rice,nuice,riceco2, nuiceco2,qsurf(:,igcm_co2,1),
+     &              rice,nuice,riceco2, nuiceco2,
+     &              qsurf_meshavg(:,igcm_co2),
      &              rstormdust,rtopdust,totstormfract,
      &              clearatm,dsords,dsotop,
@@ -1313,5 +1317,5 @@
 c               for radiative transfer
      &                      clearatm,icount,zday,zls,
-     &                      tsurf,qsurf(:,igcm_co2,iflat),igout,
+     &                      tsurf,qsurf_meshavg(:,igcm_co2),igout,
      &                      totstormfract,tauscaling,
      &                      dust_rad_adjust,IRtoVIScoef,
@@ -1377,5 +1381,5 @@
      &                zzlay,zdtsw,zdtlw,
      &                icount,zday,zls,tsurf(:,iflat),
-     &                qsurf(:,igcm_co2,iflat), 
+     &                qsurf_meshavg(:,igcm_co2), 
      &                igout,aerosol,tauscaling,dust_rad_adjust,
      &                IRtoVIScoef,totstormfract,clearatm,
@@ -1500,6 +1504,6 @@
              dwatercap(ig,:)=dwatercap(ig,:)+dwatercap_dif(ig,:)
           ENDDO
-c          call compute_meshgridavg(ngrid,nq,albedo,emis,tsurf,zdqsdif,
-c     &   albedo_meshavg,emis_meshavg,tsurf_meshavg,zdqsdif_meshavg_tmp)
+          call compute_meshgridavg(ngrid,nq,albedo,emis,tsurf,zdqsdif,
+     &   albedo_meshavg,emis_meshavg,tsurf_meshavg,zdqsdif_meshavg_tmp)
          IF (.not.turb_resolved) THEN
           DO l=1,nlayer
@@ -1539,7 +1543,7 @@
              zdqdif(1:ngrid,1,1:nq)=0.
              zdqdif(1:ngrid,1,igcm_dust_number) = 
-     .                  -zdqsdif(1:ngrid,igcm_dust_number,iflat)
+     .        -zdqsdif_meshavg_tmp(1:ngrid,igcm_dust_number)
              zdqdif(1:ngrid,1,igcm_dust_mass) = 
-     .                  -zdqsdif(1:ngrid,igcm_dust_mass,iflat)
+     .        -zdqsdif_meshavg_tmp(1:ngrid,igcm_dust_mass)
              zdqdif(1:ngrid,2:nlayer,1:nq) = 0.
              DO iq=1, nq
@@ -2605,8 +2609,7 @@
           call physdem1("restartfi.nc",nsoilmx,ngrid,nlayer,nq,
      .                ptimestep,ztime_fin,
-     .                tsurf(:,iflat),tsoil(:,:,iflat),albedo(:,:,iflat),
-     .                emis(:,iflat),
-     .                q2,qsurf(:,:,iflat),tauscaling,totcloudfrac,wstar,
-     .                watercap(:,iflat))
+     .                tsurf,tsoil,albedo,
+     .                emis,q2,qsurf,tauscaling,totcloudfrac,wstar,
+     .                watercap)
           
          ENDIF ! of IF (write_restart)
@@ -3155,4 +3158,9 @@
          call WRITEDIAGFI(ngrid,"emis","Surface emissivity","w.m-1",2,
      &                  emis(:,iflat))
+         do islope=1,nslope
+           write(str2(1:2),'(i2.2)') islope
+           call WRITEDIAGFI(ngrid,"emis_slope"//str2,
+     &    "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)
@@ -3165,10 +3173,26 @@
          call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",2,
      &                  tsurf(:,iflat))
+         do islope=1,nslope
+          write(str2(1:2),'(i2.2)') islope
+         call WRITEDIAGFI(ngrid,"tsurf_slope"//str2,
+     &            "Surface temperature","K",2,
+     &                  tsurf(:,islope))
+         ENDDO
          call WRITEDIAGFI(ngrid,"ps","surface pressure","Pa",2,ps)
          call WRITEDIAGFI(ngrid,"co2ice","co2 ice thickness"
      &                              ,"kg.m-2",2,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",2,qsurf(:,igcm_co2,islope))
+         ENDDO
          call WRITEDIAGFI(ngrid,"watercap","Water ice thickness"
      &         ,"kg.m-2",2,watercap(:,iflat))
-
+         do islope=1,nslope
+          write(str2(1:2),'(i2.2)') islope
+         call WRITEDIAGFI(ngrid,"watercap_slope"//str2,
+     &         "Water ice thickness"
+     &         ,"kg.m-2",2,watercap(:,islope))
+         ENDDO
          call WRITEDIAGFI(ngrid,"temp_layer1","temperature in layer 1",
      &                  "K",2,zt(1,1))
@@ -3177,6 +3201,18 @@
          call WRITEDIAGFI(ngrid,"fluxsurf_lw","fluxsurf_lw","W.m-2",2,
      &                  fluxsurf_lw(:,iflat))
+         do islope=1,nslope
+          write(str2(1:2),'(i2.2)') islope
+         call WRITEDIAGFI(ngrid,"fluxsurf_lw_slope"//str2,
+     &              "fluxsurf_lw","W.m-2",2,
+     &                  fluxsurf_lw(:,islope))
+         ENDDO
          call WRITEDIAGFI(ngrid,"fluxsurf_dn_sw","fluxsurf_dn_sw",
      &                  "W.m-2",2,fluxsurf_dn_sw_tot(:,iflat))
+         do islope=1,nslope
+          write(str2(1:2),'(i2.2)') islope
+         call WRITEDIAGFI(ngrid,"fluxsurf_dn_sw_slope"//str2,
+     &                  "fluxsurf_dn_sw",
+     &                  "W.m-2",2,fluxsurf_dn_sw_tot(:,islope))
+         ENDDO
          call WRITEDIAGFI(ngrid,"fluxtop_lw","fluxtop_lw","W.m-2",2,
      &                  fluxtop_lw)
@@ -3306,4 +3342,10 @@
      &                     'kg.m-2',2,
      &                     qsurf(1:ngrid,igcm_h2o_ice,iflat))
+         do islope=1,nslope
+          write(str2(1:2),'(i2.2)') islope
+          call WRITEDIAGFI(ngrid,'qsurf02_slope'//str2,
+     &         'surface tracer','kg.m-2',2,
+     &                     qsurf(1:ngrid,igcm_h2o_ice,islope))
+         ENDDO
 #endif
           call WRITEDIAGFI(ngrid,'mtot',
@@ -3400,9 +3442,21 @@
             call WRITEDIAGFI(ngrid,'h2o_ice_s',
      &                       'surface h2o_ice',
-     &                       'kg.m-2',2,qsurf(1,igcm_h2o_ice,iflat))
+     &            'kg.m-2',2,qsurf(1,igcm_h2o_ice,iflat))
+         do islope=1,nslope
+          write(str2(1:2),'(i2.2)') islope
+            call WRITEDIAGFI(ngrid,'h2o_ice_s_slope'//str2,
+     &                       'surface h2o_ice',
+     &             'kg.m-2',2,qsurf(1,igcm_h2o_ice,islope))
+         ENDDO
             if (hdo) then
             call WRITEDIAGFI(ngrid,'hdo_ice_s',
      &                       'surface hdo_ice',
-     &                       'kg.m-2',2,qsurf(1,igcm_hdo_ice,iflat))
+     &           'kg.m-2',2,qsurf(1,igcm_hdo_ice,iflat))
+         do islope=1,nslope
+          write(str2(1:2),'(i2.2)') islope
+            call WRITEDIAGFI(ngrid,'hdo_ice_s_slope'//str2,
+     &                       'surface hdo_ice',
+     &           'kg.m-2',2,qsurf(1,igcm_hdo_ice,islope))
+         ENDDO
 
                 do ig=1,ngrid
@@ -3423,4 +3477,10 @@
      &                         '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 (tifeedback) then
                  call WRITEDIAGSOIL(ngrid,"soiltemp",
@@ -3429,5 +3489,11 @@
                  call WRITEDIAGSOIL(ngrid,'soilti',
      &                       'Soil Thermal Inertia',
-     &                       'J.s-1/2.m-2.K-1',3,inertiesoil(:,:,iflat))
+     &         'J.s-1/2.m-2.K-1',3,inertiesoil(:,:,iflat))
+         do islope=1,nslope
+          write(str2(1:2),'(i2.2)') islope
+                 call WRITEDIAGSOIL(ngrid,'soilti_slope'//str2,
+     &                       'Soil Thermal Inertia',
+     &          'J.s-1/2.m-2.K-1',3,inertiesoil(:,:,islope))
+         ENDDO
               endif
 !A. Pottier
@@ -3528,4 +3594,9 @@
                call WRITEDIAGFI(ngrid,'qsurf'//str2,'qsurf',
      &                         'kg.m-2',2,qsurf(1,iq,iflat))
+         do islope=1,nslope
+          write(str2(1:2),'(i2.2)') islope
+               call WRITEDIAGFI(ngrid,'qsurf_slope'//str2,'qsurf',
+     &                         'kg.m-2',2,qsurf(1,iq,islope))
+         ENDDO
              end do
            endif ! (doubleq)
@@ -3555,7 +3626,19 @@
      &                 'stormdust injection',
      &                 'kg.m-2',2,qsurf(:,igcm_stormdust_mass,iflat))
+         do islope=1,nslope
+          write(str2(1:2),'(i2.2)') islope
+             call WRITEDIAGFI(ngrid,'qsurf_slope'//str2,
+     &                 'stormdust injection',
+     &          'kg.m-2',2,qsurf(:,igcm_stormdust_mass,islope))
+         ENDDO
              call WRITEDIAGFI(ngrid,'pdqsurf',
      &                  'tendancy stormdust mass at surface',
-     &                  'kg.m-2',2,dqsurf(:,igcm_stormdust_mass,iflat))
+     &          'kg.m-2',2,dqsurf(:,igcm_stormdust_mass,iflat))
+         do islope=1,nslope
+          write(str2(1:2),'(i2.2)') islope
+             call WRITEDIAGFI(ngrid,'pdqsurf_slope'//str2,
+     &                  'tendancy stormdust mass at surface',
+     &            'kg.m-2',2,dqsurf(:,igcm_stormdust_mass,islope))
+         ENDDO
              call WRITEDIAGFI(ngrid,'wspeed','vertical speed stormdust',
      &                        'm/s',3,wspeed(:,1:nlayer))
@@ -3620,6 +3703,18 @@
              call WRITEDIAGFI(ngrid,'surfccnq','Surf nuclei mass mr',
      &                        'kg.m-2',2,qsurf(1,igcm_ccn_mass,iflat))
+         do islope=1,nslope
+          write(str2(1:2),'(i2.2)') islope
+             call WRITEDIAGFI(ngrid,'surfccnq_slope'//str2,
+     &               'Surf nuclei mass mr',
+     &               'kg.m-2',2,qsurf(1,igcm_ccn_mass,islope))
+         ENDDO
              call WRITEDIAGFI(ngrid,'surfccnN','Surf nuclei number',
      &                        'kg.m-2',2,qsurf(1,igcm_ccn_number,iflat))
+         do islope=1,nslope
+          write(str2(1:2),'(i2.2)') islope
+             call WRITEDIAGFI(ngrid,'surfccnN_slope'//str2,
+     &                 'Surf nuclei number',
+     &                 'kg.m-2',2,qsurf(1,igcm_ccn_number,islope))
+         ENDDO
            endif ! (scavenging)
 
@@ -3772,4 +3867,10 @@
         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
          ! Write surface temperature
 !        call writediagsoil(ngrid,"tsurf","Surface temperature","K",
@@ -3802,4 +3903,10 @@
         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
@@ -3826,4 +3933,10 @@
          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)
@@ -3840,4 +3953,10 @@
          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
@@ -3878,6 +3997,18 @@
              call WRITEDIAGFI(ngrid,'dqsdifdustq','diffusion',
      &          'kg.m-2.s-1',0,zdqsdif(1,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))
+         ENDDO
              call WRITEDIAGFI(ngrid,'dqsdifrdsq','diffusion',
      &          'kg.m-2.s-1',0,zdqsdif(1,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))
+         ENDDO
              call WRITEDIAGFI(ngrid,'mstormdtot',
      &                        'total mass of stormdust only',
@@ -3903,7 +4034,19 @@
      &                 '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)
@@ -3996,5 +4139,10 @@
      &                         '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',
