Fix uninitialized variables

src/modlsm.f90

@@ -306,7 +306,7 @@ subroutine calc_liquid_reservoir
             end do
 
             ! Tendency due to interception of precipitation by vegetation
-            if (imicro == 0) then
+            if (imicro == 0 .or. imicro == 1) then
                 rainrate = 0.
             else
                !rainrate = -sed_qr(i,j,1)/rhow
@@ -1796,8 +1796,34 @@ subroutine init_heterogeneous_nc
     tile(nlu)%lveg    = .false.
     tile(nlu)%laqu    = .false.
 
+    ! initialize to 0 to avoid un-initialized values in the halo or elsewhere
+    tile(nlu)%base_frac = 0
+    tile(nlu)%z0m = 0
+    tile(nlu)%z0h = 0
+    tile(nlu)%lambda_stable = 0
+    tile(nlu)%lambda_unstable = 0
+    tile(nlu)%rs_min = 0
+    tile(nlu)%lai = 0
+    tile(nlu)%a_r = 0
+    tile(nlu)%b_r = 0
+    tile(nlu)%gD = 0
+    tile(nlu)%tskin = 0
+
     ! 2D surface fields
     do ilu=1,nlu-1
+       ! initialize to 0 to avoid un-initialized values in the halo
+      tile(ilu)%base_frac = 0
+      tile(ilu)%z0m = 0
+      tile(ilu)%z0h = 0
+      tile(ilu)%lambda_stable = 0
+      tile(ilu)%lambda_unstable = 0
+      tile(ilu)%rs_min = 0
+      tile(ilu)%lai = 0
+      tile(ilu)%a_r = 0
+      tile(ilu)%b_r = 0
+      tile(ilu)%gD = 0
+      tile(ilu)%tskin = 0
+
       write(*,*) 'reading variables for LU type: ', trim(tile(ilu)%lushort)
       ! LU cover
       call check( nf90_inq_varid( ncid, 'cover_'//trim(tile(ilu)%lushort), varid) )
@@ -1983,6 +2009,7 @@ subroutine init_heterogeneous_nc
     end do
 
     ! Calculate vegetation fraction, and limit to prevent div/0's
+    cveg = 0
     do ilu=1,nlu
       if (tile(ilu)%lveg) then
         cveg(:,:) = cveg(:,:) + tile(ilu)%base_frac(:,:)
@@ -1991,6 +2018,7 @@ subroutine init_heterogeneous_nc
     where (cveg == 0) cveg = eps1
 
     ! Calculate land fraction, and limit to prevent div/0's
+    land_frac = 0
     do ilu=1,nlu
       if (tile(ilu)%laqu) then
         land_frac(:,:) = 1.-tile(ilu)%base_frac(:,:)

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
 
@@ -444,9 +445,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

src/modtimestat.f90

@@ -591,7 +591,7 @@ subroutine timestat
       end do
     end if
 
-    if (imicro > 0) then
+    if (imicro > 1) then
       !$acc parallel loop collapse(2) default(present) reduction(+:qrintavl) &
       !$acc& private(qrint) async
       do j = 2, j1