Changeset 3109


Ignore:
Timestamp:
Nov 2, 2023, 10:23:18 AM (13 months ago)
Author:
llange
Message:

MARS PCM
Correction of a small bug in soil_settings: the interpolation on the new
grid for the thermal inertia was not used with the old thermal inertia vector
but only with the last value of the oldthermal inertia
(oldthermalinertia(ig,dimlen) instead of (oldthermalinertia(ig,1:dimlen))
Cleaning of the code
LL

Location:
trunk/LMDZ.MARS
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.MARS/changelog.txt

    r3107 r3109  
    42974297== 26/10/2023 == JBC
    42984298Few small fixes following r3098.
     4299
     4300== 02/11/2023 == LL
     4301Correction of a small bug in soil_settings: the interpolation on the new grid
     4302for the thermal inertia was not used with the old thermal inertia vector but
     4303only with the last value of the oldthermal inertia
     4304Cleaning of the code

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

    r2952 r3109  
    33implicit none
    44! nsoilmx : number of subterranean layers
    5 !EM: old soil routine:      integer, parameter :: nsoilmx = 10
    65  integer, parameter :: nsoilmx = 18
    76

  • trunk/LMDZ.MARS/libf/phymars/dyn1d/init_testphys1d_mod.F90

    r3098 r3109  
    476476    ice_depth = -1 ! default value: no ice
    477477    call getin("subsurface_ice_depth",ice_depth)
    478 
     478    write(*,*) " subsurface_ice_depth = ",ice_depth
    479479
    480480    z0(1) = z0_default ! default value for roughness
    … …  
    654654
    655655    inertiesoil(1,:,1) = inertiedat(1,:)
    656 
    657656    tsoil(:,:,1) = tsurf(1,1) ! soil temperature
    658657endif !(.not. startfiles_1D)

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

    r2951 r3109  
    8989! 1.1 Start by reading how many layers of soil there are
    9090
    91 !       ierr=nf90_inq_dimid(nid,"subsurface_layers",dimid)
    92 !       if (ierr.ne.nf90_noerr) then
    93 !        write(*,*)'soil_settings: Problem reading <subsurface_layers>'
    94 !         write(*,*)trim(nf90_strerror(ierr))
    95 !        call abort
    96 !       endif
    97 
    98 !       ierr=nf90_inquire_dimension(nid,dimid,len=dimlen)
    99 !       if (ierr.ne.nf90_noerr) then
    100 !        write(*,*)'soil_settings: Problem getting <subsurface_layers>',
    101 !     &             'length'
    102 !        write(*,*)trim(nf90_strerror(ierr))
    103 !         call abort
    104 !       endif
    10591        dimlen=inquire_dimension_length("subsurface_layers")
    10692
    … …  
    124110! 1.2 Find out the # of dimensions <inertiedat> was defined as using
    125111!     (in ye old days, thermal inertia was only given at the "surface")
    126       ! Look for thermal inertia data
    127 !      ierr=nf90_inq_varid(nid,"inertiedat",nvarid)
    128 !      if (ierr.NE.nf90_noerr) then
    129 !         write(*,*)'soil_settings: Field <inertiedat> not found!'
    130 !         write(*,*)trim(nf90_strerror(ierr))
    131 !         call abort
    132 !      endif
    133 !
    134 !      ! Read the # of dimensions <inertidat> was defined as using
    135 !      ierr=nf90_inquire_variable(nid,nvarid,ndims=ndims)
    136 !      ! if (ndims.eq.1) then we have the "old 2D-surface" format
     112
    137113      ndims=inquire_field_ndims("inertiedat")
    138114
    … …  
    148124        enddo
    149125      else ! Look for depth
    150 !        ierr=nf90_inq_varid(nid,"soildepth",nvarid)
    151 !        if (ierr.ne.nf90_noerr) then
    152 !          write(*,*)'soil_settings: Field <soildepth> not found!'
    153 !          write(*,*)trim(nf90_strerror(ierr))
    154 !         call abort
    155 !        endif
    156126        ! read <depth> coordinate
    157127        if (interpol) then !put values in oldmlayer
    158 !          ierr=nf90_get_var(nid,nvarid,oldmlayer)
    159 !          if (ierr.ne.nf90_noerr) then
    160 !           write(*,*)'soil_settings: Problem while reading <soildepth>'
    161 !           write(*,*)trim(nf90_strerror(ierr))
    162 !           call abort
    163 !          endif
    164128          call get_var("soildepth",oldmlayer,found)
    165129          if (.not.found) then
    … …  
    167131          endif
    168132        else ! put values in mlayer
    169 !          ierr=nf90_get_var(nid,nvarid,mlayer)
    170 !          if (ierr.ne.nf90_noerr) then
    171 !           write(*,*)'soil_settings: Problem while reading <soildepth>'
    172 !           write(*,*)trim(nf90_strerror(ierr))
    173 !           call abort
    174 !          endif
    175133          call get_var("soildepth",mlayer,found)
    176134          if (.not.found) then
    … …  
    184142      ! default mlayer distribution, following a power law:
    185143      !  mlayer(k)=lay1*alpha**(k-1/2)
    186         lay1=2.e-4
    187         alpha=2
    188         do iloop=0,nsoil-1
    189           mlayer(iloop)=lay1*(alpha**(iloop-0.5))
    190         enddo
     144         lay1=2.e-4
     145         alpha=2
     146         do iloop=0,nsoil-1
     147           mlayer(iloop)=lay1*(alpha**(iloop-0.5))
     148         enddo
    191149      endif
    192150
    … …  
    214172        volcapa=default_volcapa
    215173      endif
    216 ! Look for it
    217 !      ierr=NF_INQ_VARID(nid,"volcapa",nvarid)
    218 !      if (ierr.NE.nf90_noerr) then
    219 !         write(*,*)'soil_settings: Field <volcapa> not found!'
    220 !         write(*,*)'Initializing Volumetric heat capacity to ',
    221 !     &             default_volcapa
    222 !         volcapa=default_volcapa
    223 !      else
    224 !#ifdef NC_DOUBLE
    225 !       ierr = NF_GET_VAR_DOUBLE(nid,nvarid,volcapa)
    226 !#else
    227 !       ierr = NF_GET_VAR_REAL(nid,nvarid,volcapa)
    228 !#endif
    229 !        if (ierr.ne.nf90_noerr) then
    230 !         write(*,*)'soil_settings: Problem while reading <volcapa>'
    231 !         call abort
    232 !       endif
    233 !      endif
    234 
    235 ! 3. Thermal inertia for presetend climate -inertiedat- (note: it is declared in comsoil_h)
     174
     175! 3. Thermal inertia for present day climate -inertiedat- & the one given by the pem -ineritesoil-(note: it is declared in comsoil_h)
    236176! ------------------
    237 
    238 ! 3.1 Look (again) for thermal inertia data (to reset nvarid)
    239 !      ierr=nf90_inq_varid(nid,"inertiedat",nvarid)
    240 !      if (ierr.NE.nf90_noerr) then
    241 !         write(*,*)'soil_settings: Field <inertiedat> not found!'
    242 !         write(*,*)trim(nf90_strerror(ierr))
    243 !         call abort
    244 !      endif
    245 
    246 ! 3.2 Knowing the # of dimensions <inertidat> was defined as using,
     177!     Knowing the # of dimensions <inertidat> was defined as using,
    247178!     read/build thermal inertia
     179
     180! 3.1 Present day - inertiedat
    248181
    249182      if (ndims.eq.1) then ! "old 2D-surface" format
    … …  
    283216           endif
    284217         endif ! of if (.not.allocated(oldinertiedat))
    285 !         ierr=nf90_get_var(nid,nvarid,oldinertiedat)
    286 !        if (ierr.NE.nf90_noerr) then
    287 !         write(*,*)'soil_settings: Problem while reading <inertiedat>'
    288 !         write(*,*)trim(nf90_strerror(ierr))
    289 !         call abort
    290 !        endif
     218
    291219        call get_field("inertiedat",oldinertiedat,found)
    292220        if (.not.found) then
    … …  
    296224        endif
    297225       else ! put values in therm_i
    298 !        ierr=nf90_get_var(nid,nvarid,inertiedat)
    299 !        if (ierr.NE.nf90_noerr) then
    300 !         write(*,*)'soil_settings: Problem while reading <inertiedat>'
    301 !         write(*,*)trim(nf90_strerror(ierr))
    302 !         call abort
    303226         call get_field("inertiedat",inertiedat,found)
    304227         if (.not.found) then
    … …  
    307230     &        "failed loading <inertiedat>",1)
    308231         endif
    309 !        endif
    310232       endif ! of if (interpol)
    311233      endif ! of if (ndims.eq.1)
    312234
    313 
    314 ! 4. Read soil temperatures
    315 ! -------------------------
    316 
    317 !      ierr=nf90_inq_varid(nid,"tsoil",nvarid)
     235! 3.2 Inertie from the PEM
     236
    318237      ok=inquire_field("inertiesoil")
    319 !      if (ierr.ne.nf90_noerr) then
     238
    320239      if (.not.ok) then
    321240        write(*,*)'soil_settings: Field <inertiesoil> not found!'
    … …  
    326245        enddo
    327246      else ! <inertiesoil> found
    328        if (interpol) then ! put values in oldtsoil
     247       if (interpol) then ! put values in oldinertiesoil
    329248         if (.not.allocated(oldinertiesoil)) then
    330249           allocate(oldinertiesoil(ngrid,dimlen,nslope),stat=ierr)
    … …  
    336255           endif
    337256         endif
    338 !        ierr=nf90_get_var(nid,nvarid,oldtsoil)
    339 !        if (ierr.ne.nf90_noerr) then
    340 !        write(*,*)'soil_settings: Problem while reading <tsoil>'
    341 !         write(*,*)trim(nf90_strerror(ierr))
    342 !        call abort
    343 !       endif
     257
    344258         call get_field("inertiesoil",oldinertiesoil,found)
    345259         if (.not.found) then
    … …  
    348262     &          "failed loading <inertiesoil>",1)
    349263         endif
    350        else ! put values in tsoil
    351 !        corner(1)=1
    352 !        corner(2)=1
    353 !        corner(3)=indextime
    354 !        edges(1)=ngrid
    355 !        edges(2)=nsoil
    356 !        edges(3)=1
    357 !        !ierr=nf90_get_var(nid,nvarid,tsoil,corner,edges)
    358 !        ierr=nf90_get_var(nid,nvarid,tsoil)
    359 !        if (ierr.ne.nf90_noerr) then
    360 !        write(*,*)'soil_settings: Problem while reading <tsoil>'
    361 !         write(*,*)trim(nf90_strerror(ierr))
    362 !        call abort
    363 !       endif
     264       else ! put values in inertiesoil
    364265         call get_field("inertiesoil",inertiesoil,found,
    365266     &                 timeindex=indextime)
    … …  
    373274
    374275
    375 
    376 
    377      
    378 ! 5. Read soil temperatures
     276     
     277! 4. Read soil temperatures
    379278! -------------------------
    380279
    381 !      ierr=nf90_inq_varid(nid,"tsoil",nvarid)
    382280      ok=inquire_field("tsoil")
    383 !      if (ierr.ne.nf90_noerr) then
     281
    384282      if (.not.ok) then
    385283        write(*,*)'soil_settings: Field <tsoil> not found!'
    … …  
    401299           endif
    402300         endif
    403 !        ierr=nf90_get_var(nid,nvarid,oldtsoil)
    404 !        if (ierr.ne.nf90_noerr) then
    405 !        write(*,*)'soil_settings: Problem while reading <tsoil>'
    406 !         write(*,*)trim(nf90_strerror(ierr))
    407 !        call abort
    408 !       endif
    409301         call get_field("tsoil",oldtsoil,found)
    410302         if (.not.found) then
    … …  
    414306         endif
    415307       else ! put values in tsoil
    416 !        corner(1)=1
    417 !        corner(2)=1
    418 !        corner(3)=indextime
    419 !        edges(1)=ngrid
    420 !        edges(2)=nsoil
    421 !        edges(3)=1
    422 !        !ierr=nf90_get_var(nid,nvarid,tsoil,corner,edges)
    423 !        ierr=nf90_get_var(nid,nvarid,tsoil)
    424 !        if (ierr.ne.nf90_noerr) then
    425 !        write(*,*)'soil_settings: Problem while reading <tsoil>'
    426 !         write(*,*)trim(nf90_strerror(ierr))
    427 !        call abort
    428 !       endif
    429308         call get_field("tsoil",tsoil,found,timeindex=indextime)
    430309         if (.not.found) then
    … …  
    436315      endif! of if (.not.ok)
    437316
    438 ! 6. If necessary, interpolate soil temperatures and thermal inertias
     317
     318! 5. If necessary, interpolate soil temperatures and thermal inertias
    439319! -------------------------------------------------------------------
    440320
    … …  
    474354      endif
    475355
    476       ! thermal inertia - present day
     356      ! thermal inertia - present day (inertiedat)
    477357      do ig=1,ngrid
    478358          ! copy data
    479           oldval(1:dimlen)=oldinertiedat(ig,dimlen)
     359          oldval(1:dimlen)=oldinertiedat(ig,1:dimlen)
    480360          ! interpolate
    481361          call interp_line(oldmlayer,oldval,dimlen,mlayer,newval,nsoil)
    482362          !copy result in inertiedat
    483363          inertiedat(ig,:)=newval(:)
    484           enddo
    485       ! thermal inertia - general case
    486       do ig=1,ngrid
    487        do islope=1,nslope
     364      ! thermal inertia - general case with PEM (inertie soil)
     365          do islope=1,nslope
    488366          ! copy data
    489              oldval(1:dimlen)=oldinertiesoil(ig,dimlen,islope)
     367             oldval(1:dimlen)=oldinertiesoil(ig,1:dimlen,islope)
    490368          ! interpolate
    491369             call interp_line(oldmlayer,oldval,dimlen,mlayer,newval,nsoil)
    492370          !copy result in inertiedat
    493           inertiesoil(ig,:,islope)=newval(:)
    494           enddo
    495           enddo
     371             inertiesoil(ig,:,islope)=newval(:)
     372           enddo
     373        enddo ! ig
    496374
    497375       
    … …  
    510388          tsoil(ig,:,islope)=newval(:)
    511389          enddo
    512           enddo
     390         enddo
    513391       
    514392        !cleanup
    515       deallocate(oldgrid)
     393         deallocate(oldgrid)
    516394          deallocate(oldval)
    517395          deallocate(newval)
    … …  
    523401
    524402     
    525 ! 7. Load Geothermal Flux
     403! 6. Load Geothermal Flux
    526404! ----------------------------------------------------------------------
    527405
    … …  
    536414
    537415     
    538 ! 8. Report min and max values of soil temperatures and thermal inertias
     416! 7. Report min and max values of soil temperatures and thermal inertias
    539417! ----------------------------------------------------------------------
    540418
Note: See TracChangeset for help on using the changeset viewer.