modopenboundary: fix initialization

boundary(...)%uphase is used in openboundary_tend but
was set in subroutine openboundary_phasevelocity,
called in the time loop *after* openboundary_tend.

For now fixed by calling openboundary_phasevelocity near the end of
startup().  Note the placement is tricky, because it has to be after
rhobf has been initialized.

Also moved the allocation and initialization of rhointi from initopenboundary
to startup, so that it's done before the first openboundary_phasevelocity
but after rhobf.

src/modopenboundary.f90

@@ -43,8 +43,9 @@ module modopenboundary
   subroutine initopenboundary
     ! Initialisation routine for openboundaries
     use modmpi, only : myidx, myidy, nprocx, nprocy, myid
-    use modglobal, only : imax,jmax,kmax,i1,j1,k1,dx,dy,itot,jtot,solver_id,nsv,cu,cv
+    use modglobal, only : imax,jmax,kmax,i1,j1,k1,dx,dy,itot,jtot,solver_id,nsv,cu,cv,dzf
     use modboundary, only: dsv
+    use modfields, only: rhobf
     implicit none
     integer :: i
 
@@ -168,7 +169,7 @@ subroutine openboundary_initfields
     ! Variables not present in the input file are initialised to their profiles
     use modglobal, only : imax,jmax,kmax,i1,j1,cexpnr,nsv,i2,j2,k1
     use modfields, only : u0,um,v0,vm,w0,wm,thl0,thlm,qt0,qtm,e120,e12m, sv0, svm, &
-      uprof,vprof,thlprof,qtprof,e12prof,svprof
+      uprof,vprof,thlprof,qtprof,e12prof,svprof,up,vp,wp
     use modmpi, only : myidx,myidy,myid
     implicit none
     integer :: VARID,STATUS,NCID,n
@@ -177,6 +178,7 @@ subroutine openboundary_initfields
     if(.not.linithetero) return
     u0 = 0.; um = 0.; v0 = 0.; vm = 0.; w0 = 0.; wm = 0.;
     thl0 = 0.; thlm = 0.; qt0 = 0; qtm = 0; e120 = 0.; e12m = 0.
+    up = 0.; vp = 0.; wp = 0. ! do we need this?
     !--- open initfields.input.xxx.nc ---
     STATUS = NF90_OPEN('initfields.inp.'//cexpnr//'.nc', nf90_nowrite, NCID)
     if (STATUS .ne. nf90_noerr) call handle_err(STATUS)
@@ -444,9 +446,6 @@ subroutine openboundary_divcorr()
     implicit none
     real(field_r) :: sumdiv,div,divpart,divnew,divold
     integer :: i,j,k,it,icalc
-    ! Create 1/int(rho)
-    allocate(rhointi(k1))
-    rhointi = 1./(rhobf*dzf)
     ! Divergence correction
     if(myid==0) print *, "Start divergence correction boundaries"
     do it = 1,ntboundary

src/modstartup.f90

@@ -77,15 +77,15 @@ subroutine startup(path)
                                   solver_id, maxiter, maxiter_precond, tolerance, n_pre, n_post, precond_id, checknamelisterror, &
                                   loutdirs, output_prefix, &
                                   lopenbc,linithetero,lperiodic,dxint,dyint,dzint,dxturb,dyturb,taum,tauh,pbc,lsynturb,nmodes,tau,lambda,lambdas,lambdas_x,lambdas_y,lambdas_z,iturb, &
-                                  hypre_logging,rdt,rk3step,i1,j1,k1,ih,jh,lboundary,lconstexner, lstart_netcdf
+                                  hypre_logging,rdt,rk3step,i1,j1,k1,ih,jh,lboundary,lconstexner, lstart_netcdf,dzf
     use modforces,         only : lforce_user
     use modsurfdata,       only : z0,ustin,wtsurf,wqsurf,wsvsurf,ps,thls,isurf
     use modsurface,        only : initsurface
     use moddatetime,       only : initdatetime
     use modemission,       only : initemission
     use modlsm,            only : initlsm, kmax_soil
     use moddrydeposition,  only : initdrydep
-    use modfields,         only : initfields,um,vm,wm,u0,v0,w0,up,vp,wp
+    use modfields,         only : initfields,um,vm,wm,u0,v0,w0,up,vp,wp,rhobf
     use modtracers,        only : inittracers, allocate_tracers
     use modpois,           only : initpois,poisson
     use modradiation,      only : initradiation
@@ -104,7 +104,9 @@ subroutine startup(path)
     use tstep,             only : inittstep
     use modchem,           only : initchem
     use modversion,        only : git_version
-    use modopenboundary,   only : initopenboundary,openboundary_divcorr,openboundary_excjs,lbuoytop
+    use modopenboundary,   only : initopenboundary,openboundary_divcorr,openboundary_excjs,lbuoytop,&
+                                  rhointi, openboundary_phasevelocity
+
     use modchecksim,       only : chkdiv
 #if defined(_OPENACC)
     use modgpu,             only : initgpu
@@ -379,6 +381,12 @@ subroutine startup(path)
     call inittimedep !depends on modglobal,modfields, modmpi, modsurf, modradiation
     call initpois ! hypre solver needs grid and baseprofiles
     if(lopenbc) then  ! Correct boundaries and initial field for divergence
+      ! Create 1/int(rho) - must be after rhobf has been initialized
+      allocate(rhointi(k1))
+      rhointi = 1./(rhobf*dzf)
+
+      call openboundary_phasevelocity() ! needed for initialization, called late in the time loop
+
       call chkdiv
       call openboundary_divcorr ! Remove divergence from large scale input
       ! Use poisson solver to get rid of divergence in initial field, needs to
@@ -596,9 +604,9 @@ subroutine readinitfiles
 
         else if (lstart_netcdf) then
           call init_from_netcdf('init.'//cexpnr//'.nc', height, uprof, vprof, &
-                                thlprof, qtprof, e12prof, ug, vg, wfls, & 
+                                thlprof, qtprof, e12prof, ug, vg, wfls, &
                                 dqtdxls, dqtdyls, dqtdtls, thlpcar, kmax)
-          call tracer_profs_from_netcdf('tracers.'//cexpnr//'.nc', & 
+          call tracer_profs_from_netcdf('tracers.'//cexpnr//'.nc', &
                                         tracer_prop, nsv, svprof(1:kmax,:))
         else
           open (ifinput,file='prof.inp.'//cexpnr,status='old',iostat=ierr)
@@ -687,7 +695,7 @@ subroutine readinitfiles
             endif
           enddo
           ! write(*,*) 'scalar_indices: ', scalar_indices
-          
+
           close(ifinput)
 
           write (6,*) 'height   sv(1) --------- sv(nsv) '
@@ -841,10 +849,10 @@ subroutine readinitfiles
       else
         call boundary
       end if
-      
+
       call thermodynamics
       call surface
-      
+
       if ( lopenbc ) then
         call openboundary_ghost()
       else
@@ -1260,7 +1268,7 @@ subroutine writerestartfiles
     if ((timee>=tnextrestart .and. trestart > 0) .or. (timeleft==0 .and. trestart >= 0)) then
       tnextrestart = tnextrestart+itrestart
 #if defined(_OPENACC)
-      call update_host 
+      call update_host
 #endif
       call do_writerestartfiles
     end if
@@ -1775,12 +1783,12 @@ end subroutine baseprofs
   !! \param dqtdtls Tendency of the total water mixing ratio.
   !! \param dthlrad Tendency of the liquid water potential temperature due to radiative heating.
   !! \param kmax Index of highest vertical level.
-  !! 
+  !!
   !! \note Tracers are read from tracers.XXX.nc, not here.
   !! \todo Make DEPHY-compatible.
   subroutine init_from_netcdf(filename, height, uprof, vprof, thlprof, qtprof, &
                               e12prof, ug, vg, wfls, dqtdxls, dqtdyls, &
-                              dqtdtls, dthlrad, kmax) 
+                              dqtdtls, dthlrad, kmax)
     character(*),   intent(in)  :: filename
     real(field_r),  intent(out) :: height(:)
     real(field_r),  intent(out) :: uprof(:)
@@ -1831,7 +1839,7 @@ subroutine init_from_netcdf(filename, height, uprof, vprof, thlprof, qtprof, &
                        fillvalue=0._field_r)
 
     call nchandle_error(nf90_close(ncid))
-    
+
   end subroutine init_from_netcdf
 
 end module modstartup