Changeset 224 for trunk/LMDZ.MARS/libf/phymars

Timestamp:

Jul 15, 2011, 9:06:27 AM (13 years ago)

Author:

emillour

Message:

Mars GCM:

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.

Location:

trunk/LMDZ.MARS/libf/phymars

Files:

• datareadnc.F (modified) (9 diffs)
• dustdevil.F (modified) (2 diffs)
• phyetat0.F (modified) (1 diff)
• physdem1.F (modified) (3 diffs)
• physiq.F (modified) (1 diff)
• planete.h (modified) (1 diff)
• surfdat.h (modified) (2 diffs)
• tabfi.F (modified) (6 diffs)
• testphys1d.F (modified) (2 diffs)
• vdif_cd.F (modified) (4 diffs)
• vdifc.F (modified) (1 diff)

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

@@ -1 +1 @@
 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 +64 @@
       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 +96 @@
 
       CHARACTER*20 string
-      DIMENSION string(4)
+      DIMENSION string(0:4)
 
 
@@ -123 +124 @@
         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 +178 @@
 c=======================================================================
 
+      string(0) = 'z0'
       string(1) = 'albedo'
       string(2) = 'thermal'
@@ -192 +192 @@
  
 
-      DO k=1,4
+      DO k=0,4
           write(*,*) 'string',k,string(k)
           
@@ -207 +207 @@
 
       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 +217 @@
       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 +270 @@
 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)

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

@@ -30 +30 @@
 #include "comcstfi.h"      
 c#include "comconst.h"        ! TEMPORAIRE AVEC ANLDEVIL !!!!
-#include "planete.h"
+#include "surfdat.h"
 #include "comgeomfi.h"
 #include "tracer.h"
@@ -107 +107 @@
 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.

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

@@ -464 +464 @@
       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

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

@@ -67 +67 @@
       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 +178 @@
 
 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 +516 @@
       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

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

@@ -1072 +1072 @@
       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

trunk/LMDZ.MARS/libf/phymars/planete.h

@@ -2 +2 @@
 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-----------------------------------------------------------------------

trunk/LMDZ.MARS/libf/phymars/surfdat.h

@@ -4 +4 @@
      &   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 +19 @@
       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)

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

@@ -112 +112 @@
 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 +190 @@
       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 +244 @@
       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 +277 @@
       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 +317 @@
 
         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 +493 @@
       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

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

@@ -122 +122 @@
 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 +399 @@
       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

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

@@ -37 +37 @@
 
       INTEGER ngrid,nlay
-      REAL pz0
+      REAL pz0(ngrid)
       REAL pg,pz(ngrid,nlay)
       REAL pu(ngrid,nlay),pv(ngrid,nlay)
@@ -62 +62 @@
 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 +81 @@
 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 +101 @@
 
       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

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

@@ -56 +56 @@
       REAL co2ice (ngrid), ppopsk(ngrid,nlay)
       logical lecrit
-      REAL pz0
+
+      REAL pz0(ngridmx) ! surface roughness length (m)
 
 c    Traceurs :