source: trunk/LMDZ.MARS/libf/phymars/soil_settings.F @ 3109

Last change on this file since 3109 was 3109, checked in by llange, 13 months ago

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

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