Index: trunk/LMDZ.PLUTO/libf/phypluto/soil_settings.F =================================================================== --- trunk/LMDZ.PLUTO/libf/phypluto/soil_settings.F (revision 3502) +++ (revision ) @@ -1,365 +1,0 @@ - subroutine soil_settings(nid,ngrid,nsoil,tsurf,tsoil,indextime) - -! use netcdf - use comsoil_h, only: layer, mlayer, inertiedat, volcapa, - & lay1_soil, alpha_soil - use iostart, only: inquire_field_ndims, get_var, get_field, - & inquire_field, inquire_dimension_length - 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 - real,intent(in) :: tsurf(ngrid) ! surface temperature - integer,intent(in) :: indextime ! position on time axis -! output: - real,intent(out) :: tsoil(ngrid,nsoil) ! soil temperature - -!====================================================================== -! 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 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 :: oldtsoil - ! for interpolation - real,dimension(:),allocatable :: oldgrid - real,dimension(:),allocatable :: oldval - real,dimension(:),allocatable :: newval - - real malpha,mlay1_soil ! coefficients for building layers - real xmin,xmax ! to display min and max of a field - - real,parameter :: default_volcapa=1.e6 - - logical :: found,ok -!====================================================================== -! 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 - ! if dimlen doesn't match nsoil, then interpolation of - ! soil temperatures and thermal inertia will be requiered - interpol=.true. - endif - - ! allocate oldmlayer - allocate(oldmlayer(dimlen),stat=ierr) - if (ierr.ne.0) then - write(*,*) 'soil_settings: failed allocation of oldmlayer!' - stop - endif - - ! check if olmlayer distribution matches current one - call get_var("soildepth",oldmlayer,found) - if (found) then - malpha=oldmlayer(2)/oldmlayer(1) - if ((abs(malpha-alpha_soil)/alpha_soil).gt.1.e-6) then - ! alpha values are too different, intepolation needed - interpol=.true. - endif - ! check if loaded mid-layer depth value differs from current one - if (abs((oldmlayer(1)-lay1_soil*alpha_soil**(-1./2.))/ - & (lay1_soil*alpha_soil**(-1./2.))).gt.1.e-6) then - interpol=.true. - endif - endif - - if (interpol) then - write(*,*)'soil_settings: Interpolation of soil temperature ', - & 'and thermal inertia will be required!' - endif - -! 1.2 Find out the # of dimensions was defined as using - 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. - ! allocate oldmlayer - if (.not.allocated(oldmlayer)) then - allocate(oldmlayer(dimlen),stat=ierr) - if (ierr.ne.0) then - write(*,*) 'soil_settings: failed allocation of oldmlayer!' - stop - endif - endif - 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) - print*,"mlayer",mlayer - 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 - ! default mlayer distribution, following a power law: - ! mlayer(k)=lay1_soil*alpha_soil**(k-1/2) - do iloop=0,nsoil-1 - mlayer(iloop)=lay1_soil*(alpha_soil**(iloop-0.5)) - enddo - endif -! 1.5 Build layer(); following the same law as mlayer() - ! Assuming layer distribution follows mid-layer law: - ! layer(k)=lay1_soil*alpha_soil**(k-1) - mlay1_soil=sqrt(mlayer(0)*mlayer(1)) - malpha=mlayer(1)/mlayer(0) - ! Check that these values are the same as those prescibed for mlayers - if ((abs(mlay1_soil-lay1_soil)/lay1_soil).gt.1.e-6) then - write(*,*) "soil_settings error: mlay1_soil=",mlay1_soil - write(*,*) " does not match comsoil_h lay1_soil=",lay1_soil - stop - endif - if ((abs(malpha-alpha_soil)/alpha_soil).gt.1.e-6) then - write(*,*) "soil_settings error: malpha=",malpha - write(*,*) " does not match comsoil_h alpha_soil=",alpha_soil - stop - endif - do iloop=1,nsoil - layer(iloop)=mlay1_soil*(malpha**(iloop-1)) - enddo - -! 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 (note: it is declared in comsoil_h) -! ------------------ -! 3.1 Look (again) for thermal inertia data (to reset nvarid) - -! 3.2 Knowing the # of dimensions was defined as using, -! read/build thermal inertia - - 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)) - call get_field("inertiedat",surfinertia,found) - if (.not.found) then - write(*,*) "soil_settings: Failed loading " - call abort - 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!' - stop - endif - endif ! of if (.not.allocated(oldinertiedat)) - call get_field("inertiedat",oldinertiedat,found) - if (.not.found) then - write(*,*) "soil_settings: Failed loading " - call abort - endif - else ! put values in therm_i - call get_field("inertiedat",inertiedat,found) - if (.not.found) then - write(*,*) "soil_settings: Failed loading " - call abort - endif -! endif - endif ! of if (interpol) - endif ! of if (ndims.eq.1) - -! 4. Read soil temperatures -! ------------------------- -! ierr=nf90_inq_varid(nid,"tsoil",nvarid) - ok=inquire_field("tsoil") -! if (ierr.ne.nf90_noerr) then - if (.not.ok) then - write(*,*)'soil_settings: Field not found!' - write(*,*)' => Building from surface values ' - do iloop=1,nsoil - tsoil(:,iloop)=tsurf(:) - enddo - else ! found - if (interpol) then ! put values in oldtsoil - if (.not.allocated(oldtsoil)) then - allocate(oldtsoil(ngrid,dimlen),stat=ierr) - if (ierr.ne.0) then - write(*,*) 'soil_settings: failed allocation of ', - & 'oldtsoil!' - stop - endif - endif - call get_field("tsoil",oldtsoil,found) - if (.not.found) then - write(*,*) "soil_settings: Failed loading " - call abort - 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 - 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 - oldval(1)=tsurf(ig) - oldval(2:dimlen+1)=oldtsoil(ig,1:dimlen) - ! build vertical coordinate - oldgrid(1)=0. ! ground - oldgrid(2:dimlen+1)=oldmlayer(1:dimlen)* - & (inertiedat(ig,1)/volcapa) - ! interpolate - call interp_line(oldgrid,oldval,dimlen+1,mlayer,newval,nsoil) - ! copy result in tsoil - tsoil(ig,:)=newval(:) - 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 - do ig=1,ngrid - ! copy data - oldval(1:dimlen)=oldinertiedat(ig,dimlen) - ! interpolate - call interp_line(oldmlayer,oldval,dimlen,mlayer,newval,nsoil) - !copy result in inertiedat - inertiedat(ig,:)=newval(:) - enddo - - ! soil temperature - ! vertical coordinate - oldgrid(1)=0.0 - oldgrid(2:dimlen+1)=oldmlayer(1:dimlen) - do ig=1,ngrid - ! data - oldval(1)=tsurf(ig) - oldval(2:dimlen+1)=oldtsoil(ig,1:dimlen) - ! interpolate - call interp_line(oldgrid,oldval,dimlen+1,mlayer,newval,nsoil) - ! copy result in inertiedat - tsoil(ig,:)=newval(:) - enddo - - !cleanup - deallocate(oldgrid) - deallocate(oldval) - deallocate(newval) - deallocate(oldinertiedat) - deallocate(oldtsoil) - endif ! of if (interpol) - -! 6. 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 = 1.0E+20 - xmax = -1.0E+20 - xmin = MINVAL(tsoil) - xmax = MAXVAL(tsoil) - write(*,*)'Soil temperature :',xmin,xmax - end Index: trunk/LMDZ.PLUTO/libf/phypluto/soil_settings.F90 =================================================================== --- trunk/LMDZ.PLUTO/libf/phypluto/soil_settings.F90 (revision 3503) +++ trunk/LMDZ.PLUTO/libf/phypluto/soil_settings.F90 (revision 3503) @@ -0,0 +1,366 @@ + subroutine soil_settings(nid,ngrid,nsoil,tsurf,tsoil,indextime) + +! use netcdf + use comsoil_h, only: layer, mlayer, inertiedat, volcapa, & + lay1_soil, alpha_soil + use iostart, only: inquire_field_ndims, get_var, get_field, & + inquire_field, inquire_dimension_length + use mod_phys_lmdz_para, only : is_master + 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 + real,intent(in) :: tsurf(ngrid) ! surface temperature + integer,intent(in) :: indextime ! position on time axis +! output: + real,intent(out) :: tsoil(ngrid,nsoil) ! soil temperature + +!====================================================================== +! 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 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 :: oldtsoil + ! for interpolation + real,dimension(:),allocatable :: oldgrid + real,dimension(:),allocatable :: oldval + real,dimension(:),allocatable :: newval + + real malpha,mlay1_soil ! coefficients for building layers + real xmin,xmax ! to display min and max of a field + + real,parameter :: default_volcapa=1.e6 + + logical :: found,ok +!====================================================================== +! 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 + ! if dimlen doesn't match nsoil, then interpolation of + ! soil temperatures and thermal inertia will be requiered + interpol=.true. + endif + + ! allocate oldmlayer + allocate(oldmlayer(dimlen),stat=ierr) + if (ierr.ne.0) then + if (is_master) write(*,*) 'soil_settings: failed allocation of oldmlayer!' + stop + endif + + ! check if olmlayer distribution matches current one + call get_var("soildepth",oldmlayer,found) + if (found) then + malpha=oldmlayer(2)/oldmlayer(1) + if ((abs(malpha-alpha_soil)/alpha_soil).gt.1.e-6) then + ! alpha values are too different, intepolation needed + interpol=.true. + endif + ! check if loaded mid-layer depth value differs from current one + if (abs((oldmlayer(1)-lay1_soil*alpha_soil**(-1./2.))/& + (lay1_soil*alpha_soil**(-1./2.))).gt.1.e-6) then + interpol=.true. + endif + endif + + if (interpol) then + if (is_master) write(*,*)'soil_settings: Interpolation of soil temperature ', & + 'and thermal inertia will be required!' + endif + +! 1.2 Find out the # of dimensions was defined as using + 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 + if (is_master) write(*,*)'soil_settings: no field expected' + if (is_master) write(*,*)'building one mimicking old definitions' + olddepthdef=.true. + interpol=.true. + ! allocate oldmlayer + if (.not.allocated(oldmlayer)) then + allocate(oldmlayer(dimlen),stat=ierr) + if (ierr.ne.0) then + if (is_master) write(*,*) 'soil_settings: failed allocation of oldmlayer!' + stop + endif + endif + 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 + if (is_master) write(*,*)'soil_settings: Problem while reading ' + endif + else ! put values in mlayer + call get_var("soildepth",mlayer,found) + if (is_master) print*,"mlayer",mlayer + if (.not.found) then + if (is_master) 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 + ! default mlayer distribution, following a power law: + ! mlayer(k)=lay1_soil*alpha_soil**(k-1/2) + do iloop=0,nsoil-1 + mlayer(iloop)=lay1_soil*(alpha_soil**(iloop-0.5)) + enddo + endif +! 1.5 Build layer(); following the same law as mlayer() + ! Assuming layer distribution follows mid-layer law: + ! layer(k)=lay1_soil*alpha_soil**(k-1) + mlay1_soil=sqrt(mlayer(0)*mlayer(1)) + malpha=mlayer(1)/mlayer(0) + ! Check that these values are the same as those prescibed for mlayers + if ((abs(mlay1_soil-lay1_soil)/lay1_soil).gt.1.e-6) then + if (is_master) write(*,*) "soil_settings error: mlay1_soil=",mlay1_soil + if (is_master) write(*,*) " does not match comsoil_h lay1_soil=",lay1_soil + stop + endif + if ((abs(malpha-alpha_soil)/alpha_soil).gt.1.e-6) then + if (is_master) write(*,*) "soil_settings error: malpha=",malpha + if (is_master) write(*,*) " does not match comsoil_h alpha_soil=",alpha_soil + stop + endif + do iloop=1,nsoil + layer(iloop)=mlay1_soil*(malpha**(iloop-1)) + enddo + +! 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 + if (is_master) write(*,*)'soil_settings: Warning, volcapa = ',volcapa + if (is_master) write(*,*)' That doesn t seem right' + if (is_master) write(*,*)' Initializing Volumetric heat capacity to ', & + default_volcapa + volcapa=default_volcapa + endif + +! 3. Thermal inertia (note: it is declared in comsoil_h) +! ------------------ +! 3.1 Look (again) for thermal inertia data (to reset nvarid) + +! 3.2 Knowing the # of dimensions was defined as using, +! read/build thermal inertia + + if (ndims.eq.1) then ! "old 2D-surface" format + if (is_master) write(*,*)'soil_settings: Thermal inertia is only given as surfac& + e data!' + ! Read Surface thermal inertia + allocate(surfinertia(ngrid)) + call get_field("inertiedat",surfinertia,found) + if (.not.found) then + if (is_master) write(*,*) "soil_settings: Failed loading " + call abort + endif + + if (is_master) 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 + if (is_master) write(*,*) 'soil_settings: failed allocation of ', & + 'oldinertiedat!' + stop + endif + endif ! of if (.not.allocated(oldinertiedat)) + call get_field("inertiedat",oldinertiedat,found) + if (.not.found) then + if (is_master) write(*,*) "soil_settings: Failed loading " + call abort + endif + else ! put values in therm_i + call get_field("inertiedat",inertiedat,found) + if (.not.found) then + if (is_master) write(*,*) "soil_settings: Failed loading " + call abort + endif +! endif + endif ! of if (interpol) + endif ! of if (ndims.eq.1) + +! 4. Read soil temperatures +! ------------------------- +! ierr=nf90_inq_varid(nid,"tsoil",nvarid) + ok=inquire_field("tsoil") +! if (ierr.ne.nf90_noerr) then + if (.not.ok) then + if (is_master) write(*,*)'soil_settings: Field not found!' + if (is_master) write(*,*)' => Building from surface values ' + do iloop=1,nsoil + tsoil(:,iloop)=tsurf(:) + enddo + else ! found + if (interpol) then ! put values in oldtsoil + if (.not.allocated(oldtsoil)) then + allocate(oldtsoil(ngrid,dimlen),stat=ierr) + if (ierr.ne.0) then + if (is_master) write(*,*) 'soil_settings: failed allocation of ', & + 'oldtsoil!' + stop + endif + endif + call get_field("tsoil",oldtsoil,found) + if (.not.found) then + if (is_master) write(*,*) "soil_settings: Failed loading " + call abort + endif + else ! put values in tsoil + call get_field("tsoil",tsoil,found,timeindex=indextime) + if (.not.found) then + if (is_master) write(*,*) "soil_settings: Failed loading " + call abort + 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 + oldval(1)=tsurf(ig) + oldval(2:dimlen+1)=oldtsoil(ig,1:dimlen) + ! build vertical coordinate + oldgrid(1)=0. ! ground + oldgrid(2:dimlen+1)=oldmlayer(1:dimlen)*& + (inertiedat(ig,1)/volcapa) + ! interpolate + call interp_line(oldgrid,oldval,dimlen+1,mlayer,newval,nsoil) + ! copy result in tsoil + tsoil(ig,:)=newval(:) + enddo + + ! cleanup + deallocate(oldgrid) + deallocate(oldval) + deallocate(newval) + interpol=.false. ! no need for interpolation any more + endif !of if (olddepthdef) + + if (interpol) then + if (is_master) 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 + do ig=1,ngrid + ! copy data + oldval(1:dimlen)=oldinertiedat(ig,dimlen) + ! interpolate + call interp_line(oldmlayer,oldval,dimlen,mlayer,newval,nsoil) + !copy result in inertiedat + inertiedat(ig,:)=newval(:) + enddo + + ! soil temperature + ! vertical coordinate + oldgrid(1)=0.0 + oldgrid(2:dimlen+1)=oldmlayer(1:dimlen) + do ig=1,ngrid + ! data + oldval(1)=tsurf(ig) + oldval(2:dimlen+1)=oldtsoil(ig,1:dimlen) + ! interpolate + call interp_line(oldgrid,oldval,dimlen+1,mlayer,newval,nsoil) + ! copy result in inertiedat + tsoil(ig,:)=newval(:) + enddo + + !cleanup + deallocate(oldgrid) + deallocate(oldval) + deallocate(newval) + deallocate(oldinertiedat) + deallocate(oldtsoil) + endif ! of if (interpol) + +! 6. Report min and max values of soil temperatures and thermal inertias +! ---------------------------------------------------------------------- + + if (is_master) write(*,*) + if (is_master) write(*,*)'Soil volumetric heat capacity:',volcapa + + xmin = MINVAL(inertiedat) + xmax = MAXVAL(inertiedat) + if (is_master) write(*,*)'Soil thermal inertia :',xmin,xmax + + xmin = 1.0E+20 + xmax = -1.0E+20 + xmin = MINVAL(tsoil) + xmax = MAXVAL(tsoil) + if (is_master) write(*,*)'Soil temperature :',xmin,xmax + end Index: trunk/LMDZ.PLUTO/libf/phypluto/surfini.F =================================================================== --- trunk/LMDZ.PLUTO/libf/phypluto/surfini.F (revision 3502) +++ (revision ) @@ -1,76 +1,0 @@ - SUBROUTINE surfini(ngrid,nq,qsurf,albedo,albedo_bareground, - & albedo_snow_SPECTV,albedo_n2_ice_SPECTV) - - USE surfdat_h, only: albedodat - USE tracer_h, only: igcm_haze - use planetwide_mod, only: planetwide_maxval, planetwide_minval - use radinc_h, only : L_NSPECTV - use callkeys_mod, only : albedosnow, albedon2ice - - IMPLICIT NONE - - -ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc -cccccccccccccc cccccccccccccc -cccccccccccccc Spectral Albedo Initialisation - Routine modified by MT2015. cccccccccccccc -cccccccccccccc cccccccccccccc -ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc - - -c-------------------- -c Declarations: -c-------------------- - - INTEGER,INTENT(IN) :: ngrid - INTEGER,INTENT(IN) :: nq - REAL,INTENT(OUT) :: albedo(ngrid,L_NSPECTV) - REAL,INTENT(OUT) :: albedo_bareground(ngrid) - REAL,INTENT(OUT) :: albedo_snow_SPECTV(L_NSPECTV) - REAL,INTENT(OUT) :: albedo_n2_ice_SPECTV(L_NSPECTV) - REAL,INTENT(IN) :: qsurf(ngrid,nq) ! tracer on surface (kg/m2) - - INTEGER :: ig,nw - REAL :: min_albedo,max_albedo - -c======================================================================= - - ! Step 1 : Initialisation of the Spectral Albedos. - DO nw=1,L_NSPECTV - albedo_snow_SPECTV(nw)=albedosnow - albedo_n2_ice_SPECTV(nw)=albedon2ice - ENDDO - - - ! Step 2 : We get the bare ground albedo from the start files. - DO ig=1,ngrid - albedo_bareground(ig)=0.7 ! albedodat(ig) - DO nw=1,L_NSPECTV - albedo(ig,nw)=0.7 !albedo_bareground(ig) - ENDDO - ENDDO - call planetwide_minval(albedo_bareground,min_albedo) - call planetwide_maxval(albedo_bareground,max_albedo) - write(*,*) 'surfini: minimum bare ground albedo',min_albedo - write(*,*) 'surfini: maximum bare ground albedo',max_albedo - - - ! Step 3 : We modify the albedo considering some N2 at the surface. We dont take into account water ice (this is made in hydrol after the first timestep) ... - if (igcm_haze.ne.0) then - DO ig=1,ngrid - IF (qsurf(ig,igcm_haze) .GT. 1.) THEN ! This was changed by MT2015. Condition for ~1mm of N2 ice deposit. - DO nw=1,L_NSPECTV - albedo(ig,nw)=albedo_n2_ice_SPECTV(nw) - ENDDO - END IF - ENDDO - else - write(*,*) "surfini: No N2 ice tracer on surface ..." - write(*,*) " and therefore no albedo change." - endif - call planetwide_minval(albedo,min_albedo) - call planetwide_maxval(albedo,max_albedo) - write(*,*) 'surfini: minimum corrected initial albedo',min_albedo - write(*,*) 'surfini: maximum corrected initial albedo',max_albedo - - - END Index: trunk/LMDZ.PLUTO/libf/phypluto/surfini.F90 =================================================================== --- trunk/LMDZ.PLUTO/libf/phypluto/surfini.F90 (revision 3503) +++ trunk/LMDZ.PLUTO/libf/phypluto/surfini.F90 (revision 3503) @@ -0,0 +1,76 @@ + SUBROUTINE surfini(ngrid,nq,qsurf,albedo,albedo_bareground,& + albedo_snow_SPECTV,albedo_n2_ice_SPECTV) + + USE surfdat_h, only: albedodat + USE tracer_h, only: igcm_haze + use planetwide_mod, only: planetwide_maxval, planetwide_minval + use radinc_h, only : L_NSPECTV + use callkeys_mod, only : albedosnow, albedon2ice + + IMPLICIT NONE + + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!! !!!!!!! +!!!!!!! Spectral Albedo Initialisation - Routine modified by MT2015. !!!!!!! +!!!!!!! !!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + +!-------------------- +! Declarations: +!-------------------- + + INTEGER,INTENT(IN) :: ngrid + INTEGER,INTENT(IN) :: nq + REAL,INTENT(OUT) :: albedo(ngrid,L_NSPECTV) + REAL,INTENT(OUT) :: albedo_bareground(ngrid) + REAL,INTENT(OUT) :: albedo_snow_SPECTV(L_NSPECTV) + REAL,INTENT(OUT) :: albedo_n2_ice_SPECTV(L_NSPECTV) + REAL,INTENT(IN) :: qsurf(ngrid,nq) ! tracer on surface (kg/m2) + + INTEGER :: ig,nw + REAL :: min_albedo,max_albedo + +!======================================================================= + + ! Step 1 : Initialisation of the Spectral Albedos. + DO nw=1,L_NSPECTV + albedo_snow_SPECTV(nw)=albedosnow + albedo_n2_ice_SPECTV(nw)=albedon2ice + ENDDO + + + ! Step 2 : We get the bare ground albedo from the start files. + DO ig=1,ngrid + albedo_bareground(ig)=0.7 ! albedodat(ig) + DO nw=1,L_NSPECTV + albedo(ig,nw)=0.7 !albedo_bareground(ig) + ENDDO + ENDDO + call planetwide_minval(albedo_bareground,min_albedo) + call planetwide_maxval(albedo_bareground,max_albedo) + write(*,*) 'surfini: minimum bare ground albedo',min_albedo + write(*,*) 'surfini: maximum bare ground albedo',max_albedo + + + ! Step 3 : We modify the albedo considering some N2 at the surface. We dont take into account water ice (this is made in hydrol after the first timestep) ... + if (igcm_haze.ne.0) then + DO ig=1,ngrid + IF (qsurf(ig,igcm_haze) .GT. 1.) THEN ! This was changed by MT2015. Condition for ~1mm of N2 ice deposit. + DO nw=1,L_NSPECTV + albedo(ig,nw)=albedo_n2_ice_SPECTV(nw) + ENDDO + END IF + ENDDO + else + write(*,*) "surfini: No N2 ice tracer on surface ..." + write(*,*) " and therefore no albedo change." + endif + call planetwide_minval(albedo,min_albedo) + call planetwide_maxval(albedo,max_albedo) + write(*,*) 'surfini: minimum corrected initial albedo',min_albedo + write(*,*) 'surfini: maximum corrected initial albedo',max_albedo + + + END