Actual source code: ex37f90.F90

  1: #include "finclude/petscdef.h"
  2: !
  3: !   Notes:
  4: !     This uses Fortran 90 free-form, this means the lines can be up to 132 columns wide
  5: !       CHKR(ierr) is put on the same line as the call except when it violates the 132 columns
  6: !
  7: !     Eventually this file can be split into several files
  8: !
  9: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 10: !  Error handler that aborts when error is detected
 11: !
 12:       subroutine HE1(ierr,line)
 13:       use mex37f90
 14:       use mex37f90interfaces

 16:       call PetscError(ierr,line,0,'',ierr)
 17:       call MPI_Abort(PETSC_COMM_WORLD,ierr,ierr)
 18:       return
 19:       end
 20: #define CHKQ(n) if(n .ne. 0)call HE1(n,__LINE__)

 22: !  Error handler forces traceback of where error occurred
 23: !
 24:       subroutine HE2(ierr,line)
 25:       use mex37f90
 26:       use mex37f90interfaces

 28:       call PetscError(ierr,line,0,'',ierr)
 29:       return
 30:       end
 31: #define CHKR(n) if(n .ne. 0)then;call HE2(n,__LINE__);return;endif

 33: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 34: !   Utilities to map between global (linear algebra) to local (ghosted, physics) representations
 35: !
 36: !    Allocates the local work arrays and vectors
 37:       subroutine LocalFormCreate(dm,lf,ierr)
 38:       use mex37f90
 39:       use mex37f90interfaces
 40:       implicit none
 41:       DMComposite  dm
 42:       type(LocalForm) lf
 43:       PetscErrorCode ierr

 45:       call DMCompositeGetEntries(dm,lf%da,lf%np,lf%da,lf%np,ierr);CHKR(ierr) ! get the da's and sizes that define the unknowns
 46:       call DAGetLocalInfof90(lf%da,lf%dainfo,ierr);CHKR(ierr)
 47:       call DMCompositeGetLocalVectors(dm,lf%vIHX,lf%HotPool,lf%vCore,lf%ColdPool,ierr);
 48:       CHKR(ierr)
 49:       return
 50:       end

 52: !     Updates the (ghosted) IHX and Core arrays from the global vector x
 53:       subroutine LocalFormUpdate(dm,x,lf,ierr)
 54:       use mex37f90
 55:       use mex37f90interfaces
 56:       implicit none
 57:       DMComposite  dm
 58:       type(LocalForm) lf
 59:       PetscErrorCode ierr
 60:       Vec x

 62:       call DMCompositeScatter(dm,x,lf%vIHX,lf%HotPool,lf%vCore,lf%ColdPool,ierr);
 63:       CHKR(ierr)
 64:       call DAVecGetArrayf90(lf%da,lf%vIHX,lf%IHX,ierr);CHKR(ierr)
 65:       call DAVecGetArrayf90(lf%da,lf%vCore,lf%Core,ierr);CHKR(ierr)
 66:       return
 67:       end

 69: !     You should call this when you no longer need access to %IHX and %Core
 70:       subroutine LocalFormRestore(dm,lf,ierr)
 71:       use mex37f90
 72:       use mex37f90interfaces
 73:       implicit none
 74:       DMComposite  dm
 75:       type(LocalForm) lf
 76:       PetscErrorCode ierr

 78:       call DAVecRestoreArrayf90(lf%da,lf%vIHX,lf%IHX,ierr);CHKR(ierr)
 79:       call DAVecRestoreArrayf90(lf%da,lf%vCore,lf%Core,ierr);CHKR(ierr)
 80:       return
 81:       end

 83: !     Destroys the data structure
 84:       subroutine LocalFormDestroy(dm,lf,ierr)
 85:       use mex37f90
 86:       use mex37f90interfaces
 87:       implicit none
 88:       DMComposite  dm
 89:       type(LocalForm) lf
 90:       PetscErrorCode ierr

 92:       call DMCompositeRestoreLocalVectors(dm,lf%vIHX,lf%HotPool,lf%vCore,lf%ColdPool,ierr);
 93:       CHKR(ierr)
 94:       return
 95:       end

 97: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 98: !      Implements a nonlinear solver for a simple domain
 99: !     consisting of 2 zero dimensional problems linked by
100: !     2 one dimensional problems.
101: !
102: !                           Channel1
103: !                       -------------------------
104: !               Pool 1                              Pool 2
105: !                       -------------------------
106: !                           Channel2
107: !VAM
108: !VAM
109: !VAM
110: !VAM                         Hot Pool 1
111: !VAM                 |                       |
112: !VAM                 |                       |
113: !VAM                 |                       |
114: !VAM                 |                       |
115: !VAM  Core Channel 2 |                       | IHX Channel 1
116: !VAM                 |                       |
117: !VAM                 |                       |
118: !VAM                 |                       |
119: !VAM                 |                       |
120: !VAM                         Cold Pool 2
121: !VAM
122: !
123: !     For Newton's method with block Jacobi (one block for
124: !     each channel and one block for each pool) run with
125: !
126: !       -dmmg_iscoloring_type global -snes_mf_operator -pc_type lu
127: !       -help lists all options

129:       program ex37f90
130:       use mex37f90
131:  #include finclude/petscsys.h
132:  #include finclude/petscviewer.h
133:  #include finclude/petscvec.h

135:       DMMGArray        dmmg
136:       DMMG             dmmglevel
137:       PetscErrorCode   ierr
138:       DA               da
139:       DMComposite      dm
140:       type(appctx)     app
141:       external         FormInitialGuess,FormFunction,FormGraph
142:       Vec              x

144:       double precision time
145:       integer i
146: !BARRY
147:        PetscViewer        view0,view1
148: !BARRY

150:       call PetscInitialize(PETSC_NULL_CHARACTER,ierr)

152: !      Create the composite DM object that manages the grid

154:       call DMCompositeCreate(PETSC_COMM_WORLD,dm,ierr);CHKR(ierr)
155:       call DACreate1d(PETSC_COMM_WORLD,DA_NONPERIODIC,app%nxc,app%nc,1,PETSC_NULL_INTEGER,da,ierr)
156:       CHKR(ierr)
157:       call DMCompositeAddDM(dm,da,ierr);CHKR(ierr)
158:       call DADestroy(da,ierr);CHKR(ierr)
159:       call DMCompositeAddArray(dm,0,app%np,ierr);CHKR(ierr)
160:       call DMCompositeAddDM(dm,da,ierr);CHKR(ierr)
161:       call DMCompositeAddArray(dm,0,app%np,ierr);CHKR(ierr)

163: !       Create solver object and associate it with the unknowns (on the grid)

165:       call DMMGCreate(PETSC_COMM_WORLD,1,PETSC_NULL_OBJECT,dmmg,ierr);CHKR(ierr)
166:       call DMMGSetDM(dmmg,dm,ierr);CHKR(ierr)
167:       call DMCompositeDestroy(dm,ierr);CHKR(ierr)
168:       call DMMGSetUser(dmmg,0,app,ierr);CHKR(ierr)  ! currently only one level solver
169:       call DMMGSetInitialGuess(dmmg,FormInitialGuess,ierr);CHKR(ierr)
170:       call DMMGSetSNES(dmmg,FormFunction,PETSC_NULL_OBJECT,ierr);CHKR(ierr)
171:       call DMMGSetFromOptions(dmmg,ierr)
172: !BARRY
173:       call PetscViewerDrawOpen(PETSC_COMM_WORLD,PETSC_NULL_CHARACTER,'core',0,0,300,300,view0,ierr)
174:       CHKR(ierr)
175:       call PetscViewerDrawOpen(PETSC_COMM_WORLD,PETSC_NULL_CHARACTER,'ihx',320,0,300,300,view1,ierr)
176:       CHKR(ierr)
177: !BARRY

179:       call DMMGGetX(dmmg,x,ierr);CHKR(ierr)
180:       call VecDuplicate(x,app%xold,ierr);CHKR(ierr)
181:       call VecCopy(x,app%xold,ierr);CHKR(ierr)
182: !BARRY
183:        call DMMGArrayGetDMMG(dmmg,dmmglevel,ierr);CHKR(ierr)
184:        call FormGraph(dmmglevel,x,view0,view1,ierr);CHKR(ierr)
185: !BARRY

187:       time = 0.0d+0

189:       write(*,*)
190:       write(*,*)  'initial time'
191:       write(*,*)

193:       do i = 1,app%nstep

195:         time = time + app%dt

197:         write(*,*)
198:         write(*,*)  'time =',time
199:         write(*,*)
200: !
201: !  copy new to old
202: !
203: !
204: !   Solve the nonlinear system
205: !
206:         call DMMGSolve(dmmg,ierr);CHKR(ierr)

208:         call DMMGGetX(dmmg,x,ierr);CHKR(ierr)
209:         call VecCopy(x,app%xold,ierr);CHKR(ierr)
210: !BARRY
211:        call DMMGArrayGetDMMG(dmmg,dmmglevel,ierr);CHKR(ierr)
212:        call FormGraph(dmmglevel,x,view0,view1,ierr);CHKR(ierr)
213: !BARRY
214:       enddo
215: !BARRY
216:       call PetscViewerDestroy(view0,ierr);CHKR(ierr)
217:       call PetscViewerDestroy(view1,ierr);CHKR(ierr)
218: !BARRY
219:       write(*,*)
220:       write(*,*)  'final time'
221:       write(*,*)

223:       call VecDestroy(app%xold,ierr);CHKR(ierr)
224:       call DMMGDestroy(dmmg,ierr);CHKR(ierr)
225:       call PetscFinalize(ierr)
226:       end

228: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
229: !
230: !     The physics  :-)

232:       subroutine FormFunction(snes,x,f,dmmg,ierr)
233:       use mex37f90
234:       use mex37f90interfaces

236:       SNES snes
237:       Vec x,f
238:       DMMG dmmg
239:       PetscErrorCode ierr

241:       DMComposite dm
242:       Vec  fvc1,fvc2
243:       PetscInt np,i
244:       DA da
245:       type(poolfield), pointer :: fHotPool,fColdPool
246:       type(channelfield), pointer :: fIHX(:),fCore(:)
247:       type(appctx), pointer :: app
248:       PetscMPIInt rank
249:       type(LocalForm) xlf,xoldlf

251:       double precision dhhpl, mult, dhcpl, phpco, pcpio, pcpci, phpii

253:       logical debug

255:       debug = .false.

257:       call DMMGGetUser(dmmg,app,ierr);CHKR(ierr)   ! access user context

259:       call DMMGGetDM(dmmg,dm,ierr);CHKR(ierr)      ! access the grid information

261:       call LocalFormCreate(dm,xlf,ierr)            ! Access the local parts of x
262:       call LocalFormUpdate(dm,x,xlf,ierr)

264: !       Access the global (non-ghosted) parts of f
265:       call DMCompositeGetEntries(dm,da,np,da,np,ierr);CHKR(ierr) ! get the da's and sizes that define the unknowns
266:       call DMCompositeGetAccess(dm,f,fvc1,fHotPool,fvc2,fColdPool,ierr);CHKR(ierr)
267:       call DAVecGetArrayf90(da,fvc1,fIHX,ierr);CHKR(ierr)
268:       call DAVecGetArrayf90(da,fvc2,fCore,ierr);CHKR(ierr)

270: !BARRY
271: !
272: !  At this point I need to create old time values from "xold" passed in through
273: !  app for
274: !
275: !  xHotPool%vol, xHotPool%vel, xColdPool%vol, xColdPool%vel
276: !  xIHX(i)%press, xCore(i)%press
277: !
278: !  I don't know how to do this.
279: !
280: !BARRY
281:       call LocalFormCreate(dm,xoldlf,ierr)            ! Access the local parts of xold
282:       call LocalFormUpdate(dm,app%xold,xoldlf,ierr)

284:       call MPI_Comm_rank(app%comm,rank,ierr);CHKR(ierr)
285: !      fPool only lives on zeroth process, ghosted xPool lives on all processes
286:       if (rank == 0) then
287: !
288: !  compute velocity into hotpool from core
289: !
290:          dhhpl = app%dhhpl0 + ( (xlf%HotPool%vol - app%hpvol0) / app%ahp )
291:          phpco = app%P0 + ( app%rho * app%grav * (dhhpl - app%dhco) )
292:          mult = app%dt / (app%dxc * app%rho)
293:          fHotPool%vel = xlf%HotPool%vel - (mult * (xlf%Core(app%nxc-1)%press - phpco) ) + (abs(xlf%HotPool%vel) * xlf%HotPool%vel)
294: !
295: ! compute change in hot pool volume
296: !
297:          fHotPool%vol = xlf%HotPool%vol - app%hpvol0 - (app%dt * app%acore * (xlf%HotPool%vel -xlf%ColdPool%vel) )
298: !
299: !  compute velocity into coldpool from IHX
300: !
301:          dhcpl = app%dhcpl0 + ( (xlf%ColdPool%vol - app%cpvol0) / app%acp )
302:          pcpio = app%P0 + ( app%rho * app%grav * (dhcpl - app%dhio) )
303:          mult = app%dt / (app%dxc * app%rho)
304:          fColdPool%vel = xlf%ColdPool%vel-(mult*(xlf%IHX(app%nxc-1)%press-pcpio))+(abs(xlf%ColdPool%vel)*xlf%ColdPool%vel)
305: !
306: !  compute the change in cold pool volume
307: !
308:          fColdPool%vol = xlf%ColdPool%vol - app%cpvol0 - (app%dt * app%aihx * (xlf%ColdPool%vel - xlf%HotPool%vel) )
309:       endif
310: !
311: !  Compute the pressure distribution in IHX and core
312: !
313:       pcpci = app%P0 + ( app%rho * app%grav * (dhcpl - app%dhci) )
314:       phpii = app%P0 + ( app%rho * app%grav * (dhhpl - app%dhii) )

316:       do i=xlf%dainfo%xs,xlf%dainfo%xs+xlf%dainfo%xm-1

318:         fIHX(i)%press = xlf%IHX(i)%press - phpii - (app%rho * app%grav * dble(i) * app%dxi)
319:         fCore(i)%press = xlf%Core(i)%press - pcpci + (app%rho * app%grav * dble(i) * app%dxc)

321:       enddo

323:       if (debug) then
324:         write(*,*)
325:         write(*,*) 'Hot vel,vol ',xlf%HotPool%vel,xlf%HotPool%vol
326:         write(*,*) 'delta p = ', xlf%Core(app%nxc-1)%press - phpco,xlf%Core(app%nxc-1)%press,phpco
327:         write(*,*)

329:         do i=xlf%dainfo%xs,xlf%dainfo%xs+xlf%dainfo%xm-1
330:           write(*,*) 'xlf%IHX(',i,')%press = ',xlf%IHX(i)%press
331:         enddo

333:         write(*,*)
334:         write(*,*) 'Cold vel,vol ',xlf%ColdPool%vel,xlf%ColdPool%vol
335:         write(*,*) 'delta p = ',xlf%IHX(app%nxc-1)%press - pcpio,xlf%IHX(app%nxc-1)%press, pcpio
336:         write(*,*)


339:         do i=xlf%dainfo%xs,xlf%dainfo%xs+xlf%dainfo%xm-1
340:           write(*,*) 'xlf%Core(',i,')%press = ',xlf%Core(i)%press
341:         enddo

343:       endif

345:       call DAVecRestoreArrayf90(da,fvc1,fIHX,ierr);CHKR(ierr)
346:       call DAVecRestoreArrayf90(da,fvc2,fCore,ierr);CHKR(ierr)
347:       call DMCompositeRestoreAccess(dm,f,fvc1,fHotPool,fvc2,fColdPool,ierr);CHKR(ierr)

349:       call LocalFormRestore(dm,xoldlf,ierr)
350:       call LocalFormDestroy(dm,xoldlf,ierr)

352:       call LocalFormRestore(dm,xlf,ierr)
353:       call LocalFormDestroy(dm,xlf,ierr)
354:       return
355:       end

357:       subroutine FormGraph(dmmg,x,view0,view1,ierr)

359: ! --------------------------------------------------------------------------------------
360: !
361: !  FormGraph - Forms Graphics output
362: !
363: !  Input Parameter:
364: !  x - vector
365: !  time - scalar
366: !
367: !  Output Parameters:
368: !  ierr - error code
369: !
370: !  Notes:
371: !  This routine serves as a wrapper for the lower-level routine
372: !  "ApplicationXmgr", where the actual computations are
373: !  done using the standard Fortran style of treating the local
374: !  vector data as a multidimensional array over the local mesh.
375: !  This routine merely accesses the local vector data via
376: !  VecGetArray() and VecRestoreArray().
377: !
378:       use mex37f90
379:       use mex37f90interfaces

381:       Vec       x,xvc1,xvc2   !,corep,ihxp
382:       DMMG      dmmg
383:       PetscErrorCode ierr
384:       DMComposite dm
385:       PetscInt np            !,i
386:       DA da
387:       type(DALocalInfof90) dainfo
388:       type(poolfield), pointer :: xHotPool,xColdPool
389:       type(channelfield), pointer :: xIHX(:),xCore(:)
390:       type(appctx), pointer :: app

392:       PetscViewer        view0,view1

394:       integer iplotnum
395:       save iplotnum
396:       character*8 grfile

398:       data iplotnum / -1 /
399: !BARRY
400: !
401: !  This is my attempt to get the information out of vector "x" to plot
402: !  it.  I need to be able to build  xCore(i)%press and xIHX(i)%press
403: !  from the vector "x" and I do not know how to do this
404: !
405: !BARRY

407:       call DMMGGetUser(dmmg,app,ierr);CHKR(ierr)   ! access user context
408:       call DMMGGetDM(dmmg,dm,ierr);CHKR(ierr)      ! Access the grid information

410:       call DMCompositeGetEntries(dm,da,np,da,np,ierr);CHKR(ierr) ! get the da's and sizes that define the unknowns
411:       call DAGetLocalInfof90(da,dainfo,ierr);CHKR(ierr)

413:       call DMCompositeGetLocalVectors(dm,xvc1,xHotPool,xvc2,xColdPool,ierr);CHKR(ierr)
414:       call DMCompositeScatter(dm,x,xvc1,xHotPool,xvc2,xColdPool,ierr);CHKR(ierr)
415:       call DAVecGetArrayf90(da,xvc1,xIHX,ierr);CHKR(ierr)
416:       call DAVecGetArrayf90(da,xvc2,xCore,ierr);CHKR(ierr)

418: !
419: !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
420: !
421:       iplotnum = iplotnum + 1
422:       0
423: !
424: !
425: !  plot corep vector
426: !
427:       call VecView(xvc1,view0,ierr);CHKR(ierr)
428: !
429: !  make xmgr plot of corep
430: !
431:       write(grfile,4000) iplotnum
432:  4000 format('CoreP',i3.3)

434:       open (unit=44,file=grfile,status='unknown')

436:       do i = 1,app%nxc
437:         write(44,1000) dble(i)*app%dxc, xCore(i)%press
438:       enddo

440:       close(44)
441: !
442: !  plot ihxp vector
443: !
444:       call VecView(xvc2,view1,ierr);CHKR(ierr)
445: !
446: !  make xmgr plot of ihxp
447: !
448:       write(grfile,8000) iplotnum
449:  8000 format('IHXPr',i3.3)

451:       open (unit=44,file=grfile,status='unknown')

453:       do i = 1,app%nxc
454:         write(44,1000) dble(i)*app%dxi, xIHX(i)%press
455:       enddo

457:       close(44)



461:  1000 format(3(e18.12,2x))

463:       call DAVecRestoreArrayf90(da,xvc1,xIHX,ierr);CHKR(ierr)
464:       call DAVecRestoreArrayf90(da,xvc2,xCore,ierr);CHKR(ierr)
465:       call DMCompositeRestoreLocalVectors(dm,xvc1,xHotPool,xvc2,xColdPool,ierr);CHKR(ierr)

467:       return
468:       end

470: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

472:       subroutine FormInitialGuess(dmmg,v,ierr)
473:       use mex37f90
474:       use mex37f90interfaces

476:       DMMG dmmg
477:       Vec v
478:       PetscErrorCode ierr

480:       DMComposite dm
481:       Vec  vc1,vc2
482:       PetscInt np,i
483:       DA da
484:       type(DALocalInfof90) dainfo
485:       type(poolfield), pointer :: HotPool,ColdPool
486:       type(channelfield), pointer :: IHX(:),Core(:)
487:       type(appctx), pointer :: app
488:       PetscMPIInt rank

490:       double precision pcpci, phpii

492:       logical debug

494:       debug = .false.

496:       call DMMGGetUser(dmmg,app,ierr);CHKR(ierr)   ! access user context

498: !       Access the grid information

500:       call DMMGGetDM(dmmg,dm,ierr);CHKR(ierr)
501:       call DMCompositeGetEntries(dm,da,np,da,np,ierr);CHKR(ierr) ! get the da's and sizes that define the unknowns
502:       call DAGetLocalInfof90(da,dainfo,ierr);CHKR(ierr)

504: !      Acess the vector unknowns in global (nonghosted) form

506:       call DMCompositeGetAccess(dm,v,vc1,HotPool,vc2,ColdPool,ierr);CHKR(ierr)
507:       call DAVecGetArrayf90(da,vc1,IHX,ierr);CHKR(ierr)
508:       call DAVecGetArrayf90(da,vc2,Core,ierr);CHKR(ierr)

510:       call MPI_Comm_rank(app%comm,rank,ierr);CHKR(ierr)
511: !
512: !  Set the input values
513: !

515:        app%P0 = 1.0d+5
516:        app%rho = 1.0d+3
517:        app%grav = 9.8d+0

519:        app%dhhpl0 = 12.2d+0
520:        app%dhcpl0 = 10.16d+0

522:        app%lenc = 3.0d+0
523:        app%leni = 4.46d+0

525:        app%dhci = 0.83d+0
526:        app%dhco = app%dhci + app%lenc

528:        app%dhii = 7.83d+0
529:        app%dhio = app%dhii - app%leni

531:        app%dxc = app%lenc / dble(app%nxc)
532:        app%dxi = app%leni / dble(app%nxc)

534:        app%dt = 1.0d+0

536:        app%ahp = 7.0d+0
537:        app%acp = 7.0d+0

539:        app%acore = 0.8d+0
540:        app%aihx  = 5.0d-2

542:        app%hpvol0 = app%dhhpl0 * app%ahp
543:        app%cpvol0 = app%dhcpl0 * app%acp
544: !
545: !      the pools are only unknowns on process 0
546: !
547:       if (rank == 0) then
548:          HotPool%vel = -1.0d+0
549:          HotPool%vol = app%hpvol0
550:          ColdPool%vel = 1.0d+0
551:          ColdPool%vol = app%cpvol0
552:       endif
553: !
554: !  Construct and initial guess at the pressure
555: !
556:       pcpci = app%P0 + ( app%rho * app%grav * (app%dhcpl0 - app%dhci) )
557:       phpii = app%P0 + ( app%rho * app%grav * (app%dhhpl0 - app%dhii) )

559:       if (debug) then
560:         write(*,*)
561:         write(*,*) 'pcpci = ',pcpci
562:         write(*,*) 'phpii = ',phpii
563:         write(*,*) 'app%P0 = ',app%P0
564:         write(*,*) 'dhcpl0 - app%dhci ',app%dhcpl0 - app%dhci,app%dhcpl0, app%dhci
565:         write(*,*) 'dhhpl0 - app%dhii ',app%dhhpl0 - app%dhii,app%dhhpl0, app%dhii
566:         write(*,*)
567:       endif

569:       do i=dainfo%xs,dainfo%xs+dainfo%xm-1

571:         IHX(i)%press = phpii  + (app%rho * app%grav * dble(i) * app%dxi)

573:         Core(i)%press = pcpci - (app%rho * app%grav * dble(i) * app%dxc)

575:       enddo

577:       call DAVecRestoreArrayf90(da,vc1,IHX,ierr);CHKR(ierr)
578:       call DAVecRestoreArrayf90(da,vc2,Core,ierr);CHKR(ierr)
579:       call DMCompositeRestoreAccess(dm,v,vc1,HotPool,vc2,ColdPool,ierr);CHKR(ierr)

581:       CHKMEMQ
582:       return
583:       end