Commit 69e17830 authored by hofmeist's avatar hofmeist
Browse files

added NO_BAROCLINIC functionality, Fixed default real consts

parent db2b1f68
...@@ -630,7 +630,7 @@ ...@@ -630,7 +630,7 @@
<element name="depthmin" type="float" label="minimum layer height" unit="m"/> <element name="depthmin" type="float" label="minimum layer height" unit="m"/>
<element name="cNN" type="float" label="zoom. tendency stratification" minInclusive="0" maxInclusive="1"/> <element name="cNN" type="float" label="zoom. tendency stratification" minInclusive="0" maxInclusive="1"/>
<element name="cSS" type="float" label="zoom. tendency shear" minInclusive="0" maxInclusive="1"/> <element name="cSS" type="float" label="zoom. tendency shear" minInclusive="0" maxInclusive="1"/>
<element name="cNN" type="float" label="zoom. tendency boundaries" minInclusive="0" maxInclusive="1"/> <element name="cdd" type="float" label="zoom. tendency boundaries" minInclusive="0" maxInclusive="1"/>
<element name="d_vel" type="float" unit="m/s" label="vel. difference to resolve"/> <element name="d_vel" type="float" unit="m/s" label="vel. difference to resolve"/>
<element name="d_dens" type="float" unit="kg/m³" label="dens. difference to resolve"/> <element name="d_dens" type="float" unit="kg/m³" label="dens. difference to resolve"/>
<element name="dsurf" type="float" label="reference distance to boundaries" unit="m"/> <element name="dsurf" type="float" label="reference distance to boundaries" unit="m"/>
......
...@@ -57,11 +57,11 @@ endif ...@@ -57,11 +57,11 @@ endif
OBJ += \ OBJ += \
OBJ += \ OBJ += \
$(LIB)(uv_advect_3d.o) \ $(LIB)(uv_advect_3d.o) \
$(LIB)(uv_diffusion_3d.o) $(LIB)(uv_diffusion_3d.o) \
$(LIB)(adaptive_coordinates.o) \
$(LIB)(preadapt_coordinates.o)
ifneq ($(GETM_NO_BAROCLINIC),true) ifneq ($(GETM_NO_BAROCLINIC),true)
OBJ += \ OBJ += \
$(LIB)(adaptive_coordinates.o) \
$(LIB)(preadapt_coordinates.o) \
$(LIB)(ip_blumberg_mellor.o) \ $(LIB)(ip_blumberg_mellor.o) \
$(LIB)(ip_blumberg_mellor_lin.o) \ $(LIB)(ip_blumberg_mellor_lin.o) \
$(LIB)(ip_z_interpol.o) \ $(LIB)(ip_z_interpol.o) \
......
...@@ -60,8 +60,11 @@ ...@@ -60,8 +60,11 @@
! ADAPTIVE-BEGIN ! ADAPTIVE-BEGIN
use parameters, only: g,rho_0 use parameters, only: g,rho_0
use variables_3d, only: uu,vv,NN,SS use variables_3d, only: uu,vv,SS
#ifndef NO_BAROCLINIC
use variables_3d, only: NN
use variables_3d, only: rho use variables_3d, only: rho
#endif
use domain, only: ddu,ddl use domain, only: ddu,ddl
use halo_zones, only: update_3d_halo,wait_halo use halo_zones, only: update_3d_halo,wait_halo
use halo_zones, only: H_TAG,U_TAG,V_TAG use halo_zones, only: H_TAG,U_TAG,V_TAG
...@@ -100,33 +103,33 @@ ...@@ -100,33 +103,33 @@
REALTYPE, save, dimension(:,:,:), allocatable :: zpos ! new pos. of z-levels REALTYPE, save, dimension(:,:,:), allocatable :: zpos ! new pos. of z-levels
REALTYPE, save, dimension(:,:,:), allocatable :: zposo! old pos. of z-levels REALTYPE, save, dimension(:,:,:), allocatable :: zposo! old pos. of z-levels
REALTYPE, save, dimension(:,:,:), allocatable :: work2,work3 REALTYPE, save, dimension(:,:,:), allocatable :: work2,work3
REALTYPE :: faclag=0.0 ! Factor on Lagrangian coords., 0.le.faclag.le.1 REALTYPE :: faclag=_ZERO_ ! Factor on Lagrangian coords., 0.le.faclag.le.1
REALTYPE :: facdif=0.3 ! Factor on thickness filter, 0.le.faclag.le.1 REALTYPE :: facdif=3*_TENTH_ ! Factor on thickness filter, 0.le.faclag.le.1
REALTYPE :: fachor=0.1 ! Factor on position filter, 0.le.faclag.le.1 REALTYPE :: fachor=_TENTH_ ! Factor on position filter, 0.le.faclag.le.1
REALTYPE :: faciso=0.e-6 ! Factor for isopycnal tendency REALTYPE :: faciso=_ZERO_ ! Factor for isopycnal tendency
REALTYPE :: depthmin=0.2 REALTYPE :: depthmin=_ONE_/5
REALTYPE :: Ncrit=1.e-6 REALTYPE :: Ncrit=_ONE_/1000000
integer :: mhor=1 ! this number is experimental - it has to be 1 for now- integer :: mhor=1 ! this number is experimental - it has to be 1 for now-
integer :: iw=2 ! stencil for isopycnal tendency integer :: iw=2 ! stencil for isopycnal tendency
REALTYPE :: rm REALTYPE :: rm
INTEGER :: im,iii,jjj,ii INTEGER :: im,iii,jjj,ii
integer :: split=1 !splits the vertical adaption into #split steps integer :: split=1 ! splits the vertical adaption into #split steps
REALTYPE :: c1ad=0.2 ! dependence on NN REALTYPE :: c1ad=_ONE_/5 ! dependence on NN
REALTYPE :: c2ad=0.0 ! dependence on SS REALTYPE :: c2ad=_ZERO_ ! dependence on SS
REALTYPE :: c3ad=0.2 ! distance from surface and bottom REALTYPE :: c3ad=_ONE_/5 ! distance from surface and bottom
REALTYPE :: c4ad=0.6 ! background value REALTYPE :: c4ad=6*_TENTH_ ! background value
REALTYPE :: d_vel=0.1 ! Typical value of absolute shear REALTYPE :: d_vel=_TENTH_ ! Typical value of absolute shear
REALTYPE :: d_dens=0.5 ! Typical value of BVF squared REALTYPE :: d_dens=_HALF_ ! Typical value of BVF squared
REALTYPE :: dsurf=1. ! reference value for surface/bottom distance REALTYPE :: dsurf=20*_ONE_ ! reference value for surface/bottom distance
REALTYPE :: tgrid=21600. ! Time scale of grid adaptation REALTYPE :: tgrid=21600*_ONE_ ! Time scale of grid adaptation
REALTYPE :: dtgrid REALTYPE :: dtgrid
REALTYPE :: aau(0:kmax),bu(0:kmax) REALTYPE :: aau(0:kmax),bu(0:kmax)
REALTYPE :: cu(0:kmax),du(0:kmax) REALTYPE :: cu(0:kmax),du(0:kmax)
REALTYPE :: facupper=1.0 REALTYPE :: facupper=_ONE_
REALTYPE :: faclower=1.0 REALTYPE :: faclower=_ONE_
REALTYPE :: cNN,cSS,cdd,csum REALTYPE :: cNN,cSS,cdd,csum
REALTYPE :: cbg=0.6 REALTYPE :: cbg=6*_TENTH_
REALTYPE :: tfac_hor=3600.0 ! factor introducing a hor. adaption timescale = dt/tgrid_hor where tgrid_hor=tgrid/N REALTYPE :: tfac_hor=3600*_ONE_ ! factor introducing a hor. adaption timescale = dt/tgrid_hor where tgrid_hor=tgrid/N
integer :: iip integer :: iip
integer,save :: count=0 integer,save :: count=0
...@@ -195,6 +198,10 @@ STDERR 'adaptive_coordinates()' ...@@ -195,6 +198,10 @@ STDERR 'adaptive_coordinates()'
kmaxm1= _ONE_/float(kmax) kmaxm1= _ONE_/float(kmax)
LEVEL2 "allocated memory" LEVEL2 "allocated memory"
! In the case of no-baroclinic, the nnloc will never be updated,
! so we initiialize it here. Same for avn: Both just remain zero.
NNloc(:) = _ZERO_
avn(:) = _ZERO_
if (.not.hotstart) then if (.not.hotstart) then
! Dirty way to read initial distribution (as equidistant sigma coordinates): ! Dirty way to read initial distribution (as equidistant sigma coordinates):
do j=jmin-HALO,jmax+HALO do j=jmin-HALO,jmax+HALO
...@@ -232,13 +239,17 @@ STDERR 'adaptive_coordinates()' ...@@ -232,13 +239,17 @@ STDERR 'adaptive_coordinates()'
faclower=max(_ZERO_,ddl)/(max(_ZERO_,ddu)+max(_ZERO_,ddl)) faclower=max(_ZERO_,ddl)/(max(_ZERO_,ddu)+max(_ZERO_,ddl))
facupper=max(_ZERO_,ddu)/(max(_ZERO_,ddu)+max(_ZERO_,ddl)) facupper=max(_ZERO_,ddu)/(max(_ZERO_,ddu)+max(_ZERO_,ddl))
! norm factors for diffusive adaption ! norm factors for diffusive adaption
#ifdef NO_BAROCLINIC
c1ad=_ZERO_
#else
c1ad=max(_ZERO_,cNN) c1ad=max(_ZERO_,cNN)
#endif
c2ad=max(_ZERO_,cSS) c2ad=max(_ZERO_,cSS)
c3ad=max(_ZERO_,cdd) c3ad=max(_ZERO_,cdd)
if (cbg .lt. 1.0) then if (cbg .lt. _ONE_) then
c4ad=max(_ZERO_,cbg) c4ad=max(_ZERO_,cbg)
else else
c4ad=max(_ZERO_,(1.0-c1ad-c2ad-c3ad)) c4ad=max(_ZERO_,(_ONE_-c1ad-c2ad-c3ad))
end if end if
csum=c1ad+c2ad+c3ad+c4ad csum=c1ad+c2ad+c3ad+c4ad
if (csum .gt. _ONE_) then if (csum .gt. _ONE_) then
...@@ -268,7 +279,7 @@ STDERR 'adaptive_coordinates()' ...@@ -268,7 +279,7 @@ STDERR 'adaptive_coordinates()'
(ho(i+1,j ,k)-ho(i,j,k))*min(1,au(i ,j )) & (ho(i+1,j ,k)-ho(i,j,k))*min(1,au(i ,j )) &
+(ho(i-1,j ,k)-ho(i,j,k))*min(1,au(i-1,j )) & +(ho(i-1,j ,k)-ho(i,j,k))*min(1,au(i-1,j )) &
+(ho(i ,j+1,k)-ho(i,j,k))*min(1,av(i ,j )) & +(ho(i ,j+1,k)-ho(i,j,k))*min(1,av(i ,j )) &
+ (ho(i ,j-1,k)-ho(i,j,k))*min(1,av(i ,j-1)) & +(ho(i ,j-1,k)-ho(i,j,k))*min(1,av(i ,j-1)) &
)*facdif*_QUART_ *tfac_hor )*facdif*_QUART_ *tfac_hor
! Ensure smooth transition to cut-off around depth ! Ensure smooth transition to cut-off around depth
...@@ -306,10 +317,10 @@ STDERR 'adaptive_coordinates()' ...@@ -306,10 +317,10 @@ STDERR 'adaptive_coordinates()'
do j=jmin-HALO,jmax+HALO do j=jmin-HALO,jmax+HALO
do i=imin+1-HALO,imax+HALO do i=imin+1-HALO,imax+HALO
if ((zpos(i,j,k)-zpos(i,j,k-1)).lt.depthmin) then if ((zpos(i,j,k)-zpos(i,j,k-1)).lt.depthmin) then
work2(i ,j,k )=0 work2(i ,j,k )=_ZERO_
work2(i ,j,k-1)=0 work2(i ,j,k-1)=_ZERO_
work2(i-1,j,k )=0 work2(i-1,j,k )=_ZERO_
work2(i-1,j,k-1)=0 work2(i-1,j,k-1)=_ZERO_
endif endif
end do end do
end do end do
...@@ -325,10 +336,10 @@ STDERR 'adaptive_coordinates()' ...@@ -325,10 +336,10 @@ STDERR 'adaptive_coordinates()'
do j=jmin+1-HALO,jmax+HALO do j=jmin+1-HALO,jmax+HALO
do i=imin-HALO,imax+HALO do i=imin-HALO,imax+HALO
if((zpos(i,j,k)-zpos(i,j,k-1)).lt.depthmin) then if((zpos(i,j,k)-zpos(i,j,k-1)).lt.depthmin) then
work3(i,j,k)=0 work3(i,j ,k) = _ZERO_
work3(i,j,k-1)=0 work3(i,j ,k-1) = _ZERO_
work3(i,j-1,k)=0 work3(i,j-1,k) = _ZERO_
work3(i,j-1,k-1)=0 work3(i,j-1,k-1) = _ZERO_
endif endif
end do end do
end do end do
...@@ -338,8 +349,10 @@ STDERR 'adaptive_coordinates()' ...@@ -338,8 +349,10 @@ STDERR 'adaptive_coordinates()'
do k=1,kmax-1 do k=1,kmax-1
do j=jmin-HALO+2,jmax+HALO-2 do j=jmin-HALO+2,jmax+HALO-2
do i=imin-HALO+2,imax+HALO-2 do i=imin-HALO+2,imax+HALO-2
if (faciso .ge. 0.00001) then deltaiso=_ZERO_
rm=0 #ifndef NO_BAROCLINIC
if (faciso .ge. _ONE_/100000) then
rm=_ZERO_
im=0 im=0
do iii=max(imin-HALO,i-iw),min(imax+HALO,i+iw) do iii=max(imin-HALO,i-iw),min(imax+HALO,i+iw)
do jjj=max(jmin-HALO,j-iw),min(jmax+HALO,j+iw) do jjj=max(jmin-HALO,j-iw),min(jmax+HALO,j+iw)
...@@ -355,18 +368,14 @@ STDERR 'adaptive_coordinates()' ...@@ -355,18 +368,14 @@ STDERR 'adaptive_coordinates()'
deltaiso=hn(i,j,k+1)-depthmin deltaiso=hn(i,j,k+1)-depthmin
if (deltaiso.le.(-hn(i,j,k))) & if (deltaiso.le.(-hn(i,j,k))) &
deltaiso=-hn(i,j,k)+depthmin deltaiso=-hn(i,j,k)+depthmin
else end if
deltaiso=_ZERO_ #endif
end if zpos(i,j,k)= zposo(i,j,k) + deltaiso &
+_QUART_*fachor *tfac_hor*(_ZERO_ &
zpos(i,j,k)=zposo(i,j,k)+ & +(zposo(i+1,j ,k)-zposo(i ,j,k))*work2(i ,j,k) &
( & -(zposo(i, j ,k)-zposo(i-1,j,k))*work2(i-1,j,k) &
(zposo(i+1,j,k)-zposo(i,j,k))*work2(i,j,k) - & +(zposo(i, j+1,k)-zposo(i, j,k))*work3(i, j,k) &
(zposo(i,j,k)-zposo(i-1,j,k))*work2(i-1,j,k) & -(zposo(i, j, k)-zposo(i,j-1,k))*work3(i,j-1,k) )
+ (zposo(i,j+1,k)-zposo(i,j,k))*work3(i,j,k) - &
(zposo(i,j,k)-zposo(i,j-1,k))*work3(i,j-1,k) &
)*_QUART_*fachor *tfac_hor &
+ deltaiso
end do end do
end do end do
end do end do
...@@ -397,28 +406,29 @@ STDERR 'adaptive_coordinates()' ...@@ -397,28 +406,29 @@ STDERR 'adaptive_coordinates()'
do j=jmin,jmax do j=jmin,jmax
do i=imin,imax do i=imin,imax
if (az(i,j) .ge. 1) then if (az(i,j) .ge. 1) then
NNloc=NN(i,j,:) #ifndef NO_BAROCLINIC
NNloc=NN(i,j,:)
#endif
SSloc=SS(i,j,:) SSloc=SS(i,j,:)
do k=0,kmax do k=0,kmax
gaa(k)=(zpos(i,j,k)-ssen(i,j))/(ssen(i,j)+H(i,j)) gaa(k) =( zpos(i,j,k)-ssen(i,j))/(ssen(i,j)+H(i,j))
gaaold(k)=(zposo(i,j,k)-ssen(i,j))/(ssen(i,j)+H(i,j)) gaaold(k)=(zposo(i,j,k)-ssen(i,j))/(ssen(i,j)+H(i,j))
end do end do
do ii=1,split do ii=1,split
call col_interpol(kmax-1,1,gaaold,NN(i,j,:),kmax-1,gaa,NNloc)
call col_interpol(kmax-1,1,gaaold,SS(i,j,:),kmax-1,gaa,SSloc)
! Stratification ! Stratification
#ifndef NO_BAROCLINIC
call col_interpol(kmax-1,1,gaaold,NN(i,j,:),kmax-1,gaa,NNloc)
NNloc(kmax)=NNloc(kmax-1) NNloc(kmax)=NNloc(kmax-1)
NNloc(0)=NNloc(1) NNloc(0)=NNloc(1)
SSloc(kmax)=SSloc(kmax-1)
SSloc(0)=SSloc(1)
do k=1,kmax do k=1,kmax
avn(k)=min(_ONE_,max(_ZERO_,_HALF_*(NNloc(k)+NNloc(k-1)))/g & avn(k)=min(_ONE_,max(_ZERO_,_HALF_*(NNloc(k)+NNloc(k-1)))/g &
*rho_0/d_dens) *rho_0/d_dens)
end do end do
#endif
! Shear ! Shear
call col_interpol(kmax-1,1,gaaold,SS(i,j,:),kmax-1,gaa,SSloc)
SSloc(kmax)=SSloc(kmax-1)
SSloc(0)=SSloc(1)
do k=1,kmax do k=1,kmax
avs(k)=min(_ONE_,sqrt(max(_ZERO_,_HALF_* & avs(k)=min(_ONE_,sqrt(max(_ZERO_,_HALF_* &
(SSloc(k)+SSloc(k-1))))/d_vel) (SSloc(k)+SSloc(k-1))))/d_vel)
......
...@@ -120,9 +120,7 @@ STDERR 'coordinates(): hybrid_coordinates not coded yet' ...@@ -120,9 +120,7 @@ STDERR 'coordinates(): hybrid_coordinates not coded yet'
stop stop
case (_ADAPTIVE_COORDS_) ! adaptive vertical coordinates case (_ADAPTIVE_COORDS_) ! adaptive vertical coordinates
LEVEL2 'using adaptive vertical coordinates' LEVEL2 'using adaptive vertical coordinates'
#ifndef NO_BAROCLINIC
call adaptive_coordinates(.true.,hotstart) call adaptive_coordinates(.true.,hotstart)
#endif
case default case default
end select end select
first = .false. first = .false.
...@@ -136,9 +134,7 @@ stop ...@@ -136,9 +134,7 @@ stop
case (_HYBRID_COORDS_) ! hybrid vertical coordinates case (_HYBRID_COORDS_) ! hybrid vertical coordinates
call hybrid_coordinates(.false.) call hybrid_coordinates(.false.)
case (_ADAPTIVE_COORDS_) ! adaptive vertical coordinates case (_ADAPTIVE_COORDS_) ! adaptive vertical coordinates
#ifndef NO_BAROCLINIC
call adaptive_coordinates(.false.,hotstart) call adaptive_coordinates(.false.,hotstart)
#endif
case default case default
end select end select
end if ! first end if ! first
......
...@@ -289,9 +289,9 @@ ...@@ -289,9 +289,9 @@
call init_advection_3d(2) call init_advection_3d(2)
end if end if
end if end if
#endif
if (vert_cord .eq. _ADAPTIVE_COORDS_) call preadapt_coordinates(preadapt) if (vert_cord .eq. _ADAPTIVE_COORDS_) call preadapt_coordinates(preadapt)
#endif
#ifdef DEBUG #ifdef DEBUG
write(debug,*) 'Leaving init_3d()' write(debug,*) 'Leaving init_3d()'
......
...@@ -17,11 +17,13 @@ ...@@ -17,11 +17,13 @@
! !
! !USES: ! !USES:
use getm_timers, only: tic, toc, TIM_COORDS use getm_timers, only: tic, toc, TIM_COORDS
#ifndef NO_BAROCLINIC
use m3d, only: calc_salt, calc_temp use m3d, only: calc_salt, calc_temp
use salinity, only: init_salinity_field, do_salinity use salinity, only: init_salinity_field, do_salinity
use temperature, only: init_temperature_field, do_temperature use temperature, only: init_temperature_field, do_temperature
use eqstate, only: do_eqstate use eqstate, only: do_eqstate
use internal_pressure, only: do_internal_pressure use internal_pressure, only: do_internal_pressure
#endif
use variables_3d, only: SS use variables_3d, only: SS
use domain, only: vert_cord use domain, only: vert_cord
IMPLICIT NONE IMPLICIT NONE
...@@ -39,7 +41,6 @@ ...@@ -39,7 +41,6 @@
!BOC !BOC
call tic(TIM_COORDS) call tic(TIM_COORDS)
#ifndef NO_BAROCLINIC
if (N.ne.0) then if (N.ne.0) then
LEVEL1 'pre-adapting coordinates' LEVEL1 'pre-adapting coordinates'
do ii=1,N do ii=1,N
...@@ -47,14 +48,17 @@ ...@@ -47,14 +48,17 @@
SS=_ZERO_ SS=_ZERO_
call adaptive_coordinates(.false.,.false.) call adaptive_coordinates(.false.,.false.)
call ww_momentum_3d() call ww_momentum_3d()
#ifndef NO_BAROCLINIC
if(calc_salt) call do_salinity(1) if(calc_salt) call do_salinity(1)
if(calc_temp) call do_temperature(1) if(calc_temp) call do_temperature(1)
call do_eqstate() call do_eqstate()
#endif
call ss_nn() call ss_nn()
call stop_macro() call stop_macro()
if (mod(ii,10).eq._ZERO_) LEVEL3 ii if (mod(ii,10).eq._ZERO_) LEVEL3 ii
end do end do
#ifndef NO_BAROCLINIC
LEVEL2 'reinterpolating initial salinity' LEVEL2 'reinterpolating initial salinity'
if(calc_salt) then if(calc_salt) then
call init_salinity_field() call init_salinity_field()
...@@ -65,8 +69,8 @@ ...@@ -65,8 +69,8 @@
end if end if
call do_eqstate() call do_eqstate()
call do_internal_pressure() !assuming to run always in runtype>2 call do_internal_pressure() !assuming to run always in runtype>2
end if
#endif #endif
end if
call toc(TIM_COORDS) call toc(TIM_COORDS)
#ifdef DEBUG #ifdef DEBUG
write(debug,*) 'Leaving preadapt_coordinates()' write(debug,*) 'Leaving preadapt_coordinates()'
......
...@@ -374,22 +374,18 @@ ...@@ -374,22 +374,18 @@
call read_topo_file(bathy_format,bathymetry) call read_topo_file(bathy_format,bathymetry)
select case (vert_cord) select case (vert_cord)
case(1) case(_SIGMA_COORDS_)
LEVEL2 'Using sigma coordinates' LEVEL2 'Using sigma coordinates'
case(2) case(_Z_COORDS_)
LEVEL2 'Using z-level coordinates' LEVEL2 'Using z-level coordinates'
case(3) case(_GENERAL_COORDS_)
LEVEL2 'Using general vertical coordinates' LEVEL2 'Using general vertical coordinates'
case (4) ! hybrid vertical coordinates case (_HYBRID_COORDS_) ! hybrid vertical coordinates
LEVEL2 'using hybrid vertical coordinates' LEVEL2 'using hybrid vertical coordinates'
STDERR 'domain: hybrid_coordinates not coded yet' STDERR 'domain: hybrid_coordinates not coded yet'
stop stop
case (_ADAPTIVE_COORDS_) ! adaptive vertical coordinates case (_ADAPTIVE_COORDS_) ! adaptive vertical coordinates
LEVEL2 'using adaptive vertical coordinates' LEVEL2 'using adaptive vertical coordinates'
#ifdef NO_BAROCLINIC
LEVEL2 'adaptive coordinates are not working with NO_BAROCLINIC'
stop
#endif
case default case default
call getm_error("init_domain()", & call getm_error("init_domain()", &
"A non valid vertical coordinate system has been chosen"); "A non valid vertical coordinate system has been chosen");
......
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