standalone/moisture_main.F

   1 !
   2 !*** moisture_main.F
   3 !
   4 ! This is the main driver for the standalone moisture model.  To compile,
   5 ! create a make.inc file in this directory from one of the example
   6 ! make.inc.* files.  Make sure to set NETCDF to the prefix path of your
   7 ! netcdf installation.  Then, type `make moisture` to create moisture.exe.
   8 !
   9 ! This program requires one or more wrfout files to run.  It relies on
  10 ! surface variables T2, Q2, and PSFC to be computed from the WRF run,
  11 ! which generally requires a nonzero sfclay option.  The command line
  12 ! syntax for this program is:
  13 !
  14 !   ./moisture.exe [-r] wrfout1 [wrfout2 [wrfout3 [...]]]
  15 !
  16 ! -r, when present indicates a restart run.  The program will open up a
  17 ! previous output file named moisture.nc and use the last time step
  18 ! present in that file to initialize a new run.
  19 !
  20 ! The wrfout files specified on the command line must be in sequential
  21 ! order, but they can contain overlapping time steps.  This is useful
  22 ! for outputs from operational reanalysis runs.  When two files contain
  23 ! data for the same time step, this program will always use the second
  24 ! file (with time step closer to the reanalysis).
  25 !
  26 ! For every time step computed by this program, it will output surface
  27 ! fields to moisture.nc, a netcdf file with similar conventions as a
  28 ! wrfout file.
  29
  30 module moisture_util
  31 use module_fr_sfire_util, only : crash
  32 use esmf_mod
  33 use netcdf
  34 implicit none
  35
  36 ! size parameters that depend on WRF
  37 integer,parameter::TIMESTRLEN=19,nfmc=5
  38
  39 ! if true, the program will only output moisture fields
  40 ! but the outputs cannot be used for restarting
  41 logical,parameter::smalloutput=.false.
  42
  43 ! model array allocation size for compatibility with WRF
  44 integer,save:: ims,ime,jms,jme
  45
  46 ! model array domain size
  47 integer,save::its,ite,jts,jte
  48
  49 ! initially set to false to indicate that the sizes above
  50 ! have not been initialized
  51 logical,private,save::initialized=.false.
  52
  53 ! a simple data structure to group together all variables used
  54 ! by the moisture model at a particular time step
  55 type ncvars
  56     ! model variables
  57     real,pointer,dimension(:,:,:)::fmc_gc,fmc_equi,fmc_tend
  58     real,pointer,dimension(:,:)::t2,q2,psfc,rainc,rainnc, &
  59                t2_old,q2_old,psfc_old,rain_old,rh_fire
  60     ! time stamp
  61     type(ESMF_Time)::time
  62
  63     ! rain accumulation variable
  64     real,pointer,dimension(:,:)::rain_accum,rain_zero
  65 end type
  66
  67 ! a data structure that contains metadata information about a
  68 ! wrfout file
  69 type ncfile
  70     ! number of time steps in the file
  71     integer::ntime
  72
  73     ! array of time stamps contained in the file
  74     type(ESMF_Time),dimension(:),pointer::times
  75
  76     ! the start time of the simulation
  77     type(ESMF_Time)::starttime
  78
  79     ! file handle
  80     integer::ncid
  81 end type
  82
  83 contains
  84
  85 subroutine initialize_and_check(filename,wrffile)
  86 ! open a file and initialize an ncfile type
  87 ! if module is not initialized, then initialize
  88 ! otherwise check that sizes are consistent
  89
  90 ! the name of the file to open
  91 character(len=*),intent(in)::filename
  92
  93 ! ncfile type initialized on output
  94 type(ncfile),intent(out)::wrffile
  95
  96 ! local variables dealing with netcdf i/o
  97 integer::dimidx,dimidy,dimidtime,dimidfuel
  98 integer::varidtime
  99 integer::nx,ny,nfuel,i,rc
 100 character(len=TIMESTRLEN),dimension(:),allocatable::timestr
 101 character(len=TIMESTRLEN)::stimestr
 102
 103 ! open the file and collect metadata
 104 call check(nf90_open(filename,nf90_nowrite,wrffile%ncid))
 105 call check(nf90_inq_dimid(wrffile%ncid,'west_east',  dimidx))
 106 call check(nf90_inq_dimid(wrffile%ncid,'south_north',dimidy))
 107 call check(nf90_inq_dimid(wrffile%ncid,'Time',dimidtime))
 108 call check(nf90_inquire_dimension(wrffile%ncid,dimidx,len=nx))
 109 call check(nf90_inquire_dimension(wrffile%ncid,dimidy,len=ny))
 110 call check(nf90_inquire_dimension(wrffile%ncid,dimidtime,len=wrffile%ntime))
 111
 112 if(.not.initialized)then
 113     ! the ESMF time modules need to be initialized on first call
 114     call ESMF_Initialize(defaultCalKind=ESMF_CAL_GREGORIAN)
 115
 116     ! on first call set module level variables containing array sizes
 117     ims=1
 118     ime=nx
 119     jms=1
 120     jme=ny
 121     its=ims
 122     ite=ime
 123     jts=jms
 124     jte=jme
 125     initialized=.true.
 126 else
 127     ! on later calls check that the file sizes are consistant with the sizes
 128     ! as stored
 129     if( (nx.ne.ime-ims+1) .or. &
 130         (ny.ne.jme-jms+1) ) then
 131         call crash('Incompatible file: '//trim(filename))
 132     endif
 133 endif
 134
 135 ! allocate memory for ncfile data type
 136 allocate(wrffile%times(wrffile%ntime))
 137 allocate(timestr(wrffile%ntime))
 138
 139 call check(nf90_inq_varid(wrffile%ncid,'Times',varidtime))
 140 call check(nf90_get_var(wrffile%ncid,varidtime,timestr))
 141
 142 ! go through all time stamps in the file and store them as ESMF types
 143 do i=1,wrffile%ntime
 144     call parse_wrf_time(timestr(i),wrffile%times(i))
 145 enddo
 146
 147 ! get start time
 148 call check(nf90_get_att(wrffile%ncid,NF90_GLOBAL,'START_DATE',stimestr))
 149 call parse_wrf_time(stimestr,wrffile%starttime)
 150
 151 ! free local memory
 152 deallocate(timestr)
 153
 154 end subroutine initialize_and_check
 155
 156 subroutine alldone()
 157 implicit none
 158 ! called at then end of execution to free allocated memory
 159 call ESMF_Finalize()
 160 end subroutine alldone
 161
 162 subroutine destroy_file(wrffile)
 163 implicit none
 164 ! free memory allocated in a ncfile data structure
 165 type(ncfile),intent(inout)::wrffile
 166 call check(nf90_close(wrffile%ncid))
 167 deallocate(wrffile%times)
 168 wrffile%ntime=0
 169 end subroutine destroy_file
 170
 171 subroutine initialize_vars(vars)
 172 implicit none
 173 ! Allocate memory for the model arrays
 174 type(ncvars),intent(inout)::vars
 175 if(.not.initialized)then
 176   call crash('Module not initialized')
 177 endif
 178 allocate(vars%fmc_gc(ims:ime,1:nfmc,jms:jme))
 179 allocate(vars%fmc_equi(ims:ime,1:nfmc,jms:jme))
 180 allocate(vars%fmc_tend(ims:ime,1:nfmc,jms:jme))
 181 allocate(vars%t2(ims:ime,jms:jme))
 182 allocate(vars%q2(ims:ime,jms:jme))
 183 allocate(vars%psfc(ims:ime,jms:jme))
 184 allocate(vars%rainc(ims:ime,jms:jme))
 185 allocate(vars%rainnc(ims:ime,jms:jme))
 186 allocate(vars%t2_old(ims:ime,jms:jme))
 187 allocate(vars%q2_old(ims:ime,jms:jme))
 188 allocate(vars%psfc_old(ims:ime,jms:jme))
 189 allocate(vars%rain_old(ims:ime,jms:jme))
 190 allocate(vars%rh_fire(ims:ime,jms:jme))
 191 allocate(vars%rain_accum(ims:ime,jms:jme))
 192 allocate(vars%rain_zero(ims:ime,jms:jme))
 193 vars%rain_accum(:,:)=0.
 194 vars%rain_zero(:,:)=0.
 195 end subroutine initialize_vars
 196
 197 subroutine destroy_vars(vars)
 198 implicit none
 199 ! free memory from model arrays
 200 type(ncvars),intent(inout)::vars
 201 deallocate(vars%fmc_gc)
 202 deallocate(vars%fmc_equi)
 203 deallocate(vars%fmc_tend)
 204 deallocate(vars%t2)
 205 deallocate(vars%q2)
 206 deallocate(vars%psfc)
 207 deallocate(vars%rainc)
 208 deallocate(vars%rainnc)
 209 deallocate(vars%t2_old)
 210 deallocate(vars%q2_old)
 211 deallocate(vars%psfc_old)
 212 deallocate(vars%rain_old)
 213 deallocate(vars%rh_fire)
 214 deallocate(vars%rain_accum)
 215 deallocate(vars%rain_zero)
 216 end subroutine destroy_vars
 217
 218 subroutine create_output(filename)
 219 implicit none
 220 ! create an output netcdf file with metadata initialized
 221 ! for the model being run.... does not write any data
 222 character(len=*),intent(in)::filename
 223
 224 integer::ncid,nxid,nyid,timeid,fuelid,dateid
 225 integer::tmp
 226
 227 ! create dataset
 228 call check(nf90_create(filename,nf90_clobber,ncid))
 229
 230 ! create dimensions
 231 call check(nf90_def_dim(ncid,'west_east',ime-ims+1,nxid))
 232 call check(nf90_def_dim(ncid,'south_north',jme-jms+1,nyid))
 233 call check(nf90_def_dim(ncid,'Time',nf90_unlimited,timeid))
 234 call check(nf90_def_dim(ncid,'fuel_moisture_classes_stag',nfmc,fuelid))
 235 call check(nf90_def_dim(ncid,'DateStrLen',TIMESTRLEN,dateid))
 236
 237 ! create variables
 238 call check(nf90_def_var(ncid,'Times',NF90_CHAR,(/dateid,timeid/),tmp))
 239 call check(nf90_def_var(ncid,'FMC_GC',NF90_REAL,(/nxid,nyid,fuelid,timeid/),tmp))
 240 if(.not.smalloutput)then
 241     call check(nf90_def_var(ncid,'FMC_EQUI',NF90_REAL,(/nxid,nyid,fuelid,timeid/),tmp))
 242     call check(nf90_def_var(ncid,'FMC_TEND',NF90_REAL,(/nxid,nyid,fuelid,timeid/),tmp))
 243     call check(nf90_def_var(ncid,'T2',NF90_REAL,(/nxid,nyid,timeid/),tmp))
 244     call check(nf90_def_var(ncid,'Q2',NF90_REAL,(/nxid,nyid,timeid/),tmp))
 245     call check(nf90_def_var(ncid,'PSFC',NF90_REAL,(/nxid,nyid,timeid/),tmp))
 246     call check(nf90_def_var(ncid,'RAIN',NF90_REAL,(/nxid,nyid,timeid/),tmp))
 247     call check(nf90_def_var(ncid,'RH_FIRE',NF90_REAL,(/nxid,nyid,timeid/),tmp))
 248 endif
 249 call check(nf90_close(ncid))
 250 end subroutine create_output
 251
 252 subroutine write_output(filename,vars)
 253 implicit none
 254 ! write arrays from the model into a preinitialized output file
 255 character(len=*),intent(in)::filename
 256 type(ncvars),intent(in)::vars
 257
 258 integer::ncid,varid,dimid,itime
 259 character(len=TIMESTRLEN)::timestr
 260 real,dimension(:,:,:),allocatable::tmp
 261
 262 ! allocate temporary memory for transposing 3D arrays
 263 allocate(tmp(its:ite,jts:jte,1:nfmc))
 264
 265 print*,'writing output to '//trim(filename)
 266 call check(nf90_open(filename,NF90_WRITE,ncid))
 267
 268 ! get the number of time frames already present and set the
 269 ! output time frame to the end of the file
 270 call check(nf90_inq_dimid(ncid,'Time',dimid))
 271 call check(nf90_inquire_dimension(ncid,dimid,len=itime))
 272 itime=itime+1
 273
 274 ! write time string
 275 call check(nf90_inq_varid(ncid,'Times',varid))
 276 call ESMF_TimeGet(vars%time,timeString=timestr)
 277 call check(nf90_put_var(ncid,varid,timestr, &
 278              start=(/1,itime/),count=(/TIMESTRLEN,1/)))
 279
 280 ! write FMC_GC
 281 call check(nf90_inq_varid(ncid,'FMC_GC',varid))
 282 call transpose_var(vars%fmc_gc,tmp)
 283 call check(nf90_put_var(ncid,varid,tmp, &
 284              start=(/its,jts,1,itime/),count=(/ite,jte,nfmc,1/)))
 285
 286 ! the following variables are diagnostic only, but useful for restarting
 287 if(.not.smalloutput)then
 288
 289     call check(nf90_inq_varid(ncid,'FMC_EQUI',varid))
 290     call transpose_var(vars%fmc_equi,tmp)
 291     call check(nf90_put_var(ncid,varid,tmp, &
 292                  start=(/its,jts,1,itime/),count=(/ite,jte,nfmc,1/)))
 293
 294     call check(nf90_inq_varid(ncid,'FMC_TEND',varid))
 295     call transpose_var(vars%fmc_tend,tmp)
 296     call check(nf90_put_var(ncid,varid,tmp, &
 297                  start=(/its,jts,1,itime/),count=(/ite,jte,nfmc,1/)))
 298
 299     call check(nf90_inq_varid(ncid,'T2',varid))
 300     call check(nf90_put_var(ncid,varid,vars%t2, &
 301                  start=(/its,jts,itime/),count=(/ite,jte,1/)))
 302
 303     call check(nf90_inq_varid(ncid,'Q2',varid))
 304     call check(nf90_put_var(ncid,varid,vars%q2, &
 305                  start=(/its,jts,itime/),count=(/ite,jte,1/)))
 306
 307     call check(nf90_inq_varid(ncid,'PSFC',varid))
 308     call check(nf90_put_var(ncid,varid,vars%psfc, &
 309                  start=(/its,jts,itime/),count=(/ite,jte,1/)))
 310
 311     call check(nf90_inq_varid(ncid,'RAIN',varid))
 312     call check(nf90_put_var(ncid,varid,vars%rain_accum, &
 313                  start=(/its,jts,itime/),count=(/ite,jte,1/)))
 314
 315     call check(nf90_inq_varid(ncid,'RH_FIRE',varid))
 316     call check(nf90_put_var(ncid,varid,vars%rh_fire, &
 317                  start=(/its,jts,itime/),count=(/ite,jte,1/)))
 318
 319 endif
 320
 321 ! free allocated memory
 322 call check(nf90_close(ncid))
 323 deallocate(tmp)
 324 end subroutine write_output
 325
 326 subroutine read_file(wrffile,atime,vars,ierr,restart)
 327 implicit none
 328 ! read array data for the model from a file
 329 type(ncfile),intent(in)::wrffile ! initialized ncfile object
 330 integer,intent(in)::atime        ! the time step in the file to read
 331 type(ncvars),intent(inout)::vars ! on return contains data from the file
 332 integer,intent(out)::ierr        ! error flag, 0 on success, otherwise read error
 333 logical,optional,intent(in)::restart ! if this a restart, we will read fmc_gc
 334
 335 integer::varid,i,j,k,itime
 336 logical::r
 337 real,dimension(:,:,:),allocatable::tmp
 338 real,dimension(ims:ime,jms:jme)::rainold
 339
 340 r=.false.
 341 if(present(restart))then
 342     r=restart
 343 endif
 344
 345 ! if the time step is unspecified (<= 0) then
 346 ! just read the last time step
 347 if(atime.gt.0)then
 348     itime=atime
 349 else
 350     itime=size(wrffile%times)
 351 endif
 352 vars%time=wrffile%times(itime)
 353
 354 ! read FMC_GC if it is present (generally, this test will
 355 ! only pass when we are reading for a restart)
 356 if(r.and.nf90_inq_varid(wrffile%ncid,'FMC_GC',varid).eq.0)then
 357     allocate(tmp(its:ite,jts:jte,nfmc))
 358     call check(nf90_get_var(wrffile%ncid,varid,tmp, &
 359             start=(/its,jts,1,itime/), count=(/ite,jte,nfmc,1/)))
 360     call transpose_var(tmp,vars%fmc_gc)
 361     deallocate(tmp)
 362 endif
 363
 364 ! read variables that should always be present
 365 call check(nf90_inq_varid(wrffile%ncid,'T2',varid))
 366 call check(nf90_get_var(wrffile%ncid,varid,vars%t2, &
 367         start=(/its,jts,itime/), count=(/ite,jte,1/)))
 368
 369 call check(nf90_inq_varid(wrffile%ncid,'Q2',varid))
 370 call check(nf90_get_var(wrffile%ncid,varid,vars%q2, &
 371         start=(/its,jts,itime/), count=(/ite,jte,1/)))
 372
 373 call check(nf90_inq_varid(wrffile%ncid,'PSFC',varid))
 374 call check(nf90_get_var(wrffile%ncid,varid,vars%psfc, &
 375         start=(/its,jts,itime/), count=(/ite,jte,1/)))
 376
 377 if(nf90_inq_varid(wrffile%ncid,'RAINC',varid).eq.0)then
 378     call check(nf90_get_var(wrffile%ncid,varid,vars%rainc, &
 379             start=(/its,jts,itime/), count=(/ite,jte,1/)))
 380
 381     if(itime.gt.1)then
 382         call check(nf90_get_var(wrffile%ncid,varid,rainold, &
 383                start=(/its,jts,itime-1/),count=(/ite,jte,1/)))
 384         vars%rainc(:,:)=vars%rainc(:,:)-rainold(:,:)
 385     endif
 386
 387     call check(nf90_inq_varid(wrffile%ncid,'RAINNC',varid))
 388     call check(nf90_get_var(wrffile%ncid,varid,vars%rainnc, &
 389             start=(/its,jts,itime/), count=(/ite,jte,1/)))
 390
 391     if(itime.gt.1)then
 392         call check(nf90_get_var(wrffile%ncid,varid,rainold, &
 393                start=(/its,jts,itime-1/),count=(/ite,jte,1/)))
 394         vars%rainnc(:,:)=vars%rainnc(:,:)-rainold(:,:)
 395     endif
 396 else
 397     call check(nf90_inq_varid(wrffile%ncid,'RAIN',varid))
 398     call check(nf90_get_var(wrffile%ncid,varid,vars%rainc, &
 399             start=(/its,jts,itime/), count=(/ite,jte,1/)))
 400     vars%rainnc(:,:)=0.
 401
 402     if(itime.gt.1)then
 403         call check(nf90_get_var(wrffile%ncid,varid,rainold, &
 404                start=(/its,jts,itime-1/),count=(/ite,jte,1/)))
 405         vars%rainc(:,:)=vars%rainc(:,:)-rainold(:,:)
 406     endif
 407 endif
 408
 409 ! calculate total rain accumulation for this moisture run
 410 vars%rain_accum(:,:)=vars%rain_accum(:,:) + vars%rainc(:,:) + vars%rainnc(:,:)
 411
 412 ! here we preempt any model errors by detecting if the data is correct
 413 ! if not, the main code will skip this time step and move on to the
 414 ! next
 415 if(any(vars%t2.le.0).or.any(vars%q2.lt.0).or.any(vars%psfc.le.0))then
 416     ierr=1
 417 else
 418     ierr=0
 419 endif
 420 end subroutine read_file
 421
 422 subroutine transpose_var(A,B)
 423 implicit none
 424 ! do a transpose on a 3D array to account for different
 425 ! storage conventions between wrf output and runtime
 426 ! B(i,k,j)=A(i,j,k)
 427 real,dimension(:,:,:),intent(in)::A
 428 real,dimension(:,:,:),intent(out)::B
 429
 430 integer::i,j,k
 431 if( size(A,1).ne.size(B,1) .or. &
 432     size(A,2).ne.size(B,3) .or. &
 433     size(A,3).ne.size(B,2) )then
 434   call crash('Invalid array sizes in transpose')
 435 endif
 436 do j=1,size(A,2)
 437     do k=1,size(A,3)
 438         do i=1,size(A,1)
 439             B(i,k,j)=A(i,j,k)
 440         enddo
 441     enddo
 442 enddo
 443 end subroutine transpose_var
 444
 445 subroutine parse_wrf_time(wrfstr,time)
 446 ! parse a wrf time string into an esmf time class
 447 ! example wrf time string: "2011-01-01_00:00:00"
 448 implicit none
 449 character(len=TIMESTRLEN),intent(in)::wrfstr
 450 type(ESMF_Time),intent(out)::time
 451
 452 integer::year,month,day,hour,minute,second
 453 character(len=1)::a1,a2,a3,a4,a5
 454
 455 read(wrfstr,'(I04,A1,I02,A1,I02,A1,I02,A1,I02,A1,I02)') &
 456        year,a1,month,a2,day,a3,hour,a4,minute,a5,second
 457 call ESMF_TimeSet(time,YY=year,MM=month,DD=day,H=hour,M=minute,S=second)
 458 end subroutine parse_wrf_time
 459
 460 subroutine get_next_timestep(files,ifile,istep)
 461 ! choose a file and time step within the file to use
 462 ! for the next call to the moisture model
 463 !
 464 ! This assumes files are in sequential order, and
 465 ! allows for reanalysis runs where there may be
 466 ! time step overlaps from one file to the next.
 467 ! Here, we prefer to move on to the next file
 468 ! rather than the next time step in the old file
 469 ! to use model arrays closer to reanalysis.
 470 ! If there is no overlap, this routine will
 471 ! increment from one file to the next as expected.
 472 !
 473 ! This subroutine will return istep=ifile=0, when
 474 ! there are no more steps to process.
 475 implicit none
 476 type(ncfile),dimension(:),intent(in)::files
 477 integer,intent(inout)::ifile,istep
 478
 479 integer::nfile,nstep
 480 type(ESMF_Time)::time1,time2,time3
 481 nfile=size(files)
 482 if(ifile.lt.1)then
 483     ! for initialization only
 484     ifile=1
 485     istep=1
 486     return
 487 endif
 488 nstep=size(files(ifile)%times)
 489 time1=files(ifile)%times(istep)
 490
 491 if(istep+1.le.size(files(ifile)%times))then
 492     time2=files(ifile)%times(istep+1)
 493 else
 494     ifile=ifile+1
 495     istep=1
 496     goto 990
 497 endif
 498
 499 if(ifile+1.le.nfile)then
 500     ! if there is another file after the current one
 501     time3=files(ifile+1)%times(1)
 502 else
 503     ! if there is not another file pick next time step in current file
 504     istep=istep+1
 505     goto 990
 506 endif
 507
 508 if(time2.ge.time3)then
 509     ! if next file has a time step closer to the current time
 510     ! or if there are no more time steps in the current file,
 511     ! use the next file
 512     ifile=ifile+1
 513     istep=1
 514 else
 515     ! otherwise use the next step in the current file
 516     istep=istep+1
 517 endif
 518
 519  990 continue
 520 if(ifile.gt.nfile.or.istep.gt.size(files(ifile)%times))then
 521     ! make sure the current time step is valid
 522     ! end if not
 523     istep=0
 524     ifile=0
 525 endif
 526 end subroutine get_next_timestep
 527
 528 subroutine check(status)
 529 integer,intent(in)::status
 530 ! check a netcdf status code for error
 531 if(status.ne.nf90_noerr) then
 532     call crash(trim(nf90_strerror(status)))
 533 endif
 534 end subroutine check
 535
 536 end module moisture_util
 537
 538 program moisture_main
 539 use moisture_util
 540 use module_fr_sfire_phys, only: advance_moisture, init_fuel_cats, moisture_classes
 541 use module_fr_sfire_util, only: fire_print_msg
 542 implicit none
 543
 544 logical::initialize
 545 integer::fmoist_init=2
 546 real::moisture_dt,dt
 547 integer::istep,i,iglobal
 548 logical::restart
 549
 550 integer :: numarg,ifile,ii,nfile,s,sn,sd,ierr
 551 character(len=120) :: wrfoutfname
 552 character(len=*),parameter :: outputfname='moisture.nc'
 553 character(len=64)::timediffstr
 554
 555 type(ncfile),dimension(:),allocatable::wrffile
 556 type(ncfile)::wrfoutfile
 557 type(ncvars)::vars
 558 type(ESMF_Time)::oldtime,timenow,starttime
 559 type(ESMF_TimeInterval)::timedt
 560
 561 integer,parameter::iprint=100,jprint=100
 562 integer::k
 563 real::hours,t,q,p,rain
 564 real,dimension(5)::fmc_equi,fmc_tend,fmc_gc
 565 logical,parameter::debugprint=.false.
 566 integer,parameter::debugfile=8
 567
 568 #ifndef GFORTRAN_IARGC
 569 integer,external::iargc
 570 #endif
 571
 572 1 format(a4,a8,  a6,  a9,  a6  ,a8, a6, 5(a5,i1,a8,i1,a5,i1))
 573 2 format(a4,a8,  a6,  a9,  a6  ,a8, a6, 5(a6,a9,a6))
 574 3 format(i4,f8.2,f6.1,f9.1,f6.3,f8.1,f6.3,5(f6.3,f9.1,f6.3))
 575
 576 if(debugprint)then
 577     open(debugfile,file='moisture_output.txt',form='formatted',status='unknown')
 578     write(debugfile,1)'Step','Time', 'T','P','Q','RAIN','RH',('EQUI',k,'TLAG',k,'FMC',k,k=1,moisture_classes)
 579     write(debugfile,2)' ','hours','K','Pa','kg/kg','mm','1',('kg/kg','hours','kg/kg',k=1,moisture_classes)
 580 endif
 581 hours=0
 582
 583 ! parse commandline argument for wrfout file name
 584 numarg=iargc()
 585
 586 if(numarg .lt. 1)then
 587     print*, 'usage: ./moisture.exe [-r] wrfout1 [wrfout2 [wrfout3 [...]]]'
 588     call abort()
 589 endif
 590
 591 ! check for restart flag
 592 call getarg(1,wrfoutfname)
 593
 594 ! allocate local storage for wrfout metadata
 595 if(trim(wrfoutfname).eq.'-r')then
 596     allocate(wrffile(numarg-1))
 597 else
 598     allocate(wrffile(numarg))
 599 endif
 600
 601 ! open files and populate ncfile objects for each
 602 ! check that the sizes are consistent
 603 nfile=0
 604 restart=.false.
 605 do i=1,numarg
 606     call getarg(i,wrfoutfname)
 607     if(trim(wrfoutfname).eq.'-r')then
 608         restart=.true.
 609         call initialize_and_check(outputfname,wrfoutfile)
 610     else
 611         print*,'using: '//trim(wrfoutfname)
 612         nfile=nfile+1
 613         call initialize_and_check(wrfoutfname,wrffile(nfile))
 614     endif
 615 enddo
 616
 617 ! allocate memore for model arrays
 618 call initialize_vars(vars)
 619
 620 oldtime=wrffile(1)%times(1)
 621 timenow=oldtime
 622
 623 ! if this isn't a restart run, we need to create
 624 ! a new output file, clobber anything that might
 625 ! be there, otherwise we append data to the restart
 626 if(.not.restart) &
 627 call create_output(outputfname)
 628
 629 ! model initialization calls
 630 fire_print_msg = 2
 631 call init_fuel_cats(.true.)
 632 print *,moisture_classes,' moisture classes'
 633
 634 ! initialize various looping indices
 635 ifile=0
 636 ii=0
 637 istep=0
 638
 639 ! get the first file/time step we will read
 640 call get_next_timestep(wrffile,ifile,ii)
 641 if(restart)then
 642     starttime=wrfoutfile%times(size(wrfoutfile%times))
 643 else
 644     starttime=wrffile(1)%times(1)
 645 endif
 646
 647 ! for restarts, we have to find the first time step in the
 648 ! wrfout's after the the time in the restart
 649 do while(restart.and.wrffile(ifile)%times(ii).le.starttime)
 650     call get_next_timestep(wrffile,ifile,ii)
 651     if(ifile.eq.0)goto 1000
 652 enddo
 653  1000 continue
 654
 655 ! after all the initialization code, we crash here if there are no
 656 ! new time steps to read
 657 if(ifile.eq.0)then
 658     call crash('no new time steps available')
 659 endif
 660
 661 ! let's tell the user what we are doing
 662 print*,'starting at file ',ifile,' step ',ii
 663 call ESMF_TimeGet(wrffile(ifile)%times(ii),TimeString=timediffstr)
 664 print*,'wrffile time=',trim(timediffstr)
 665 call ESMF_TimeGet(starttime,TimeString=timediffstr)
 666 print*,'(re)start time=',trim(timediffstr)
 667
 668 ! get the next target time step
 669 call get_next_timestep(wrffile,ifile,ii)
 670
 671 iglobal=0
 672
 673 ! main program loop, go until get_next_timestep returns ifile==0
 674 do while(ifile.gt.0)
 675
 676     ! keep track of the number of loops we have made
 677     istep=istep+1
 678     iglobal=iglobal+1
 679
 680     ! read initialization data
 681     if(restart)then
 682         call read_file(wrfoutfile,-1,vars,ierr,restart=.true.)
 683         if(ierr.ne.0)then
 684             call crash('invalid restart file')
 685         endif
 686         fmoist_init=3
 687     else
 688         call read_file(wrffile(ifile),ii,vars,ierr)
 689         if(ierr.ne.0)then
 690             print*,'invalid data: skipping file ',ifile,' step ',ii
 691             goto 10001
 692         endif
 693     endif
 694
 695     ! increment time keeping
 696     oldtime=timenow
 697     timenow=wrffile(ifile)%times(ii)
 698
 699     ! let's tell the user what we are doing
 700     print*,'reading file ',ifile,' of ',nfile
 701     print*,'time step ',ii,' of ',wrffile(ifile)%ntime
 702     call ESMF_TimeGet(timenow,timeString=timediffstr)
 703     print*,'time: ',trim(timediffstr)
 704
 705     if(iglobal.eq.1)then
 706         dt=0.  ! required on initialization
 707         initialize=.true.
 708     else
 709         if(timenow.lt.oldtime)then
 710             call crash('Negative time step detected.')
 711         endif
 712         ! get the number of seconds that the model needs to advance
 713         timedt=timenow-oldtime
 714         call ESMF_TimeIntervalGet(timedt,S=s,sn=sn,sd=sd,TimeString=timediffstr)
 715         if(sd.gt.0)then
 716             dt=dble(s)+dble(sn)/dble(sd)
 717         else
 718             dt=dble(s)
 719         endif
 720         initialize=.false.
 721     endif
 722
 723
 724     call advance_moisture(    &
 725         initialize,                 & ! initialize timestepping. true on the first call at time 0, then false
 726         ims,ime,  jms,jme,          & ! memory dimensions
 727         its,ite,  jts,jte,          & ! tile dimensions
 728         nfmc,                       & ! dimension of moisture fields
 729         dt,                & ! timestep = time step time elapsed from the last call
 730         vars%rain_accum, vars%rain_zero,              & ! accumulated rain
 731         vars%t2, vars%q2, vars%psfc,               & ! temperature (K), vapor contents (kg/kg), pressure (Pa) at the surface
 732         vars%rain_old,                   & ! previous value of accumulated rain
 733         vars%t2_old, vars%q2_old, vars%psfc_old,   & ! previous values of the atmospheric state at surface
 734         vars%rh_fire,                    & ! relative humidity diagnostics
 735         vars%fmc_gc,                     & ! fuel moisture by class, updated
 736         vars%fmc_equi,                   & ! fuel moisture equilibrium by class, for diagnostics only
 737         vars%fmc_tend                    & ! fuel moisture tendency by classe, for diagnostics only
 738         )
 739
 740     if(debugprint)then
 741         hours=hours+dt/3600.
 742         t=vars%t2(iprint,jprint)
 743         q=vars%q2(iprint,jprint)
 744         p=vars%psfc(iprint,jprint)
 745         rain=vars%rain_accum(iprint,jprint) + vars%rain_zero(iprint,jprint)
 746         fmc_equi(:)=vars%fmc_equi(iprint,:,jprint)
 747         fmc_tend(:)=vars%fmc_tend(iprint,:,jprint)
 748         fmc_gc(:)=vars%fmc_gc(iprint,:,jprint)
 749         write(debugfile,3)iglobal,hours,t,p,q,rain, &
 750                   vars%rh_fire(iprint,jprint),   &
 751                (fmc_equi(k),fmc_tend(k),fmc_gc(k), &
 752                 k=1,moisture_classes)
 753     endif
 754     ! write the results to the output file
 755     call write_output(outputfname,vars)
 756
 757     10001 continue ! we skip to here on invalid data
 758
 759     ! after the first pass through the loop, nothing special needs to be done
 760     ! for restart runs, so turn of the restart flag here
 761     if(restart)then
 762         restart=.false.
 763     else
 764         call get_next_timestep(wrffile,ifile,ii)
 765     endif
 766 enddo ! end main loop
 767
 768 ! free up allocated memory
 769 call destroy_vars(vars)
 770 do i=1,nfile
 771     call destroy_file(wrffile(i))
 772 enddo
 773 call alldone()
 774 end program moisture_main