subroutine soil_settings(nid,ngrid,nsoil,nqsoil,tsurf,tsoil, & qsoil,indextime) ! use netcdf use comsoil_h, only: layer, mlayer, inertiedat, inertiesoil, & volcapa,flux_geo,adsorption_soil, & igcm_h2o_vap_soil,igcm_h2o_ice_soil, & igcm_h2o_vap_ads use iostart, only: inquire_field_ndims, get_var, get_field, & inquire_field, inquire_dimension_length USE comslope_mod, ONLY: nslope implicit none !====================================================================== ! Author: Ehouarn Millour (07/2006) ! ! Purpose: Read and/or initialise soil depths and properties ! ! Modifications: Aug.2010 EM : use NetCDF90 to load variables (enables using ! r4 or r8 restarts independently of having compiled ! the GCM in r4 or r8) ! June 2013 TN : Possibility to read files with a time axis ! ! ! This subroutine reads from a NetCDF file (opened by the caller) ! of "startfi.nc" format. ! The various actions and variable read/initialized are: ! 1. Check out the number of soil layers (in datafile); if it isn't equal ! to nsoil, then some interpolation will be required ! Also check if data in file "startfi.nc" is in older format (ie: ! thermal inertia was depth-independent; and there was no "depth" ! coordinate. ! Read/build layer (and midlayer) depths ! 2. Read volumetric specific heat (or initialise it to default value) ! 3. Read Thermal inertia ! 4. Read soil temperatures ! 5. Interpolate thermal inertia and temperature on the new grid, if ! necessary !====================================================================== !====================================================================== ! arguments ! --------- ! inputs: integer,intent(in) :: nid ! Input Netcdf file ID integer,intent(in) :: ngrid ! # of horizontal grid points integer,intent(in) :: nsoil ! # of soil layers integer,intent(in) :: nqsoil ! # of tracers in the soil real,intent(in) :: tsurf(ngrid,nslope) ! surface temperature [K] integer,intent(in) :: indextime ! position on time axis ! output: real,intent(out) :: tsoil(ngrid,nsoil,nslope) ! soil temperature [K] real,intent(out) :: qsoil(ngrid,nsoil,nqsoil,nslope) ! Tracers in the subsurface [kg/kg] !====================================================================== ! local variables: integer ierr ! status (returned by NetCDF functions) integer nvarid ! ID of NetCDF variable integer dimid ! ID of NetCDF dimension integer dimlen ! length along the "depth" dimension integer ndims ! # of dimensions of read data integer ig,iloop ! loop counters integer islope integer edges(3),corner(3) ! to read a specific time logical :: olddepthdef=.false. ! flag logical :: interpol=.false. ! flag: true if interpolation will be requiered ! to store "old" values real,dimension(:),allocatable :: surfinertia !surface thermal inertia real,dimension(:),allocatable :: oldmlayer real,dimension(:,:),allocatable :: oldinertiedat real,dimension(:,:,:),allocatable :: oldinertiesoil real,dimension(:,:,:),allocatable :: oldtsoil ! for interpolation real,dimension(:),allocatable :: oldgrid real,dimension(:),allocatable :: oldval real,dimension(:),allocatable :: newval real alpha,lay1 ! coefficients for building layers real xmin,xmax ! to display min and max of a field real,parameter :: default_volcapa=1.e6 logical :: found,ok character (len=30):: txt !====================================================================== ! 1. Depth coordinate ! ------------------- ! 1.1 Start by reading how many layers of soil there are dimlen=inquire_dimension_length("subsurface_layers") if (dimlen.ne.nsoil) then write(*,*)'soil_settings: Interpolation of soil temperature ', & 'and thermal inertia will be required!' ! if dimlen doesn't match nsoil, then interpolation of ! soil temperatures and thermal inertia will be requiered interpol=.true. ! allocate oldmlayer if (.not.allocated(oldmlayer)) then allocate(oldmlayer(dimlen),stat=ierr) if (ierr.ne.0) then write(*,*) 'soil_settings: failed allocation of oldmlayer!' call abort_physic("soil_settings", & "failed oldmlayer allocation",1) endif endif endif ! 1.2 Find out the # of dimensions was defined as using ! (in ye old days, thermal inertia was only given at the "surface") ndims=inquire_field_ndims("inertiedat") ! 1.3 Read depths values or set olddepthdef flag and values if (ndims.eq.1) then ! we know that there is none write(*,*)'soil_settings: no field expected' write(*,*)'building one mimicking old definitions' olddepthdef=.true. interpol=.true. do iloop=1,dimlen oldmlayer(iloop)=sqrt(887.75/3.14)*((2.**(iloop-0.5))-1.) enddo else ! Look for depth ! read coordinate if (interpol) then !put values in oldmlayer call get_var("soildepth",oldmlayer,found) if (.not.found) then write(*,*)'soil_settings: Problem while reading ' endif else ! put values in mlayer call get_var("soildepth",mlayer,found) if (.not.found) then write(*,*)'soil_settings: Problem while reading ' endif endif !of if (interpol) endif !of if (ndims.eq.1) ! 1.4 Build mlayer(), if necessary if (interpol) then ! mlayer(k)=lay1*(1+k^(2.9)*(1-exp(-k/20)), PYM grid lay1=2.e-4 do iloop=0,nsoil-1 mlayer(iloop) = lay1*(1+iloop**2.9*(1-exp(-real(iloop)/20.))) enddo endif ! 1.5 Build layer(); do iloop=1,nsoil-1 layer(iloop)=(mlayer(iloop)+mlayer(iloop-1))/2 enddo layer(nsoil)=2*mlayer(nsoil-1) - mlayer(nsoil-2) ! 2. Volumetric heat capacity (note: it is declared in comsoil_h) ! --------------------------- ! "volcapa" is (so far) 0D and written in "controle" table of startfi file ! volcapa is read or set when "controle" is read (see tabfi.F) ! Just in case, we check here that it is not zero. If it is, we ! set it to "default_volcapa" if (volcapa.le.0.0) then write(*,*)'soil_settings: Warning, volcapa = ',volcapa write(*,*)' That doesn t seem right' write(*,*)' Initializing Volumetric heat capacity to ', & default_volcapa volcapa=default_volcapa endif ! 3. Thermal inertia for present day climate -inertiedat- & the one given by the pem -ineritesoil-(note: it is declared in comsoil_h) ! ------------------ ! Knowing the # of dimensions was defined as using, ! read/build thermal inertia ! 3.1 Present day - inertiedat if (ndims.eq.1) then ! "old 2D-surface" format write(*,*)'soil_settings: Thermal inertia is only given as surfac &e data!' ! Read Surface thermal inertia allocate(surfinertia(ngrid)) ! ierr=nf90_get_var(nid,nvarid,surfinertia) ! if (ierr.NE.nf90_noerr) then ! write(*,*)'soil_settings: Problem while reading ' ! write(*,*)trim(nf90_strerror(ierr)) ! call abort ! endif call get_field("inertiedat",surfinertia,found) if (.not.found) then write(*,*) "soil_settings: Failed loading " call abort_physic("soil_settings", & "failed loading ",1) endif write(*,*)' => Building soil thermal inertia (using reference sur &face thermal inertia)' do iloop=1,nsoil inertiedat(:,iloop)=surfinertia(:) enddo deallocate(surfinertia) else ! "3D surface+depth" format if (interpol) then ! put values in oldinertiedat if (.not.allocated(oldinertiedat)) then allocate(oldinertiedat(ngrid,dimlen),stat=ierr) if (ierr.ne.0) then write(*,*) 'soil_settings: failed allocation of ', & 'oldinertiedat!' call abort_physic("soil_settings", & "failed allocation of oldinertiedat",1) endif endif ! of if (.not.allocated(oldinertiedat)) call get_field("inertiedat",oldinertiedat,found) if (.not.found) then write(*,*) "soil_settings: Failed loading " call abort_physic("soil_settings", & "failed loading ",1) endif else ! put values in therm_i call get_field("inertiedat",inertiedat,found) if (.not.found) then write(*,*) "soil_settings: Failed loading " call abort_physic("soil_settings", & "failed loading ",1) endif endif ! of if (interpol) endif ! of if (ndims.eq.1) ! 3.2 Inertie from the PEM ok=inquire_field("inertiesoil") if (.not.ok) then write(*,*)'soil_settings: Field not found!' write(*,*)' => Building from surface values'// & ' ' ! Case 1: No interpolation needed, we just copy past inertiedat if(.not.(interpol)) then do islope=1,nslope inertiesoil(:,:,islope)=inertiedat(:,:) enddo else ! Case 2: Interpolation needed: we copy past old value from inertiedat if (.not.allocated(oldinertiesoil)) then allocate(oldinertiesoil(ngrid,dimlen,nslope),stat=ierr) endif do islope=1,nslope oldinertiesoil(:,:,islope)=oldinertiedat(:,:) enddo endif else ! found if (interpol) then ! put values in oldinertiesoil if (.not.allocated(oldinertiesoil)) then allocate(oldinertiesoil(ngrid,dimlen,nslope),stat=ierr) if (ierr.ne.0) then write(*,*) 'soil_settings: failed allocation of ', & 'oldinertiesoil!' call abort_physic("soil_settings", & "failed allocation of oldinertiesoil",1) endif endif call get_field("inertiesoil",oldinertiesoil,found) if (.not.found) then write(*,*) "soil_settings: Failed loading " call abort_physic("soil_settings", & "failed loading ",1) endif else ! put values in inertiesoil call get_field("inertiesoil",inertiesoil,found, & timeindex=indextime) if (.not.found) then write(*,*) "soil_settings: Failed loading " call abort_physic("soil_settings", & "failed loading ",1) endif endif ! of if (interpol) endif! of if (.not.ok) ! 4. Read soil temperatures ! ------------------------- ok=inquire_field("tsoil") if (.not.ok) then write(*,*)'soil_settings: Field not found!' write(*,*)' => Building from surface values ' ! Case 1: No interpolation needed, we just copy past inertiedat if(.not.(interpol)) then do iloop=1,nsoil do islope=1,nslope tsoil(:,iloop,islope)=tsurf(:,islope) enddo enddo else ! Case 2: Interpolation needed: we copy past old value from inertiedat if (.not.allocated(oldtsoil)) then allocate(oldtsoil(ngrid,dimlen,nslope),stat=ierr) endif do iloop=1,dimlen do islope=1,nslope oldtsoil(:,iloop,islope)=tsurf(:,islope) enddo enddo endif else ! found if (interpol) then ! put values in oldtsoil if (.not.allocated(oldtsoil)) then allocate(oldtsoil(ngrid,dimlen,nslope),stat=ierr) if (ierr.ne.0) then write(*,*) 'soil_settings: failed allocation of ', & 'oldtsoil!' call abort_physic("soil_settings", & "failed allocation of oldtsoil",1) endif endif call get_field("tsoil",oldtsoil,found) if (.not.found) then write(*,*) "soil_settings: Failed loading " call abort_physic("soil_settings", & "failed loading ",1) endif else ! put values in tsoil call get_field("tsoil",tsoil,found,timeindex=indextime) if (.not.found) then write(*,*) "soil_settings: Failed loading " call abort_physic("soil_settings", & "failed loading ",1) endif endif ! of if (interpol) endif! of if (.not.ok) ! 5. If necessary, interpolate soil temperatures and thermal inertias ! ------------------------------------------------------------------- if (olddepthdef) then ! tsoil needs to be interpolated, but not therm_i allocate(oldgrid(dimlen+1)) allocate(oldval(dimlen+1)) allocate(newval(nsoil)) do ig=1,ngrid ! copy values 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 oldgrid(2:dimlen+1)=oldmlayer(1:dimlen)* & (inertiesoil(ig,1,islope)/volcapa) ! interpolate call interp_line(oldgrid,oldval,dimlen+1,mlayer,newval,nsoil) ! copy result in tsoil tsoil(ig,:,islope)=newval(:) enddo enddo ! cleanup deallocate(oldgrid) deallocate(oldval) deallocate(newval) interpol=.false. ! no need for interpolation any more endif !of if (olddepthdef) if (interpol) then write(*,*)'soil_settings: Vertical interpolation along new grid' ! interpolate soil temperatures and thermal inertias if (.not.allocated(oldgrid)) then allocate(oldgrid(dimlen+1)) allocate(oldval(dimlen+1)) allocate(newval(nsoil)) endif ! thermal inertia - present day (inertiedat) do ig=1,ngrid ! copy data oldval(1:dimlen)=oldinertiedat(ig,1:dimlen) ! interpolate call interp_line(oldmlayer,oldval,dimlen,mlayer,newval,nsoil) !copy result in inertiedat inertiedat(ig,:)=newval(:) ! thermal inertia - general case with PEM (inertie soil) do islope=1,nslope ! copy data oldval(1:dimlen)=oldinertiesoil(ig,1:dimlen,islope) ! interpolate call interp_line(oldmlayer,oldval,dimlen,mlayer,newval,nsoil) !copy result in inertiedat inertiesoil(ig,:,islope)=newval(:) enddo enddo ! ig ! soil temperature ! vertical coordinate oldgrid(1)=0.0 oldgrid(2:dimlen+1)=oldmlayer(1:dimlen) do ig=1,ngrid do islope=1,nslope ! data 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,:,islope)=newval(:) enddo enddo !cleanup deallocate(oldgrid) deallocate(oldval) deallocate(newval) deallocate(oldinertiedat) deallocate(oldinertiesoil) deallocate(oldtsoil) endif ! of if (interpol) ! 6. Load Geothermal Flux ! ---------------------------------------------------------------------- call get_field("flux_geo",flux_geo,found,timeindex=indextime) if (.not.found) then write(*,*) "soil_settings: Failed loading " write(*,*) "soil_settings: Will put to 0." flux_geo(:,:) = 0. endif ! 7. Adsorption ! ---------------------------------------------------------------------- if (adsorption_soil) then ! Subsurface water vapor txt="h2o_vap_soil" write(*,*) 'phyetat0: loading subsurface tracer', & ' h2o_vap_soil' call get_field(txt,qsoil(:,:,igcm_h2o_vap_soil,:),found, & indextime) if (.not.found) then write(*,*) "phyetat0: Failed loading <",trim(txt),">" write(*,*) " ",trim(txt)," is set to zero" qsoil(:,:,igcm_h2o_vap_soil,:)= 0. else write(*,*) "phyetat0: suburface tracer <",trim(txt), & "> range:", minval(qsoil(:,:,igcm_h2o_vap_soil,:)), & maxval(qsoil(:,:,igcm_h2o_vap_soil,:)) endif ! Subsurface ice txt="h2o_ice_soil" write(*,*) 'phyetat0: loading subsurface tracer', & ' h2o_ice_soil' call get_field(txt,qsoil(:,:,igcm_h2o_ice_soil,:),found, & indextime) if (.not.found) then write(*,*) "phyetat0: Failed loading <",trim(txt),">" write(*,*) " ",trim(txt)," is set to zero" qsoil(:,:,igcm_h2o_ice_soil,:)= 0. else write(*,*) "phyetat0: suburface tracer <",trim(txt), & "> range:", & minval(qsoil(:,:,igcm_h2o_ice_soil,:)), & maxval(qsoil(:,:,igcm_h2o_ice_soil,:)) endif ! Adsorbed water txt="h2o_vap_ads" write(*,*) 'phyetat0: loading subsurface tracer', & ' h2o_vap_ads' call get_field(txt,qsoil(:,:,igcm_h2o_vap_ads,:),found, & indextime) if (.not.found) then write(*,*) "phyetat0: Failed loading <",trim(txt),">" write(*,*) " ",trim(txt)," is set to zero" qsoil(:,:,igcm_h2o_vap_ads,:)= 0. else write(*,*) "phyetat0: suburface tracer <",trim(txt),"> & range:", & minval(qsoil(:,:,igcm_h2o_vap_ads,:)), & maxval(qsoil(:,:,igcm_h2o_vap_ads,:)) endif endif ! of adsorption_soil ! 8. Report min and max values of soil temperatures and thermal inertias ! ---------------------------------------------------------------------- write(*,*) write(*,*)'Soil volumetric heat capacity:',volcapa xmin = MINVAL(inertiedat) xmax = MAXVAL(inertiedat) write(*,*)'Soil thermal inertia :',xmin,xmax xmin = MINVAL(inertiesoil) xmax = MAXVAL(inertiesoil) write(*,*)'Soil thermal inertia :',xmin,xmax xmin = 1.0E+20 xmax = -1.0E+20 xmin = MINVAL(tsoil) xmax = MAXVAL(tsoil) write(*,*)'Soil temperature :',xmin,xmax end