Changeset 2942 for trunk/LMDZ.MARS


Ignore:
Timestamp:
Apr 14, 2023, 6:25:51 PM (20 months ago)
Author:
llange
Message:

Mars PCM
Add the possibility to use a different thermal inertia (field
'inertiesoil') than inertiedat in the PCM (for paleoclimate studies). By defaut, if there is not
inertiesoil, inertiedat is used. Soil_tifeedback still work with
inertiedat
Newstart adapted, start2archive will be modified by Romain
Vandemeulebrouck.
LL

Location:
trunk/LMDZ.MARS/libf
Files:
7 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.MARS/libf/dynphy_lonlat/phymars/newstart.F

    r2910 r2942  
    3131     &                     emis, hmons, summit, base, watercap,
    3232     &               ini_surfdat_h_slope_var,end_surfdat_h_slope_var
    33       use comsoil_h, only: inertiedat, layer, mlayer, nsoilmx, tsoil,
    34      & ini_comsoil_h_slope_var, end_comsoil_h_slope_var
     33      use comsoil_h, only: inertiedat, inertiesoil,layer, mlayer,
     34     & nsoilmx, tsoil,ini_comsoil_h_slope_var, end_comsoil_h_slope_var
    3535      use control_mod, only: day_step, iphysiq, anneeref, planet_type
    3636      use geometry_mod, only: longitude,latitude,cell_area
    … …  
    13871387            inertiedat(1:ngridmx,isoil)=surfithfi(1:ngridmx)
    13881388          enddo
     1389      do islope = 1,nslope
     1390        inertiesoil(:,:,islope) =  inertiedat(:,:)
     1391      enddo
    13891392          write(*,*)'OK: Soil thermal inertia has been reset to referenc
    13901393     &e surface values'

  • trunk/LMDZ.MARS/libf/phymars/co2condens_mod.F

    r2934 r2942  
    736736           ! NB: Updated surface pressure, at grid point 'ngrid', is
    737737           !     ps(ngrid)=pplev(ngrid,1)+pdpsrf(ngrid)*ptimestep
    738              write(*,*) "co2condens: South pole: latitude(ngrid)=",
    739      &                                           latitude(ngrid)
    740738             ztcondsol(ngrid)=
    741739     &          1./(bcond-acond*log(.01*vmr_co2(ngrid,1)*

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

    r2919 r2942  
    88  real,save,allocatable,dimension(:) :: layer      ! soil layer depths
    99  real,save,allocatable,dimension(:) :: mlayer     ! soil mid-layer depths
    10   real,save,allocatable,dimension(:,:) :: inertiedat ! soil thermal inertia
     10  real,save,allocatable,dimension(:,:) :: inertiedat ! soil thermal inertia for present climate
     11  real,save,allocatable,dimension(:,:,:) :: inertiesoil ! soil thermal inertia
    1112  real,save :: volcapa    ! soil volumetric heat capacity
    1213       ! NB: volcapa is read fromn control(35) from physicq start file
    … …  
    3839    allocate(layer(nsoilmx)) !soil layer depths
    3940    allocate(mlayer(0:nsoilmx-1)) ! soil mid-layer depths
    40     allocate(inertiedat(ngrid,nsoilmx)) ! soil thermal inertia
    41  
     41    allocate(inertiedat(ngrid,nsoilmx)) ! soil thermal inertia for present climate
     42    allocate(inertiesoil(ngrid,nsoilmx,nslope)) ! soil thermal inertia
    4243    allocate(tsoil(ngrid,nsoilmx,nslope)) ! soil temperatures
    43 
    4444    allocate(mthermdiff(ngrid,0:nsoilmx-1,nslope))
    4545    allocate(thermdiff(ngrid,nsoilmx-1,nslope))
    … …  
    6060    if (allocated(mlayer)) deallocate(mlayer)
    6161    if (allocated(inertiedat)) deallocate(inertiedat)
     62    if (allocated(inertiesoil)) deallocate(inertiesoil)
    6263    if (allocated(tsoil)) deallocate(tsoil)
    6364    if (allocated(mthermdiff)) deallocate(mthermdiff)
    … …  
    7778 
    7879    allocate(tsoil(ngrid,nsoilmx,nslope)) ! soil temperatures
     80    allocate(inertiesoil(ngrid,nsoilmx,nslope)) ! soil thermal inertia
    7981    allocate(mthermdiff(ngrid,0:nsoilmx-1,nslope))
    8082    allocate(thermdiff(ngrid,nsoilmx-1,nslope))
    … …  
    9193
    9294    if (allocated(tsoil)) deallocate(tsoil)
     95    if (allocated(inertiesoil)) deallocate(inertiesoil)
    9396    if (allocated(mthermdiff)) deallocate(mthermdiff)
    9497    if (allocated(thermdiff)) deallocate(thermdiff)

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

    r2936 r2942  
    66      use mod_grid_phy_lmdz, only : regular_lonlat
    77      use infotrac, only: nqtot, tname, nqperes,nqfils
    8       use comsoil_h, only: volcapa, layer, mlayer, inertiedat, nsoilmx,
    9      &                     flux_geo
     8      use comsoil_h, only: volcapa, layer, mlayer, inertiedat,
     9     &                     inertiesoil,nsoilmx,flux_geo
    1010      use comgeomfi_h, only: sinlat, ini_fillgeom
    1111      use surfdat_h, only: albedodat, z0_default, emissiv, emisice,

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

    r2913 r2942  
    3939  real,intent(in) :: airefi(ngrid)
    4040  real,intent(in) :: alb(ngrid)
    41   real,intent(in) :: ith(ngrid,nsoil)
     41  real,intent(in) :: ith(ngrid,nsoil) ! thermal inertia for present day
    4242  real, intent(in) :: def_slope(nslope+1) !boundaries for bining of the slopes
    4343  real, intent(in) :: subslope_dist(ngrid,nslope) !undermesh statistics
    … …  
    134134     
    135135  ! Soil thermal inertia
    136   call put_field("inertiedat","Soil thermal inertia",ith)
     136  call put_field("inertiedat","Soil thermal inertia - present day",ith)
    137137 
    138138  ! Surface roughness length
    … …  
    159159
    160160subroutine physdem1(filename,nsoil,ngrid,nlay,nq, &
    161                     phystep,time,tsurf,tsoil,albedo,emis,q2,qsurf,&
     161                    phystep,time,tsurf,tsoil,inertiesoil, &
     162                    albedo,emis,q2,qsurf,&
    162163                    tauscaling,totcloudfrac,wstar, &
    163164                    watercap)
    … …  
    185186  real,intent(in) :: tsurf(ngrid,nslope)
    186187  real,intent(in) :: tsoil(ngrid,nsoil,nslope)
     188  real,intent(in) :: inertiesoil(ngrid,nsoil,nslope)
    187189  real,intent(in) :: albedo(ngrid,2,nslope)
    188190  real,intent(in) :: emis(ngrid,nslope)
    … …  
    215217  ! Surface temperature
    216218  call put_field("tsurf","Surface temperature",tsurf,time)
     219
     220  ! Soil temperature
     221  call put_field("inertiesoil","Soil thermal inertia",inertiesoil,time)
    217222 
    218223  ! Soil temperature

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

    r2934 r2942  
    4141     &                      igcm_topdust_number,
    4242     &                      qperemin
    43       use comsoil_h, only: inertiedat, ! soil thermal inertia
     43      use comsoil_h, only: inertiedat, inertiesoil,! dat: soil thermal inertia for present climate, inertiesoil is the TI read in the start
    4444     &                     tsoil, nsoilmx,!number of subsurface layers
    4545     &                     mlayer,layer ! soil mid layer depths
    … …  
    272272      real nuiceco2(ngrid,nlayer) ! Estimated effective variance of the
    273273                                  !   size distribution
    274       REAL inertiesoil(ngrid,nsoilmx,nslope) ! Time varying subsurface
     274      REAL inertiesoil_tifeedback(ngrid,nsoilmx,nslope) ! Time varying subsurface
    275275                                      ! thermal inertia (J.s-1/2.m-2.K-1)
    276276                                      ! (used only when tifeedback=.true.)
    … …  
    540540      real :: qsurf_meshavg(ngrid,nq) ! surface tracer mesh averaged [kg/m^2]
    541541      real :: qsurf_tmp(ngrid,nq) ! temporary qsurf for chimie
    542       real :: inertiedat_tmp(ngrid,nsoilmx,nslope) ! temporary inertiedat for soil.F
    543542      integer :: islope
    544543      real :: zdqsdif_meshavg_tmp(ngrid,nq) ! temporary for dust lifting
    … …  
    586585     &         q2,qsurf,tauscaling,totcloudfrac,wstar,
    587586     &         watercap,def_slope,def_slope_mean,subslope_dist)
    588 
    589        DO islope=1,nslope
    590        sky_slope(islope) = (1.+cos(pi*def_slope_mean(islope)/180.))/2.
    591        END DO
    592587
    593588!     Sky view:
    … …  
    614609      PRINT*,'check: tsurf ',tsurf(1,:),tsurf(ngrid,:)
    615610      PRINT*,'check: tsoil ',tsoil(1,1,:),tsoil(ngrid,nsoilmx,:)
    616       PRINT*,'check: inert ',inertiedat(1,1),inertiedat(ngrid,nsoilmx)
     611      PRINT*,'check: inert ',inertiedat(1,1),inertiedat(ngrid,nsoilmx)     
    617612      PRINT*,'check: midlayer,layer ', mlayer(:),layer(:)
    618613      PRINT*,'check: tracernames ', noms
    … …  
    628623      tauscaling(:)=1.0 !! probably important
    629624      totcloudfrac(:)=1.0
    630       albedo(:,1)=albedodat(:)
    631       albedo(:,2)=albedo(:,1)
    632       watercap(:)=0.0
     625      DO islope = 1,nslope
     626      albedo(:,1,islope)=albedodat(:)
     627      albedo(:,2,islope)=albedo(:,1,islope)
     628      inertiesoil(:,:,islope) = inertiedat(:,:)
     629      watercap(:,:)=0.0
     630      ENDDO
    633631#endif
    634632#ifndef MESOSCALE
    … …  
    640638              !  mlayer(k)=lay1*alpha**(k-1/2)
    641639              lay1=2.e-4
    642               alpha=2
     640                  alpha=2
    643641              do iloop=0,nsoilmx-1
    644                  mlayer(iloop)=lay1*(alpha**(iloop-0.5))
     642                  mlayer(iloop)=lay1*(alpha**(iloop-0.5))
    645643              enddo
    646644              lay1=sqrt(mlayer(0)*mlayer(1))
    … …  
    660658           write(*,*) "Physiq: initializing inertiedat !!"
    661659           inertiedat(:,:)=400.
     660           inertiesoil(:,:,:)=400.
    662661           write(*,*) "Physiq: initializing day_ini to pdat !"
    663662           day_ini=pday
    664663        endif
    665        DO islope = 1,nslope
    666          inertiedat_tmp(:,:,islope) = inertiedat(:,:)
    667        ENDDO
    668664#endif
    669665         if (pday.ne.day_ini) then
    … …  
    695691             DO islope = 1,nslope
    696692                CALL soil_tifeedback(ngrid,nsoilmx,
    697      s               qsurf(:,:,islope),inertiesoil(:,:,islope))
     693     s               qsurf(:,:,islope),
     694     s               inertiesoil_tifeedback(:,:,islope))
    698695             ENDDO
    699                 CALL soil(ngrid,nsoilmx,firstcall,inertiesoil,
     696                CALL soil(ngrid,nsoilmx,firstcall,
     697     s             inertiesoil_tifeedback,
    700698     s             ptimestep,tsurf,tsoil,capcal,fluxgrd)
    701699            ELSE
    702                 CALL soil(ngrid,nsoilmx,firstcall,inertiedat_tmp,
     700                CALL soil(ngrid,nsoilmx,firstcall,inertiesoil,
    703701     s             ptimestep,tsurf,tsoil,capcal,fluxgrd)
    704702            ENDIF ! of IF (tifeedback)
    … …  
    11391137              sl_the = theta_sl(ig)
    11401138              IF (sl_the .ne. 0.) THEN
    1141                 ztim1=fluxsurf_dn_sw(ig,1,1)+fluxsurf_dn_sw(ig,2,1)
     1139                ztim1=fluxsurf_dn_sw(ig,1,iflat)
     1140     &               +fluxsurf_dn_sw(ig,2,iflat)
    11421141                DO l=1,2
    11431142                 sl_lct = ptime*24. + 180.*longitude(ig)/pi/15.
    … …  
    11451144                 sl_lat = 180.*latitude(ig)/pi
    11461145                 sl_tau = tau(ig,1) !il faudrait iaerdust(iaer)
    1147                  sl_alb = albedo(ig,l,1)
     1146                 sl_alb = albedo(ig,l,iflat)
    11481147                 sl_psi = psi_sl(ig)
    1149                  sl_fl0 = fluxsurf_dn_sw(ig,l,1)
     1148                 sl_fl0 = fluxsurf_dn_sw(ig,l,iflat)
    11501149                 sl_di0 = 0.
    11511150                 if ((mu0(ig) .gt. 0.).and.(ztim1.gt.0.)) then
    1152                   sl_di0 = mu0(ig)*(exp(-sl_tau/mu0(ig)))
     1151                          sl_di0 = mu0(ig)*(exp(-sl_tau/mu0(ig)))
    11531152                  sl_di0 = sl_di0*flux_1AU/dist_sol/dist_sol
    11541153                  sl_di0 = sl_di0/ztim1
    1155                   sl_di0 = fluxsurf_dn_sw(ig,l,1)*sl_di0
     1154                  sl_di0 = fluxsurf_dn_sw(ig,l,iflat)*sl_di0
    11561155                 endif
    11571156                 ! you never know (roundup concern...)
    … …  
    11691168              ENDIF
    11701169            ENDDO
    1171 c          ELSE        ! not calling subslope, nslope might be > 1
    1172 c          DO islope = 1,nslope
    1173 c            IF(def_slope_mean(islope).ge.0.) THEN
    1174 c              sl_psi = 0. !Northward slope
    1175 c            ELSE
    1176 c             sl_psi = 180. !Southward slope
    1177 c            ENDIF
    1178 c            sl_the=abs(def_slope_mean(islope))         
    1179 c            IF (sl_the .gt. 1e-6) THEN
    1180 c              DO ig=1,ngrid 
    1181 c                ztim1=fluxsurf_dn_sw(ig,1,islope)
    1182 c     s                +fluxsurf_dn_sw(ig,2,islope)
    1183 c                DO l=1,2
    1184 c                 sl_lct = ptime*24. + 180.*longitude(ig)/pi/15.
    1185 c                 sl_ra  = pi*(1.0-sl_lct/12.)
    1186 c                 sl_lat = 180.*latitude(ig)/pi
    1187 c                 sl_tau = tau(ig,1) !il faudrait iaerdust(iaer)
    1188 c                 sl_alb = albedo(ig,l,islope)
    1189 c                 sl_psi = psi_sl(ig)
    1190 c                 sl_fl0 = fluxsurf_dn_sw(ig,l,islope)
    1191 c                 sl_di0 = 0.
    1192 c                 if (mu0(ig) .gt. 0.) then
    1193 c                         sl_di0 = mu0(ig)*(exp(-sl_tau/mu0(ig)))
    1194 c                  sl_di0 = sl_di0*flux_1AU/dist_sol/dist_sol
    1195 c                  sl_di0 = sl_di0/ztim1
    1196 c                  sl_di0 = fluxsurf_dn_sw(ig,l,islope)*sl_di0
    1197 c                 endif
    1198 c                 ! you never know (roundup concern...)
    1199 c                 if (sl_fl0 .lt. sl_di0) sl_di0=sl_fl0
    1200 c                 !!!!!!!!!!!!!!!!!!!!!!!!!!
    1201 c                 CALL param_slope( mu0(ig), declin, sl_ra, sl_lat,
    1202 c     &                             sl_tau, sl_alb, sl_the, sl_psi,
    1203 c     &                             sl_di0, sl_fl0, sl_flu )
    1204 c                 !!!!!!!!!!!!!!!!!!!!!!!!!!
    1205 c                 fluxsurf_dn_sw(ig,l,islope) = sl_flu
    1206 c                ENDDO
    1207 c              !!! compute correction on IR flux as well
    1208 c
    1209 c              fluxsurf_lw(ig,islope)= fluxsurf_lw(ig,islope)
    1210 c     &                                *sky_slope(islope)
    1211 c              ENDDO
    1212 c            ENDIF
    1213 c          ENDDO ! islope = 1,nslope
     1170          ELSE        ! not calling subslope, nslope might be > 1
     1171          DO islope = 1,nslope
     1172            IF(def_slope_mean(islope).ge.0.) THEN
     1173              sl_psi = 0. !Northward slope
     1174            ELSE
     1175             sl_psi = 180. !Southward slope
     1176            ENDIF
     1177            sl_the=abs(def_slope_mean(islope)) 
     1178            IF (sl_the .gt. 1e-6) THEN
     1179              DO ig=1,ngrid 
     1180                ztim1=fluxsurf_dn_sw(ig,1,islope)
     1181     s                +fluxsurf_dn_sw(ig,2,islope)
     1182                DO l=1,2
     1183                 sl_lct = ptime*24. + 180.*longitude(ig)/pi/15.
     1184                 sl_ra  = pi*(1.0-sl_lct/12.)
     1185                 sl_lat = 180.*latitude(ig)/pi
     1186                 sl_tau = tau(ig,1) !il faudrait iaerdust(iaer)
     1187                 sl_alb = albedo(ig,l,islope)
     1188                 sl_psi = psi_sl(ig)
     1189                 sl_fl0 = fluxsurf_dn_sw(ig,l,islope)
     1190                 sl_di0 = 0.
     1191                 if (mu0(ig) .gt. 0.) then
     1192                          sl_di0 = mu0(ig)*(exp(-sl_tau/mu0(ig)))
     1193                  sl_di0 = sl_di0*flux_1AU/dist_sol/dist_sol
     1194                  sl_di0 = sl_di0/ztim1
     1195                  sl_di0 = fluxsurf_dn_sw(ig,l,islope)*sl_di0
     1196                 endif
     1197                 ! you never know (roundup concern...)
     1198                 if (sl_fl0 .lt. sl_di0) sl_di0=sl_fl0
     1199                 !!!!!!!!!!!!!!!!!!!!!!!!!!
     1200                 CALL param_slope( mu0(ig), declin, sl_ra, sl_lat,
     1201     &                             sl_tau, sl_alb, sl_the, sl_psi,
     1202     &                             sl_di0, sl_fl0, sl_flu )
     1203                 !!!!!!!!!!!!!!!!!!!!!!!!!!
     1204                 fluxsurf_dn_sw(ig,l,islope) = sl_flu
     1205                ENDDO
     1206              !!! compute correction on IR flux as well
     1207
     1208              fluxsurf_lw(ig,islope)= fluxsurf_lw(ig,islope)
     1209     &                                *sky_slope(islope)
     1210              ENDDO
     1211            ENDIF ! sub grid is not flat
     1212          ENDDO ! islope = 1,nslope
    12141213          ENDIF ! callslope
    12151214
    … …  
    23132312         DO islope = 1,nslope
    23142313           CALL soil_tifeedback(ngrid,nsoilmx,
    2315      s               qsurf(:,:,islope),inertiesoil(:,:,islope))
     2314     s               qsurf(:,:,islope),
     2315     s               inertiesoil_tifeedback(:,:,islope))
    23162316         ENDDO
    2317          CALL soil(ngrid,nsoilmx,.false.,inertiesoil,
     2317         CALL soil(ngrid,nsoilmx,.false.,inertiesoil_tifeedback,
    23182318     s          ptimestep,tsurf,tsoil,capcal,fluxgrd)
    23192319        ELSE
    2320          CALL soil(ngrid,nsoilmx,.false.,inertiedat_tmp,
     2320         CALL soil(ngrid,nsoilmx,.false.,inertiesoil,
    23212321     s          ptimestep,tsurf,tsoil,capcal,fluxgrd)
    23222322        ENDIF
    … …  
    25962596          call physdem1("restartfi.nc",nsoilmx,ngrid,nlayer,nq,
    25972597     .                ptimestep,ztime_fin,
    2598      .                tsurf,tsoil,albedo,
     2598     .                tsurf,tsoil,inertiesoil,albedo,
    25992599     .                emis,q2,qsurf,tauscaling,totcloudfrac,wstar,
    26002600     .                watercap)
    … …  
    34623462                 call write_output('soilti',
    34633463     &                       'Soil Thermal Inertia',
    3464      &         'J.s-1/2.m-2.K-1',inertiesoil(:,:,iflat))
     3464     &         'J.s-1/2.m-2.K-1',inertiesoil_tifeedback(:,:,iflat))
    34653465         do islope=1,nslope
    34663466          write(str2(1:2),'(i2.2)') islope
    34673467                 call write_output('soilti_slope'//str2,
    34683468     &                       'Soil Thermal Inertia',
    3469      &          'J.s-1/2.m-2.K-1',inertiesoil(:,:,islope))
     3469     &          'J.s-1/2.m-2.K-1',inertiesoil_tifeedback(:,:,islope))
    34703470         ENDDO
    34713471              endif

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

    r2919 r2942  
    22
    33!      use netcdf
    4       use comsoil_h, only: layer, mlayer, inertiedat, volcapa,flux_geo
     4      use comsoil_h, only: layer, mlayer, inertiedat, inertiesoil,
     5     &                     volcapa,flux_geo
    56      use iostart, only: inquire_field_ndims, get_var, get_field,
    67     &                   inquire_field, inquire_dimension_length
    … …  
    6667      real,dimension(:),allocatable :: oldmlayer
    6768      real,dimension(:,:),allocatable :: oldinertiedat
     69      real,dimension(:,:,:),allocatable :: oldinertiesoil
    6870      real,dimension(:,:,:),allocatable :: oldtsoil
    6971     
    … …  
    231233!      endif
    232234
    233 ! 3. Thermal inertia (note: it is declared in comsoil_h)
     235! 3. Thermal inertia for presetend climate -inertiedat- (note: it is declared in comsoil_h)
    234236! ------------------
    235237
    … …  
    308310       endif ! of if (interpol)
    309311      endif ! of if (ndims.eq.1)
     312
     313
     314! 4. Read soil temperatures
     315! -------------------------
     316
     317!      ierr=nf90_inq_varid(nid,"tsoil",nvarid)
     318      ok=inquire_field("inertiesoil")
     319!      if (ierr.ne.nf90_noerr) then
     320      if (.not.ok) then
     321        write(*,*)'soil_settings: Field <inertiesoil> not found!'
     322        write(*,*)' => Building <inertiesoil> from surface values <inertiedat>'
     323        do iloop=1,nsoil
     324      do islope=1,nslope
     325                inertiesoil(:,:,islope)=inertiedat(:,:)
     326       enddo
     327        enddo
     328      else ! <inertiesoil> found
     329       if (interpol) then ! put values in oldtsoil
     330         if (.not.allocated(oldinertiesoil)) then
     331           allocate(oldinertiesoil(ngrid,dimlen,nslope),stat=ierr)
     332           if (ierr.ne.0) then
     333             write(*,*) 'soil_settings: failed allocation of ',
     334     &                  'oldinertiesoil!'
     335             call abort_physic("soil_settings",
     336     &        "failed allocation of oldinertiesoil",1)
     337           endif
     338         endif
     339!        ierr=nf90_get_var(nid,nvarid,oldtsoil)
     340!        if (ierr.ne.nf90_noerr) then
     341!        write(*,*)'soil_settings: Problem while reading <tsoil>'
     342!         write(*,*)trim(nf90_strerror(ierr))
     343!        call abort
     344!       endif
     345         call get_field("inertiesoil",oldinertiesoil,found)
     346         if (.not.found) then
     347           write(*,*) "soil_settings: Failed loading <inertiesoil>"
     348           call abort_physic("soil_settings",
     349     &          "failed loading <inertiesoil>",1)
     350         endif
     351       else ! put values in tsoil
     352!        corner(1)=1
     353!        corner(2)=1
     354!        corner(3)=indextime
     355!        edges(1)=ngrid
     356!        edges(2)=nsoil
     357!        edges(3)=1
     358!        !ierr=nf90_get_var(nid,nvarid,tsoil,corner,edges)
     359!        ierr=nf90_get_var(nid,nvarid,tsoil)
     360!        if (ierr.ne.nf90_noerr) then
     361!        write(*,*)'soil_settings: Problem while reading <tsoil>'
     362!         write(*,*)trim(nf90_strerror(ierr))
     363!        call abort
     364!       endif
     365         call get_field("inertiesoil",inertiesoil,found,
     366     &                 timeindex=indextime)
     367         if (.not.found) then
     368           write(*,*) "soil_settings: Failed loading <inertiesoil>"
     369           call abort_physic("soil_settings",
     370     &          "failed loading <inertiesoil>",1)
     371         endif
     372       endif ! of if (interpol)
     373      endif! of if (.not.ok)
     374
     375
     376
     377
    310378     
    311 ! 4. Read soil temperatures
     379! 5. Read soil temperatures
    312380! -------------------------
    313381
    … …  
    369437      endif! of if (.not.ok)
    370438
    371 ! 5. If necessary, interpolate soil temperatures and thermal inertias
     439! 6. If necessary, interpolate soil temperatures and thermal inertias
    372440! -------------------------------------------------------------------
    373441
    … …  
    376444        allocate(oldgrid(dimlen+1))
    377445        allocate(oldval(dimlen+1))
    378         allocate(newval(nsoil))
    379 
     446            allocate(newval(nsoil))
    380447        do ig=1,ngrid
    381448          ! copy values
    … …  
    386453          oldgrid(1)=0. ! ground
    387454          oldgrid(2:dimlen+1)=oldmlayer(1:dimlen)*
    388      &                (inertiedat(ig,1)/volcapa)
     455     &                (inertiesoil(ig,1,islope)/volcapa)
    389456          ! interpolate
    390457          call interp_line(oldgrid,oldval,dimlen+1,mlayer,newval,nsoil)
    … …  
    392459          tsoil(ig,:,islope)=newval(:)
    393460          enddo
    394         enddo
    395 
     461          enddo
    396462        ! cleanup
    397         deallocate(oldgrid)
    398         deallocate(oldval)
    399         deallocate(newval)
    400         interpol=.false. ! no need for interpolation any more     
     463          deallocate(oldgrid)
     464          deallocate(oldval)
     465          deallocate(newval)
     466          interpol=.false. ! no need for interpolation any more     
    401467      endif !of if (olddepthdef)
    402 
    403468      if (interpol) then
    404469      write(*,*)'soil_settings: Vertical interpolation along new grid'
    405470      ! interpolate soil temperatures and thermal inertias
    406         if (.not.allocated(oldgrid)) then
     471      if (.not.allocated(oldgrid)) then
    407472          allocate(oldgrid(dimlen+1))
    408473          allocate(oldval(dimlen+1))
    409474          allocate(newval(nsoil))
    410         endif
    411 
    412       ! thermal inertia
    413         do ig=1,ngrid
     475      endif
     476
     477      ! thermal inertia - present day
     478      do ig=1,ngrid
    414479          ! copy data
    415480          oldval(1:dimlen)=oldinertiedat(ig,dimlen)
    … …  
    418483          !copy result in inertiedat
    419484          inertiedat(ig,:)=newval(:)
    420         enddo
     485          enddo
     486      ! thermal inertia - general case
     487      do ig=1,ngrid
     488       do islope=1,nslope
     489          ! copy data
     490             oldval(1:dimlen)=oldinertiesoil(ig,dimlen,islope)
     491          ! interpolate
     492             call interp_line(oldmlayer,oldval,dimlen,mlayer,newval,nsoil)
     493          !copy result in inertiedat
     494          inertiesoil(ig,:,islope)=newval(:)
     495          enddo
     496          enddo
     497
    421498       
    422499      ! soil temperature
    423500        ! vertical coordinate
    424         oldgrid(1)=0.0
    425         oldgrid(2:dimlen+1)=oldmlayer(1:dimlen)
     501          oldgrid(1)=0.0
     502          oldgrid(2:dimlen+1)=oldmlayer(1:dimlen)
    426503        do ig=1,ngrid
    427504          do islope=1,nslope
    … …  
    434511          tsoil(ig,:,islope)=newval(:)
    435512          enddo
    436         enddo
     513          enddo
    437514       
    438515        !cleanup
    439         deallocate(oldgrid)
    440         deallocate(oldval)
    441         deallocate(newval)
    442         deallocate(oldinertiedat)
    443         deallocate(oldtsoil)
     516      deallocate(oldgrid)
     517          deallocate(oldval)
     518          deallocate(newval)
     519          deallocate(oldinertiedat)
     520          deallocate(oldinertiesoil)
     521          deallocate(oldtsoil)
    444522      endif ! of if (interpol)
    445523
    446524
    447525     
    448 ! 6. Load Geothermal Flux
     526! 7. Load Geothermal Flux
    449527! ----------------------------------------------------------------------
    450528
    … …  
    459537
    460538     
    461 ! 7. Report min and max values of soil temperatures and thermal inertias
     539! 8. Report min and max values of soil temperatures and thermal inertias
    462540! ----------------------------------------------------------------------
    463541
    … …  
    468546      xmax = MAXVAL(inertiedat)
    469547      write(*,*)'Soil thermal inertia <inertiedat>:',xmin,xmax
     548
     549      xmin = MINVAL(inertiesoil)
     550      xmax = MAXVAL(inertiesoil)
     551      write(*,*)'Soil thermal inertia <inertiesoil>:',xmin,xmax
    470552
    471553      xmin = 1.0E+20
Note: See TracChangeset for help on using the changeset viewer.