Changeset 1216 for trunk/LMDZ.GENERIC/libf/phystd/soil_settings.F
- Timestamp:
- Apr 3, 2014, 9:09:47 AM (11 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.GENERIC/libf/phystd/soil_settings.F
r787 r1216 1 subroutine soil_settings(nid,ngrid,nsoil,tsurf,tsoil) 2 3 use comsoil_h 4 1 subroutine soil_settings(nid,ngrid,nsoil,tsurf,tsoil,indextime) 2 3 ! use netcdf 4 use comsoil_h, only: layer, mlayer, inertiedat, volcapa 5 use iostart, only: inquire_field_ndims, get_var, get_field, 6 & inquire_field, inquire_dimension_length 5 7 implicit none 6 8 … … 9 11 ! 10 12 ! Purpose: Read and/or initialise soil depths and properties 13 ! 14 ! Modifications: Aug.2010 EM : use NetCDF90 to load variables (enables using 15 ! r4 or r8 restarts independently of having compiled 16 ! the GCM in r4 or r8) 17 ! June 2013 TN : Possibility to read files with a time axis 18 ! 11 19 ! 12 20 ! This subroutine reads from a NetCDF file (opened by the caller) … … 26 34 !====================================================================== 27 35 28 #include "dimensions.h"29 #include "dimphys.h"30 31 #include "netcdf.inc"32 36 !====================================================================== 33 37 ! arguments 34 38 ! --------- 35 39 ! inputs: 36 integer nid ! Input Netcdf file ID37 integer ngrid ! # of horizontal grid points38 integer nsoil ! # of soil layers39 integer sta ! # at which reading starts40 real tsurf(ngrid) ! surface temperature40 integer,intent(in) :: nid ! Input Netcdf file ID 41 integer,intent(in) :: ngrid ! # of horizontal grid points 42 integer,intent(in) :: nsoil ! # of soil layers 43 real,intent(in) :: tsurf(ngrid) ! surface temperature 44 integer,intent(in) :: indextime ! position on time axis 41 45 ! output: 42 real tsoil(ngrid,nsoilmx) ! soil temperature46 real,intent(out) :: tsoil(ngrid,nsoil) ! soil temperature 43 47 44 48 !====================================================================== … … 50 54 integer ndims ! # of dimensions of read <inertiedat> data 51 55 integer ig,iloop ! loop counters 56 57 integer edges(3),corner(3) ! to read a specific time 52 58 53 59 logical :: olddepthdef=.false. ! flag … … 70 76 real,parameter :: default_volcapa=1.e6 71 77 72 integer, dimension(2) :: sta2d 73 integer, dimension(2) :: siz2d 74 75 !====================================================================== 76 77 !! this is defined in comsoil_h 78 IF (.not.ALLOCATED(layer)) 79 . ALLOCATE(layer(nsoilmx)) 80 IF (.not.ALLOCATED(mlayer)) 81 . ALLOCATE(mlayer(0:nsoilmx-1)) 82 IF (.not.ALLOCATED(inertiedat)) 83 . ALLOCATE(inertiedat(ngrid,nsoilmx)) 84 85 !! this is defined in dimphys.h 86 sta = cursor 78 logical :: found,ok 79 80 !====================================================================== 87 81 88 82 ! 1. Depth coordinate … … 91 85 ! 1.1 Start by reading how many layers of soil there are 92 86 93 ierr=NF_INQ_DIMID(nid,"subsurface_layers",dimid) 94 if (ierr.ne.NF_NOERR) then 95 write(*,*)'soil_settings: Problem reading <subsurface_layers>' 96 call abort 97 endif 98 99 ierr=NF_INQ_DIMLEN(nid,dimid,dimlen) 100 if (ierr.ne.NF_NOERR) then 101 write(*,*)'soil_settings: Problem getting <subsurface_layers>', 102 & 'length' 103 call abort 104 endif 87 dimlen=inquire_dimension_length("subsurface_layers") 105 88 106 89 if (dimlen.ne.nsoil) then … … 112 95 endif 113 96 114 sta2d = (/sta,1/)115 siz2d = (/ngrid,dimlen/)116 117 118 97 ! 1.2 Find out the # of dimensions <inertiedat> was defined as using 119 ! (in ye old days, thermal inertia was only given at the "surface") 120 ! Look for thermal inertia data 121 ierr=NF_INQ_VARID(nid, "inertiedat", nvarid) 122 if (ierr.NE.NF_NOERR) then 123 write(*,*)'soil_settings: Field <inertiedat> not found!' 124 call abort 125 endif 126 127 ! Read the # of dimensions <inertidat> was defined as using 128 ierr=NF_INQ_VARNDIMS(nid,nvarid,ndims) 129 ! if (ndims.eq.1) then we have the "old 2D-surface" format 98 99 ndims=inquire_field_ndims("inertiedat") 130 100 131 101 ! 1.3 Read depths values or set olddepthdef flag and values … … 147 117 enddo 148 118 else ! Look for depth 149 ierr=NF_INQ_VARID(nid,"soildepth",nvarid)150 if (ierr.ne.NF_NOERR) then151 write(*,*)'soil_settings: Field <soildepth> not found!'152 call abort153 endif154 119 ! read <depth> coordinate 155 120 if (interpol) then !put values in oldmlayer 156 if (.not.allocated(oldmlayer)) then 157 allocate(oldmlayer(dimlen),stat=ierr) 158 endif 159 #ifdef NC_DOUBLE 160 ierr = NF_GET_VAR_DOUBLE(nid,nvarid,oldmlayer) 161 #else 162 ierr = NF_GET_VAR_REAL(nid,nvarid,oldmlayer) 163 #endif 164 if (ierr.ne.NF_NOERR) then 165 write(*,*)'soil_settings: Problem while reading <soildepth>' 166 call abort 121 call get_var("soildepth",oldmlayer,found) 122 if (.not.found) then 123 write(*,*)'soil_settings: Problem while reading <soildepth>' 167 124 endif 168 125 else ! put values in mlayer 169 #ifdef NC_DOUBLE 170 ierr = NF_GET_VAR_DOUBLE(nid,nvarid,mlayer) 171 #else 172 ierr = NF_GET_VAR_REAL(nid,nvarid,mlayer) 173 #endif 174 if (ierr.ne.NF_NOERR) then 175 write(*,*)'soil_settings: Problem while reading <soildepth>' 176 call abort 126 call get_var("soildepth",mlayer,found) 127 if (.not.found) then 128 write(*,*)'soil_settings: Problem while reading <soildepth>' 177 129 endif 178 130 endif !of if (interpol) … … 183 135 ! default mlayer distribution, following a power law: 184 136 ! mlayer(k)=lay1*alpha**(k-1/2) 185 lay1=2.e-4 !mars case (nsoilmx=18) 186 ! lay1=3.e-2 !earth case (nsoilmx=13) 137 lay1=2.e-4 187 138 alpha=2 188 139 do iloop=0,nsoil-1 … … 200 151 enddo 201 152 202 ! 2. Volumetric heat capacity (note: it is declared in comsoil .h)153 ! 2. Volumetric heat capacity (note: it is declared in comsoil_h) 203 154 ! --------------------------- 204 155 ! "volcapa" is (so far) 0D and written in "controle" table of startfi file … … 214 165 volcapa=default_volcapa 215 166 endif 216 ! Look for it 217 ! ierr=NF_INQ_VARID(nid,"volcapa",nvarid) 218 ! if (ierr.NE.NF_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.NF_NOERR) then 230 ! write(*,*)'soil_settings: Problem while reading <volcapa>' 231 ! call abort 232 ! endif 233 ! endif 234 235 ! 3. Thermal inertia (note: it is declared in comsoil.h) 167 168 ! 3. Thermal inertia (note: it is declared in comsoil_h) 236 169 ! ------------------ 237 170 238 171 ! 3.1 Look (again) for thermal inertia data (to reset nvarid) 239 ierr=NF_INQ_VARID(nid, "inertiedat", nvarid)240 if (ierr.NE.NF_NOERR) then241 write(*,*)'soil_settings: Field <inertiedat> not found!'242 call abort243 endif244 172 245 173 ! 3.2 Knowing the # of dimensions <inertidat> was defined as using, … … 251 179 ! Read Surface thermal inertia 252 180 allocate(surfinertia(ngrid)) 253 #ifdef NC_DOUBLE 254 ierr = NF_GET_VARA_DOUBLE(nid, nvarid, sta, ngrid, surfinertia) 255 #else 256 ierr = NF_GET_VARA_REAL(nid, nvarid, sta, ngrid, surfinertia) 257 #endif 258 if (ierr.NE.NF_NOERR) then 259 write(*,*)'soil_settings: Problem while reading <inertiedat>' 181 call get_field("inertiedat",surfinertia,found) 182 if (.not.found) then 183 write(*,*) "soil_settings: Failed loading <inertiedat>" 260 184 call abort 261 185 endif 262 186 263 187 write(*,*)' => Building soil thermal inertia (using reference sur … … 277 201 stop 278 202 endif 279 endif 280 #ifdef NC_DOUBLE 281 ierr = NF_GET_VARA_DOUBLE(nid,nvarid,sta2d,siz2d,oldinertiedat) 282 #else 283 ierr = NF_GET_VARA_REAL(nid,nvarid,sta2d,siz2d,oldinertiedat) 284 #endif 285 if (ierr.NE.NF_NOERR) then 286 write(*,*)'soil_settings: Problem while reading <inertiedat>' 287 call abort 203 endif ! of if (.not.allocated(oldinertiedat)) 204 call get_field("inertiedat",oldinertiedat,found) 205 if (.not.found) then 206 write(*,*) "soil_settings: Failed loading <inertiedat>" 207 call abort 288 208 endif 289 209 else ! put values in therm_i 290 #ifdef NC_DOUBLE 291 ierr = NF_GET_VARA_DOUBLE(nid,nvarid,sta2d,siz2d,inertiedat) 292 #else 293 ierr = NF_GET_VARA_REAL(nid,nvarid,sta2d,siz2d,inertiedat) 294 #endif 295 if (ierr.NE.NF_NOERR) then 296 write(*,*)'soil_settings: Problem while reading <inertiedat>' 297 call abort 298 endif 210 call get_field("inertiedat",inertiedat,found) 211 if (.not.found) then 212 write(*,*) "soil_settings: Failed loading <inertiedat>" 213 call abort 214 endif 215 ! endif 299 216 endif ! of if (interpol) 300 217 endif ! of if (ndims.eq.1) … … 303 220 ! ------------------------- 304 221 305 ierr=NF_INQ_VARID(nid,"tsoil",nvarid) 306 if (ierr.ne.NF_NOERR) then 222 ! ierr=nf90_inq_varid(nid,"tsoil",nvarid) 223 ok=inquire_field("tsoil") 224 ! if (ierr.ne.nf90_noerr) then 225 if (.not.ok) then 307 226 write(*,*)'soil_settings: Field <tsoil> not found!' 308 227 write(*,*)' => Building <tsoil> from surface values <tsurf>' … … 320 239 endif 321 240 endif 322 #ifdef NC_DOUBLE 323 ierr=NF_GET_VARA_DOUBLE(nid,nvarid,sta2d,siz2d,oldtsoil) 324 #else 325 ierr=NF_GET_VARA_REAL(nid,nvarid,sta2d,siz2d,oldtsoil) 326 #endif 327 if (ierr.ne.NF_NOERR) then 328 write(*,*)'soil_settings: Problem while reading <tsoil>' 329 call abort 330 endif 241 call get_field("tsoil",oldtsoil,found) 242 if (.not.found) then 243 write(*,*) "soil_settings: Failed loading <tsoil>" 244 call abort 245 endif 331 246 else ! put values in tsoil 332 #ifdef NC_DOUBLE 333 ierr=NF_GET_VARA_DOUBLE(nid,nvarid,sta2d,siz2d,tsoil) 334 #else 335 ierr=NF_GET_VARA_REAL(nid,nvarid,sta2d,siz2d,tsoil) 336 #endif 337 if (ierr.ne.NF_NOERR) then 338 write(*,*)'soil_settings: Problem while reading <tsoil>' 339 call abort 340 endif 247 call get_field("tsoil",tsoil,found,timeindex=indextime) 248 if (.not.found) then 249 write(*,*) "soil_settings: Failed loading <tsoil>" 250 call abort 251 endif 341 252 endif ! of if (interpol) 342 endif 253 endif! of if (.not.ok) 343 254 344 255 ! 5. If necessary, interpolate soil temperatures and thermal inertias
Note: See TracChangeset
for help on using the changeset viewer.