Changeset 2900 for trunk/LMDZ.MARS/libf

Timestamp:

Feb 21, 2023, 5:00:24 PM (21 months ago)

Author:

romain.vande

Message:

Mars PCM:
Following r2896, 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!
Utilitaries are not adapted yet.
Dimension of some variables (30 listed below) is changed to add the nslope dimension.
Outputs (diagfi and restartfi) are not changed yet.
nslope is now read and written in the startfi.nc

Changed variables:
_ In surfdat_h:

• qsurf
• tsurf
• watercap
• emis
• capcal
• fluxgrd

_ In comsoil_h

• tsoil
• mthermdiff
• thermdiff
• coefd
• alph
• beta

_ In dimradmars_mod

• albedo
• fluxrad_sky
• fluxrad

_ In physiq_mod

• inertiesoil
• fluxsurf_lw
• fluxsurf_dn_sw
• dqsurf
• zdtsurf
• zdtsdif
• zdtsurfc
• zdqsdif
• zdqsc
• dwatercap
• dwatercap_dif
• zflubid
• fluxsurf_dn_sw_tot
• inertiedat_tmp

New variables called VAR_meshavg correspond to the mesh average over all the subslope distribution of the variable

Location:

trunk/LMDZ.MARS/libf/phymars

Files:

• comslope_mod.F90 (modified) (2 diffs)
• comsoil_h.F90 (modified) (2 diffs)
• dimradmars_mod.F90 (modified) (4 diffs)
• phys_state_var_init_mod.F90 (modified) (4 diffs)
• physiq_mod.F (modified) (87 diffs)
• soil.F (modified) (9 diffs)
• surfdat_h.F90 (modified) (3 diffs)
• writediagfi.F (modified) (1 diff)

trunk/LMDZ.MARS/libf/phymars/comslope_mod.F90

@@ -13 +13 @@
 IMPLICIT NONE
 
-INTEGER, parameter :: nslope = 1                            ! Number of slopes for the statistic  (1)
+INTEGER, SAVE :: nslope                                ! Number of slopes for the statistic  (1)
 INTEGER, SAVE :: iflat                                      ! Number of slopes for the statistic  (1)
 REAL,SAVE,ALLOCATABLE,DIMENSION(:) :: def_slope             ! bound of the slope statistics / repartitions (°)
@@ -20 +20 @@
 REAL,SAVE,ALLOCATABLE,DIMENSION(:,:) :: subslope_dist       ! distribution  of the slopes (1)
 INTEGER,SAVE,ALLOCATABLE,DIMENSION(:) :: major_slope        ! Index of the subslope that occupies most of themesh (1)
-!$OMP THREADPRIVATE(def_slope,def_slope_mean,sky_slope,subslope_dist,iflat,major_slope)
+!$OMP THREADPRIVATE(nslope,def_slope,def_slope_mean,sky_slope,subslope_dist,iflat,major_slope)
 
 CONTAINS 

trunk/LMDZ.MARS/libf/phymars/comsoil_h.F90

@@ -17 +17 @@
 
   ! variables (FC: built in firstcall in soil.F)
-  REAL,SAVE,ALLOCATABLE :: tsoil(:,:)       ! sub-surface temperatures (K)
-  real,save,allocatable :: mthermdiff(:,:)  ! (FC) mid-layer thermal diffusivity
-  real,save,allocatable :: thermdiff(:,:)   ! (FC) inter-layer thermal diffusivity
+  REAL,SAVE,ALLOCATABLE :: tsoil(:,:,:)       ! sub-surface temperatures (K)
+  real,save,allocatable :: mthermdiff(:,:,:)  ! (FC) mid-layer thermal diffusivity
+  real,save,allocatable :: thermdiff(:,:,:)   ! (FC) inter-layer thermal diffusivity
   real,save,allocatable :: coefq(:)         ! (FC) q_{k+1/2} coefficients
-  real,save,allocatable :: coefd(:,:)       ! (FC) d_k coefficients
-  real,save,allocatable :: alph(:,:)        ! (FC) alpha_k coefficients
-  real,save,allocatable :: beta(:,:)        ! beta_k coefficients
+  real,save,allocatable :: coefd(:,:,:)       ! (FC) d_k coefficients
+  real,save,allocatable :: alph(:,:,:)        ! (FC) alpha_k coefficients
+  real,save,allocatable :: beta(:,:,:)        ! beta_k coefficients
   real,save :: mu
   real,save,allocatable :: flux_geo(:)       ! Geothermal Flux (W/m^2)
@@ -31 +31 @@
 contains
 
-  subroutine ini_comsoil_h(ngrid)
+  subroutine ini_comsoil_h(ngrid,nslope)
   
   implicit none
   integer,intent(in) :: ngrid ! number of atmospheric columns
-  
+  integer,intent(in) :: nslope ! number of sub grid slopes
     allocate(layer(nsoilmx)) !soil layer depths
     allocate(mlayer(0:nsoilmx-1)) ! soil mid-layer depths
     allocate(inertiedat(ngrid,nsoilmx)) ! soil thermal inertia
  
-    allocate(tsoil(ngrid,nsoilmx)) ! soil temperatures
+    allocate(tsoil(ngrid,nsoilmx,nslope)) ! soil temperatures
 
-    allocate(mthermdiff(ngrid,0:nsoilmx-1))
-    allocate(thermdiff(ngrid,nsoilmx-1))
+    allocate(mthermdiff(ngrid,0:nsoilmx-1,nslope))
+    allocate(thermdiff(ngrid,nsoilmx-1,nslope))
     allocate(coefq(0:nsoilmx-1))
-    allocate(coefd(ngrid,nsoilmx-1))
-    allocate(alph(ngrid,nsoilmx-1))
-    allocate(beta(ngrid,nsoilmx-1))
+    allocate(coefd(ngrid,nsoilmx-1,nslope))
+    allocate(alph(ngrid,nsoilmx-1,nslope))
+    allocate(beta(ngrid,nsoilmx-1,nslope))
     allocate(flux_geo(ngrid))
  

trunk/LMDZ.MARS/libf/phymars/dimradmars_mod.F90

@@ -117 +117 @@
 !!
   REAL,SAVE,ALLOCATABLE :: dtrad(:,:) ! Net atm. radiative heating rate (K.s-1)
-  REAL,SAVE,ALLOCATABLE :: fluxrad_sky(:) ! rad. flux from sky absorbed by surface (W.m-2)
-  REAL,SAVE,ALLOCATABLE :: fluxrad(:) ! Net radiative surface flux (W.m-2)
-  REAL,SAVE,ALLOCATABLE :: albedo(:,:) ! Surface albedo in each solar band
+  REAL,SAVE,ALLOCATABLE :: fluxrad_sky(:,:) ! rad. flux from sky absorbed by surface (W.m-2)
+  REAL,SAVE,ALLOCATABLE :: fluxrad(:,:) ! Net radiative surface flux (W.m-2)
+  REAL,SAVE,ALLOCATABLE :: albedo(:,:,:) ! Surface albedo in each solar band
   REAL,SAVE,ALLOCATABLE :: totcloudfrac(:) ! total cloud fraction over the column
 ! aerosol (dust or ice) extinction optical depth  at reference wavelength 
@@ -181 +181 @@
 contains
 
-  subroutine ini_dimradmars_mod(ngrid,nlayer)
+  subroutine ini_dimradmars_mod(ngrid,nlayer,nslope)
   
   implicit none
@@ -187 +187 @@
   integer,intent(in) :: ngrid ! number of atmospheric columns
   integer,intent(in) :: nlayer ! number of atmospheric layers
-
+  integer,intent(in) :: nslope ! number of subgrid scale slopes
    nflev=nlayer
 !  ndomainsz=ngrid
@@ -195 +195 @@
    ndlo2=ndlon
 
-   allocate(albedo(ngrid,2))
+   allocate(albedo(ngrid,2,nslope))
    allocate(dtrad(ngrid,nlayer))
-   allocate(fluxrad_sky(ngrid))
-   allocate(fluxrad(ngrid))
+   allocate(fluxrad_sky(ngrid,nslope))
+   allocate(fluxrad(ngrid,nslope))
    allocate(nueffdust(ngrid,nlayer))
    allocate(totcloudfrac(ngrid))

trunk/LMDZ.MARS/libf/phymars/phys_state_var_init_mod.F90

@@ -61 +61 @@
       use dust_rad_adjust_mod, only: ini_dust_rad_adjust_mod, &
                                      end_dust_rad_adjust_mod
+      use comslope_mod, ONLY: nslope
+      use netcdf
+      USE mod_phys_lmdz_para, ONLY: is_master,bcast
+
       IMPLICIT NONE
       
@@ -71 +75 @@
       INTEGER,INTENT(in) :: dyn_nqperes
       INTEGER,INTENT(in) :: dyn_nqfils(nq)
+      character(len=10) :: filename  ! name of the startfi.nc
+      integer :: ncid, status, nslope_dim_id
+      integer :: nslope_read
+
+      filename="startfi.nc"
+      if(is_master) then
+        status = nf90_open(filename, nf90_nowrite, ncid)
+        if (status /= nf90_noerr) then
+          nslope=1
+        else
+          status = nf90_inq_dimid(ncid, "nslope", nslope_dim_id)
+          if (status /= nf90_noerr) then
+            nslope=1
+          else
+            status = nf90_inquire_dimension(ncid, nslope_dim_id, len = nslope_read)
+            if (status /= nf90_noerr) then
+              call abort_physic("phys_state_var_init","nslope present but not readable",1)
+            else
+              nslope=nslope_read
+            endif
+          endif
+        endif
+      endif
+      call bcast(nslope)
 
       ! set dimension and allocate arrays in tracer_mod
       call end_tracer_mod
       call ini_tracer_mod(nq,tname,dyn_nqperes,dyn_nqfils)! MVals: variables isotopes
-
 
       ! initialize "print_control" constants/flags ("prt_level","lunout", etc.)
@@ -101 +128 @@
       ! allocate "surfdat_h" arrays
       call end_surfdat_h
-      call ini_surfdat_h(ngrid,nq)
+      call ini_surfdat_h(ngrid,nq,nslope)
 
       ! allocate "comgeomfi_h" arrays
@@ -109 +136 @@
       ! allocate "comsoil_h" arrays
       call end_comsoil_h
-      call ini_comsoil_h(ngrid)
+      call ini_comsoil_h(ngrid,nslope)
 
       ! set some variables in "dimradmars_mod"
       call end_dimradmars_mod
-      call ini_dimradmars_mod(ngrid,nlayer)
+      call ini_dimradmars_mod(ngrid,nlayer,nslope)
 
       ! allocate arrays in "yomlw_h"

trunk/LMDZ.MARS/libf/phymars/physiq_mod.F

@@ -271 +271 @@
       real nuiceco2(ngrid,nlayer) ! Estimated effective variance of the 
                                   !   size distribution
-      REAL inertiesoil(ngrid,nsoilmx) ! Time varying subsurface
+      REAL inertiesoil(ngrid,nsoilmx,nslope) ! Time varying subsurface
                                       ! thermal inertia (J.s-1/2.m-2.K-1)
                                       ! (used only when tifeedback=.true.)
@@ -301 +301 @@
       INTEGER l,ig,ierr,igout,iq,tapphys,isoil
 
-      REAL fluxsurf_lw(ngrid)      !incident LW (IR) surface flux (W.m-2)
-      REAL fluxsurf_dn_sw(ngrid,2) ! Incident SW (solar) surface flux (W.m-2)
+      REAL fluxsurf_lw(ngrid,nslope)      !incident LW (IR) surface flux (W.m-2)
+      REAL fluxsurf_dn_sw(ngrid,2,nslope) ! Incident SW (solar) surface flux (W.m-2)
       REAL fluxsurf_up_sw(ngrid,2) ! Reflected SW (solar) surface flux (W.m-2)
       REAL fluxtop_lw(ngrid)       !Outgoing LW (IR) flux to space (W.m-2)
@@ -333 +333 @@
 
 c     Tendancies due to various processes:
-      REAL dqsurf(ngrid,nq)  ! tendency for tracers on surface (Kg/m2/s)
+      REAL dqsurf(ngrid,nq,nslope)  ! tendency for tracers on surface (Kg/m2/s)
       REAL zdtlw(ngrid,nlayer)     ! (K/s)
       REAL zdtsw(ngrid,nlayer)     ! (K/s)
@@ -340 +340 @@
       REAL zdtnirco2(ngrid,nlayer) ! (K/s)
       REAL zdtnlte(ngrid,nlayer)   ! (K/s)
-      REAL zdtsurf(ngrid)            ! (K/s)
+      REAL zdtsurf(ngrid,nslope)            ! (K/s)
       REAL zdtcloud(ngrid,nlayer),zdtcloudco2(ngrid,nlayer)
       REAL zdvdif(ngrid,nlayer),zdudif(ngrid,nlayer)  ! (m.s-2)
-      REAL zdhdif(ngrid,nlayer), zdtsdif(ngrid)         ! (K/s)
+      REAL zdhdif(ngrid,nlayer), zdtsdif(ngrid,nslope)         ! (K/s)
       REAL zdvadj(ngrid,nlayer),zduadj(ngrid,nlayer)  ! (m.s-2)
       REAL zdhadj(ngrid,nlayer)                           ! (K/s)
       REAL zdtgw(ngrid,nlayer)                            ! (K/s)
       REAL zdugw(ngrid,nlayer),zdvgw(ngrid,nlayer)    ! (m.s-2)
-      REAL zdtc(ngrid,nlayer),zdtsurfc(ngrid)
+      REAL zdtc(ngrid,nlayer),zdtsurfc(ngrid,nslope)
       REAL zdvc(ngrid,nlayer),zduc(ngrid,nlayer)
 
-      REAL zdqdif(ngrid,nlayer,nq), zdqsdif(ngrid,nq)
+      REAL zdqdif(ngrid,nlayer,nq), zdqsdif(ngrid,nq,nslope)
       REAL zdqsed(ngrid,nlayer,nq), zdqssed(ngrid,nq)
       REAL zdqdev(ngrid,nlayer,nq), zdqsdev(ngrid,nq)
@@ -357 +357 @@
       REAL zdqc(ngrid,nlayer,nq)
       REAL zdqcloudco2(ngrid,nlayer,nq)
-      REAL zdqsc(ngrid,nq)
+      REAL zdqsc(ngrid,nq,nslope)
 
       REAL zdteuv(ngrid,nlayer)    ! (K/s)
@@ -366 +366 @@
       real*8 PhiEscH,PhiEscH2,PhiEscD
 
-      REAL dwatercap(ngrid), dwatercap_dif(ngrid)     ! (kg/m-2)
+      REAL dwatercap(ngrid,nslope), dwatercap_dif(ngrid,nslope)     ! (kg/m-2)
 
 c     Local variable for local intermediate calcul:
-      REAL zflubid(ngrid)
+      REAL zflubid(ngrid,nslope)
       REAL zplanck(ngrid),zpopsk(ngrid,nlayer)
       REAL zdum1(ngrid,nlayer)
@@ -393 +393 @@
 
       REAL fluxtop_dn_sw_tot(ngrid), fluxtop_up_sw_tot(ngrid)
-      REAL fluxsurf_dn_sw_tot(ngrid), fluxsurf_up_sw_tot(ngrid)
+      REAL fluxsurf_dn_sw_tot(ngrid,nslope), fluxsurf_up_sw_tot(ngrid)
       character*2 str2
 !      character*5 str5
@@ -534 +534 @@
 
 c Sub-grid scale slopes
+      real :: tsurf_meshavg(ngrid) ! Surface temperature grid box averaged [K]
+      real :: albedo_meshavg(ngrid,2) ! albedo temperature grid box averaged [1]
+      real :: emis_meshavg(ngrid,2) ! emis temperature grid box averaged [1]
+      real :: qsurf_meshavg(ngrid,nq) ! surface tracer mesh averaged [kg/m^2]
+      real :: qsurf_tmp(ngrid,nq) ! temporary qsurf for chimie
+      real :: inertiedat_tmp(ngrid,nsoilmx,nslope) ! temporary inertiedat for soil.F
       integer :: islope
+      real :: zdqsdif_meshavg_tmp(ngrid,nq) ! temporary for dust lifting
 
       logical :: write_restart
@@ -557 +564 @@
 
 #ifndef MESOSCALE
-         fluxrad(:)=0
+         fluxrad(:,:)=0
          wstar(:)=0.
 #endif
@@ -600 +607 @@
        PRINT*,'corresponding criterium = ',def_slope_mean(iflat)
 
+       DO islope = 1,nslope
+         inertiedat_tmp(:,:,islope) = inertiedat(:,:)
+       ENDDO
+
 #else
 ! MESOSCALE. Supposedly everything is already set in modules.
@@ -606 +617 @@
       print*,"check: rad,cpp,g,r,rcp,daysec"
       print*,rad,cpp,g,r,rcp,daysec
-      PRINT*,'check: tsurf ',tsurf(1),tsurf(ngrid)
-      PRINT*,'check: tsoil ',tsoil(1,1),tsoil(ngrid,nsoilmx)
+      PRINT*,'check: tsurf ',tsurf(1,:),tsurf(ngrid,:)
+      PRINT*,'check: tsoil ',tsoil(1,1,:),tsoil(ngrid,nsoilmx,:)
       PRINT*,'check: inert ',inertiedat(1,1),inertiedat(ngrid,nsoilmx)
       PRINT*,'check: midlayer,layer ', mlayer(:),layer(:)
       PRINT*,'check: tracernames ', noms
-      PRINT*,'check: emis ',emis(1),emis(ngrid)
+      PRINT*,'check: emis ',emis(1,:),emis(ngrid,:)
       PRINT*,'check: q2 ',q2(1,1),q2(ngrid,nlayer+1)
-      PRINT*,'check: qsurf ',qsurf(1,1),qsurf(ngrid,nq)
-      PRINT*,'check: co2ice ',qsurf(1,igcm_co2),qsurf(ngrid,igcm_co2)
+      PRINT*,'check: qsurf ',qsurf(1,1,:),qsurf(ngrid,nq,:)
+      PRINT*,'check: co2ice ',qsurf(1,igcm_co2,:),qsurf(ngrid,igcm_co2,:)
       !!!
       day_ini = pday
@@ -645 +656 @@
            endif
            ! additionnal "academic" initialization of physics
-           write(*,*) "Physiq: initializing tsurf(:) to pt(:,1) !!"
-           tsurf(:)=pt(:,1)
+           do islope = 1,nslope
+             tsurf(:,islope)=pt(:,1)
+           enddo
            write(*,*) "Physiq: initializing tsoil(:) to pt(:,1) !!"
            do isoil=1,nsoilmx
-             tsoil(1:ngrid,isoil)=tsurf(1:ngrid)
+             tsoil(1:ngrid,isoil,:)=tsurf(1:ngrid,:)
            enddo
            write(*,*) "Physiq: initializing inertiedat !!"
@@ -683 +695 @@
 c           Thermal inertia feedback:
             IF (tifeedback) THEN
-                CALL soil_tifeedback(ngrid,nsoilmx,qsurf,inertiesoil)
+             DO islope = 1,nslope
+                CALL soil_tifeedback(ngrid,nsoilmx,
+     s               qsurf(:,:,islope),inertiesoil(:,:,islope))
+             ENDDO
                 CALL soil(ngrid,nsoilmx,firstcall,inertiesoil,
      s             ptimestep,tsurf,tsoil,capcal,fluxgrd)
             ELSE
-                CALL soil(ngrid,nsoilmx,firstcall,inertiedat,
+                CALL soil(ngrid,nsoilmx,firstcall,inertiedat_tmp,
      s             ptimestep,tsurf,tsoil,capcal,fluxgrd)
             ENDIF ! of IF (tifeedback)
@@ -694 +709 @@
      &     'PHYSIQ WARNING! Thermal conduction in the soil turned off'
             DO ig=1,ngrid
-               capcal(ig)=1.e5
-               fluxgrd(ig)=0.
+               capcal(ig,:)=1.e5
+               fluxgrd(ig,:)=0.
             ENDDO
          ENDIF
@@ -793 +808 @@
       pdq(:,:,:)=0
       pdpsrf(:)=0
-      zflubid(:)=0
-      zdtsurf(:)=0
-      dqsurf(:,:)=0
+      zflubid(:,:)=0
+      zdtsurf(:,:)=0
+      dqsurf(:,:,:)=0
       dsodust(:,:)=0.
       dsords(:,:)=0.
       dsotop(:,:)=0.
-      dwatercap(:)=0
+      dwatercap(:,:)=0
       
       ! Dust scenario conversion coefficient from IRabs to VISext
@@ -899 +914 @@
      &        +(pdq(ig,l,igcm_co2)+pdq(ig,l,igcm_co2_ice))*ptimestep)
           end do
-          co2totA = co2totA + qsurf(ig,igcm_co2)
+          do islope = 1,nslope
+          co2totA = co2totA + qsurf(ig,igcm_co2,islope)*
+     &    subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)
+          enddo
         end do
       else
@@ -909 +927 @@
      &        +pdq(ig,l,igcm_co2)*ptimestep)
           end do
-          co2totA = co2totA + qsurf(ig,igcm_co2)
+          do islope = 1,nslope
+          co2totA = co2totA + qsurf(ig,igcm_co2,islope)*
+     &    subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)
+          enddo
         end do
       endif ! of if (igcm_co2_ice.ne.0)
@@ -1007 +1028 @@
            ! callradite for the part with clouds
            clearsky=.false. ! part with clouds for both cases CLFvarying true/false
-           CALL callradite(icount,ngrid,nlayer,nq,zday,zls,pq,albedo,
-     &     emis,mu0,zplev,zplay,pt,tsurf,fract,dist_sol,igout,
-     &     zdtlw,zdtsw,fluxsurf_lw,fluxsurf_dn_sw,fluxsurf_up_sw,
+           CALL callradite(icount,ngrid,nlayer,nq,zday,zls,pq,
+     &     albedo(:,:,1),
+     &     emis(:,1),mu0,zplev,zplay,pt,tsurf(:,1),fract,dist_sol,igout,
+     &     zdtlw,zdtsw,fluxsurf_lw(:,iflat),fluxsurf_dn_sw(:,:,iflat),
+     &     fluxsurf_up_sw,
      &     fluxtop_lw,fluxtop_dn_sw,fluxtop_up_sw,
      &     tau_pref_scenario,tau_pref_gcm,
      &     tau,aerosol,dsodust,tauscaling,dust_rad_adjust,IRtoVIScoef,
      &     taucloudtes,rdust,rice,nuice,riceco2,nuiceco2,
-     &     qsurf(:,igcm_co2),rstormdust,rtopdust,totstormfract,
+     &     qsurf(:,igcm_co2,1),rstormdust,rtopdust,totstormfract,
      &     clearatm,dsords,dsotop,nohmons,clearsky,totcloudfrac)
+
+           DO islope=1,nslope
+             fluxsurf_lw(:,islope) =fluxsurf_lw(:,iflat)
+             fluxsurf_dn_sw(:,:,islope) =fluxsurf_dn_sw(:,:,iflat)
+           ENDDO
 
            ! case of sub-grid water ice clouds: callradite for the clear case
@@ -1025 +1053 @@
                clearsky=.true.  
                CALL callradite(icount,ngrid,nlayer,nq,zday,zls,pq,
-     &              albedo,emis,mu0,zplev,zplay,pt,tsurf,fract,
+     &              albedo(:,:,1),emis(:,1),mu0,zplev,zplay,pt,
+     &              tsurf(:,1),fract,
      &              dist_sol,igout,zdtlwclf,zdtswclf,
      &              fluxsurf_lwclf,fluxsurf_dn_swclf,fluxsurf_up_swclf,
@@ -1032 +1061 @@
      &              dsodust,tauscaling,dust_rad_adjust,IRtoVIScoef,
      &              taucloudtesclf,rdust,
-     &              rice,nuice,riceco2, nuiceco2,qsurf(:,igcm_co2),
+     &              rice,nuice,riceco2, nuiceco2,qsurf(:,igcm_co2,1),
      &              rstormdust,rtopdust,totstormfract,
      &              clearatm,dsords,dsotop,
@@ -1042 +1071 @@
                   tf_clf=totcloudfrac(ig)
                   ntf_clf=1.-tf_clf
-                  fluxsurf_lw(ig) = ntf_clf*fluxsurf_lwclf(ig) 
-     &                                      + tf_clf*fluxsurf_lw(ig)
-                  fluxsurf_dn_sw(ig,1:2) = 
+                  DO islope=1,nslope
+                  fluxsurf_lw(ig,islope) = ntf_clf*fluxsurf_lwclf(ig) 
+     &                          + tf_clf*fluxsurf_lw(ig,islope)
+                  fluxsurf_dn_sw(ig,1:2,islope) = 
      &                           ntf_clf*fluxsurf_dn_swclf(ig,1:2) 
-     &                          + tf_clf*fluxsurf_dn_sw(ig,1:2)
+     &                          + tf_clf*fluxsurf_dn_sw(ig,1:2,islope)
+                  ENDDO
                   fluxsurf_up_sw(ig,1:2) = 
      &                           ntf_clf*fluxsurf_up_swclf(ig,1:2) 
@@ -1099 +1130 @@
 c          ---------------------------------------------------------
            IF(callslope) THEN
+           ! assume that in this case, nslope = 1
+            if(nslope.ne.1) then
+              call abort_physic(
+     &      "physiq","callslope=true but nslope.ne.1",1)
+            endif
             print *, 'Slope scheme is on and computing...'
             DO ig=1,ngrid  
               sl_the = theta_sl(ig)
               IF (sl_the .ne. 0.) THEN
-                ztim1=fluxsurf_dn_sw(ig,1)+fluxsurf_dn_sw(ig,2)
+                ztim1=fluxsurf_dn_sw(ig,1,1)+fluxsurf_dn_sw(ig,2,1)
                 DO l=1,2
                  sl_lct = ptime*24. + 180.*longitude(ig)/pi/15.
@@ -1109 +1145 @@
                  sl_lat = 180.*latitude(ig)/pi
                  sl_tau = tau(ig,1) !il faudrait iaerdust(iaer)
-                 sl_alb = albedo(ig,l)
+                 sl_alb = albedo(ig,l,1)
                  sl_psi = psi_sl(ig)
-                 sl_fl0 = fluxsurf_dn_sw(ig,l)
+                 sl_fl0 = fluxsurf_dn_sw(ig,l,1)
                  sl_di0 = 0.
                  if ((mu0(ig) .gt. 0.).and.(ztim1.gt.0.)) then
@@ -1117 +1153 @@
                   sl_di0 = sl_di0*flux_1AU/dist_sol/dist_sol
                   sl_di0 = sl_di0/ztim1
-                  sl_di0 = fluxsurf_dn_sw(ig,l)*sl_di0
+                  sl_di0 = fluxsurf_dn_sw(ig,l,1)*sl_di0
                  endif
                  ! you never know (roundup concern...)
@@ -1126 +1162 @@
      &                             sl_di0, sl_fl0, sl_flu )
                  !!!!!!!!!!!!!!!!!!!!!!!!!!
-                 fluxsurf_dn_sw(ig,l) = sl_flu
+                 fluxsurf_dn_sw(ig,l,1) = sl_flu
                 ENDDO
               !!! compute correction on IR flux as well
               sky= (1.+cos(pi*theta_sl(ig)/180.))/2.
-              fluxsurf_lw(ig)= fluxsurf_lw(ig)*sky
+              fluxsurf_lw(ig,:)= fluxsurf_lw(ig,:)*sky
               ENDIF
             ENDDO
-           ENDIF
+c          ELSE        ! not calling subslope, nslope might be > 1
+c          DO islope = 1,nslope
+c            IF(def_slope_mean(islope).ge.0.) THEN
+c              sl_psi = 0. !Northward slope
+c            ELSE
+c             sl_psi = 180. !Southward slope
+c            ENDIF
+c            sl_the=abs(def_slope_mean(islope))         
+c            IF (sl_the .gt. 1e-6) THEN
+c              DO ig=1,ngrid  
+c                ztim1=fluxsurf_dn_sw(ig,1,islope)
+c     s                +fluxsurf_dn_sw(ig,2,islope)
+c                DO l=1,2
+c                 sl_lct = ptime*24. + 180.*longitude(ig)/pi/15.
+c                 sl_ra  = pi*(1.0-sl_lct/12.)
+c                 sl_lat = 180.*latitude(ig)/pi
+c                 sl_tau = tau(ig,1) !il faudrait iaerdust(iaer)
+c                 sl_alb = albedo(ig,l,islope)
+c                 sl_psi = psi_sl(ig)
+c                 sl_fl0 = fluxsurf_dn_sw(ig,l,islope)
+c                 sl_di0 = 0.
+c                 if (mu0(ig) .gt. 0.) then
+c                         sl_di0 = mu0(ig)*(exp(-sl_tau/mu0(ig)))
+c                  sl_di0 = sl_di0*flux_1AU/dist_sol/dist_sol
+c                  sl_di0 = sl_di0/ztim1
+c                  sl_di0 = fluxsurf_dn_sw(ig,l,islope)*sl_di0
+c                 endif
+c                 ! you never know (roundup concern...)
+c                 if (sl_fl0 .lt. sl_di0) sl_di0=sl_fl0
+c                 !!!!!!!!!!!!!!!!!!!!!!!!!!
+c                 CALL param_slope( mu0(ig), declin, sl_ra, sl_lat, 
+c     &                             sl_tau, sl_alb, sl_the, sl_psi,
+c     &                             sl_di0, sl_fl0, sl_flu )
+c                 !!!!!!!!!!!!!!!!!!!!!!!!!!
+c                 fluxsurf_dn_sw(ig,l,islope) = sl_flu
+c                ENDDO
+c              !!! compute correction on IR flux as well
+c
+c              fluxsurf_lw(ig,islope)= fluxsurf_lw(ig,islope)
+c     &                                *sky_slope(islope)
+c              ENDDO
+c            ENDIF
+c          ENDDO ! islope = 1,nslope
+          ENDIF ! callslope
 
 c          CO2 near infrared absorption
@@ -1146 +1225 @@
 c          ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
            DO ig=1,ngrid
-               fluxrad_sky(ig)=emis(ig)*fluxsurf_lw(ig)
-     $         +fluxsurf_dn_sw(ig,1)*(1.-albedo(ig,1))
-     $         +fluxsurf_dn_sw(ig,2)*(1.-albedo(ig,2))
+             DO islope = 1,nslope
+               fluxrad_sky(ig,islope) =
+     $           emis(ig,islope)*fluxsurf_lw(ig,islope)
+     $         +fluxsurf_dn_sw(ig,1,islope)*(1.-albedo(ig,1,islope))         
+     $         +fluxsurf_dn_sw(ig,2,islope)*(1.-albedo(ig,2,islope))
+             ENDDO
            ENDDO
 
@@ -1175 +1257 @@
 c
            DO ig=1,ngrid
-               zplanck(ig)=tsurf(ig)*tsurf(ig)
-               zplanck(ig)=emis(ig)*
+             DO islope = 1,nslope
+               zplanck(ig)=tsurf(ig,islope)*tsurf(ig,islope)
+               zplanck(ig)=emis(ig,islope)*
      $         stephan*zplanck(ig)*zplanck(ig)
-               fluxrad(ig)=fluxrad_sky(ig)-zplanck(ig)
+               fluxrad(ig,islope)=fluxrad_sky(ig,islope)-zplanck(ig)
                IF(callslope) THEN
                  sky= (1.+cos(pi*theta_sl(ig)/180.))/2.
-                 fluxrad(ig)=fluxrad(ig)+(1.-sky)*zplanck(ig)
+                 fluxrad(ig,nslope)=fluxrad(ig,nslope)+
+     $                             (1.-sky)*zplanck(ig)
+               ELSE 
+                 fluxrad(ig,islope)=fluxrad(ig,islope) +
+     $              (1.-sky_slope(iflat))*emis(ig,iflat)*
+     $              stephan*tsurf(ig,iflat)**4  
                ENDIF
            ENDDO
+         ENDDO
 
          DO l=1,nlayer
@@ -1224 +1313 @@
 c               for radiative transfer
      &                      clearatm,icount,zday,zls,
-     &                      tsurf,qsurf(:,igcm_co2),igout,
+     &                      tsurf,qsurf(:,igcm_co2,iflat),igout,
      &                      totstormfract,tauscaling,
      &                      dust_rad_adjust,IRtoVIScoef,
@@ -1287 +1376 @@
      &                pq,pdq,pt,pdt,zplev,zplay,zzlev,
      &                zzlay,zdtsw,zdtlw,
-     &                icount,zday,zls,tsurf,qsurf(:,igcm_co2), 
+     &                icount,zday,zls,tsurf(:,iflat),
+     &                qsurf(:,igcm_co2,iflat), 
      &                igout,aerosol,tauscaling,dust_rad_adjust,
      &                IRtoVIScoef,totstormfract,clearatm,
@@ -1349 +1439 @@
       IF (calldifv) THEN
          DO ig=1,ngrid
-            zflubid(ig)=fluxrad(ig)+fluxgrd(ig)
+           DO islope = 1,nslope
+            zflubid(ig,islope)=fluxrad(ig,islope)
+     &                        +fluxgrd(ig,islope)
+           ENDDO
          ENDDO
-
          zdum1(:,:)=0
          zdum2(:,:)=0
@@ -1371 +1463 @@
 
           DO ig=1, ngrid
-             IF (zh(ig,1) .lt. tsurf(ig)) THEN
+             IF (zh(ig,1) .lt. tsurf(ig,1)) THEN
                wstar(ig)=1.
                hfmax_th(ig)=0.2
@@ -1390 +1482 @@
 
 c        Calling vdif (Martian version WITH CO2 condensation)
-         dwatercap_dif(:) = 0.
+         dwatercap_dif(:,:) = 0.
          zcdh(:) = 0.
          zcdv(:) = 0.
@@ -1405 +1497 @@
 
           DO ig=1,ngrid
-             zdtsurf(ig)=zdtsurf(ig)+zdtsdif(ig)
-             dwatercap(ig)=dwatercap(ig)+dwatercap_dif(ig)
+             zdtsurf(ig,:)=zdtsurf(ig,:)+zdtsdif(ig,:)
+             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)
          IF (.not.turb_resolved) THEN
           DO l=1,nlayer
@@ -1429 +1522 @@
            DO iq=1, nq
               DO ig=1,ngrid
-                 dqsurf(ig,iq)=dqsurf(ig,iq) + zdqsdif(ig,iq)
+                 dqsurf(ig,iq,:)=dqsurf(ig,iq,:) + zdqsdif(ig,iq,:)
               ENDDO
            ENDDO
@@ -1446 +1539 @@
              zdqdif(1:ngrid,1,1:nq)=0.
              zdqdif(1:ngrid,1,igcm_dust_number) = 
-     .                  -zdqsdif(1:ngrid,igcm_dust_number)
+     .                  -zdqsdif(1:ngrid,igcm_dust_number,iflat)
              zdqdif(1:ngrid,1,igcm_dust_mass) = 
-     .                  -zdqsdif(1:ngrid,igcm_dust_mass)
+     .                  -zdqsdif(1:ngrid,igcm_dust_mass,iflat)
              zdqdif(1:ngrid,2:nlayer,1:nq) = 0.
              DO iq=1, nq
                IF ((iq .ne. igcm_dust_mass)
      &          .and. (iq .ne. igcm_dust_number)) THEN
-                 zdqsdif(:,iq)=0.
+                 zdqsdif(:,iq,:)=0.
                ENDIF
              ENDDO
            ELSE
              zdqdif(1:ngrid,1:nlayer,1:nq) = 0.
-             zdqsdif(1:ngrid,1:nq) = 0.
+             zdqsdif(1:ngrid,1:nq,1:nslope) = 0.
            ENDIF
          ENDIF
       ELSE    
          DO ig=1,ngrid
-            zdtsurf(ig)=zdtsurf(ig)+
-     s        (fluxrad(ig)+fluxgrd(ig))/capcal(ig)
+           DO islope=1,nslope
+             zdtsurf(ig,islope)=zdtsurf(ig,islope)+
+     s   (fluxrad(ig,islope)+fluxgrd(ig,islope))/capcal(ig,islope)
+           ENDDO
          ENDDO
          IF (turb_resolved) THEN
@@ -1707 +1802 @@
           DO iq=1, nq
             DO ig=1,ngrid
-              dqsurf(ig,iq)=dqsurf(ig,iq)+zdqssed_ccn(ig,iq)
+             DO islope = 1,nslope
+              dqsurf(ig,iq,islope)=dqsurf(ig,iq,islope)+
+     &         zdqssed_ccn(ig,iq)*cos(pi*def_slope_mean(islope)/180.)
+             ENDDO !(islope)
             ENDDO  ! (ig)
-           ENDDO    ! (iq)
+           ENDDO    ! (iq)q)
 c Temperature variation due to latent heat release
                pdt(1:ngrid,1:nlayer) =
@@ -1884 +1982 @@
               do iq=1,nq
                  DO ig=1,ngrid
-                    dqsurf(ig,iq)= dqsurf(ig,iq) + zdqsdev(ig,iq)
+                   DO islope = 1,nslope
+                    dqsurf(ig,iq,islope)= dqsurf(ig,iq,islope) +
+     &          zdqsdev(ig,iq)*cos(pi*def_slope_mean(islope)/180.)
+                   ENDDO
                  ENDDO
               enddo
@@ -1933 +2034 @@
            DO iq=1, nq
              DO ig=1,ngrid
-                dqsurf(ig,iq)= dqsurf(ig,iq) + zdqssed(ig,iq)
+               DO islope = 1,nslope
+                    dqsurf(ig,iq,islope)= dqsurf(ig,iq,islope) +
+     &          zdqssed(ig,iq)*cos(pi*def_slope_mean(islope)/180.)
+               ENDDO
              ENDDO
            ENDDO
@@ -1950 +2054 @@
            DO iq=1, nq
               DO ig=1,ngrid
-                 dqsurf(ig,iq)=dqsurf(ig,iq) + zdqsdif(ig,iq)
+                 dqsurf(ig,iq,:)=dqsurf(ig,iq,:) + zdqsdif(ig,iq,:)
               ENDDO
            ENDDO
@@ -1974 +2078 @@
      $                       surfdust, surfice)
 !           call photochemistry
+            DO ig = 1,ngrid
+              qsurf_tmp(ig,:) = qsurf(ig,:,major_slope(ig))
+            ENDDO
             call calchim(ngrid,nlayer,nq,
      &                   ptimestep,zplay,zplev,pt,pdt,dist_sol,mu0,
      $                   zzlev,zzlay,zday,pq,pdq,zdqchim,zdqschim,
-     $                   zdqcloud,zdqscloud,tau(:,1),qsurf(:,igcm_co2),
+     $                   zdqcloud,zdqscloud,tau(:,1),
+     $                   qsurf_tmp(:,igcm_co2),
      $                   pu,pdu,pv,pdv,surfdust,surfice)
 
@@ -2006 +2114 @@
                       ! tracers is zero anyways
              DO ig=1,ngrid
-               dqsurf(ig,iq)=dqsurf(ig,iq)+zdqschim(ig,iq)
+              DO islope = 1,nslope
+               dqsurf(ig,iq,islope)=dqsurf(ig,iq,islope) +
+     &           zdqschim(ig,iq)*cos(pi*def_slope_mean(islope)/180.)
+              ENDDO
              ENDDO
            ENDDO ! of DO iq=1,nq
@@ -2013 +2124 @@
            if (igcm_h2o2.ne.0) then
              DO ig=1,ngrid
-               dqsurf(ig,igcm_h2o2)=dqsurf(ig,igcm_h2o2)
-     &                                +zdqscloud(ig,igcm_h2o2)
+              DO islope = 1,nslope
+              dqsurf(ig,igcm_h2o2,islope)=dqsurf(ig,igcm_h2o2,islope)+
+     &      zdqscloud(ig,igcm_h2o2)*cos(pi*def_slope_mean(islope)/180.)
+              ENDDO
              ENDDO
            endif
@@ -2029 +2142 @@
       if (callthermos) then
         call thermosphere(ngrid,nlayer,nq,zplev,zplay,dist_sol,
-     $     mu0,ptimestep,ptime,zday,tsurf,zzlev,zzlay,
+     $     mu0,ptimestep,ptime,zday,tsurf(:,iflat),zzlev,zzlay,
      &     pt,pq,pu,pv,pdt,pdq,
      $     zdteuv,zdtconduc,zdumolvis,zdvmolvis,zdqmoldiff,
@@ -2056 +2169 @@
          CALL geticecover(ngrid, 180.*zls/pi,
      .                  180.*longitude/pi, 180.*latitude/pi,
-     .                  qsurf(:,igcm_co2) )
-         qsurf(:,igcm_co2) = qsurf(:,igcm_co2) * 10000.
+     .                  qsurf_tmp(:,igcm_co2) )
+         qsurf_tmp(:,igcm_co2) = qsurf_tmp(:,igcm_co2) * 10000.
       ENDIF
       
@@ -2063 +2176 @@
       IF (callcond) THEN
          zdtc(:,:) = 0.
-         zdtsurfc(:) = 0.
+         zdtsurfc(:,:) = 0.
          zduc(:,:) = 0.
          zdvc(:,:) = 0.
          zdqc(:,:,:) = 0.
-         zdqsc(:,:) = 0.
+         zdqsc(:,:,:) = 0.
          CALL co2condens(ngrid,nlayer,nq,ptimestep,
      $              capcal,zplay,zplev,tsurf,pt,
      $              pphi,pdt,pdu,pdv,zdtsurf,pu,pv,pq,pdq,
-     $              qsurf(:,igcm_co2),albedo,emis,rdust,
+     $              qsurf(:,igcm_co2,:),albedo,emis,rdust,
      $              zdtc,zdtsurfc,pdpsrf,zduc,zdvc,zdqc,
      $              fluxsurf_dn_sw,zls,
@@ -2079 +2192 @@
          DO iq=1, nq
            DO ig=1,ngrid
-              dqsurf(ig,iq)=dqsurf(ig,iq)+zdqsc(ig,iq)
+              dqsurf(ig,iq,:)=dqsurf(ig,iq,:)+zdqsc(ig,iq,:)
            ENDDO  ! (ig)
          ENDDO    ! (iq)
@@ -2090 +2203 @@
          ENDDO
          DO ig=1,ngrid
-           zdtsurf(ig) = zdtsurf(ig) + zdtsurfc(ig)
+           zdtsurf(ig,:) = zdtsurf(ig,:) + zdtsurfc(ig,:)
          ENDDO
 
@@ -2130 +2243 @@
         DO iq=1, nq
           DO ig=1,ngrid
-
-            qsurf(ig,iq)=qsurf(ig,iq)+ptimestep*dqsurf(ig,iq)
-
+           DO islope = 1,nslope
+            qsurf(ig,iq,islope)=qsurf(ig,iq,islope)+
+     &                ptimestep*dqsurf(ig,iq,islope)
+           ENDDO 
           ENDDO  ! (ig)
         ENDDO    ! (iq)
@@ -2144 +2258 @@
 
       DO ig=1,ngrid
-         tsurf(ig)=tsurf(ig)+ptimestep*zdtsurf(ig) 
+        DO islope = 1,nslope
+            tsurf(ig,islope)=tsurf(ig,islope)+
+     &            ptimestep*zdtsurf(ig,islope) 
+       ENDDO
       ENDDO
 
@@ -2169 +2286 @@
 c       -------------------------------------------------------------
          do ig=1,ngrid
-           if ((qsurf(ig,igcm_co2).eq.0).and.
-     &        (qsurf(ig,igcm_h2o_ice).gt.frost_albedo_threshold)) then
+          do islope = 1,nslope
+           if ((qsurf(ig,igcm_co2,islope).eq.0).and.
+     &        (qsurf(ig,igcm_h2o_ice,islope)
+     &       .gt.frost_albedo_threshold)) then
              if ((watercaptag(ig)).and.(cst_cap_albedo)) then
-               albedo(ig,1) = albedo_h2o_cap
-               albedo(ig,2) = albedo_h2o_cap
+               albedo(ig,1,islope) = albedo_h2o_cap
+               albedo(ig,2,islope) = albedo_h2o_cap
              else
-               albedo(ig,1) = albedo_h2o_frost
-               albedo(ig,2) = albedo_h2o_frost
+               albedo(ig,1,islope) = albedo_h2o_frost
+               albedo(ig,2,islope) = albedo_h2o_frost
              endif !((watercaptag(ig)).and.(cst_cap_albedo)) then
 c              write(*,*) "frost thickness", qsurf(ig,igcm_h2o_ice)
@@ -2182 +2301 @@
 c     &        ,latitude(ig)*180./pi, longitude(ig)*180./pi
            endif
+          enddo ! islope
          enddo  ! of do ig=1,ngrid
       ENDIF  ! of IF (water)
@@ -2191 +2311 @@
 c       Thermal inertia feedback
         IF (tifeedback) THEN
-         CALL soil_tifeedback(ngrid,nsoilmx,qsurf,inertiesoil)
+         DO islope = 1,nslope
+           CALL soil_tifeedback(ngrid,nsoilmx,
+     s               qsurf(:,:,islope),inertiesoil(:,:,islope))
+         ENDDO
          CALL soil(ngrid,nsoilmx,.false.,inertiesoil,
      s          ptimestep,tsurf,tsoil,capcal,fluxgrd)
         ELSE
-         CALL soil(ngrid,nsoilmx,.false.,inertiedat,
+         CALL soil(ngrid,nsoilmx,.false.,inertiedat_tmp,
      s          ptimestep,tsurf,tsoil,capcal,fluxgrd)
         ENDIF
@@ -2242 +2365 @@
 
       DO ig=1,ngrid
-         watercap(ig)=watercap(ig)+ptimestep*dwatercap(ig)
+        DO islope = 1,nslope
+         watercap(ig,islope)=watercap(ig,islope)+
+     s   ptimestep*dwatercap(ig,islope)
+        ENDDO
       ENDDO
 
@@ -2248 +2374 @@
 
         DO ig=1,ngrid
-           if (watercaptag(ig).and.
-     &        (qsurf(ig,igcm_h2o_ice).gt.frost_metam_threshold)) then
-
-                watercap(ig)=watercap(ig)+qsurf(ig,igcm_h2o_ice)
-     &                  - frost_metam_threshold
-                qsurf(ig,igcm_h2o_ice) = frost_metam_threshold
+         DO islope = 1,nslope
+           if (watercaptag(ig).and. (qsurf(ig,igcm_h2o_ice,islope)
+     &        .gt.frost_metam_threshold)) then
+
+                watercap(ig,islope)=watercap(ig,islope)
+     &          +qsurf(ig,igcm_h2o_ice,islope)
+     &         - frost_metam_threshold
+                qsurf(ig,igcm_h2o_ice,islope) = frost_metam_threshold
            endif ! (watercaptag(ig).and.
+          ENDDO
         ENDDO
 
@@ -2263 +2392 @@
 c  13. Write output files
 c  ----------------------
+
+      call compute_meshgridavg(ngrid,nq,albedo,emis,tsurf,qsurf,
+     &       albedo_meshavg,emis_meshavg,tsurf_meshavg,qsurf_meshavg)
 
 c    -------------------------------
@@ -2312 +2444 @@
       fluxtop_up_sw_tot(1:ngrid)=fluxtop_up_sw(1:ngrid,1) +
      &                           fluxtop_up_sw(1:ngrid,2)
-      fluxsurf_dn_sw_tot(1:ngrid)=fluxsurf_dn_sw(1:ngrid,1) +
-     &                            fluxsurf_dn_sw(1:ngrid,2)
+      fluxsurf_dn_sw_tot(1:ngrid,1:nslope)=
+     &                           fluxsurf_dn_sw(1:ngrid,1,1:nslope) +
+     &                           fluxsurf_dn_sw(1:ngrid,2,1:nslope)
       fluxsurf_up_sw_tot(1:ngrid)=fluxsurf_up_sw(1:ngrid,1) +
      &                            fluxsurf_up_sw(1:ngrid,2)
@@ -2344 +2477 @@
          WRITE(*,'(a6,a10,2a15)') 'Ts','dTs','ps','dps'
          WRITE(*,'(2f10.5,2f15.5)')
-     s   tsurf(igout),zdtsurf(igout)*ptimestep,
+     s   tsurf(igout,:),zdtsurf(igout,:)*ptimestep,
      s   zplev(igout,1),pdpsrf(igout)*ptimestep
          WRITE(*,'(a4,a6,5a10)') 'l','u','du','v','dv','T','dT'
@@ -2374 +2507 @@
 
            call pbl_parameters(ngrid,nlayer,ps,zplay,z0,
-     & g,zzlay,zzlev,zu,zv,wstar,hfmax_th,zmax_th,tsurf,zh,z_out,n_out,
-     & T_out,u_out,ustar,tstar,L_mo,vhf,vvv)
+     & g,zzlay,zzlev,zu,zv,wstar,hfmax_th,zmax_th,tsurf(:,iflat),
+     & zh,z_out,n_out,T_out,u_out,ustar,tstar,L_mo,vhf,vvv)
                    ! pourquoi ustar recalcule ici? fait dans vdifc.
 
@@ -2472 +2605 @@
           call physdem1("restartfi.nc",nsoilmx,ngrid,nlayer,nq,
      .                ptimestep,ztime_fin,
-     .                tsurf,tsoil,albedo,emis,
-     .                q2,qsurf,tauscaling,totcloudfrac,wstar,
-     .                watercap)
+     .                tsurf(:,iflat),tsoil(:,:,iflat),albedo(:,:,iflat),
+     .                emis(:,iflat),
+     .                q2,qsurf(:,:,iflat),tauscaling,totcloudfrac,wstar,
+     .                watercap(:,iflat))
           
          ENDIF ! of IF (write_restart)
@@ -2677 +2811 @@
          
         call wstats(ngrid,"ps","Surface pressure","Pa",2,ps)
-        call wstats(ngrid,"tsurf","Surface temperature","K",2,tsurf)
+        call wstats(ngrid,"tsurf","Surface temperature","K",2
+     &      ,tsurf(:,iflat))
         call wstats(ngrid,"co2ice","CO2 ice cover",
-     &                "kg.m-2",2,qsurf(:,igcm_co2))
+     &                "kg.m-2",2,qsurf(:,igcm_co2,iflat))
         call wstats(ngrid,"watercap","H2O ice cover",
-     &                "kg.m-2",2,watercap)
+     &                "kg.m-2",2,watercap(:,iflat))
         call wstats(ngrid,"tau_pref_scenario",
      &              "prescribed visible dod at 610 Pa","NU",
@@ -2690 +2825 @@
         call wstats(ngrid,"fluxsurf_lw",
      &                "Thermal IR radiative flux to surface","W.m-2",2,
-     &                fluxsurf_lw)
+     &                fluxsurf_lw(:,iflat))
         call wstats(ngrid,"fluxsurf_dn_sw",
      &        "Incoming Solar radiative flux to surface","W.m-2",2,
-     &                fluxsurf_dn_sw_tot)
+     &                fluxsurf_dn_sw_tot(:,iflat))
         call wstats(ngrid,"fluxsurf_up_sw",
      &        "Reflected Solar radiative flux from surface","W.m-2",2,
@@ -2718 +2853 @@
      &                "m2.s-2",3,q2)
           call wstats(ngrid,"emis","Surface emissivity","w.m-1",2,
-     &                emis)
+     &                emis(:,iflat))
 c          call wstats(ngrid,"ssurf","Surface stress","N.m-2",
 c    &                2,zstress)
@@ -2764 +2899 @@
                call wstats(ngrid,"h2o_ice_s",
      &                    "surface h2o_ice","kg/m2",
-     &                    2,qsurf(1,igcm_h2o_ice))
+     &                    2,qsurf(1,igcm_h2o_ice,iflat))
                call wstats(ngrid,'albedo',
      &                         'albedo',
-     &                         '',2,albedo(1,1))
+     &                         '',2,albedo(1,1,iflat))
                call wstats(ngrid,"mtot",
      &                    "total mass of water vapor","kg/m2",
@@ -2852 +2987 @@
      &                         'kg/kg',3,zq(1,1,iq))
                  call wstats(ngrid,'qsurf'//str2,'qsurf',
-     &                         'kg.m-2',2,qsurf(1,iq))
+     &                         'kg.m-2',2,qsurf(1,iq,iflat))
                end do
              endif ! (doubleq)
@@ -3019 +3154 @@
 c        for any variables !!
          call WRITEDIAGFI(ngrid,"emis","Surface emissivity","w.m-1",2,
-     &                  emis)
+     &                  emis(:,iflat))
 c        call WRITEDIAGFI(ngrid,"pplay","Pressure","Pa",3,zplay)
 c        call WRITEDIAGFI(ngrid,"pplev","Pressure","Pa",3,zplev)
@@ -3029 +3164 @@
      &                    "m2s-2",2,phisfi)
          call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",2,
-     &                  tsurf)
+     &                  tsurf(:,iflat))
          call WRITEDIAGFI(ngrid,"ps","surface pressure","Pa",2,ps)
          call WRITEDIAGFI(ngrid,"co2ice","co2 ice thickness"
-     &                              ,"kg.m-2",2,qsurf(:,igcm_co2))
+     &                              ,"kg.m-2",2,qsurf(:,igcm_co2,iflat))
          call WRITEDIAGFI(ngrid,"watercap","Water ice thickness"
-     &                                         ,"kg.m-2",2,watercap)
+     &         ,"kg.m-2",2,watercap(:,iflat))
 
          call WRITEDIAGFI(ngrid,"temp_layer1","temperature in layer 1",
@@ -3041 +3176 @@
      &                  "K",2,zt(1,7))
          call WRITEDIAGFI(ngrid,"fluxsurf_lw","fluxsurf_lw","W.m-2",2,
-     &                  fluxsurf_lw)
+     &                  fluxsurf_lw(:,iflat))
          call WRITEDIAGFI(ngrid,"fluxsurf_dn_sw","fluxsurf_dn_sw",
-     &                  "W.m-2",2,fluxsurf_dn_sw_tot)
+     &                  "W.m-2",2,fluxsurf_dn_sw_tot(:,iflat))
          call WRITEDIAGFI(ngrid,"fluxtop_lw","fluxtop_lw","W.m-2",2,
      &                  fluxtop_lw)
@@ -3079 +3214 @@
       !!! this is to ensure correct initialisation of mesoscale model
       call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",2,
-     &                 tsurf)
+     &                 tsurf(:,iflat))
       call WRITEDIAGFI(ngrid,"ps","surface pressure","Pa",2,ps)
       call WRITEDIAGFI(ngrid,"co2ice","co2 ice thickness","kg.m-2",2,
-     &                 qsurf(:,igcm_co2))
+     &                 qsurf(:,igcm_co2,iflat))
       call WRITEDIAGFI(ngrid,"temp","temperature","K",3,zt)
       call WRITEDIAGFI(ngrid,"u","Zonal wind","m.s-1",3,zu)
@@ -3089 +3224 @@
      &                 emis)
       call WRITEDIAGFI(ngrid,"tsoil","Soil temperature",
-     &                 "K",3,tsoil)
+     &                 "K",3,tsoil(:,:,iflat))
       call WRITEDIAGFI(ngrid,"inertiedat","Soil inertia",
      &                 "K",3,inertiedat)
@@ -3170 +3305 @@
           call WRITEDIAGFI(ngrid,'qsurf02','surface tracer',
      &                     'kg.m-2',2,
-     &                     qsurf(1:ngrid,igcm_h2o_ice))
+     &                     qsurf(1:ngrid,igcm_h2o_ice,iflat))
 #endif
           call WRITEDIAGFI(ngrid,'mtot',
@@ -3265 +3400 @@
             call WRITEDIAGFI(ngrid,'h2o_ice_s',
      &                       'surface h2o_ice',
-     &                       'kg.m-2',2,qsurf(1,igcm_h2o_ice))
+     &                       'kg.m-2',2,qsurf(1,igcm_h2o_ice,iflat))
             if (hdo) then
             call WRITEDIAGFI(ngrid,'hdo_ice_s',
      &                       'surface hdo_ice',
-     &                       'kg.m-2',2,qsurf(1,igcm_hdo_ice))
+     &                       'kg.m-2',2,qsurf(1,igcm_hdo_ice,iflat))
 
                 do ig=1,ngrid
-                if (qsurf(ig,igcm_h2o_ice).gt.qperemin) then
-                    DoH_surf(ig) = 0.5*( qsurf(ig,igcm_hdo_ice)/
-     &                  qsurf(ig,igcm_h2o_ice) )/155.76e-6
+                if (qsurf_meshavg(ig,igcm_h2o_ice).gt.qperemin) then
+           DoH_surf(ig) = 0.5*( qsurf_meshavg(ig,igcm_hdo_ice)/
+     &          qsurf_meshavg(ig,igcm_h2o_ice) )/155.76e-6
                 else
                     DoH_surf(ig) = 0.
@@ -3287 +3422 @@
             CALL WRITEDIAGFI(ngrid,'albedo',
      &                         'albedo',
-     &                         '',2,albedo(1,1))
+     &                         '',2,albedo(1,1,iflat))
               if (tifeedback) then
                  call WRITEDIAGSOIL(ngrid,"soiltemp",
@@ -3294 +3429 @@
                  call WRITEDIAGSOIL(ngrid,'soilti',
      &                       'Soil Thermal Inertia',
-     &                       'J.s-1/2.m-2.K-1',3,inertiesoil)
+     &                       'J.s-1/2.m-2.K-1',3,inertiesoil(:,:,iflat))
               endif
 !A. Pottier
@@ -3392 +3527 @@
      &                         'kg/kg',3,zq(1,1,iq))
                call WRITEDIAGFI(ngrid,'qsurf'//str2,'qsurf',
-     &                         'kg.m-2',2,qsurf(1,iq))
+     &                         'kg.m-2',2,qsurf(1,iq,iflat))
              end do
            endif ! (doubleq)
@@ -3419 +3554 @@
              call WRITEDIAGFI(ngrid,'qsurf',
      &                 'stormdust injection',
-     &                 'kg.m-2',2,qsurf(:,igcm_stormdust_mass))
+     &                 'kg.m-2',2,qsurf(:,igcm_stormdust_mass,iflat))
              call WRITEDIAGFI(ngrid,'pdqsurf',
      &                  'tendancy stormdust mass at surface',
-     &                  'kg.m-2',2,dqsurf(:,igcm_stormdust_mass))
+     &                  'kg.m-2',2,dqsurf(:,igcm_stormdust_mass,iflat))
              call WRITEDIAGFI(ngrid,'wspeed','vertical speed stormdust',
      &                        'm/s',3,wspeed(:,1:nlayer))
@@ -3484 +3619 @@
      &                        'part/kg',3,nccn)
              call WRITEDIAGFI(ngrid,'surfccnq','Surf nuclei mass mr',
-     &                        'kg.m-2',2,qsurf(1,igcm_ccn_mass))
+     &                        'kg.m-2',2,qsurf(1,igcm_ccn_mass,iflat))
              call WRITEDIAGFI(ngrid,'surfccnN','Surf nuclei number',
-     &                        'kg.m-2',2,qsurf(1,igcm_ccn_number))
+     &                        'kg.m-2',2,qsurf(1,igcm_ccn_number,iflat))
            endif ! (scavenging)
 
@@ -3636 +3771 @@
          ! Write soil temperature
         call writediagsoil(ngrid,"soiltemp","Soil temperature","K",
-     &                     3,tsoil)
+     &                     3,tsoil(:,:,iflat))
          ! Write surface temperature
 !        call writediagsoil(ngrid,"tsurf","Surface temperature","K",
@@ -3666 +3801 @@
 c        CALL writeg1d(ngrid,nlayer,pw,'w','m.s-1')
         call writediagsoil(ngrid,"soiltemp","Soil temperature","K",
-     &                     3,tsoil)
+     &                     3,tsoil(:,:,iflat))
 
 ! THERMALS STUFF 1D
@@ -3690 +3825 @@
          call WRITEDIAGFI(ngrid,"q2","q2","kg.m-3",1,q2)
          call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",0,
-     &                  tsurf)
+     &                  tsurf(:,iflat))
          call WRITEDIAGFI(ngrid,"u","u wind","m/s",1,zu)
          call WRITEDIAGFI(ngrid,"v","v wind","m/s",1,zv)
@@ -3704 +3839 @@
      &                   'w.m-2',1,zdtlw)
          call WRITEDIAGFI(ngrid,"co2ice","co2 ice thickness"
-     &                                   ,"kg.m-2",0,qsurf(:,igcm_co2))
+     &       ,"kg.m-2",0,qsurf(:,igcm_co2,iflat))
 
          if (igcm_co2.ne.0) then
@@ -3742 +3877 @@
      &                        ,1,aerosol(:,:,iaer_stormdust_doubleq))
              call WRITEDIAGFI(ngrid,'dqsdifdustq','diffusion',
-     &                       'kg.m-2.s-1',0,zdqsdif(1,igcm_dust_mass))
+     &          'kg.m-2.s-1',0,zdqsdif(1,igcm_dust_mass,iflat))
              call WRITEDIAGFI(ngrid,'dqsdifrdsq','diffusion',
-     &                  'kg.m-2.s-1',0,zdqsdif(1,igcm_stormdust_mass))
+     &          'kg.m-2.s-1',0,zdqsdif(1,igcm_stormdust_mass,iflat))
              call WRITEDIAGFI(ngrid,'mstormdtot',
      &                        'total mass of stormdust only',
@@ -3767 +3902 @@
              call WRITEDIAGFI(ngrid,'rdsqsurf',
      &                 'stormdust at surface',
-     &                 'kg.m-2',0,qsurf(:,igcm_stormdust_mass))
+     &                 'kg.m-2',0,qsurf(:,igcm_stormdust_mass,iflat))
              call WRITEDIAGFI(ngrid,'qsurf',
      &                  'dust mass at surface',
-     &                  'kg.m-2',0,qsurf(:,igcm_dust_mass))
+     &                  'kg.m-2',0,qsurf(:,igcm_dust_mass,iflat))
              call WRITEDIAGFI(ngrid,'wspeed','vertical speed stormdust',
      &                        'm/s',1,wspeed)
@@ -3824 +3959 @@
            mtot = 0
            icetot = 0
-           h2otot = qsurf(1,igcm_h2o_ice)
+           h2otot = qsurf(1,igcm_h2o_ice,iflat)
            if (hdo) THEN
            mtotD = 0
            icetotD = 0
-           hdotot = qsurf(1,igcm_hdo_ice)
+           hdotot = qsurf(1,igcm_hdo_ice,iflat)
            ENDIF !hdo
 
@@ -3860 +3995 @@
              CALL WRITEDIAGFI(ngrid,'albedo',
      &                         'albedo',
-     &                         '',2,albedo(1,1))
+     &                         '',2,albedo(1,1,iflat))
 
              IF (hdo) THEN
@@ -3970 +4105 @@
            do iq=1,nq
              call WRITEDIAGFI(ngrid,trim(noms(iq))//'_s',
-     &            trim(noms(iq))//'_s','kg/kg',0,qsurf(1,iq))
+     &       trim(noms(iq))//'_s','kg/kg',0,qsurf(1,iq,iflat))
            end do
      
@@ -4030 +4165 @@
      &        +(pdq(ig,l,igcm_co2)+pdq(ig,l,igcm_co2_ice))*ptimestep)
           enddo
-          co2totB = co2totB + qsurf(ig,igcm_co2)
+          co2totB = co2totB + qsurf(ig,igcm_co2,iflat)
         enddo
       else
@@ -4039 +4174 @@
      &             (pq(ig,l,igcm_co2)+pdq(ig,l,igcm_co2)*ptimestep)
           enddo
-          co2totB = co2totB + qsurf(ig,igcm_co2)
+          co2totB = co2totB + qsurf(ig,igcm_co2,iflat)
         enddo
       endif ! of if (igcm_co2_ice.ne.0)

trunk/LMDZ.MARS/libf/phymars/soil.F

@@ -8 +8 @@
      &                     coefd, alph, beta, mu,flux_geo
       use surfdat_h, only: watercaptag, inert_h2o_ice
-
+      use comslope_mod, ONLY: nslope
       implicit none
 
@@ -29 +29 @@
       integer nsoil     ! number of soil layers 
       logical firstcall ! identifier for initialization call 
-      real therm_i(ngrid,nsoil) ! thermal inertia
+      real therm_i(ngrid,nsoil,nslope) ! thermal inertia
       real timestep         ! time step
-      real tsurf(ngrid)   ! surface temperature
+      real tsurf(ngrid,nslope)   ! surface temperature
 ! outputs:
-      real tsoil(ngrid,nsoil) ! soil (mid-layer) temperature
-      real capcal(ngrid) ! surface specific heat
-      real fluxgrd(ngrid) ! surface diffusive heat flux
+      real tsoil(ngrid,nsoil,nslope) ! soil (mid-layer) temperature
+      real capcal(ngrid,nslope) ! surface specific heat
+      real fluxgrd(ngrid,nslope) ! surface diffusive heat flux
 
 ! local variables:
-      integer ig,ik
+      integer ig,ik,islope
 
 ! 0. Initialisations and preprocessing step
@@ -46 +46 @@
 ! 0.1 Build mthermdiff(:), the mid-layer thermal diffusivities
       do ig=1,ngrid
+       do islope = 1,nslope
         if (watercaptag(ig)) then
           do ik=0,nsoil-1
 ! If we have permanent ice, we use the water ice thermal inertia from ground to last layer.
-              mthermdiff(ig,ik)=inert_h2o_ice*inert_h2o_ice/volcapa
+              mthermdiff(ig,ik,islope)=inert_h2o_ice*
+     &            inert_h2o_ice/volcapa
           enddo
         else
           do ik=0,nsoil-1
-           mthermdiff(ig,ik)=therm_i(ig,ik+1)*therm_i(ig,ik+1)/volcapa
+           mthermdiff(ig,ik,islope)=therm_i(ig,ik+1,islope)*
+     &         therm_i(ig,ik+1,islope)/volcapa
           enddo
         endif
       enddo
+      enddo
 
 #ifdef MESOSCALE
       do ig=1,ngrid
-        if ( therm_i(ig,1) .ge. inert_h2o_ice ) then
-          print *, "limit max TI ", therm_i(ig,1), inert_h2o_ice
+       do islope = 1,nslope
+        if ( therm_i(ig,1,islope) .ge. inert_h2o_ice ) then
+          print *, "limit max TI ", therm_i(ig,1,islope), inert_h2o_ice
           do ik=0,nsoil-1
-               mthermdiff(ig,ik)=inert_h2o_ice*inert_h2o_ice/volcapa
+               mthermdiff(ig,ik,islope)=inert_h2o_ice*
+     &    inert_h2o_ice/volcapa
           enddo
         endif
+       enddo
       enddo
 #endif
@@ -71 +78 @@
 ! 0.2 Build thermdiff(:), the "interlayer" thermal diffusivities
       do ig=1,ngrid
+       do islope = 1,nslope
         do ik=1,nsoil-1
-      thermdiff(ig,ik)=((layer(ik)-mlayer(ik-1))*mthermdiff(ig,ik)
-     &                +(mlayer(ik)-layer(ik))*mthermdiff(ig,ik-1))
+      thermdiff(ig,ik,islope)=((layer(ik)-mlayer(ik-1))
+     &                        *mthermdiff(ig,ik,islope)
+     &                +(mlayer(ik)-layer(ik))
+     &                *mthermdiff(ig,ik-1,islope))
      &                    /(mlayer(ik)-mlayer(ik-1))
 !       write(*,*),'soil: ik: ',ik,' thermdiff:',thermdiff(ig,ik)
         enddo
+       enddo
       enddo
 
@@ -92 +103 @@
 
       do ig=1,ngrid
+       do islope = 1,nslope
         ! d_k
         do ik=1,nsoil-1
-          coefd(ig,ik)=thermdiff(ig,ik)/(mlayer(ik)-mlayer(ik-1))
+          coefd(ig,ik,islope)=thermdiff(ig,ik,islope)
+     &        /(mlayer(ik)-mlayer(ik-1))
         enddo
         
         ! alph_{N-1}
-        alph(ig,nsoil-1)=coefd(ig,nsoil-1)/
-     &                  (coefq(nsoil-1)+coefd(ig,nsoil-1))
+        alph(ig,nsoil-1,islope)=coefd(ig,nsoil-1,islope)/
+     &                  (coefq(nsoil-1)+coefd(ig,nsoil-1,islope))
         ! alph_k
         do ik=nsoil-2,1,-1
-          alph(ig,ik)=coefd(ig,ik)/(coefq(ik)+coefd(ig,ik+1)*
-     &                              (1.-alph(ig,ik+1))+coefd(ig,ik))
+          alph(ig,ik,islope)=coefd(ig,ik,islope)/
+     &         (coefq(ik)+coefd(ig,ik+1,islope)*
+     &          (1.-alph(ig,ik+1,islope))+coefd(ig,ik,islope))
         enddo
 
         ! capcal
 ! Cstar
-        capcal(ig)=volcapa*layer(1)+
-     &              (thermdiff(ig,1)/(mlayer(1)-mlayer(0)))*
-     &              (timestep*(1.-alph(ig,1)))
+        capcal(ig,islope)=volcapa*layer(1)+
+     &         (thermdiff(ig,1,islope)/(mlayer(1)-mlayer(0)))*
+     &         (timestep*(1.-alph(ig,1,islope)))
 ! Cs
-        capcal(ig)=capcal(ig)/(1.+mu*(1.0-alph(ig,1))*
-     &                         thermdiff(ig,1)/mthermdiff(ig,0))
+        capcal(ig,islope)=capcal(ig,islope)/
+     &            (1.+mu*(1.0-alph(ig,1,islope))*
+     &            thermdiff(ig,1,islope)/mthermdiff(ig,0,islope))
 !      write(*,*)'soil: ig=',ig,' capcal(ig)=',capcal(ig)
+      enddo ! islope
       enddo ! of do ig=1,ngrid
             
@@ -122 +138 @@
       IF (.not.firstcall) THEN
 ! First layer:
-      do ig=1,ngrid
-        tsoil(ig,1)=(tsurf(ig)+mu*beta(ig,1)*
-     &                         thermdiff(ig,1)/mthermdiff(ig,0))/
-     &              (1.+mu*(1.0-alph(ig,1))*
-     &               thermdiff(ig,1)/mthermdiff(ig,0))
+      do islope = 1,nslope
+      do ig=1,ngrid
+        tsoil(ig,1,islope)=(tsurf(ig,islope)+mu*beta(ig,1,islope)*
+     &      thermdiff(ig,1,islope)/mthermdiff(ig,0,islope))/
+     &      (1.+mu*(1.0-alph(ig,1,islope))*
+     &      thermdiff(ig,1,islope)/mthermdiff(ig,0,islope))
       enddo
 ! Other layers:
       do ik=1,nsoil-1
         do ig=1,ngrid
-          tsoil(ig,ik+1)=alph(ig,ik)*tsoil(ig,ik)+beta(ig,ik)
-        enddo
-      enddo
-      
+          tsoil(ig,ik+1,islope)=alph(ig,ik,islope)*
+     &    tsoil(ig,ik,islope)+beta(ig,ik,islope)
+        enddo
+      enddo
+       enddo ! islope
       ENDIF! of if (.not.firstcall)
 
 !  2. Compute beta coefficients (preprocessing for next time step)
 ! Bottom layer, beta_{N-1}
-      do ig=1,ngrid
-        beta(ig,nsoil-1)=coefq(nsoil-1)*tsoil(ig,nsoil)
-     &            /(coefq(nsoil-1)+coefd(ig,nsoil-1))
-     &            +flux_geo(ig)/(coefq(nsoil-1) + coefd(ig,nsoil-1))
+       do islope = 1,nslope
+      do ig=1,ngrid
+        beta(ig,nsoil-1,islope)=coefq(nsoil-1)*tsoil(ig,nsoil,islope)
+     &      /(coefq(nsoil-1)+coefd(ig,nsoil-1,islope))
+     &      +flux_geo(ig)/(coefq(nsoil-1) + coefd(ig,nsoil-1,islope))
       enddo
    
@@ -148 +167 @@
       do ik=nsoil-2,1,-1
         do ig=1,ngrid
-          beta(ig,ik)=(coefq(ik)*tsoil(ig,ik+1)+
-     &                 coefd(ig,ik+1)*beta(ig,ik+1))/
-     &                 (coefq(ik)+coefd(ig,ik+1)*(1.0-alph(ig,ik+1))
-     &                  +coefd(ig,ik))
+          beta(ig,ik,islope)=(coefq(ik)*tsoil(ig,ik+1,islope)+
+     &                 coefd(ig,ik+1,islope)*beta(ig,ik+1,islope))/
+     &                 (coefq(ik)+coefd(ig,ik+1,islope)*
+     &                 (1.0-alph(ig,ik+1,islope))
+     &                  +coefd(ig,ik,islope))
         enddo
       enddo
@@ -162 +182 @@
 !     &              (timestep*(1.-alph(ig,1)))
 ! Fstar
-        fluxgrd(ig)=(thermdiff(ig,1)/(mlayer(1)-mlayer(0)))*
-     &              (beta(ig,1)+(alph(ig,1)-1.0)*tsoil(ig,1))
+        fluxgrd(ig,islope)=(thermdiff(ig,1,islope)/
+     &              (mlayer(1)-mlayer(0)))*
+     &              (beta(ig,1,islope)+(alph(ig,1,islope)-1.0)
+     &              *tsoil(ig,1,islope))
 
 !        mu=mlayer(ig,0)/(mlayer(ig,1)-mlayer(ig,0))
@@ -169 +191 @@
 !     &                         thermdiff(ig,1)/mthermdiff(ig,0))
 ! Fs
-        fluxgrd(ig)=fluxgrd(ig)+(capcal(ig)/timestep)*
-     &              (tsoil(ig,1)*(1.+mu*(1.0-alph(ig,1))*
-     &                         thermdiff(ig,1)/mthermdiff(ig,0))
-     &               -tsurf(ig)-mu*beta(ig,1)*
-     &                          thermdiff(ig,1)/mthermdiff(ig,0))
-      enddo
-
+        fluxgrd(ig,islope)=fluxgrd(ig,islope)+(capcal(ig,islope)
+     &              /timestep)*
+     &              (tsoil(ig,1,islope)*
+     &              (1.+mu*(1.0-alph(ig,1,islope))*
+     &               thermdiff(ig,1,islope)/mthermdiff(ig,0,islope))
+     &               -tsurf(ig,islope)-mu*beta(ig,1,islope)*
+     &                thermdiff(ig,1,islope)/mthermdiff(ig,0,islope))
+      enddo
+      enddo ! islope
       end
 

trunk/LMDZ.MARS/libf/phymars/surfdat_h.F90

@@ -44 +44 @@
 
   !! variables
-  REAL,SAVE,ALLOCATABLE :: tsurf(:)   ! Surface temperature (K)
-  REAL,SAVE,ALLOCATABLE :: emis(:)    ! Thermal IR surface emissivity
-  REAL,SAVE,ALLOCATABLE :: capcal(:) ! surface heat capacity (J m-2 K-1)
-  REAL,SAVE,ALLOCATABLE :: fluxgrd(:) ! surface conduction flux (W.m-2)
-  REAL,ALLOCATABLE,SAVE :: qsurf(:,:) ! tracer on surface (e.g. kg.m-2)
-  REAL,SAVE,ALLOCATABLE :: watercap(:) ! Surface water ice (kg.m-2)
+  REAL,SAVE,ALLOCATABLE :: tsurf(:,:)   ! Surface temperature (K)
+  REAL,SAVE,ALLOCATABLE :: emis(:,:)    ! Thermal IR surface emissivity
+  REAL,SAVE,ALLOCATABLE :: capcal(:,:) ! surface heat capacity (J m-2 K-1)
+  REAL,SAVE,ALLOCATABLE :: fluxgrd(:,:) ! surface conduction flux (W.m-2)
+  REAL,ALLOCATABLE,SAVE :: qsurf(:,:,:) ! tracer on surface (e.g. kg.m-2)
+  REAL,SAVE,ALLOCATABLE :: watercap(:,:) ! Surface water ice (kg.m-2)
 
 !$OMP THREADPRIVATE(tsurf,emis,capcal,fluxgrd,qsurf,watercap)
@@ -55 +55 @@
 contains
 
-  subroutine ini_surfdat_h(ngrid,nq)
+  subroutine ini_surfdat_h(ngrid,nq,nslope)
   
   implicit none
   integer,intent(in) :: ngrid ! number of atmospheric columns
   integer,intent(in) :: nq ! number of tracers  
-
+  integer,intent(in) :: nslope ! number of sub-grid scale slope  
     allocate(albedodat(ngrid))
     allocate(phisfi(ngrid))
@@ -71 +71 @@
     allocate(zthe(ngrid))
     allocate(z0(ngrid))
-    allocate(qsurf(ngrid,nq))
-    allocate(tsurf(ngrid))
-    allocate(watercap(ngrid))
-    allocate(emis(ngrid))
-    allocate(capcal(ngrid))
-    allocate(fluxgrd(ngrid))
+    allocate(qsurf(ngrid,nq,nslope))
+    allocate(tsurf(ngrid,nslope))
+    allocate(watercap(ngrid,nslope))
+    allocate(emis(ngrid,nslope))
+    allocate(capcal(ngrid,nslope))
+    allocate(fluxgrd(ngrid,nslope))
     allocate(hmons(ngrid))
     allocate(summit(ngrid))

trunk/LMDZ.MARS/libf/phymars/writediagfi.F

@@ -78 +78 @@
 
       integer,save :: zitau=0
-      character(len=20),save :: firstnom='1234567890'
+      character(len=27),save :: firstnom='1234567890'
 !$OMP THREADPRIVATE(zitau,firstnom)