From 71b3b06143a85239b4b8ff007654b9b9d2e41808 Mon Sep 17 00:00:00 2001 From: Knut Date: Mon, 11 Apr 2016 15:14:44 +0200 Subject: [PATCH] save_fluxes: flexible output of velx,vely,velx3d,vely3d --- src/2d/Makefile | 1 + src/2d/dynamic_allocations_2d.h | 6 +++ src/2d/dynamic_declarations_2d.h | 1 + src/2d/m2d.F90 | 14 +++++ src/2d/static_2d.h | 1 + src/2d/variables_2d.F90 | 1 + src/2d/velocity_update.F90 | 83 +++++++++++++++++++++++++++++ src/3d/Makefile | 3 +- src/3d/dynamic_allocations_3d.h | 9 ++++ src/3d/dynamic_declarations_3d.h | 1 + src/3d/m3d.F90 | 4 ++ src/3d/static_3d.h | 2 + src/3d/variables_3d.F90 | 1 + src/3d/velocity_update_3d.F90 | 59 ++++++++++++++++++++ src/CMakeLists.txt | 2 + src/getm/register_all_variables.F90 | 5 ++ src/ncdf/save_2d_ncdf.F90 | 15 ++---- src/ncdf/save_3d_ncdf.F90 | 29 ++-------- 18 files changed, 199 insertions(+), 38 deletions(-) create mode 100644 src/2d/velocity_update.F90 create mode 100644 src/3d/velocity_update_3d.F90 diff --git a/src/2d/Makefile b/src/2d/Makefile index d22ed800..3b42f7eb 100644 --- a/src/2d/Makefile +++ b/src/2d/Makefile @@ -22,6 +22,7 @@ ${LIB}(cfl_check.o) \ ${LIB}(bottom_friction.o) \ ${LIB}(depth_update.o) \ ${LIB}(momentum.o) \ +${LIB}(velocity_update.o) \ ${LIB}(sealevel.o) \ ${LIB}(uv_advect.o) \ ${LIB}(uv_diffusion.o) \ diff --git a/src/2d/dynamic_allocations_2d.h b/src/2d/dynamic_allocations_2d.h index 8932be34..e5bf8b4b 100644 --- a/src/2d/dynamic_allocations_2d.h +++ b/src/2d/dynamic_allocations_2d.h @@ -30,6 +30,12 @@ allocate(V(E2DFIELD),stat=rc) if (rc /= 0) stop 'init_2d: Error allocating memory (V)' + allocate(velx(E2DFIELD),stat=rc) + if (rc /= 0) stop 'init_2d: Error allocating memory (velx)' + + allocate(vely(E2DFIELD),stat=rc) + if (rc /= 0) stop 'init_2d: Error allocating memory (vely)' + allocate(UEx(E2DFIELD),stat=rc) if (rc /= 0) stop 'init_2d: Error allocating memory (UEx)' diff --git a/src/2d/dynamic_declarations_2d.h b/src/2d/dynamic_declarations_2d.h index db541c43..8179455e 100644 --- a/src/2d/dynamic_declarations_2d.h +++ b/src/2d/dynamic_declarations_2d.h @@ -5,6 +5,7 @@ REALTYPE,dimension(:,:),allocatable,target :: t_zo,t_z,D,Dvel,DU,DV REALTYPE,dimension(:,:),allocatable :: U,V REALTYPE,dimension(:,:),allocatable :: UEx,VEx + REALTYPE,dimension(:,:),allocatable,target :: velx,vely REALTYPE,dimension(:,:),allocatable :: ru,rv REALTYPE,dimension(:,:),allocatable :: Uint,Vint REALTYPE,dimension(:,:),allocatable :: Uinto,Vinto diff --git a/src/2d/m2d.F90 b/src/2d/m2d.F90 index 572cf4cf..2409e13f 100644 --- a/src/2d/m2d.F90 +++ b/src/2d/m2d.F90 @@ -43,6 +43,17 @@ REALTYPE,dimension(E2DFIELD),intent(out) :: D,Dvel,DU,DV end subroutine depth_update + subroutine velocity_update(dt,z,zo,Dvel,U,DU,V,DV,wwm,wwp,missing, & + velx,vely) + use domain, only: imin,imax,jmin,jmax + IMPLICIT NONE + REALTYPE,intent(in) :: dt + REALTYPE,dimension(E2DFIELD),intent(in) :: z,zo,Dvel,U,DU,V,DV + REALTYPE,dimension(E2DFIELD),target,intent(in),optional :: wwm,wwp + REALTYPE,intent(in),optional :: missing + REALTYPE,dimension(E2DFIELD),intent(out) :: velx,vely + end subroutine velocity_update + subroutine uv_advect(U,V,D,Dvel,DU,DV) use domain, only: imin,imax,jmin,jmax IMPLICIT NONE @@ -412,6 +423,8 @@ end if + call velocity_update(dtm,z,zo,Dvel,U,DU,V,DV,velx=velx,vely=vely) + return end subroutine postinit_2d !EOC @@ -488,6 +501,7 @@ if (have_boundaries) call update_2d_bdy(loop,bdy2d_ramp) call sealevel(loop) call depth_update(zo,z,D,Dvel,DU,DV) + call velocity_update(dtm,z,zo,Dvel,U,DU,V,DV,velx=velx,vely=vely) if(residual .gt. 0) then call tic(TIM_INTEGR2D) diff --git a/src/2d/static_2d.h b/src/2d/static_2d.h index d6dd02c5..27da815d 100644 --- a/src/2d/static_2d.h +++ b/src/2d/static_2d.h @@ -8,6 +8,7 @@ REALTYPE,dimension(E2DFIELD),target :: t_zo,t_z,D,Dvel,DU,DV REALTYPE U(E2DFIELD) REALTYPE V(E2DFIELD) + REALTYPE,dimension(E2DFIELD),target :: velx,vely REALTYPE UEx(E2DFIELD) REALTYPE VEx(E2DFIELD) REALTYPE ru(E2DFIELD) diff --git a/src/2d/variables_2d.F90 b/src/2d/variables_2d.F90 index 26954a38..a2937cd4 100644 --- a/src/2d/variables_2d.F90 +++ b/src/2d/variables_2d.F90 @@ -125,6 +125,7 @@ D = _ZERO_ ; Dvel = _ZERO_ U = _ZERO_; DU = _ZERO_; Uint = _ZERO_; UEx = _ZERO_ V = _ZERO_; DV = _ZERO_; Vint = _ZERO_; VEx = _ZERO_ + velx = -9999.0 ; vely = -9999.0 ru = _ZERO_; ruu=_ZERO_; Uinto=_ZERO_ rv = _ZERO_; rvv=_ZERO_; Vinto=_ZERO_ diff --git a/src/2d/velocity_update.F90 b/src/2d/velocity_update.F90 new file mode 100644 index 00000000..cc4d9143 --- /dev/null +++ b/src/2d/velocity_update.F90 @@ -0,0 +1,83 @@ +#include "cppdefs.h" +!----------------------------------------------------------------------- +!BOP +! +! !ROUTINE: velocity_update - calculate new velocities. +! +! !INTERFACE: + subroutine velocity_update(dt,z,zo,Dvel,U,DU,V,DV,wwm,wwp,missing, & + velx,vely) + +! Note (KK): keep in sync with interface in m2d.F90 +! +! !DESCRIPTION: +! +! !USES: + use domain + IMPLICIT NONE +! +! !INPUT PARAMETERS: + REALTYPE,intent(in) :: dt + REALTYPE,dimension(E2DFIELD),intent(in) :: z,zo,Dvel,U,DU,V,DV + REALTYPE,dimension(E2DFIELD),target,intent(in),optional :: wwm,wwp + REALTYPE,intent(in),target,optional :: missing +! +! !OUTPUT PARAMETERS: + REALTYPE,dimension(E2DFIELD),intent(out) :: velx,vely +! +! !OUTPUT PARAMETERS: +! +! !REVISION HISTORY: +! Original author(s): Knut Klingbeil +! +! !LOCAL VARIABLES: + REALTYPE,dimension(E2DFIELD),target :: zeros + REALTYPE,dimension(:,:),pointer :: p_wwm,p_wwp + REALTYPE :: missval + REALTYPE,parameter :: vel_missing=-9999.0 +!EOP +!----------------------------------------------------------------------- +!BOC +#ifdef DEBUG + integer, save :: Ncall = 0 + Ncall = Ncall+1 + write(debug,*) 'velocity_update() # ',Ncall +#endif + + zeros = _ZERO_ + if (present(wwm)) then + p_wwm => wwm + else + p_wwm => zeros + end if + if (present(wwp)) then + p_wwp => wwp + else + p_wwp => zeros + end if + if (present(missing)) then + missval = missing + else + missval = vel_missing + end if + + call to_u(imin,jmin,imax,jmax,az, & + dt,grid_type, & + dxv,dyu,arcd1, & + xc,xu,xv,z,zo,Dvel,U,DU,V,DV,p_wwm,p_wwp,missval,velx) + call to_v(imin,jmin,imax,jmax,az, & + dt,grid_type, & + dxv,dyu,arcd1, & + yc,yu,yv,z,zo,Dvel,U,DU,V,DV,p_wwm,p_wwp,missval,vely) + +#ifdef DEBUG + write(debug,*) 'Leaving velocity_update()' + write(debug,*) +#endif + return + end subroutine velocity_update +!EOC + +!----------------------------------------------------------------------- +! Copyright (C) 2016 - Knut Klingbeil ! +!----------------------------------------------------------------------- diff --git a/src/3d/Makefile b/src/3d/Makefile index 63780a49..4a89a7cc 100644 --- a/src/3d/Makefile +++ b/src/3d/Makefile @@ -51,7 +51,8 @@ $(LIB)(hybrid_coordinates.o) \ $(LIB)(bottom_friction_3d.o) \ $(LIB)(uu_momentum_3d.o) \ $(LIB)(vv_momentum_3d.o) \ -$(LIB)(ww_momentum_3d.o) +$(LIB)(ww_momentum_3d.o) \ +$(LIB)(velocity_update_3d.o) ifeq ($(GETM_STRUCTURE_FRICTION),true) OBJ += \ $(LIB)(structure_friction_3d.o) diff --git a/src/3d/dynamic_allocations_3d.h b/src/3d/dynamic_allocations_3d.h index 1574be49..768372e7 100644 --- a/src/3d/dynamic_allocations_3d.h +++ b/src/3d/dynamic_allocations_3d.h @@ -7,6 +7,15 @@ allocate(ww(I3DFIELD),stat=rc) ! 3D field for w-velocity if (rc /= 0) stop 'init_3d: Error allocating memory (ww)' + allocate(velx3d(I3DFIELD),stat=rc) ! 3D field for u-velocity + if (rc /= 0) stop 'init_3d: Error allocating memory (velx3d)' + + allocate(vely3d(I3DFIELD),stat=rc) ! 3D field for v-velocity + if (rc /= 0) stop 'init_3d: Error allocating memory (vely3d)' + + allocate(w(I3DFIELD),stat=rc) ! 3D field for w-velocity + if (rc /= 0) stop 'init_3d: Error allocating memory (w)' + #ifdef _MOMENTUM_TERMS_ allocate(tdv_u(I3DFIELD),stat=rc) ! 3D field for tdv_u if (rc /= 0) stop 'init_3d: Error allocating memory (tdv_u)' diff --git a/src/3d/dynamic_declarations_3d.h b/src/3d/dynamic_declarations_3d.h index 18c5afd5..5019986f 100644 --- a/src/3d/dynamic_declarations_3d.h +++ b/src/3d/dynamic_declarations_3d.h @@ -7,6 +7,7 @@ REALTYPE, dimension(:,:,:), allocatable :: uu,vv REALTYPE, dimension(:,:,:), allocatable, target :: ww + REALTYPE, dimension(:,:,:), allocatable, target :: velx3d,vely3d,w #ifdef _MOMENTUM_TERMS_ REALTYPE, dimension(:,:,:), allocatable, target :: tdv_u,adv_u,vsd_u diff --git a/src/3d/m3d.F90 b/src/3d/m3d.F90 index 2699e419..3254e7f0 100644 --- a/src/3d/m3d.F90 +++ b/src/3d/m3d.F90 @@ -415,6 +415,8 @@ end if + call velocity_update_3d() + return end subroutine postinit_3d !EOC @@ -632,6 +634,8 @@ end if + call velocity_update_3d() + call slow_terms() #endif call tic(TIM_INTEGR3D) diff --git a/src/3d/static_3d.h b/src/3d/static_3d.h index 1289899c..c5f9510f 100644 --- a/src/3d/static_3d.h +++ b/src/3d/static_3d.h @@ -15,6 +15,8 @@ REALTYPE :: uu(I3DFIELD) REALTYPE :: vv(I3DFIELD) REALTYPE, target :: ww(I3DFIELD) + REALTYPE,dimension(I3DFIELD),target :: velx3d,vely3d,w + #ifdef _MOMENTUM_TERMS_ REALTYPE :: tdv_u(I3DFIELD) REALTYPE :: adv_u(I3DFIELD) diff --git a/src/3d/variables_3d.F90 b/src/3d/variables_3d.F90 index 3bf1dca6..ce7dd56f 100644 --- a/src/3d/variables_3d.F90 +++ b/src/3d/variables_3d.F90 @@ -194,6 +194,7 @@ hn = _ZERO_ ; hvel = _ZERO_ ; hun = _ZERO_ ; hvn = _ZERO_ uu = _ZERO_ ; vv = _ZERO_ ; ww = _ZERO_ + velx3d = -9999.0 ; vely3d = -9999.0 ; w = -9999.0 #ifdef _MOMENTUM_TERMS_ tdv_u = _ZERO_ ; adv_u = _ZERO_ ; vsd_u = _ZERO_ ; hsd_u = _ZERO_ diff --git a/src/3d/velocity_update_3d.F90 b/src/3d/velocity_update_3d.F90 new file mode 100644 index 00000000..5585682f --- /dev/null +++ b/src/3d/velocity_update_3d.F90 @@ -0,0 +1,59 @@ +#include "cppdefs.h" +!----------------------------------------------------------------------- +!BOP +! +! !ROUTINE: velocity_update_3d - calculate new 3D velocities. +! +! !INTERFACE: + subroutine velocity_update_3d() +! +! !DESCRIPTION: +! +! !USES: + use domain + use m2d, only: velocity_update + use variables_3d, only: kmin + use variables_3d, only: hn,ho,uu,vv,ww,hun,hvn,hvel,velx3d,vely3d,w,dt + IMPLICIT NONE +! +! !INPUT PARAMETERS: +! +! !OUTPUT PARAMETERS: +! +! !REVISION HISTORY: +! Original author(s): Knut Klingbeil +! +! !LOCAL VARIABLES: + integer :: k + REALTYPE,parameter :: vel_missing=-9999.0 +!EOP +!----------------------------------------------------------------------- +!BOC +#ifdef DEBUG + integer, save :: Ncall = 0 + Ncall = Ncall+1 + write(debug,*) 'velocity_update_3d() # ',Ncall +#endif + + do k=1,kmax + call velocity_update(dt,hn(:,:,k),ho(:,:,k),hvel(:,:,k), & + uu(:,:,k),hun(:,:,k),vv(:,:,k),hvn(:,:,k), & + wwm=ww(:,:,k-1),wwp=ww(:,:,k), & + velx=velx3d(:,:,k),vely=vely3d(:,:,k)) + end do + call to_w(imin,jmin,imax,jmax,kmin,kmax,az, & + dt, & + dxv,dyu,arcd1, & + H,HU,HV,hn,ho,hvel,uu,hun,vv,hvn,ww,vel_missing,w) + +#ifdef DEBUG + write(debug,*) 'Leaving velocity_update_3d()' + write(debug,*) +#endif + return + end subroutine velocity_update_3d +!EOC + +!----------------------------------------------------------------------- +! Copyright (C) 2016 - Knut Klingbeil ! +!----------------------------------------------------------------------- diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index d5d482a9..fd47bf48 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -184,6 +184,7 @@ add_library(2d OBJECT 2d/uv_diff_2dh.F90 2d/uv_diffusion.F90 2d/variables_2d.F90 + 2d/velocity_update.F90 ) set_property(TARGET 2d APPEND PROPERTY INCLUDE_DIRECTORIES "${PROJECT_SOURCE_DIR}/2d") @@ -231,6 +232,7 @@ add_library(3d OBJECT 3d/uv_advect_3d.F90 3d/uv_diffusion_3d.F90 3d/variables_3d.F90 + 3d/velocity_update_3d.F90 3d/vv_momentum_3d.F90 3d/ww_momentum_3d.F90 ) diff --git a/src/getm/register_all_variables.F90 b/src/getm/register_all_variables.F90 index 23d4a42c..6590c686 100644 --- a/src/getm/register_all_variables.F90 +++ b/src/getm/register_all_variables.F90 @@ -315,6 +315,8 @@ call fm%register('D', 'm', 'water depth', standard_name='water depth', fill_value=0.0_rk, data2d=D(_2D_W_), category="2d") call fm%register('U', 'm2/s', 'transport in local x-direction', standard_name='', data2d=U(_2D_W_), category='2d', output_level=output_level_debug) call fm%register('V', 'm2/s', 'transport in local y-direction', standard_name='', data2d=V(_2D_W_), category='2d', output_level=output_level_debug) + call fm%register('velx', 'm/s', 'velocity in global x-direction', standard_name='', data2d=velx(_2D_W_), category='2d', fill_value=-9999.0_rk, output_level=output_level_debug) + call fm%register('vely', 'm/s', 'velocity in global y-direction', standard_name='', data2d=vely(_2D_W_), category='2d', fill_value=-9999.0_rk, output_level=output_level_debug) return end subroutine register_2d_variables @@ -448,6 +450,9 @@ call fm%register('ssvn', 'm', 'elevtion at V-points (3D)', standard_name='', data2d=ssvn(_2D_W_), category='3d', output_level=output_level_debug) call fm%register('uu', 'm2/s', 'transport in local x-direction (3D)', standard_name='', dimensions=(/id_dim_z/), data3d=uu(_3D_W_), category='3d', output_level=output_level_debug) call fm%register('vv', 'm2/s', 'transport in local y-direction (3D)', standard_name='', dimensions=(/id_dim_z/), data3d=vv(_3D_W_), category='3d', output_level=output_level_debug) + call fm%register('velx3d', 'm/s', 'velocity in global x-direction (3D)', standard_name='', dimensions=(/id_dim_z/), data3d=velx3d(_3D_W_), category='3d', fill_value=-9999.0_rk, output_level=output_level_debug) + call fm%register('vely3d', 'm/s', 'velocity in global y-direction (3D)', standard_name='', dimensions=(/id_dim_z/), data3d=vely3d(_3D_W_), category='3d', fill_value=-9999.0_rk, output_level=output_level_debug) + call fm%register('w', 'm2/s', 'vertical velocity', standard_name='', dimensions=(/id_dim_z/), data3d=w(_3D_W_), category='3d', fill_value=-9999.0_rk, output_level=output_level_debug) call fm%register('SS', 's-2', 'shear frequency squared', standard_name='', dimensions=(/id_dim_z/), data3d=SS(_3D_W_), category='3d', output_level=output_level_debug) end if diff --git a/src/ncdf/save_2d_ncdf.F90 b/src/ncdf/save_2d_ncdf.F90 index f3d753cf..4e213e95 100644 --- a/src/ncdf/save_2d_ncdf.F90 +++ b/src/ncdf/save_2d_ncdf.F90 @@ -21,6 +21,7 @@ use domain, only: dxv,dyu,arcd1 use m2d, only: dtm use variables_2d, only: zo,z,D,Dvel,U,DU,V,DV,res_u,res_v + use variables_2d, only: velx,vely #ifdef USE_BREAKS use variables_2d, only: break_stat #endif @@ -93,21 +94,11 @@ ! velocites if (u_id .ne. -1) then - wrk = _ZERO_ - call to_u(imin,jmin,imax,jmax,az, & - dtm,grid_type, & - dxv,dyu,arcd1, & - xc,xu,xv,z,zo,Dvel,U,DU,V,DV,wrk,wrk,vel_missing,ws) - err = nf90_put_var(ncid,u_id,ws(_2D_W_),start,edges) + err = nf90_put_var(ncid,u_id,velx(_2D_W_),start,edges) if (err .NE. NF90_NOERR) go to 10 end if if (v_id .ne. -1) then - wrk = _ZERO_ - call to_v(imin,jmin,imax,jmax,az, & - dtm,grid_type, & - dxv,dyu,arcd1, & - yc,yu,yv,z,zo,Dvel,U,DU,V,DV,wrk,wrk,vel_missing,ws) - err = nf90_put_var(ncid,v_id,ws(_2D_W_),start,edges) + err = nf90_put_var(ncid,v_id,vely(_2D_W_),start,edges) if (err .NE. NF90_NOERR) go to 10 end if diff --git a/src/ncdf/save_3d_ncdf.F90 b/src/ncdf/save_3d_ncdf.F90 index f7d2c0a1..b81585a9 100644 --- a/src/ncdf/save_3d_ncdf.F90 +++ b/src/ncdf/save_3d_ncdf.F90 @@ -22,6 +22,7 @@ use variables_2d, only: Uinto,Vinto use variables_3d, only: sseo,ssen,Dn,Dveln,Dun,Dvn use variables_3d, only: dt,kmin,ho,hn,hvel,uu,hun,vv,hvn,ww,hcc,SS + use variables_3d, only: velx3d,vely3d,w use variables_3d, only: taubx,tauby #ifdef _MOMENTUM_TERMS_ use variables_3d, only: tdv_u,adv_u,vsd_u,hsd_u,cor_u,epg_u,ipg_u @@ -210,37 +211,15 @@ ! velocites if (uu_id .ne. -1) then - ws(:,:,0) = vel_missing - do k=1,kmax - call to_u(imin,jmin,imax,jmax,az, & - dt,grid_type, & - dxv,dyu,arcd1, & - xc,xu,xv,hn(:,:,k),ho(:,:,k),hvel(:,:,k), & - uu(:,:,k),hun(:,:,k),vv(:,:,k),hvn(:,:,k), & - ww(:,:,k-1),ww(:,:,k),vel_missing,ws(:,:,k)) - end do - err = nf90_put_var(ncid,uu_id,ws(_3D_W_),start,edges) + err = nf90_put_var(ncid,uu_id,velx3d(_3D_W_),start,edges) if (err .NE. NF90_NOERR) go to 10 end if if (vv_id .ne. -1) then - ws(:,:,0) = vel_missing - do k=1,kmax - call to_v(imin,jmin,imax,jmax,az, & - dt,grid_type, & - dxv,dyu,arcd1, & - yc,yu,yv,hn(:,:,k),ho(:,:,k),hvel(:,:,k), & - uu(:,:,k),hun(:,:,k),vv(:,:,k),hvn(:,:,k), & - ww(:,:,k-1),ww(:,:,k),vel_missing,ws(:,:,k)) - end do - err = nf90_put_var(ncid,vv_id,ws(_3D_W_),start,edges) + err = nf90_put_var(ncid,vv_id,vely3d(_3D_W_),start,edges) if (err .NE. NF90_NOERR) go to 10 end if if (w_id .ne. -1) then - call to_w(imin,jmin,imax,jmax,kmin,kmax,az, & - dt, & - dxv,dyu,arcd1, & - H,HU,HV,hn,ho,hvel,uu,hun,vv,hvn,ww,vel_missing,ws) - err = nf90_put_var(ncid,w_id,ws(_3D_W_),start,edges) + err = nf90_put_var(ncid,w_id,w(_3D_W_),start,edges) if (err .NE. NF90_NOERR) go to 10 end if -- GitLab