diff --git a/src/modopenboundary.f90 b/src/modopenboundary.f90 index 8cf6926f..a42b7fcc 100644 --- a/src/modopenboundary.f90 +++ b/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 diff --git a/src/modstartup.f90 b/src/modstartup.f90 index c44f8cc4..44f601a4 100644 --- a/src/modstartup.f90 +++ b/src/modstartup.f90 @@ -77,7 +77,7 @@ 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 @@ -85,7 +85,7 @@ subroutine startup(path) 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