Index: /trunk/LMDZ.MARS/README
===================================================================
--- /trunk/LMDZ.MARS/README	(revision 223)
+++ /trunk/LMDZ.MARS/README	(revision 224)
@@ -778,2 +778,24 @@
 - Tidying up dust properties in DATAFILE for better consistency (cf. suaer.F90)
 - Cosmetic changes in aeropacity.F (changed comment and put a print inside a water flag)
+
+== 15/07/2011 == EM
+>> Implemented using 'z0' roughness length map (important: 'z0' reference
+   field is in datafile surface.nc, which has also been updated).
+- made z0 a z0(ngridmx) array and moved 'z0' from 'planete.h' to 'surfdat.h';
+  added a 'z0_default' (common in surfdat.h) corresponding to the 'control'
+  array value (contole(19) in startfi.nc).
+- adapted 'tabfi.F' to use 'z0_default'.
+- adapted 'phyetat0.F' to look for a 'z0' field in startfi.nc. If 'z0' is not
+  found in the startfi.nc file, then the uniform default value (z0_default)
+  is used.
+- modified 'physdem1.F' to write 'z0' field to restart.nc
+- adapted use of z0() in 'physiq.F' (diagnostic computation of surface stress),
+  'vdifc.F' and 'vdif_cd.F'.
+- adapted 'dustdevil.F' to use 'z0_default'.
+- 'testphys1d.F' now uses 'z0_default', and the value to use can be set
+  in run.def (with "z0=TheValueYouWant").
+- modified 'datareadnc.F' to load reference map of 'z0' from surface.nc, and
+  added a 'z0' option in 'newstart.F' to force a uniform value of z0. Note
+  that the use of the z0 map is automatic when using newstart, but only when
+  it loads a start_archive.nc file.
+
Index: /trunk/LMDZ.MARS/libf/dyn3d/newstart.F
===================================================================
--- /trunk/LMDZ.MARS/libf/dyn3d/newstart.F	(revision 223)
+++ /trunk/LMDZ.MARS/libf/dyn3d/newstart.F	(revision 224)
@@ -90,4 +90,5 @@
       real zmeaS(iip1,jjp1),zstdS(iip1,jjp1)
       real zsigS(iip1,jjp1),zgamS(iip1,jjp1),ztheS(iip1,jjp1)
+      real z0S(iip1,jjp1)
 
 c variable physique
@@ -365,6 +366,6 @@
       call getin("datadir",datafile) ! in case user specified another path
 
-      CALL datareadnc(relief,phis,alb,surfith,zmeaS,zstdS,zsigS,zgamS,
-     .          ztheS)
+      CALL datareadnc(relief,phis,alb,surfith,z0S,
+     &          zmeaS,zstdS,zsigS,zgamS,ztheS)
 
       CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,phisfi)
@@ -372,4 +373,5 @@
       CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,surfith,surfithfi)
       CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
+      CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,z0S,z0)
       CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zmeaS,zmea)
       CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,zstdS,zstd)
@@ -470,4 +472,5 @@
       write(*,*) 'flat : no topography ("aquaplanet")'
       write(*,*) 'bilball : uniform albedo and thermal inertia'
+      write(*,*) 'z0 : set a uniform surface roughness length'
       write(*,*) 'coldspole : cold subsurface and high albedo at S.pole'
       write(*,*) 'qname : change tracer name'
@@ -567,4 +570,19 @@
           CALL gr_dyn_fi(nsoilmx,iip1,jjp1,ngridmx,ith,ithfi)
           CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,alb,albfi)
+        
+         ! also reset surface roughness length to default value
+         write(*,*) 'surface roughness length set to:',z0_default,' m'
+         z0(:)=z0_default
+
+!       z0 : set surface roughness length to a constant value
+!       -----------------------------------------------------
+        else if (trim(modif) .eq. 'z0') then
+          write(*,*) 'set a uniform surface roughness length'
+          write(*,*) ' value for z0_default (ex: ',z0_default,')?'
+          ierr=1
+          do while (ierr.ne.0)
+            read(*,*,iostat=ierr) z0_default
+          enddo
+          z0(:)=z0_default
 
 c       coldspole : sous-sol de la calotte sud toujours froid
Index: /trunk/LMDZ.MARS/libf/phymars/datareadnc.F
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/datareadnc.F	(revision 223)
+++ /trunk/LMDZ.MARS/libf/phymars/datareadnc.F	(revision 224)
@@ -1,5 +1,5 @@
 c=======================================================================
-      SUBROUTINE datareadnc(relief,phisinit,alb,ith,
-     .                    zmea,zstd,zsig,zgam,zthe)
+      SUBROUTINE datareadnc(relief,phisinit,alb,ith,z0,
+     &                    zmea,zstd,zsig,zgam,zthe)
 c=======================================================================
 c
@@ -64,20 +64,21 @@
       parameter    (iimp1=iim+1-1/iim)
 
-      CHARACTER relief*3
-
+! Arguments:
+      CHARACTER(len=3),intent(inout) :: relief
+      REAL,intent(out) :: phisinit(iimp1*jjp1)
+      REAL,intent(out) :: alb(iimp1*jjp1)
+      REAL,intent(out) :: ith(iimp1*jjp1)
+      REAL,intent(out) :: z0(iimp1*jjp1)
+      REAL,intent(out) :: zmea(imdp1*jmdp1)
+      REAL,intent(out) :: zstd(imdp1*jmdp1)
+      REAL,intent(out) :: zsig(imdp1*jmdp1)
+      REAL,intent(out) :: zgam(imdp1*jmdp1)
+      REAL,intent(out) :: zthe(imdp1*jmdp1)
+
+! Local variables:
       REAL        zdata(imd*jmdp1)
       REAL        zdataS(imdp1*jmdp1)
       REAL        pfield(iimp1*jjp1)
 
-      REAL        alb(iimp1*jjp1)
-      REAL        ith(iimp1*jjp1)
-      REAL        phisinit(iimp1*jjp1)
-
-      REAL        zmea(imdp1*jmdp1)
-      REAL        zstd(imdp1*jmdp1)
-      REAL        zsig(imdp1*jmdp1)
-      REAL        zgam(imdp1*jmdp1)
-      REAL        zthe(imdp1*jmdp1)
-
       INTEGER     ierr
 
@@ -95,5 +96,5 @@
 
       CHARACTER*20 string
-      DIMENSION string(4)
+      DIMENSION string(0:4)
 
 
@@ -123,6 +124,4 @@
         write(*,*)'1) You can set this path in the callphys.def file:'
         write(*,*)'   datadir=/path/to/the/datafiles'
-        write(*,*)'   OR specify the path to use in callphys.def, as:'
-        write(*,*)'   datadir=/path/to/the/directory'
         write(*,*)'2) If necessary, surface.nc (and other datafiles)'
         write(*,*)'   can be obtained online on:'
@@ -179,4 +178,5 @@
 c=======================================================================
 
+      string(0) = 'z0'
       string(1) = 'albedo'
       string(2) = 'thermal'
@@ -192,5 +192,5 @@
  
 
-      DO k=1,4
+      DO k=0,4
           write(*,*) 'string',k,string(k)
           
@@ -207,4 +207,9 @@
 
       ierr = NF_INQ_VARID (unit, string(k), nvarid)
+      if (ierr.ne.nf_noerr) then
+        write(*,*) 'datareadnc error, cannot find ',trim(string(k))
+        write(*,*) ' in file ',trim(datafile),'/surface.nc'
+        stop
+      endif
 #ifdef NC_DOUBLE
       ierr = NF_GET_VAR_DOUBLE(unit, nvarid, zdata)
@@ -212,4 +217,8 @@
       ierr = NF_GET_VAR_REAL(unit, nvarid, zdata)
 #endif
+      if (ierr.ne.nf_noerr) then
+        write(*,*) 'datareadnc: error failed loading ',trim(string(k))
+        stop
+      endif
 
 c-----------------------------------------------------------------------
@@ -261,5 +270,9 @@
 c-----------------------------------------------------------------------
 
-      if (k.eq.1) then                    ! albedo
+      if (k.eq.0) then                    ! z0
+         z0(1:iimp1*jjp1)=pfield(1:iimp1*jjp1)*.01
+         ! multiplied by 0.01 to have z0 in m
+         CALL dump2d(iimp1,jjp1,z0,'z0 in m')
+      elseif (k.eq.1) then                    ! albedo
          do i=1,iimp1*jjp1
               alb(i) = pfield(i)
Index: /trunk/LMDZ.MARS/libf/phymars/dustdevil.F
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/dustdevil.F	(revision 223)
+++ /trunk/LMDZ.MARS/libf/phymars/dustdevil.F	(revision 224)
@@ -30,5 +30,5 @@
 #include "comcstfi.h"      
 c#include "comconst.h"        ! TEMPORAIRE AVEC ANLDEVIL !!!!
-#include "planete.h"
+#include "surfdat.h"
 #include "comgeomfi.h"
 #include "tracer.h"
@@ -107,5 +107,5 @@
 c 
          z1= -0.5*13.e3*log(pplev(1,2)/pplev(1,1))
-         Cd = (0.4/log(z1/z0))**2
+         Cd = (0.4/log(z1/z0_default))**2
 
          firstcall=.false.
Index: /trunk/LMDZ.MARS/libf/phymars/phyetat0.F
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/phyetat0.F	(revision 223)
+++ /trunk/LMDZ.MARS/libf/phymars/phyetat0.F	(revision 224)
@@ -464,4 +464,28 @@
       xmax = MAXVAL(emis)
       PRINT*,'Surface emissivity <emis>:', xmin, xmax
+
+!
+! surface roughness length (NB: z0 is a common in surfdat.h)
+!
+      ierr=nf90_inq_varid(nid,"z0",nvarid)
+      if (ierr.ne.nf90_noerr) then
+        write(*,*) 'phyetat0: did not find z0 field!'
+        write(*,*) 'will use constant value of z0_default instead'
+        z0(:)=z0_default
+      else
+        ierr=nf90_get_var(nid,nvarid,z0)
+        if (ierr.ne.nf90_noerr) then
+          write(*,*) 'phyetat0: Failed loading <z0>'
+          write(*,*)trim(nf90_strerror(ierr))
+          call abort
+        endif
+        xmin=1.0E+20
+        xmax=-1.0E+20
+        xmin=minval(z0)
+        xmax=maxval(z0)
+        write(*,*)'Surface roughness <z0>:',xmin,xmax
+      endif
+
+
 c
 c pbl wind variance
Index: /trunk/LMDZ.MARS/libf/phymars/physdem1.F
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/physdem1.F	(revision 223)
+++ /trunk/LMDZ.MARS/libf/phymars/physdem1.F	(revision 224)
@@ -67,6 +67,4 @@
       real pzsig(ngridmx),pzgam(ngridmx),pzthe(ngridmx)
       integer ig
-      INTEGER lnblnk
-      EXTERNAL lnblnk
 
 ! flag which identifies if we are using old tracer names (qsurf01,...)
@@ -180,5 +178,5 @@
 
 c Boundary layer and turbulence
-      tab_cntrl(19) = z0        ! surface roughness (m) ~0.01
+      tab_cntrl(19) = z0_default   ! default surface roughness (m) ~0.01
       tab_cntrl(20) = lmixmin   ! mixing length ~100
       tab_cntrl(21) = emin_turb ! minimal energy ~1.e-8
@@ -518,4 +516,22 @@
       ierr = NF_PUT_VAR_REAL (nid,nvarid,emis)
 #endif
+
+! surface roughness length (z0 is a common in surfdat.h)
+
+      ierr = NF_REDEF (nid)
+#ifdef NC_DOUBLE
+      ierr = NF_DEF_VAR (nid, "z0", NF_DOUBLE, 1, idim2,nvarid)
+#else
+      ierr = NF_DEF_VAR (nid, "z0", NF_FLOAT, 1, idim2,nvarid)
+#endif
+      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 24,
+     .                        "Surface roughness length")
+      ierr = NF_ENDDEF(nid)
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,z0)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,z0)
+#endif
+
 
 c planetary boundary layer
Index: /trunk/LMDZ.MARS/libf/phymars/physiq.F
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/physiq.F	(revision 223)
+++ /trunk/LMDZ.MARS/libf/phymars/physiq.F	(revision 224)
@@ -1072,7 +1072,7 @@
       ENDDO
 
-c    Compute surface stress : (NB: z0 is a common in planete.h)
+c    Compute surface stress : (NB: z0 is a common in surfdat.h)
 c     DO ig=1,ngrid
-c        cd = (0.4/log(zzlay(ig,1)/z0))**2
+c        cd = (0.4/log(zzlay(ig,1)/z0(ig)))**2
 c        zstress(ig) = rho(ig,1)*cd*(zu(ig,1)**2 + zv(ig,1)**2)
 c     ENDDO
Index: /trunk/LMDZ.MARS/libf/phymars/planete.h
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/planete.h	(revision 223)
+++ /trunk/LMDZ.MARS/libf/phymars/planete.h	(revision 224)
@@ -2,13 +2,13 @@
 c INCLUDE planet.h
 
-      COMMON/planet/aphelie,periheli,year_day,peri_day,
-     $       obliquit,
-     $       z0,lmixmin,emin_turb,coefvis,coefir,
-     $       timeperi,e_elips,p_elips,unitastr
+      COMMON/planet/aphelie,periheli,year_day,peri_day,                 &
+     &       obliquit,                                                  &
+     &       lmixmin,emin_turb,coefvis,coefir,                          &
+     &       timeperi,e_elips,p_elips,unitastr
 
-      REAL aphelie,periheli,year_day,peri_day,
-     $     obliquit,
-     $     z0,lmixmin,emin_turb,coefvis,coefir,
-     $       timeperi,e_elips,p_elips,unitastr
+      REAL aphelie,periheli,year_day,peri_day,                          &
+     &     obliquit,                                                    &
+     &     lmixmin,emin_turb,coefvis,coefir,                            &
+     &       timeperi,e_elips,p_elips,unitastr
 
 c-----------------------------------------------------------------------
Index: /trunk/LMDZ.MARS/libf/phymars/surfdat.h
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/surfdat.h	(revision 223)
+++ /trunk/LMDZ.MARS/libf/phymars/surfdat.h	(revision 224)
@@ -4,5 +4,6 @@
      &   iceradius(2) , dtemisice(2),                                   &
      &   zmea(ngridmx),zstd(ngridmx),                                   &
-     &   zsig(ngridmx),zgam(ngridmx),zthe(ngridmx)
+     &   zsig(ngridmx),zgam(ngridmx),zthe(ngridmx),                     &
+     &   z0(ngridmx),z0_default
 
       COMMON/surfdatl/TESicealbedo
@@ -18,2 +19,4 @@
       real iceradius , dtemisice
       real zmea,zstd,zsig,zgam,zthe
+      real z0 ! surface roughness lenght (m)
+      real z0_default ! default (constant over planet) surface roughness (m)
Index: /trunk/LMDZ.MARS/libf/phymars/tabfi.F
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/tabfi.F	(revision 223)
+++ /trunk/LMDZ.MARS/libf/phymars/tabfi.F	(revision 224)
@@ -112,5 +112,5 @@
 c Boundary layer and turbulence
 c----------------------------
-      z0 =  1.e-2               ! surface roughness (m) ~0.01
+      z0_default =  1.e-2       ! surface roughness (m) ~0.01
       emin_turb = 1.e-6         ! minimal energy ~1.e-8
       lmixmin = 30              ! mixing length ~100
@@ -190,5 +190,5 @@
       obliquit = tab_cntrl(tab0+18)
 c boundary layer and turbeulence
-      z0 = tab_cntrl(tab0+19)
+      z0_default = tab_cntrl(tab0+19)
       lmixmin = tab_cntrl(tab0+20)
       emin_turb = tab_cntrl(tab0+21)
@@ -244,5 +244,5 @@
       write(*,5) '(18)       obliquit',tab_cntrl(tab0+18),obliquit
 
-      write(*,6) '(19)             z0',tab_cntrl(tab0+19),z0
+      write(*,6) '(19)     z0_default',tab_cntrl(tab0+19),z0_default
       write(*,6) '(21)      emin_turb',tab_cntrl(tab0+21),emin_turb
       write(*,5) '(20)        lmixmin',tab_cntrl(tab0+20),lmixmin
@@ -277,5 +277,5 @@
       write(*,*)
       write(*,*) '(3)          day_ini : Initial day (=0 at Ls=0)'
-      write(*,*) '(19)              z0 :  surface roughness (m)'
+      write(*,*) '(19)              z0 : default surface roughness (m)'
       write(*,*) '(21)       emin_turb :  minimal energy (PBL)'
       write(*,*) '(20)         lmixmin : mixing length (PBL)'
@@ -317,10 +317,10 @@
 
         else if (modif(1:lnblnk(modif)) .eq. 'z0') then
-          write(*,*) 'current value:',z0
-          write(*,*) 'enter new value:'
- 102      read(*,*,iostat=ierr) z0
+          write(*,*) 'current value (m):',z0_default
+          write(*,*) 'enter new value (m):'
+ 102      read(*,*,iostat=ierr) z0_default
           if(ierr.ne.0) goto 102
           write(*,*) ' '
-          write(*,*) ' z0 (new value):',z0
+          write(*,*) ' z0 (new value):',z0_default
 
         else if (modif(1:lnblnk(modif)) .eq. 'emin_turb') then
@@ -493,5 +493,5 @@
       write(*,5) '(18)       obliquit',tab_cntrl(tab0+18),obliquit
  
-      write(*,6) '(19)             z0',tab_cntrl(tab0+19),z0
+      write(*,6) '(19)     z0_default',tab_cntrl(tab0+19),z0_default
       write(*,6) '(21)      emin_turb',tab_cntrl(tab0+21),emin_turb
       write(*,5) '(20)        lmixmin',tab_cntrl(tab0+20),lmixmin
Index: /trunk/LMDZ.MARS/libf/phymars/testphys1d.F
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/testphys1d.F	(revision 223)
+++ /trunk/LMDZ.MARS/libf/phymars/testphys1d.F	(revision 224)
@@ -122,5 +122,5 @@
 c     Parametres Couche limite et Turbulence 
 c     --------------------------------------
-      z0 =  1.e-2                ! surface roughness (m) ~0.01 
+      z0_default =  1.e-2                ! surface roughness (m) ~0.01 
       emin_turb = 1.e-6          ! energie minimale ~1.e-8
       lmixmin = 30               ! longueur de melange ~100
@@ -399,4 +399,9 @@
       call getin("inertia",inertiedat(1,1))
       write(*,*) " inertia = ",inertiedat(1,1)
+
+      z0(1)=z0_default ! default value for roughness
+      write(*,*) 'Surface roughness length z0 (m)?'
+      call getin("z0",z0(1))
+      write(*,*) " z0 = ",z0(1)
 c
 c  pour le schema d'ondes de gravite
Index: /trunk/LMDZ.MARS/libf/phymars/vdif_cd.F
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/vdif_cd.F	(revision 223)
+++ /trunk/LMDZ.MARS/libf/phymars/vdif_cd.F	(revision 224)
@@ -37,5 +37,5 @@
 
       INTEGER ngrid,nlay
-      REAL pz0
+      REAL pz0(ngrid)
       REAL pg,pz(ngrid,nlay)
       REAL pu(ngrid,nlay),pv(ngrid,nlay)
@@ -62,5 +62,5 @@
 c     DO ig=1,ngrid
 c        zu2=pu(ig)*pu(ig)+pv(ig)*pv(ig)+umin2
-c        pcdv(ig)=pz0*(1.+sqrt(zu2))
+c        pcdv(ig)=pz0(ig)*(1.+sqrt(zu2))
 c        pcdh(ig)=pcdv(ig)
 c     ENDDO
@@ -81,5 +81,5 @@
 c            zri=0.E+0
 c._.
-c         z1=1.+pz(ig)/pz0
+c         z1=1.+pz(ig)/pz0(ig)
 c         zcd0=karman/log(z1)
 c._.         zcd0=zcd0*zcd0*sqrt(zu2)
@@ -101,5 +101,5 @@
 
       DO ig=1,ngrid
-         z1=1.E+0 + pz(ig,1)/pz0
+         z1=1.E+0 + pz(ig,1)/pz0(ig)
          zcd0=karman/log(z1)
          zcd0=zcd0*zcd0
Index: /trunk/LMDZ.MARS/libf/phymars/vdifc.F
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/vdifc.F	(revision 223)
+++ /trunk/LMDZ.MARS/libf/phymars/vdifc.F	(revision 224)
@@ -56,5 +56,6 @@
       REAL co2ice (ngrid), ppopsk(ngrid,nlay)
       logical lecrit
-      REAL pz0
+
+      REAL pz0(ngridmx) ! surface roughness length (m)
 
 c    Traceurs :
