source: trunk/LMDZ.GENERIC/libf/phystd/soil_settings.F @ 3562

Last change on this file since 3562 was 3552, checked in by emillour, 10 days ago

Generic PCM

Modify iostart routines to take netcdf index as an argument in order to
be used to read/write in other files than startfi.nc (preparing 1D
restart)

MM

File size: 13.0 KB
RevLine 
[1216]1      subroutine soil_settings(nid,ngrid,nsoil,tsurf,tsoil,indextime)
[787]2
[1216]3!      use netcdf
[1516]4      use comsoil_h, only: layer, mlayer, inertiedat, volcapa,
[1538]5     &                     lay1_soil, alpha_soil
[1216]6      use iostart, only: inquire_field_ndims, get_var, get_field,
7     &                   inquire_field, inquire_dimension_length
[135]8      implicit none
9
10!======================================================================
11!  Author: Ehouarn Millour (07/2006)
12!
13!  Purpose: Read and/or initialise soil depths and properties
14!
[1216]15! Modifications: Aug.2010 EM : use NetCDF90 to load variables (enables using
16!                      r4 or r8 restarts independently of having compiled
17!                      the GCM in r4 or r8)
18!                June 2013 TN : Possibility to read files with a time axis
19!
20!
[135]21!  This subroutine reads from a NetCDF file (opened by the caller)
22!  of "startfi.nc" format.
23!  The various actions and variable read/initialized are:
24!  1. Check out the number of soil layers (in datafile); if it isn't equal
25!     to nsoil, then some interpolation will be required
26!     Also check if data in file "startfi.nc" is in older format (ie:
27!     thermal inertia was depth-independent; and there was no "depth"
28!     coordinate.
29!     Read/build layer (and midlayer) depths
30!  2. Read volumetric specific heat (or initialise it to default value)
31!  3. Read Thermal inertia
32!  4. Read soil temperatures
33!  5. Interpolate thermal inertia and temperature on the new grid, if
34!     necessary
35!======================================================================
36
37!======================================================================
38!  arguments
39!  ---------
40!  inputs:
[1216]41      integer,intent(in) :: nid ! Input Netcdf file ID
42      integer,intent(in) :: ngrid       ! # of horizontal grid points
43      integer,intent(in) :: nsoil       ! # of soil layers
44      real,intent(in) :: tsurf(ngrid)   ! surface temperature
45      integer,intent(in) :: indextime   ! position on time axis
[135]46!  output:
[1216]47      real,intent(out) :: tsoil(ngrid,nsoil)    ! soil temperature
[135]48
49!======================================================================
50! local variables:
51      integer ierr      ! status (returned by NetCDF functions)
52      integer nvarid    ! ID of NetCDF variable
53      integer dimid     ! ID of NetCDF dimension
54      integer dimlen    ! length along the "depth" dimension
55      integer ndims     ! # of dimensions of read <inertiedat> data
56      integer ig,iloop  ! loop counters
[1216]57     
58      integer edges(3),corner(3) ! to read a specific time
[135]59
60      logical :: olddepthdef=.false. ! flag
61      logical :: interpol=.false. ! flag: true if interpolation will be requiered
62
63      ! to store "old" values
64      real,dimension(:),allocatable :: surfinertia !surface thermal inertia
65      real,dimension(:),allocatable :: oldmlayer
66      real,dimension(:,:),allocatable :: oldinertiedat
67      real,dimension(:,:),allocatable :: oldtsoil
68      ! for interpolation
69      real,dimension(:),allocatable :: oldgrid
70      real,dimension(:),allocatable :: oldval
71      real,dimension(:),allocatable :: newval
72
[1538]73      real malpha,mlay1_soil ! coefficients for building layers
[135]74      real xmin,xmax ! to display min and max of a field
75
76      real,parameter :: default_volcapa=1.e6
77
[1216]78      logical :: found,ok
[135]79!======================================================================
80! 1. Depth coordinate
81! -------------------
82! 1.1 Start by reading how many layers of soil there are
[3552]83        dimlen=inquire_dimension_length(nid,"subsurface_layers")
[135]84
85        if (dimlen.ne.nsoil) then
86        ! if dimlen doesn't match nsoil, then interpolation of
87        ! soil temperatures and thermal inertia will be requiered
88          interpol=.true.
[1538]89        endif
90
91        ! allocate oldmlayer 
92        allocate(oldmlayer(dimlen),stat=ierr)
93        if (ierr.ne.0) then
94          write(*,*) 'soil_settings: failed allocation of oldmlayer!'
95          stop
96        endif
97
98        ! check if olmlayer distribution matches current one
[3552]99        call get_var(nid,"soildepth",oldmlayer,found)
[1538]100        if (found) then
101          malpha=oldmlayer(2)/oldmlayer(1)
102          if ((abs(malpha-alpha_soil)/alpha_soil).gt.1.e-6) then
103            ! alpha values are too different, intepolation needed
104            interpol=.true.
[1350]105          endif
[1538]106          ! check if loaded mid-layer depth value differs from current one
107          if (abs((oldmlayer(1)-lay1_soil*alpha_soil**(-1./2.))/
108     &             (lay1_soil*alpha_soil**(-1./2.))).gt.1.e-6) then
109            interpol=.true.
110          endif
111        endif
112
113        if (interpol) then
114          write(*,*)'soil_settings: Interpolation of soil temperature ',
115     &              'and thermal inertia will be required!'
116        endif
117       
[135]118! 1.2 Find out the # of dimensions <inertiedat> was defined as using
[3552]119      ndims=inquire_field_ndims(nid,"inertiedat")
[135]120! 1.3 Read depths values or set olddepthdef flag and values
121      if (ndims.eq.1) then ! we know that there is none
122        write(*,*)'soil_settings: no <soildepth> field expected'
123        write(*,*)'building one mimicking old definitions'
124        olddepthdef=.true.
125        interpol=.true.
126        ! allocate oldmlayer
127        if (.not.allocated(oldmlayer)) then
128          allocate(oldmlayer(dimlen),stat=ierr)
129          if (ierr.ne.0) then
130            write(*,*) 'soil_settings: failed allocation of oldmlayer!'
131            stop
132          endif
133        endif
134        do iloop=1,dimlen
135          oldmlayer(iloop)=sqrt(887.75/3.14)*((2.**(iloop-0.5))-1.)
136        enddo
137      else ! Look for depth
138        ! read <depth> coordinate
139        if (interpol) then !put values in oldmlayer
[3552]140          call get_var(nid,"soildepth",oldmlayer,found)
[1216]141          if (.not.found) then
142            write(*,*)'soil_settings: Problem while reading <soildepth>'
[135]143          endif
144        else ! put values in mlayer
[3552]145          call get_var(nid,"soildepth",mlayer,found)
[1297]146          print*,"mlayer",mlayer
147          if (.not.found) then
[1216]148            write(*,*)'soil_settings: Problem while reading <soildepth>'
[135]149          endif
150        endif !of if (interpol)
151      endif !of if (ndims.eq.1)
152! 1.4 Build mlayer(), if necessary
153      if (interpol) then
154      ! default mlayer distribution, following a power law:
[1538]155      !  mlayer(k)=lay1_soil*alpha_soil**(k-1/2)
[135]156        do iloop=0,nsoil-1
[1538]157          mlayer(iloop)=lay1_soil*(alpha_soil**(iloop-0.5))
[1297]158        enddo
[135]159      endif
160! 1.5 Build layer(); following the same law as mlayer()
161      ! Assuming layer distribution follows mid-layer law:
[1538]162      ! layer(k)=lay1_soil*alpha_soil**(k-1)
163      mlay1_soil=sqrt(mlayer(0)*mlayer(1))
[1516]164      malpha=mlayer(1)/mlayer(0)
165      ! Check that these values are the same as those prescibed for mlayers
[1538]166      if ((abs(mlay1_soil-lay1_soil)/lay1_soil).gt.1.e-6) then
167         write(*,*) "soil_settings error: mlay1_soil=",mlay1_soil
168         write(*,*) " does not match comsoil_h lay1_soil=",lay1_soil
[1516]169         stop
170      endif
[1538]171      if ((abs(malpha-alpha_soil)/alpha_soil).gt.1.e-6) then
[1516]172         write(*,*) "soil_settings error: malpha=",malpha
[1538]173         write(*,*) " does not match comsoil_h alpha_soil=",alpha_soil
[1516]174         stop
175      endif
[135]176      do iloop=1,nsoil
[1538]177        layer(iloop)=mlay1_soil*(malpha**(iloop-1))
[135]178      enddo
[1516]179
[1216]180! 2. Volumetric heat capacity (note: it is declared in comsoil_h)
[135]181! ---------------------------
182! "volcapa" is (so far) 0D and written in "controle" table of startfi file
183! volcapa is read or set when "controle" is read (see tabfi.F)
184! Just in case, we check here that it is not zero. If it is, we
185! set it to "default_volcapa"
186
187      if (volcapa.le.0.0) then
188        write(*,*)'soil_settings: Warning, volcapa = ',volcapa
189        write(*,*)'               That doesn t seem right'
190        write(*,*)'        Initializing Volumetric heat capacity to ',
191     &             default_volcapa
192        volcapa=default_volcapa
193      endif
194
[1216]195! 3. Thermal inertia (note: it is declared in comsoil_h)
[135]196! ------------------
197! 3.1 Look (again) for thermal inertia data (to reset nvarid)
198
199! 3.2 Knowing the # of dimensions <inertidat> was defined as using,
200!     read/build thermal inertia
201
202      if (ndims.eq.1) then ! "old 2D-surface" format
203       write(*,*)'soil_settings: Thermal inertia is only given as surfac
204     &e data!'
205       ! Read Surface thermal inertia
206       allocate(surfinertia(ngrid))
[3552]207       call get_field(nid,"inertiedat",surfinertia,found)
[1216]208       if (.not.found) then
209         write(*,*) "soil_settings: Failed loading <inertiedat>"
[135]210         call abort
[1216]211       endif
[135]212       
213       write(*,*)' => Building soil thermal inertia (using reference sur
214     &face thermal inertia)'
215       do iloop=1,nsoil
216         inertiedat(:,iloop)=surfinertia(:)
217       enddo
218       deallocate(surfinertia)
219
220      else ! "3D surface+depth" format
221       if (interpol) then ! put values in oldinertiedat
222         if (.not.allocated(oldinertiedat)) then
223           allocate(oldinertiedat(ngrid,dimlen),stat=ierr)
224           if (ierr.ne.0) then
225            write(*,*) 'soil_settings: failed allocation of ',
226     &                 'oldinertiedat!'
227            stop
228           endif
[1216]229         endif ! of if (.not.allocated(oldinertiedat))
[3552]230        call get_field(nid,"inertiedat",oldinertiedat,found)
[1216]231        if (.not.found) then
232          write(*,*) "soil_settings: Failed loading <inertiedat>"
233          call abort
[135]234        endif
235       else ! put values in therm_i
[3552]236         call get_field(nid,"inertiedat",inertiedat,found)
[1216]237         if (.not.found) then
238           write(*,*) "soil_settings: Failed loading <inertiedat>"
239           call abort
240         endif
241!        endif
[135]242       endif ! of if (interpol)
243      endif ! of if (ndims.eq.1)
244     
245! 4. Read soil temperatures
246! -------------------------
[1216]247!      ierr=nf90_inq_varid(nid,"tsoil",nvarid)
[3552]248      ok=inquire_field(nid,"tsoil")
[1216]249!      if (ierr.ne.nf90_noerr) then
250      if (.not.ok) then
[135]251        write(*,*)'soil_settings: Field <tsoil> not found!'
252        write(*,*)' => Building <tsoil> from surface values <tsurf>'
253        do iloop=1,nsoil
254          tsoil(:,iloop)=tsurf(:)
255        enddo
256      else ! <tsoil> found
257       if (interpol) then ! put values in oldtsoil
258         if (.not.allocated(oldtsoil)) then
259           allocate(oldtsoil(ngrid,dimlen),stat=ierr)
260           if (ierr.ne.0) then
261             write(*,*) 'soil_settings: failed allocation of ',
262     &                  'oldtsoil!'
263             stop
264           endif
265         endif
[3552]266         call get_field(nid,"tsoil",oldtsoil,found)
[1216]267         if (.not.found) then
268           write(*,*) "soil_settings: Failed loading <tsoil>"
269           call abort
270         endif
[135]271       else ! put values in tsoil
[3552]272         call get_field(nid,"tsoil",tsoil,found,
273     &                  timeindex=indextime)
[1216]274         if (.not.found) then
275           write(*,*) "soil_settings: Failed loading <tsoil>"
276           call abort
277         endif
[135]278       endif ! of if (interpol)
[1216]279      endif! of if (.not.ok)
[135]280
281! 5. If necessary, interpolate soil temperatures and thermal inertias
282! -------------------------------------------------------------------
283      if (olddepthdef) then
284      ! tsoil needs to be interpolated, but not therm_i
285        allocate(oldgrid(dimlen+1))
286        allocate(oldval(dimlen+1))
287        allocate(newval(nsoil))
288
289        do ig=1,ngrid
290          ! copy values
291          oldval(1)=tsurf(ig)
292          oldval(2:dimlen+1)=oldtsoil(ig,1:dimlen)
293          ! build vertical coordinate
294          oldgrid(1)=0. ! ground
295          oldgrid(2:dimlen+1)=oldmlayer(1:dimlen)*
296     &                (inertiedat(ig,1)/volcapa)
297          ! interpolate
298          call interp_line(oldgrid,oldval,dimlen+1,mlayer,newval,nsoil)
299          ! copy result in tsoil
300          tsoil(ig,:)=newval(:)
301        enddo
302
303        ! cleanup
304        deallocate(oldgrid)
305        deallocate(oldval)
306        deallocate(newval)
307        interpol=.false. ! no need for interpolation any more     
308      endif !of if (olddepthdef)
309
310      if (interpol) then
311      write(*,*)'soil_settings: Vertical interpolation along new grid'
312      ! interpolate soil temperatures and thermal inertias
313        if (.not.allocated(oldgrid)) then
314          allocate(oldgrid(dimlen+1))
315          allocate(oldval(dimlen+1))
316          allocate(newval(nsoil))
317        endif
318
319      ! thermal inertia
320        do ig=1,ngrid
321          ! copy data
322          oldval(1:dimlen)=oldinertiedat(ig,dimlen)
323          ! interpolate
324          call interp_line(oldmlayer,oldval,dimlen,mlayer,newval,nsoil)
325          !copy result in inertiedat
326          inertiedat(ig,:)=newval(:)
327        enddo
328       
329      ! soil temperature
330        ! vertical coordinate
331        oldgrid(1)=0.0
332        oldgrid(2:dimlen+1)=oldmlayer(1:dimlen)
333        do ig=1,ngrid
334          ! data
335          oldval(1)=tsurf(ig)
336          oldval(2:dimlen+1)=oldtsoil(ig,1:dimlen)
337          ! interpolate
338          call interp_line(oldgrid,oldval,dimlen+1,mlayer,newval,nsoil)
339          ! copy result in inertiedat
340          tsoil(ig,:)=newval(:)
341        enddo
342       
343        !cleanup
344        deallocate(oldgrid)
345        deallocate(oldval)
346        deallocate(newval)
347        deallocate(oldinertiedat)
348        deallocate(oldtsoil)
349      endif ! of if (interpol)
350     
351! 6. Report min and max values of soil temperatures and thermal inertias
352! ----------------------------------------------------------------------
353
354      write(*,*)
355      write(*,*)'Soil volumetric heat capacity:',volcapa
356
357      xmin = MINVAL(inertiedat)
358      xmax = MAXVAL(inertiedat)
359      write(*,*)'Soil thermal inertia <inertiedat>:',xmin,xmax
360
361      xmin = 1.0E+20
362      xmax = -1.0E+20
363      xmin = MINVAL(tsoil)
364      xmax = MAXVAL(tsoil)
365      write(*,*)'Soil temperature <tsoil>:',xmin,xmax
366      end
Note: See TracBrowser for help on using the repository browser.