Commit 03ad1fb8 authored by Hans Burchard's avatar Hans Burchard

output of momentum terms added

parent 7bea26b6
......@@ -7,6 +7,52 @@
allocate(ww(I3DFIELD),stat=rc) ! 3D field for w-velocity
if (rc /= 0) stop 'init_3d: Error allocating memory (ww)'
#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)'
allocate(adv_u(I3DFIELD),stat=rc) ! 3D field for adv_u
if (rc /= 0) stop 'init_3d: Error allocating memory (adv_u)'
allocate(vsd_u(I3DFIELD),stat=rc) ! 3D field for vsd_u
if (rc /= 0) stop 'init_3d: Error allocating memory (vsd_u)'
allocate(hsd_u(I3DFIELD),stat=rc) ! 3D field for hsd_u
if (rc /= 0) stop 'init_3d: Error allocating memory (hsd_u)'
allocate(cor_u(I3DFIELD),stat=rc) ! 3D field for cor_u
if (rc /= 0) stop 'init_3d: Error allocating memory (cor_u)'
allocate(epg_u(I3DFIELD),stat=rc) ! 3D field for epg_u
if (rc /= 0) stop 'init_3d: Error allocating memory (epg_u)'
allocate(ipg_u(I3DFIELD),stat=rc) ! 3D field for ipg_u
if (rc /= 0) stop 'init_3d: Error allocating memory (ipg_u)'
allocate(tdv_v(I3DFIELD),stat=rc) ! 3D field for tdv_v
if (rc /= 0) stop 'init_3d: Error allocating memory (tdv_v)'
allocate(adv_v(I3DFIELD),stat=rc) ! 3D field for adv_v
if (rc /= 0) stop 'init_3d: Error allocating memory (adv_v)'
allocate(vsd_v(I3DFIELD),stat=rc) ! 3D field for vsd_v
if (rc /= 0) stop 'init_3d: Error allocating memory (vsd_v)'
allocate(hsd_v(I3DFIELD),stat=rc) ! 3D field for hsd_v
if (rc /= 0) stop 'init_3d: Error allocating memory (hsd_v)'
allocate(cor_v(I3DFIELD),stat=rc) ! 3D field for cor_v
if (rc /= 0) stop 'init_3d: Error allocating memory (cor_v)'
allocate(epg_v(I3DFIELD),stat=rc) ! 3D field for epg_v
if (rc /= 0) stop 'init_3d: Error allocating memory (epg_v)'
allocate(ipg_v(I3DFIELD),stat=rc) ! 3D field for ipg_v
if (rc /= 0) stop 'init_3d: Error allocating memory (ipg_v)'
#endif
#ifdef STRUCTURE_FRICTION
allocate(sf(I3DFIELD),stat=rc) ! 3D field for velocity in T-points
if (rc /= 0) stop 'init_3d: Error allocating memory (sf)'
......
......@@ -307,7 +307,7 @@
do k=1,kmax ! Calculating v-interface high-order fluxes !
do j=jmin-1,jmax
do i=imin,imax
uuu=vv(i,j,k)*dt/delyv(i,j)
uuu=vv(i,j,k)/hvn(i,j,k)*dt/delyv(i,j)
vvv=0.25*( &
uu(i-1,j,k)/hun(i-1,j,k) &
+uu(i-1,j+1,k)/hun(i-1,j+1,k) &
......
......@@ -23,6 +23,23 @@
REALTYPE :: uu(I3DFIELD)
REALTYPE :: vv(I3DFIELD)
REALTYPE, target :: ww(I3DFIELD)
#ifdef _MOMENTUM_TERMS_
REALTYPE :: tdv_u(I3DFIELD)
REALTYPE :: adv_u(I3DFIELD)
REALTYPE :: vsd_u(I3DFIELD)
REALTYPE :: hsd_u(I3DFIELD)
REALTYPE :: cor_u(I3DFIELD)
REALTYPE :: epg_u(I3DFIELD)
REALTYPE :: ipg_u(I3DFIELD)
REALTYPE :: tdv_v(I3DFIELD)
REALTYPE :: adv_v(I3DFIELD)
REALTYPE :: vsd_v(I3DFIELD)
REALTYPE :: hsd_v(I3DFIELD)
REALTYPE :: cor_v(I3DFIELD)
REALTYPE :: epg_v(I3DFIELD)
REALTYPE :: ipg_v(I3DFIELD)
#endif
#ifdef STRUCTURE_FRICTION
REALTYPE :: sf(I3DFIELD)
#endif
......
......@@ -61,6 +61,9 @@
use variables_3d, only: dt,cnpar,kumin,uu,vv,huo,hun,hvo,uuEx,ww,hvn
use variables_3d, only: num,nuh,sseo,ssun,rru
use variables_3d, only: ssuo
#ifdef _MOMENTUM_TERMS_
use variables_3d, only: tdv_u,cor_u,ipg_u,epg_u,vsd_u
#endif
#ifdef STRUCTURE_FRICTION
use variables_3d, only: sf
#endif
......@@ -168,10 +171,16 @@
#else
ex(k)=coru(i,j)*Vloc
#endif
#ifdef _MOMENTUM_TERMS_
cor_u(i,j,k)=-dry_u(i,j)*ex(k)
#endif
#ifdef NO_BAROCLINIC
ex(k)=dry_u(i,j)*(ex(k)-uuEx(i,j,k))
#else
ex(k)=dry_u(i,j)*(ex(k)-uuEx(i,j,k)+ip_fac*idpdx(i,j,k))
#ifdef _MOMENTUM_TERMS_
ipg_u(i,j,k)=-dry_u(i,j)*ip_fac*idpdx(i,j,k)
#endif
#endif
end do
ex(kmax)=ex(kmax) &
......@@ -249,10 +258,27 @@
#endif
do k=kumin(i,j),kmax
#ifdef _MOMENTUM_TERMS_
tdv_u(i,j,k)=uu(i,j,k)
epg_u(i,j,k)=_HALF_*(huo(i,j,k)+hun(i,j,k))*g*zx &
-hun(i,j,k)*Diff/dt
vsd_u(i,j,k)=auxo(k)*(uu(i,j,k+1)/huo(i,j,k+1) &
-uu(i,j,k)/huo(i,j,k)) &
-auxo(k-1)*(uu(i,j,k)/huo(i,j,k) &
-uu(i,j,k-1)/huo(i,j,k-1))
#endif
#ifndef NO_BAROTROPIC
uu(i,j,k)=Res(k) +hun(i,j,k)*Diff
#else
uu(i,j,k)=Res(k)
#endif
#ifdef _MOMENTUM_TERMS_
tdv_u(i,j,k)=(uu(i,j,k)-tdv_u(i,j,k))/dt
vsd_u(i,j,k)=-(vsd_u(i,j,k) &
+auxn(k)*(uu(i,j,k+1)/hun(i,j,k+1) &
-uu(i,j,k)/hun(i,j,k)) &
-auxn(k-1)*(uu(i,j,k)/hun(i,j,k) &
-uu(i,j,k-1)/hun(i,j,k-1)))/dt
#endif
end do
else ! if (kmax .eq. kumin(i,j))
......@@ -260,7 +286,7 @@
end if
end if
end do
end do
end do
!$OMP END DO
#ifdef SLICE_MODEL
......
......@@ -267,6 +267,9 @@
use variables_3d, only: dt,kumin,kvmin,uu,vv,ww,hun,hvn,huo,hvo,uuEx,vvEx
#ifdef UV_TVD
use variables_3d, only: uadv,vadv,wadv,huadv,hvadv,hoadv,hnadv
#endif
#ifdef _MOMENTUM_TERMS_
use variables_3d, only: adv_u,adv_v
#endif
use advection_3d, only: do_advection_3d
use halo_zones, only: update_3d_halo,wait_halo,U_TAG,V_TAG
......@@ -300,7 +303,9 @@
REALTYPE :: dti,dxdyi
REALTYPE :: vel2(I3DFIELD)
REALTYPE :: vel2o(I3DFIELD)
REALTYPE :: vel(I3DFIELD)
REALTYPE :: numdiss(I3DFIELD)
REALTYPE :: momtoto=0.,momtot=0.
#endif
!EOP
!-----------------------------------------------------------------------
......@@ -331,6 +336,7 @@
do j=jmin-HALO,jmax+HALO
do i=imin-HALO,imax+HALO-1
uadv(i,j,k)=_HALF_*(uu(i+1,j,k)+uu(i,j,k))
if ((au(i+1,j).eq.0).or.(au(i,j).eq.0)) uadv(i,j,k)=_ZERO_
vadv(i,j,k)=_HALF_*(vv(i+1,j,k)+vv(i,j,k))
wadv(i,j,k)=_HALF_*(ww(i+1,j,k)+ww(i,j,k))
huadv(i,j,k)=_HALF_*(hun(i+1,j,k)+hun(i,j,k))
......@@ -342,11 +348,12 @@
!$OMP END DO NOWAIT
end do
uuEx=_ZERO_
do k=1,kmax ! uuEx is here the velocity to be transported.
!$OMP DO SCHEDULE(RUNTIME)
do j=jmin,jmax
do i=imin,imax
uuEx(i,j,k)=uu(i,j,k)/huo(i,j,k)
uuEx(i,j,k)=uu(i,j,k)/hoadv(i,j,k)
end do
end do
!$OMP END DO NOWAIT
......@@ -374,6 +381,13 @@
dyvadv(i,j)=dy
area_inv(i,j)=dxdyi
#endif
do k=1,kmax
hnadv(i,j,k)=hoadv(i,j,k)-dt*((uadv(i,j,k)*dyuadv(i,j) &
-uadv(i-1,j,k)*dyuadv(i-1,j)) &
+(vadv(i,j,k)*dxvadv(i,j) &
-vadv(i,j-1,k)*dxvadv(i,j-1)))*area_inv(i,j) &
-dt*(wadv(i,j,k)-wadv(i,j,k-1))
end do
end do
end do
......@@ -382,17 +396,30 @@
call wait_halo(U_TAG)
call toc(TIM_UVADV3DH)
vel=_ONE_
if (do_numerical_analyses) then
do k=1,kmax ! calculate square of u-velocity before advection step
do j=jmin,jmax
do i=imin,imax
vel2(i,j,k)=uuEx(i,j,k)**2
if (au(i,j).eq.0) vel(i,j,k)=0.
! momtoto=momtoto+uuEx(i,j,k)*hoadv(i,j,k)
momtoto=momtoto+vel(i,j,k)*hoadv(i,j,k)
end do
end do
end do
vel2o=vel2
end if
call do_advection_3d(dt,vel,uadv,vadv,wadv,huadv,hvadv,hoadv,hnadv,&
dxuadv,dxvadv,dyuadv,dyvadv,area_inv, &
azadv,auadv,avadv,hor_adv,ver_adv,adv_split,AH)
i=50
j=15
k=30
write(101,*) uuEx(i,j,k),'old'
call do_advection_3d(dt,uuEx,uadv,vadv,wadv,huadv,hvadv,hoadv,hnadv,&
dxuadv,dxvadv,dyuadv,dyvadv,area_inv, &
azadv,auadv,avadv,hor_adv,ver_adv,adv_split,AH)
......@@ -402,13 +429,20 @@
dxuadv,dxvadv,dyuadv,dyvadv,area_inv, &
azadv,auadv,avadv,hor_adv,ver_adv,adv_split,AH)
i=50
j=15
k=30
write(101,*) uuEx(i,j,k),'new'
do k=1,kmax ! calculate kinetic energy dissipaion rate for u-velocity
do j=jmin,jmax
do i=imin,imax
numdiss(i,j,k)=(vel2(i,j,k)-uuEx(i,j,k)**2)/dt
! momtot= momtot + uuEx(i,j,k)*hnadv(i,j,k)
momtot= momtot + vel(i,j,k)*hnadv(i,j,k)
end do
end do
end do
write(100,*) momtot-momtoto,momtot
numdis2d=_ZERO_
do k=1,kmax ! calculate kinetic energy dissipaion rate for u-velocity
......@@ -433,6 +467,9 @@
do j=jmin-HALO,jmax+HALO
do i=imin-HALO,imax+HALO
uuEx(i,j,k)=dti*(uu(i,j,k)-uuEx(i,j,k)*hun(i,j,k))
#ifdef _MOMENTUM_TERMS_
adv_u(i,j,k)=uuEx(i,j,k)
#endif
end do
end do
!$OMP END DO NOWAIT
......@@ -543,6 +580,9 @@
do j=jmin-HALO,jmax+HALO
do i=imin-HALO,imax+HALO
vvEx(i,j,k)=dti*(vv(i,j,k)-vvEx(i,j,k)*hvn(i,j,k))
#ifdef _MOMENTUM_TERMS_
adv_v(i,j,k)=vvEx(i,j,k)
#endif
end do
end do
!$OMP END DO NOWAIT
......
......@@ -119,6 +119,9 @@
#endif
use variables_3d, only: kumin,kvmin,uu,vv,ww,hn,hun,hvn,uuEx,vvEx
use getm_timers, only: tic, toc, TIM_UVDIFF3D
#ifdef _MOMENTUM_TERMS_
use variables_3d, only: hsd_u,hsd_v
#endif
!$ use omp_lib
IMPLICIT NONE
!
......@@ -168,6 +171,9 @@
if (au(i,j) .ge. 1) then
if (k .ge. kumin(i,j)) then
uuEx(i,j,k)=uuEx(i,j,k)-(PP(i+1,j,k)-PP(i,j,k))*ARUD1
#ifdef _MOMENTUM_TERMS_
hsd_u(i,j,k)=-(PP(i+1,j,k)-PP(i,j,k))*ARUD1
#endif
end if
end if
end do
......@@ -202,6 +208,9 @@
if (au(i,j) .ge. 1) then
if (k .ge. kumin(i,j)) then
uuEx(i,j,k)=uuEx(i,j,k)-(PP(i,j,k)-PP(i,j-1,k))*ARUD1
#ifdef _MOMENTUM_TERMS_
hsd_u(i,j,k)=hsd_u(i,j,k)-(PP(i+1,j,k)-PP(i,j,k))*ARUD1
#endif
end if
end if
end do
......@@ -236,6 +245,9 @@
if (av(i,j) .ge. 1) then
if (k .ge. kvmin(i,j)) then
vvEx(i,j,k)=vvEx(i,j,k)-(PP(i,j,k)-PP(i-1,j,k))*ARVD1
#ifdef _MOMENTUM_TERMS_
hsd_v(i,j,k)=-(PP(i,j,k)-PP(i-1,j,k))*ARVD1
#endif
end if
end if
end do
......@@ -269,6 +281,9 @@
if (av(i,j) .ge. 1) then
if (k .ge. kvmin(i,j)) then
vvEx(i,j,k)=(vvEx(i,j,k)-(PP(i,j+1,k)-PP(i,j,k))*ARVD1)
#ifdef _MOMENTUM_TERMS_
hsd_v(i,j,k)=hsd_v(i,j,k)-(PP(i,j+1,k)-PP(i,j,k))*ARVD1
#endif
end if
end if
end do
......
......@@ -187,6 +187,12 @@
hn = _ZERO_ ; hun = _ZERO_ ; hvn = _ZERO_
uu = _ZERO_ ; vv = _ZERO_ ; ww = _ZERO_
#ifdef _MOMENTUM_TERMS_
tdv_u = _ZERO_ ; adv_u = _ZERO_ ; vsd_u = _ZERO_ ; hsd_u = _ZERO_
cor_u = _ZERO_ ; epg_u = _ZERO_ ; ipg_u = _ZERO_
tdv_v = _ZERO_ ; adv_v = _ZERO_ ; vsd_v = _ZERO_ ; hsd_v = _ZERO_
cor_v = _ZERO_ ; epg_v = _ZERO_ ; ipg_v = _ZERO_
#endif
ssen = _ZERO_ ; ssun = _ZERO_ ; ssvn = _ZERO_
rru= _ZERO_ ; rrv= _ZERO_
uuEx= _ZERO_ ; vvEx= _ZERO_
......@@ -204,6 +210,8 @@
huadv = _ZERO_ ; hvadv = _ZERO_
#endif
#ifndef NO_BAROCLINIC
idpdx=_ZERO_
idpdy=_ZERO_
......
......@@ -66,6 +66,9 @@
use variables_3d, only: dt,cnpar,kvmin,uu,vv,huo,hvo,hvn,vvEx,ww,hun
use variables_3d, only: num,nuh,sseo,ssvn,rrv
use variables_3d, only: ssvo
#ifdef _MOMENTUM_TERMS_
use variables_3d, only: tdv_v,cor_v,ipg_v,epg_v,vsd_v
#endif
#ifdef XZ_PLUME_TEST
use variables_3d, only: buoy
#endif
......@@ -183,6 +186,9 @@
#else
ex(k)=-corv(i,j)*Uloc
#endif
#ifdef _MOMENTUM_TERMS_
cor_v(i,j,k)=-dry_v(i,j)*ex(k)
#endif
#ifdef NO_BAROCLINIC
ex(k)=dry_v(i,j)*(ex(k)-vvEx(i,j,k))
#else
......@@ -191,6 +197,9 @@
#else
ex(k)=dry_v(i,j)*(ex(k)-vvEx(i,j,k)+ip_fac*idpdy(i,j,k))
#endif
#ifdef _MOMENTUM_TERMS_
ipg_v(i,j,k)=-dry_v(i,j)*ip_fac*idpdy(i,j,k)
#endif
#endif
end do
ex(kmax)=ex(kmax) &
......@@ -271,10 +280,27 @@
do k=kvmin(i,j),kmax
#ifdef _MOMENTUM_TERMS_
tdv_v(i,j,k)=vv(i,j,k)
epg_v(i,j,k)=_HALF_*(hvo(i,j,k)+hvn(i,j,k))*g*zy &
-hvn(i,j,k)*Diff/dt
vsd_v(i,j,k)=-(auxo(k)*(vv(i,j,k+1)/hvo(i,j,k+1) &
-vv(i,j,k)/hvo(i,j,k)) &
-auxo(k-1)*(vv(i,j,k)/hvo(i,j,k) &
-vv(i,j,k-1)/hvo(i,j,k-1)))
#endif
#ifndef NO_BAROTROPIC
vv(i,j,k)=Res(k)+hvn(i,j,k)*Diff
#else
vv(i,j,k)=Res(k)
#endif
#ifdef _MOMENTUM_TERMS_
tdv_v(i,j,k)=(vv(i,j,k)-tdv_v(i,j,k))/dt
vsd_v(i,j,k)=(vsd_v(i,j,k) &
+auxn(k)*(vv(i,j,k+1)/hvn(i,j,k+1) &
-vv(i,j,k)/hvn(i,j,k)) &
-auxn(k-1)*(vv(i,j,k)/hvn(i,j,k) &
-vv(i,j,k-1)/hvn(i,j,k-1)))/dt
#endif
end do
else ! (kmax .eq. kvmin(i,j))
......
......@@ -24,6 +24,7 @@ endif
# The file ../compilers/compiler.$(FORTRAN_COMPILER) must exist
DEFINES=-D$(FORTRAN_COMPILER)
DEFINES+=-D_MOMENTUM_TERMS_
include $(GETMDIR)/compilers/compiler.$(FORTRAN_COMPILER)
# The compilation mode is obtained from $COMPILATION_MODE
......
......@@ -192,6 +192,79 @@
call set_attributes(ncid,w_id,long_name='vertical vel.',units='m/s', &
FillValue=fv,missing_value=mv,valid_range=vr)
#ifdef _MOMENTUM_TERMS_
err = nf90_def_var(ncid,'tdv_u',NCDF_FLOAT_PRECISION,f4_dims,tdv_u_id)
if (err .NE. NF90_NOERR) go to 10
call set_attributes(ncid,tdv_u_id,long_name='vertical vel.',units='m/s', &
FillValue=fv,missing_value=mv,valid_range=vr)
err = nf90_def_var(ncid,'adv_u',NCDF_FLOAT_PRECISION,f4_dims,adv_u_id)
if (err .NE. NF90_NOERR) go to 10
call set_attributes(ncid,adv_u_id,long_name='vertical vel.',units='m/s', &
FillValue=fv,missing_value=mv,valid_range=vr)
err = nf90_def_var(ncid,'vsd_u',NCDF_FLOAT_PRECISION,f4_dims,vsd_u_id)
if (err .NE. NF90_NOERR) go to 10
call set_attributes(ncid,vsd_u_id,long_name='vertical vel.',units='m/s', &
FillValue=fv,missing_value=mv,valid_range=vr)
err = nf90_def_var(ncid,'hsd_u',NCDF_FLOAT_PRECISION,f4_dims,hsd_u_id)
if (err .NE. NF90_NOERR) go to 10
call set_attributes(ncid,hsd_u_id,long_name='vertical vel.',units='m/s', &
FillValue=fv,missing_value=mv,valid_range=vr)
err = nf90_def_var(ncid,'cor_u',NCDF_FLOAT_PRECISION,f4_dims,cor_u_id)
if (err .NE. NF90_NOERR) go to 10
call set_attributes(ncid,cor_u_id,long_name='vertical vel.',units='m/s', &
FillValue=fv,missing_value=mv,valid_range=vr)
err = nf90_def_var(ncid,'epg_u',NCDF_FLOAT_PRECISION,f4_dims,epg_u_id)
if (err .NE. NF90_NOERR) go to 10
call set_attributes(ncid,epg_u_id,long_name='vertical vel.',units='m/s', &
FillValue=fv,missing_value=mv,valid_range=vr)
err = nf90_def_var(ncid,'ipg_u',NCDF_FLOAT_PRECISION,f4_dims,ipg_u_id)
if (err .NE. NF90_NOERR) go to 10
call set_attributes(ncid,ipg_u_id,long_name='vertical vel.',units='m/s', &
FillValue=fv,missing_value=mv,valid_range=vr)
err = nf90_def_var(ncid,'tdv_v',NCDF_FLOAT_PRECISION,f4_dims,tdv_v_id)
if (err .NE. NF90_NOERR) go to 10
call set_attributes(ncid,tdv_v_id,long_name='vertical vel.',units='m/s', &
FillValue=fv,missing_value=mv,valid_range=vr)
err = nf90_def_var(ncid,'adv_v',NCDF_FLOAT_PRECISION,f4_dims,adv_v_id)
if (err .NE. NF90_NOERR) go to 10
call set_attributes(ncid,adv_v_id,long_name='vertical vel.',units='m/s', &
FillValue=fv,missing_value=mv,valid_range=vr)
err = nf90_def_var(ncid,'vsd_v',NCDF_FLOAT_PRECISION,f4_dims,vsd_v_id)
if (err .NE. NF90_NOERR) go to 10
call set_attributes(ncid,vsd_v_id,long_name='vertical vel.',units='m/s', &
FillValue=fv,missing_value=mv,valid_range=vr)
err = nf90_def_var(ncid,'hsd_v',NCDF_FLOAT_PRECISION,f4_dims,hsd_v_id)
if (err .NE. NF90_NOERR) go to 10
call set_attributes(ncid,hsd_v_id,long_name='vertical vel.',units='m/s', &
FillValue=fv,missing_value=mv,valid_range=vr)
err = nf90_def_var(ncid,'cor_v',NCDF_FLOAT_PRECISION,f4_dims,cor_v_id)
if (err .NE. NF90_NOERR) go to 10
call set_attributes(ncid,cor_v_id,long_name='vertical vel.',units='m/s', &
FillValue=fv,missing_value=mv,valid_range=vr)
err = nf90_def_var(ncid,'epg_v',NCDF_FLOAT_PRECISION,f4_dims,epg_v_id)
if (err .NE. NF90_NOERR) go to 10
call set_attributes(ncid,epg_v_id,long_name='vertical vel.',units='m/s', &
FillValue=fv,missing_value=mv,valid_range=vr)
err = nf90_def_var(ncid,'ipg_v',NCDF_FLOAT_PRECISION,f4_dims,ipg_v_id)
if (err .NE. NF90_NOERR) go to 10
call set_attributes(ncid,ipg_v_id,long_name='vertical vel.',units='m/s', &
FillValue=fv,missing_value=mv,valid_range=vr)
#endif
#if defined(CURVILINEAR)
! rotated zonal velocity
err = nf90_def_var(ncid,'uurot',NCDF_FLOAT_PRECISION,f4_dims,uurot_id)
......
......@@ -24,6 +24,23 @@
integer :: elev_id,u_id,v_id
integer :: taubx_id,tauby_id
integer :: uu_id,vv_id,w_id
#ifdef _MOMENTUM_TERMS_
integer :: tdv_u_id
integer :: adv_u_id
integer :: vsd_u_id
integer :: hsd_u_id
integer :: cor_u_id
integer :: epg_u_id
integer :: ipg_u_id
integer :: tdv_v_id
integer :: adv_v_id
integer :: vsd_v_id
integer :: hsd_v_id
integer :: cor_v_id
integer :: epg_v_id
integer :: ipg_v_id
#endif
#if defined(CURVILINEAR)
integer :: uurot_id,vvrot_id
#endif
......
......@@ -27,6 +27,10 @@
use variables_2d, only: U,V,DU,DV
use variables_3d, only: dt,kmin,ho,hn,uu,hun,vv,hvn,ww,hcc,SS
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
use variables_3d, only: tdv_v,adv_v,vsd_v,hsd_v,cor_v,epg_v,ipg_v
#endif
#ifndef NO_BAROCLINIC
use variables_3d, only: S,T,rho,rad,NN
use variables_3d, only: nummix3d_S,nummix3d_T,phymix3d_S,phymix3d_T
......@@ -208,6 +212,50 @@
err = nf90_put_var(ncid,w_id,ws(_3D_W_),start,edges)
if (err .NE. NF90_NOERR) go to 10
#ifdef _MOMENTUM_TERMS_
err = nf90_put_var(ncid,tdv_u_id,tdv_u(_3D_W_),start,edges)
if (err .NE. NF90_NOERR) go to 10
err = nf90_put_var(ncid,adv_u_id,adv_u(_3D_W_),start,edges)
if (err .NE. NF90_NOERR) go to 10
err = nf90_put_var(ncid,vsd_u_id,vsd_u(_3D_W_),start,edges)
if (err .NE. NF90_NOERR) go to 10
err = nf90_put_var(ncid,hsd_u_id,hsd_u(_3D_W_),start,edges)
if (err .NE. NF90_NOERR) go to 10
err = nf90_put_var(ncid,cor_u_id,cor_u(_3D_W_),start,edges)
if (err .NE. NF90_NOERR) go to 10
err = nf90_put_var(ncid,epg_u_id,epg_u(_3D_W_),start,edges)
if (err .NE. NF90_NOERR) go to 10
err = nf90_put_var(ncid,ipg_u_id,ipg_u(_3D_W_),start,edges)
if (err .NE. NF90_NOERR) go to 10
err = nf90_put_var(ncid,tdv_v_id,tdv_v(_3D_W_),start,edges)
if (err .NE. NF90_NOERR) go to 10
err = nf90_put_var(ncid,adv_v_id,adv_v(_3D_W_),start,edges)
if (err .NE. NF90_NOERR) go to 10
err = nf90_put_var(ncid,vsd_v_id,vsd_v(_3D_W_),start,edges)
if (err .NE. NF90_NOERR) go to 10
err = nf90_put_var(ncid,hsd_v_id,hsd_v(_3D_W_),start,edges)
if (err .NE. NF90_NOERR) go to 10
err = nf90_put_var(ncid,cor_v_id,cor_v(_3D_W_),start,edges)
if (err .NE. NF90_NOERR) go to 10
err = nf90_put_var(ncid,epg_v_id,epg_v(_3D_W_),start,edges)
if (err .NE. NF90_NOERR) go to 10
err = nf90_put_var(ncid,ipg_v_id,ipg_v(_3D_W_),start,edges)
if (err .NE. NF90_NOERR) go to 10
#endif
#if defined(CURVILINEAR)
! rotated zonal and meridional velocities
do j=jmin,jmax
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment