Changeset 2913 for trunk/LMDZ.MARS/libf

Timestamp:

Mar 14, 2023, 10:07:33 AM (21 months ago)

Author:

romain.vande

Message:

Mars PCM:
Adapt start2archive.F to the subslope parametrisation.
Small correction for some dimensions of variables.
RV

Location:

trunk/LMDZ.MARS/libf

Files:

• dynphy_lonlat/phymars/ini_archive.F (modified) (6 diffs)
• dynphy_lonlat/phymars/start2archive.F (modified) (11 diffs)
• dynphy_lonlat/phymars/write_archive.F (modified) (6 diffs)
• phymars/iostart.F90 (modified) (6 diffs)
• phymars/phyetat0_mod.F90 (modified) (14 diffs)
• phymars/phyredem.F90 (modified) (4 diffs)
• phymars/soil_settings.F (modified) (8 diffs)

trunk/LMDZ.MARS/libf/dynphy_lonlat/phymars/ini_archive.F

@@ -1 +1 @@
 c=======================================================================
-      subroutine ini_archive(nid,idayref,phis,ith,tab_cntrl_fi)
+      subroutine ini_archive(nid,idayref,phis,ith,tab_cntrl_fi
+     &                      ,def_slope,subslope_dist)
 c=======================================================================
 c
@@ -40 +41 @@
       USE serre_mod, ONLY: clon,clat,grossismx,grossismy,dzoomx,dzoomy
       USE ener_mod, ONLY: etot0,ptot0,ztot0,stot0,ang0
+      use comslope_mod, ONLY: nslope
       implicit none
 
@@ -66 +68 @@
       REAL              phis(ip1jmp1)
       real ith(ip1jmp1,nsoilmx)
+      real subslope_dist(ip1jmp1,nslope)
+      real def_slope(nslope+1)
       REAL              tab_cntrl_fi(length)
 
@@ -77 +81 @@
       INTEGER idim_tim
       INTEGER idim_nsoilmx ! "subsurface_layers" dimension ID #
+      INTEGER idim_nslope, idim_nslope_p1
       INTEGER nid,nvarid
       real sig_s(llm),s(llm)
@@ -163 +168 @@
       ierr = NF_DEF_DIM (nid,"interlayer", llmp1, idim_llmp1)
       ierr = NF_DEF_DIM (nid,"Time", NF_UNLIMITED, idim_tim)
+      ierr = NF_DEF_DIM (nid,"nslope", nslope, idim_nslope)
+      ierr = NF_DEF_DIM (nid,"nslope_plus_1",nslope+1,idim_nslope_p1)
 
 c
@@ -513 +520 @@
 #endif
 
+c Put subslope dist
+      dims3(1)=idim_rlonv
+      dims3(2)=idim_rlatu
+      dims3(3)=idim_nslope
+      ierr = NF_REDEF (nid)
+#ifdef NC_DOUBLE
+      ierr = NF_DEF_VAR (nid, "subslope_dist", NF_DOUBLE, 3,
+     .       dims3,nvarid)
+#else
+      ierr = NF_DEF_VAR (nid, "subslope_dist", NF_FLOAT, 3,
+     .       dims3,nvarid)
+#endif
+      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"long_name",13,
+     .    "subslope_dist")
+
+      ierr = NF_ENDDEF(nid)
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,subslope_dist)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,subslope_dist)
+#endif
+
+c Put def_slope
+
+      ierr = NF_REDEF (nid)
+#ifdef NC_DOUBLE
+      ierr = NF_DEF_VAR (nid, "def_slope", NF_DOUBLE, 1,
+     .       [idim_nslope_p1],nvarid)
+#else
+      ierr = NF_DEF_VAR (nid, "def_slope", NF_FLOAT, 1,
+     .       [idim_nslope_p1],nvarid)
+#endif
+      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"long_name",7,"def_slope")
+      ierr = NF_ENDDEF(nid)
+#ifdef NC_DOUBLE
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,def_slope)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,def_slope)
+#endif
+
       PRINT*,'iim,jjm,llm,idayref',iim,jjm,llm,idayref
       PRINT*,'rad,omeg,g,mugaz,kappa',

trunk/LMDZ.MARS/libf/dynphy_lonlat/phymars/start2archive.F

@@ -33 +33 @@
       USE phyetat0_mod, ONLY: phyetat0
       USE exner_hyb_m, ONLY: exner_hyb
+      use comslope_mod, ONLY: nslope,def_slope,def_slope_mean,
+     &                        subslope_dist
+      USE comcstfi_h, only: pi
       implicit none
 
@@ -61 +64 @@
 c Variable Physiques (grille physique)
 c ------------------------------------
-      REAL tsurf(ngridmx)        ! Surface temperature
-      REAL tsoil(ngridmx,nsoilmx) ! Soil temperature
-      REAL watercap(ngridmx)        ! h2o ice layer
-      REAL tauscaling(ngridmx) ! dust conversion factor
-      REAL totcloudfrac(ngridmx) ! sub-grid cloud fraction
+      REAL,ALLOCATABLE :: tsurf(:,:)        ! Surface temperature
+      REAL,ALLOCATABLE :: tsoil(:,:,:) ! Soil temperature
+      REAL,ALLOCATABLE :: watercap(:,:)        ! h2o ice layer
+      REAL :: tauscaling(ngridmx) ! dust conversion factor
+      REAL:: totcloudfrac(ngridmx) ! sub-grid cloud fraction
       REAL q2(ngridmx,llm+1)
-      REAL emis(ngridmx)
-      REAL albedo(ngridmx,2)
+      REAL,ALLOCATABLE :: emis(:,:)
+      REAL,ALLOCATABLE :: albedo(:,:,:)
       REAL wstar(ngridmx)
       INTEGER start,length
@@ -78 +81 @@
 c ------------------------------------
       REAL T(ip1jmp1,llm),us(ip1jmp1,llm),vs(ip1jmp1,llm)
-      REAL tsurfS(ip1jmp1)
-      REAL tsoilS(ip1jmp1,nsoilmx)
+      REAL,ALLOCATABLE :: tsurfS(:,:)
+      REAL,ALLOCATABLE :: tsoilS(:,:,:)
       REAL ithS(ip1jmp1,nsoilmx) ! Soil Thermal Inertia
-      REAL watercapS(ip1jmp1)
-      REAL tauscalingS(ip1jmp1)
-      REAL totcloudfracS(ip1jmp1)
+      REAL,ALLOCATABLE :: watercapS(:,:)
+      REAL :: tauscalingS(ip1jmp1)
+      REAL :: totcloudfracS(ip1jmp1)
       REAL q2S(ip1jmp1,llm+1)
-      REAL,ALLOCATABLE :: qsurfS(:,:)
-      REAL emisS(ip1jmp1)
-      REAL albedoS(ip1jmp1)
+      REAL,ALLOCATABLE :: qsurfS(:,:,:)
+      REAL,ALLOCATABLE :: emisS(:,:)
+      REAL,ALLOCATABLE :: albedoS(:,:)
+      REAL, ALLOCATABLE :: subslope_distS(:,:)
 
 c Variables intermediaires : vent naturel, mais pas coord scalaire
@@ -106 +110 @@
       data  fichier /'startfi'/
 
-      INTEGER ij, l,i,j,isoil,iq
+      INTEGER ij, l,i,j,isoil,iq,islope
       character*80      fichnom
       integer :: ierr,ntime
@@ -133 +137 @@
 
 ! allocate arrays:
-      allocate(q(ip1jmp1,llm,nqtot))
-      allocate(qsurfS(ip1jmp1,nqtot))
-      
+      allocate(q(ip1jmp1,llm,nqtot))      
 
       fichnom = 'start.nc'
@@ -161 +163 @@
       Lmodif=0
 
+      allocate(tsurf(ngridmx,nslope))
+      allocate(tsoil(ngridmx,nsoilmx,nslope))
+      allocate(watercap(ngridmx,nslope))
+      allocate(emis(ngridmx,nslope))
+      allocate(albedo(ngridmx,2,nslope))
+
+      allocate(qsurfS(ip1jmp1,nqtot,nslope))
+      allocate(tsurfS(ip1jmp1,nslope))
+      allocate(tsoilS(ip1jmp1,nsoilmx,nslope))
+      allocate(watercapS(ip1jmp1,nslope))
+      allocate(emisS(ip1jmp1,nslope))
+      allocate(albedoS(ip1jmp1,nslope))
+      allocate(subslope_distS(ip1jmp1,nslope))
+
       CALL phyetat0 (fichnom,0,Lmodif,nsoilmx,ngridmx,llm,nqtot,
      &      day_ini_fi,timefi,tsurf,tsoil,albedo,emis,q2,qsurf,
-     &      tauscaling,totcloudfrac,wstar,watercap)
+     &      tauscaling,totcloudfrac,wstar,watercap,def_slope,
+     &      def_slope_mean,subslope_dist)
 
        ierr = NF_OPEN (fichnom, NF_NOWRITE,nid1)
@@ -257 +274 @@
 c-----------------------------------------------------------------------
 
-      call gr_fi_dyn(1,ngridmx,iip1,jjp1,tsurf,tsurfS)
-      call gr_fi_dyn(1,ngridmx,iip1,jjp1,watercap,watercapS)
-      call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,tsoil,tsoilS)
+      do islope=1,nslope
+      call gr_fi_dyn(1,ngridmx,iip1,jjp1,tsurf(:,islope),
+     &    tsurfS(:,islope))
+      call gr_fi_dyn(1,ngridmx,iip1,jjp1,watercap(:,islope),
+     &    watercapS(:,islope))
+      call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,tsoil(:,:,islope),
+     &    tsoilS(:,:,islope))
       ! Note: thermal inertia "inertiedat" is in comsoil.h
+      call gr_fi_dyn(1,ngridmx,iip1,jjp1,emis(:,islope),
+     &     emisS(:,islope))
+      call gr_fi_dyn(1,ngridmx,iip1,jjp1,albedo(:,1,islope),
+     &   albedoS(:,islope))
+      call gr_fi_dyn(nqtot,ngridmx,iip1,jjp1,qsurf(:,:,islope),
+     &   qsurfS(:,:,islope))
+      call gr_fi_dyn(1,ngridmx,iip1,jjp1,subslope_dist(:,islope),
+     &    subslope_distS(:,islope))
+      enddo
       call gr_fi_dyn(nsoilmx,ngridmx,iip1,jjp1,inertiedat,ithS)
-      call gr_fi_dyn(1,ngridmx,iip1,jjp1,emis,emisS)
-      call gr_fi_dyn(1,ngridmx,iip1,jjp1,albedo(:,1),albedoS)
       call gr_fi_dyn(llm+1,ngridmx,iip1,jjp1,q2,q2S)
-      call gr_fi_dyn(nqtot,ngridmx,iip1,jjp1,qsurf,qsurfS)
       call gr_fi_dyn(1,ngridmx,iip1,jjp1,tauscaling,tauscalingS)
       call gr_fi_dyn(1,ngridmx,iip1,jjp1,totcloudfrac,totcloudfracS)
@@ -283 +310 @@
       DO j=1,jjp1
          DO i=1,iim
+           DO islope=1,nslope
            ptotal=ptotal+aire(i+(iim+1)*(j-1))*ps(i+(iim+1)*(j-1))/g
            co2icetotal = co2icetotal + 
-     &            qsurfS(i+(iim+1)*(j-1),igcm_co2)*aire(i+(iim+1)*(j-1))
+     &        qsurfS(i+(iim+1)*(j-1),igcm_co2,islope)*
+     &        aire(i+(iim+1)*(j-1))*
+     &    subslope_distS(i+(iim+1)*(j-1),islope)/
+     &    cos(pi*def_slope_mean(islope))
+           ENDDO
          ENDDO
       ENDDO
@@ -324 +356 @@
          ierr = NF_CREATE('start_archive.nc', 
      &  IOR(NF_CLOBBER,NF_64BIT_OFFSET), nid)
-         call ini_archive(nid,day_ini,phis,ithS,tab_cntrl_fi)
+         call ini_archive(nid,day_ini,phis,ithS,tab_cntrl_fi,
+     &         def_slope,subslope_distS)
       endif
 
@@ -398 +431 @@
         txt=trim(tname(iq))//"_surf"
         call write_archive(nid,ntime,txt,'Tracer on surface',
-     &  'kg.m-2',2,qsurfS(1,iq))
+     &  'kg.m-2',2,qsurfS(:,iq,:))
       enddo
 
@@ -417 +450 @@
 ! Write soil temperatures tsoil
       call write_archive(nid,ntime,'tsoil','Soil temperature',
-     &     'K',-3,tsoilS)
+     &     'K',-3,tsoilS(:,:,:))
 
 ! Write soil thermal inertia

trunk/LMDZ.MARS/libf/dynphy_lonlat/phymars/write_archive.F

@@ -33 +33 @@
 
       use comsoil_h, only: nsoilmx
+      use comslope_mod, ONLY: nslope
       implicit none
 
@@ -58 +59 @@
 
 c local
-      integer, dimension(4) :: edges,corner,id
+      integer, dimension(5) :: edges,corner,id
       integer :: varid,i,j,l
 c-----------------------------------------------------------------------
@@ -124 +125 @@
         ! get variables' ID, if it exists
         ierr=NF_INQ_VARID(nid,nom,varid)
+
+        if(nom.eq."tsoil") then
+
+         if (ierr.ne.NF_NOERR) then ! variable not defined yet
+          ! build related coordinates
+          ierr=NF_INQ_DIMID(nid,"longitude",id(1))
+          ierr=NF_INQ_DIMID(nid,"latitude",id(2))
+          ierr=NF_INQ_DIMID(nid,"subsurface_layers",id(3))
+          if (ierr.ne.NF_NOERR) then
+           write(*,*)"write_archive: dimension <subsurface_layers>",
+     &               " is missing !!!"
+           call abort
+          endif
+          ierr=NF_INQ_DIMID(nid,"nslope",id(4))
+          if (ierr.ne.NF_NOERR) then
+           write(*,*)"write_archive: dimension <nslope>",
+     &               " is missing !!!"
+           call abort
+          endif
+          ierr=NF_INQ_DIMID(nid,"Time",id(5))
+          
+          ! define the variable
+          write(*,*)"====================="
+          write(*,*)"defining ",nom
+          call def_var(nid,nom,titre,unite,5,id,varid,ierr)
+          
+         endif
+
+        ! build cedges and corners
+        corner(1)=1
+        corner(2)=1
+        corner(3)=1
+        corner(4)=1
+        corner(5)=ntime
+
+        edges(1)=iip1
+        edges(2)=jjp1
+        edges(3)=nsoilmx
+        edges(4)=nslope
+        edges(5)=1
+
+#ifdef NC_DOUBLE
+           ierr= NF_PUT_VARA_DOUBLE(nid,varid,corner,edges,px)
+#else
+           ierr= NF_PUT_VARA_REAL(nid,varid,corner,edges,px)
+#endif
+
+        else
         
          if (ierr.ne.NF_NOERR) then ! variable not defined yet
@@ -160 +209 @@
 #endif
 
+       endif
+
 
 ! For a surface 2D Variable
@@ -169 +220 @@
 
            ierr= NF_INQ_VARID(nid,nom,varid)
+
+           if(nom.eq."tauscaling" .or. nom.eq."totcloudfrac" .or. 
+     &        nom.eq."ps" .or. nom.eq."q2surf") then
+
            if (ierr /= NF_NOERR) then
 !  choix du nom des coordonnees
@@ -203 +258 @@
               call abort
            endif
+ 
+       else
+
+           if (ierr /= NF_NOERR) then
+!  choix du nom des coordonnees
+              ierr= NF_INQ_DIMID(nid,"longitude",id(1))
+              ierr= NF_INQ_DIMID(nid,"latitude",id(2))
+              ierr= NF_INQ_DIMID(nid,"nslope",id(3))
+              ierr= NF_INQ_DIMID(nid,"Time",id(4))
+
+! Creation de la variable si elle n'existait pas
+
+              write (*,*) "====================="
+              write (*,*) "creation de ",nom
+
+              call def_var(nid,nom,titre,unite,4,id,varid,ierr)
+
+           endif
+
+           corner(1)=1
+           corner(2)=1
+           corner(3)=1
+           corner(4)=ntime
+
+           edges(1)=iip1
+           edges(2)=jjp1
+           edges(3)=nslope
+           edges(4)=1
+
+
+#ifdef NC_DOUBLE
+           ierr = NF_PUT_VARA_DOUBLE (nid,varid,corner,edges,px)
+#else         
+           ierr= NF_PUT_VARA_REAL(nid,varid,corner,edges,px)
+#endif     
+
+           if (ierr.ne.NF_NOERR) then
+              write(*,*) "***** PUT_VAR matter in write_archive"
+              write(*,*) "***** with ",nom,nf_STRERROR(ierr)
+              call abort
+           endif
+
+       endif
 
 

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

@@ -15 +15 @@
     INTEGER,SAVE :: idim7 ! "Time" dimension
     INTEGER,SAVE :: idim8 ! "nslope" dimension
-    INTEGER,SAVE :: idim9 ! "inter slope" dimension
+    INTEGER,SAVE :: idim9 ! "nslope_plus_1" dimension
     INTEGER,SAVE :: timeindex ! current time index (for time-dependent fields)
     INTEGER,PARAMETER :: length=100 ! size of tab_cntrl array
@@ -564 +564 @@
       ENDIF
 
-      ierr=NF90_DEF_DIM(nid_restart,"inter slope",nslope+1,idim9)
+      ierr=NF90_DEF_DIM(nid_restart,"inter_slope",nslope+1,idim9)
       IF (ierr/=NF90_NOERR) THEN
         write(*,*)'phyredem: problem defining inter slope dimension'
@@ -660 +660 @@
   
   REAL                           :: field_glo(klon_glo,field_size)
+  REAL                           :: field_glo_reshape(klon_glo,nsoilmx,nslope,timeindex)
   INTEGER                        :: ierr
   INTEGER                        :: nvarid
@@ -804 +805 @@
         endif ! of if (.not.present(timeindex))
 
-      ELSE IF (field_size==nsoilmx) THEN
+      ELSE IF (field_size==nsoilmx .or. field_size==nsoilmx*nslope) THEN
         ! input is a 2D "subsurface field" array
         if (.not.present(time)) then ! for a time-independent field
@@ -828 +829 @@
             ierr=NF90_REDEF(nid_restart)
 #ifdef NC_DOUBLE
+            if(field_name.eq. "tsoil") then
+            ierr=NF90_DEF_VAR(nid_restart,field_name,NF90_DOUBLE,&
+                              (/idim2,idim3,idim8,idim7/),nvarid)
+            else
             ierr=NF90_DEF_VAR(nid_restart,field_name,NF90_DOUBLE,&
                               (/idim2,idim3,idim7/),nvarid)
-#else
+            endif
+#else
+            if(field_name.eq. "tsoil") then
+            ierr=NF90_DEF_VAR(nid_restart,field_name,NF90_FLOAT,&
+                              (/idim2,idim3,idim8,idim7/),nvarid)
+            else
             ierr=NF90_DEF_VAR(nid_restart,field_name,NF90_FLOAT,&
                               (/idim2,idim3,idim7/),nvarid)
+            endif
 #endif
            if (ierr.ne.NF90_NOERR) then
@@ -842 +853 @@
           endif
           ! Write the variable
+          if(field_name.eq. "tsoil") then
+
+          field_glo_reshape=RESHAPE(field_glo, (/klon_glo, nsoilmx, nslope, timeindex/))
+
+          ierr=NF90_PUT_VAR(nid_restart,nvarid,field_glo_reshape,&
+                            start=(/1,1,1,timeindex/))
+          else
           ierr=NF90_PUT_VAR(nid_restart,nvarid,field_glo,&
                             start=(/1,1,timeindex/))
-
+          endif
         endif ! of if (.not.present(time))
 

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

@@ -60 +60 @@
 
 !  outputs:
-  real,intent(out) :: tsurf(ngrid) ! surface temperature
-  real,intent(out) :: tsoil(ngrid,nsoil) ! soil temperature
-  real,intent(out) :: albedo(ngrid,2) ! surface albedo
-  real,intent(out) :: emis(ngrid) ! surface emissivity
+  real,intent(out) :: tsurf(ngrid,nslope) ! surface temperature
+  real,intent(out) :: tsoil(ngrid,nsoil,nslope) ! soil temperature
+  real,intent(out) :: albedo(ngrid,2,nslope) ! surface albedo
+  real,intent(out) :: emis(ngrid,nslope) ! surface emissivity
   real,intent(out) :: q2(ngrid,nlay+1) ! 
-  real,intent(out) :: qsurf(ngrid,nq) ! tracers on surface
+  real,intent(out) :: qsurf(ngrid,nq,nslope) ! tracers on surface
   real,intent(out) :: tauscaling(ngrid) ! dust conversion factor
   real,intent(out) :: totcloudfrac(ngrid) ! total cloud fraction
   real,intent(out) :: wstar(ngrid) ! Max vertical velocity in thermals (m/s)
-  real,intent(out) :: watercap(ngrid) ! h2o_ice_cover
+  real,intent(out) :: watercap(ngrid,nslope) ! h2o_ice_cover
   real, intent(out) :: def_slope(nslope+1) !boundaries for bining of the slopes
   real, intent(out) :: def_slope_mean(nslope)
@@ -143 +143 @@
       write(*,*)'We take default def_slope and subslope dist'
       allocate(default_def_slope(nslope+1))
-      default_def_slope(1) = 0.
-      default_def_slope(2) = 0.
+      default_def_slope(1) = -90.
+      default_def_slope(2) = 90.
       subslope_dist(:,:)=1.
+      def_slope(:)=default_def_slope(:)
      else
        write(*,*)'Variable def_slope is not present in the start and'
@@ -168 +169 @@
       default_def_slope(2) = 0.
       subslope_dist(:,:)=1.
+      def_slope(:)=default_def_slope(:)
     else
      write(*,*)'Without starfi, lets first run with nslope=1'
@@ -516 +518 @@
    endif
 else
-  tsurf(:)=0. ! will be updated afterwards in physiq !
+  tsurf(:,:)=0. ! will be updated afterwards in physiq !
 endif ! of if (startphy_file)
 write(*,*) "phyetat0: Surface temperature <tsurf> range:", &
@@ -523 +525 @@
 ! Surface albedo
 if (startphy_file) then
-   call get_field("albedo",albedo(:,1),found,indextime)
+   call get_field("albedo",albedo(:,1,:),found,indextime)
    if (.not.found) then
      write(*,*) "phyetat0: Failed loading <albedo>"
-     albedo(:,1)=albedodat(:)
-   endif
-else
-   albedo(:,1)=albedodat(:)
+     do islope=1,nslope
+       albedo(:,1,islope)=albedodat(:)
+     enddo
+   endif
+else
+   do islope=1,nslope
+     albedo(:,1,islope)=albedodat(:)
+   enddo
 endif ! of if (startphy_file)
 write(*,*) "phyetat0: Surface albedo <albedo> range:", &
-            minval(albedo(:,1)), maxval(albedo(:,1))
-albedo(:,2)=albedo(:,1)
+            minval(albedo(:,1,:)), maxval(albedo(:,1,:))
+albedo(:,2,:)=albedo(:,1,:)
 
 ! Surface emissivity
@@ -547 +553 @@
   call getin_p("surfemis_without_startfi",surfemis)
   print*,"surfemis_without_startfi",surfemis
-  emis(:)=surfemis
+  emis(:,:)=surfemis
 endif ! of if (startphy_file)
 write(*,*) "phyetat0: Surface emissivity <emis> range:", &
@@ -633 +639 @@
         ! We first check if co2ice exist in the startfi file (old way)
         ! CO2 ice cover
-        call get_field("co2ice",qsurf(:,iq),found,indextime)
+        call get_field("co2ice",qsurf(:,iq,:),found,indextime)
         ! If not, we load the corresponding tracer in qsurf
         if (.not.found) then
-          call get_field(txt,qsurf(:,iq),found,indextime)
+          call get_field(txt,qsurf(:,iq,:),found,indextime)
           if (.not.found) then
             call abort_physic(modname, &
@@ -643 +649 @@
         endif
       else ! (not the co2 tracer)
-        call get_field(txt,qsurf(:,iq),found,indextime)
+        call get_field(txt,qsurf(:,iq,:),found,indextime)
         if (.not.found) then
           write(*,*) "phyetat0: Failed loading <",trim(txt),">"
           write(*,*) "         ",trim(txt)," is set to zero"
-          qsurf(:,iq)=0.
+          qsurf(:,iq,:)=0.
         endif
       endif !endif co2
     else !(not startphy_file)
-      qsurf(:,iq)=0.
+      qsurf(:,iq,:)=0.
     endif ! of if (startphy_file)
     write(*,*) "phyetat0: Surface tracer <",trim(txt),"> range:", &
-                 minval(qsurf(:,iq)), maxval(qsurf(:,iq))
+                 minval(qsurf(:,iq,:)), maxval(qsurf(:,iq,:))
   enddo ! of do iq=1,nq
 
@@ -661 +667 @@
       ! CO2 ice cover
       if (startphy_file) then
-        call get_field("co2ice",qsurf(:,iq),found,indextime)
+        call get_field("co2ice",qsurf(:,iq,:),found,indextime)
       ! If not, we load the corresponding tracer in qsurf
         if (.not.found) then
-          call get_field(txt,qsurf(:,iq),found,indextime)
+          call get_field(txt,qsurf(:,iq,:),found,indextime)
           if (.not.found) then
             call abort_physic(modname, &
@@ -672 +678 @@
       else
        ! If we run without startfile, co2ice is set to 0
-        qsurf(:,iq)=0.
+        qsurf(:,iq,:)=0.
       endif !if (startphy_file)
         write(*,*) "phyetat0: CO2 ice cover <co2ice> range:", &
-            minval(qsurf(:,iq)), maxval(qsurf(:,iq))
+            minval(qsurf(:,iq,:)), maxval(qsurf(:,iq,:))
     endif
 
@@ -685 +691 @@
      write(*,*) "phyetat0: Failed loading <watercap> : ", &
                           "<watercap> is set to zero"
-     watercap(:)=0
+     watercap(:,:)=0
      write(*,*) 'Now transfer negative surface water ice to', &
                 ' watercap !'
@@ -693 +699 @@
        if (txt.eq."h2o_ice") then
          do ig=1,ngrid
-          if (qsurf(ig,iq).lt.0.0) then
-             watercap(ig) = qsurf(ig,iq)
-             qsurf(ig,iq) = 0.0
-          end if
+          do islope=1,nslope
+            if (qsurf(ig,iq,islope).lt.0.0) then
+              watercap(ig,islope) = qsurf(ig,iq,islope)
+              qsurf(ig,iq,islope) = 0.0
+            end if
+          enddo
          end do
        endif
@@ -702 +710 @@
        if (txt.eq."hdo_ice") then
          do ig=1,ngrid
-          if (qsurf(ig,iq).lt.0.0) then
-             qsurf(ig,iq) = 0.0
-          end if
+          do islope=1,nslope
+            if (qsurf(ig,iq,islope).lt.0.0) then
+              qsurf(ig,iq,islope) = 0.0
+            end if
+          enddo
          end do
        endif
@@ -712 +722 @@
    endif ! of if (.not.found)
 else
-   watercap(:)=0
+   watercap(:,:)=0
 endif ! of if (startphy_file) 
 write(*,*) "phyetat0: Surface water ice <watercap> range:", &

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

@@ -171 +171 @@
   use dust_param_mod, only: dustscaling_mode
   use comsoil_h,only: flux_geo
+  use comslope_mod, ONLY: nslope
   implicit none
   
@@ -182 +183 @@
   real,intent(in) :: phystep
   real,intent(in) :: time
-  real,intent(in) :: tsurf(ngrid)
-  real,intent(in) :: tsoil(ngrid,nsoil)
-  real,intent(in) :: albedo(ngrid,2)
-  real,intent(in) :: emis(ngrid)
+  real,intent(in) :: tsurf(ngrid,nslope)
+  real,intent(in) :: tsoil(ngrid,nsoil,nslope)
+  real,intent(in) :: albedo(ngrid,2,nslope)
+  real,intent(in) :: emis(ngrid,nslope)
   real,intent(in) :: q2(ngrid,nlay+1)
-  real,intent(in) :: qsurf(ngrid,nq)
+  real,intent(in) :: qsurf(ngrid,nq,nslope)
   real,intent(in) :: tauscaling(ngrid)
   real,intent(in) :: totcloudfrac(ngrid)
   real,intent(in) :: wstar(ngrid)
-  real,intent(in) :: watercap(ngrid)
+  real,intent(in) :: watercap(ngrid,nslope)
   
   integer :: iq
@@ -219 +220 @@
   
   ! Albedo of the surface
-  call put_field("albedo","Surface albedo",albedo(:,1),time)
+  call put_field("albedo","Surface albedo",albedo(:,1,:),time)
   
   ! Emissivity of the surface
@@ -312 +313 @@
       end if     
 
-      call put_field(trim(txt),"tracer on surface",qsurf(:,iq),time)
+      call put_field(trim(txt),"tracer on surface",qsurf(:,iq,:),time)
     enddo
   endif

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

@@ -5 +5 @@
       use iostart, only: inquire_field_ndims, get_var, get_field,
      &                   inquire_field, inquire_dimension_length
+      USE comslope_mod, ONLY: nslope
       implicit none
 
@@ -41 +42 @@
       integer,intent(in) :: ngrid       ! # of horizontal grid points
       integer,intent(in) :: nsoil       ! # of soil layers
-      real,intent(in) :: tsurf(ngrid)   ! surface temperature
+      real,intent(in) :: tsurf(ngrid,nslope)   ! surface temperature
       integer,intent(in) :: indextime   ! position on time axis
 !  output:
-      real,intent(out) :: tsoil(ngrid,nsoil)    ! soil temperature
+      real,intent(out) :: tsoil(ngrid,nsoil,nslope)     ! soil temperature
 
 !======================================================================
@@ -54 +55 @@
       integer ndims     ! # of dimensions of read <inertiedat> data
       integer ig,iloop  ! loop counters
+      integer islope
       
       integer edges(3),corner(3) ! to read a specific time
@@ -64 +66 @@
       real,dimension(:),allocatable :: oldmlayer
       real,dimension(:,:),allocatable :: oldinertiedat
-      real,dimension(:,:),allocatable :: oldtsoil
+      real,dimension(:,:,:),allocatable :: oldtsoil
       
       ! for interpolation
@@ -317 +319 @@
         write(*,*)' => Building <tsoil> from surface values <tsurf>'
         do iloop=1,nsoil
-          tsoil(:,iloop)=tsurf(:)
+          do islope=1,nslope
+            tsoil(:,iloop,islope)=tsurf(:,islope)
+          enddo
         enddo
       else ! <tsoil> found
        if (interpol) then ! put values in oldtsoil
          if (.not.allocated(oldtsoil)) then
-           allocate(oldtsoil(ngrid,dimlen),stat=ierr)
+           allocate(oldtsoil(ngrid,dimlen,nslope),stat=ierr)
            if (ierr.ne.0) then
              write(*,*) 'soil_settings: failed allocation of ',
@@ -376 +380 @@
         do ig=1,ngrid
           ! copy values
-          oldval(1)=tsurf(ig)
-          oldval(2:dimlen+1)=oldtsoil(ig,1:dimlen)
+          do islope=1,nslope
+          oldval(1)=tsurf(ig,islope)
+          oldval(2:dimlen+1)=oldtsoil(ig,1:dimlen,islope)
           ! build vertical coordinate
           oldgrid(1)=0. ! ground
@@ -385 +390 @@
           call interp_line(oldgrid,oldval,dimlen+1,mlayer,newval,nsoil)
           ! copy result in tsoil
-          tsoil(ig,:)=newval(:)
+          tsoil(ig,:,islope)=newval(:)
+          enddo
         enddo
 
@@ -419 +425 @@
         oldgrid(2:dimlen+1)=oldmlayer(1:dimlen)
         do ig=1,ngrid
+          do islope=1,nslope
           ! data
-          oldval(1)=tsurf(ig)
-          oldval(2:dimlen+1)=oldtsoil(ig,1:dimlen)
+          oldval(1)=tsurf(ig,islope)
+          oldval(2:dimlen+1)=oldtsoil(ig,1:dimlen,islope)
           ! interpolate
           call interp_line(oldgrid,oldval,dimlen+1,mlayer,newval,nsoil)
           ! copy result in inertiedat
-          tsoil(ig,:)=newval(:)
+          tsoil(ig,:,islope)=newval(:)
+          enddo
         enddo