Index: trunk/LMDZ.MARS/libf/phymars/getslopes.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/getslopes.F90	(revision 900)
+++ trunk/LMDZ.MARS/libf/phymars/getslopes.F90	(revision 900)
@@ -0,0 +1,97 @@
+subroutine getslopes(geopot)
+    
+implicit none
+
+#include "dimensions.h"
+#include "dimphys.h"
+#include "slope.h"
+#include "comgeomfi.h"
+#include "comcstfi.h"
+
+
+! This routine computes slope inclination and orientation for the GCM (callslope=.true. in callphys.def)
+! It works fine with a non-regular grid for zoomed simulations.
+! slope inclination angle (deg)  0 == horizontal, 90 == vertical
+! slope orientation angle (deg)  0 == Northward,  90 == Eastward, 180 == Southward, 270 == Westward
+! TN 04/1013
+
+
+real geopot(ngridmx)     ! geopotential on phy grid
+real topogrid(iim,jjm+1) ! topography on lat/lon grid with poles and only one -180/180 point
+real latigrid(iim,jjm+1),longgrid(iim,jjm+1) ! meshgrid of latitude and longitude values (radians)
+real theta_val ! slope inclination
+real psi_val   ! slope orientation
+real gradx(iim,jjm+1) ! x: latitude-wise topography gradient,  increasing northward
+real grady(iim,jjm+1) ! y: longitude-wise topography gradient, increasing westward
+integer i,j,ig0
+
+
+! rearrange topography on a 2d array
+do j=2,jjm
+   ig0= 1+(j-2)*iim
+   do i=1,iim
+      topogrid(i,j)=geopot(ig0+i)/g
+      latigrid(i,j)=lati(ig0+i)
+      longgrid(i,j)=long(ig0+i)
+   enddo
+enddo
+!poles :
+topogrid(:,1) = geopot(1)/g
+latigrid(:,1) = lati(1)
+longgrid(:,1) = long(1)
+topogrid(:,jjm+1) = geopot(ngridmx)/g
+latigrid(:,jjm+1) = lati(ngridmx)
+longgrid(:,jjm+1) = long(ngridmx)
+
+
+
+! compute topography gradient
+! topogrid and rad are both in meters
+do j=2,jjm
+   do i=1,iim
+     gradx(i,j) = (topogrid(i,j+1) - topogrid(i,j-1)) / (latigrid(i,j+1)-latigrid(i,j-1))
+     gradx(i,j) = gradx(i,j) / rad
+   enddo
+   grady(1,j) = (topogrid(2,j) - topogrid(iim,j)) / (2*pi+longgrid(2,j)-longgrid(iim,j)) 
+   grady(1,j) = grady(1,j) / rad
+   grady(iim,j) = (topogrid(1,j) - topogrid(iim-1,j)) / (2*pi+longgrid(1,j)-longgrid(iim-1,j)) 
+   grady(iim,j) = grady(iim,j) / rad
+   do i=2,iim-1
+     grady(i,j) = (topogrid(i+1,j) - topogrid(i-1,j)) / (longgrid(i+1,j)-longgrid(i-1,j)) 
+     grady(i,j) = grady(i,j) / rad
+   enddo
+enddo
+! poles :
+gradx(:,1) = 0.
+grady(:,1) = 0.
+gradx(:,jjm+1) = 0.
+grady(:,jjm+1) = 0.
+
+
+
+! compute slope inclination and orientation :
+theta_sl(:) = 0.
+psi_sl(:)   = 0.
+do j=2,jjm
+   do i=1,iim
+   
+     ig0= 1+(j-2)*iim
+   
+     theta_val=atan(sqrt( (gradx(i,j))**2 + (grady(i,j))**2 ))
+     
+     psi_val=0.
+     if (gradx(i,j) .ne. 0.) psi_val= -pi/2. - atan(grady(i,j)/gradx(i,j))
+     if (gradx(i,j) .ge. 0.) psi_val= psi_val - pi
+     psi_val = 3*pi/2. - psi_val
+     psi_val = psi_val*180./pi
+     psi_val = MODULO(psi_val,360.)
+     
+     theta_sl(ig0+i) = theta_val
+     psi_sl(ig0+i)   = psi_val
+          
+   enddo
+enddo
+
+
+
+end subroutine getslopes
Index: trunk/LMDZ.MARS/libf/phymars/physiq.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/physiq.F	(revision 899)
+++ trunk/LMDZ.MARS/libf/phymars/physiq.F	(revision 900)
@@ -475,4 +475,8 @@
      .                  albedo_h2o_ice
         ENDIF
+
+#ifndef MESOSCALE
+         if (callslope) call getslopes(phisfi)
+#endif
                    
       ENDIF        !  (end of "if firstcall")