Commit c83fdb44 authored by hb's avatar hb
Browse files

Source code documentation extended

parent 29307f62
#$Id: Makefile,v 1.9 2006-02-07 07:16:22 kbk Exp $ #$Id: Makefile,v 1.10 2006-02-10 22:41:56 hb Exp $
# #
# Makefile to build the 3D specific library - libm3d.a # Makefile to build the 3D specific library - libm3d.a
# #
...@@ -13,6 +13,8 @@ MODSRC = m3d.F90 variables_3d.F90 advection_3d.F90 eqstate.F90 \ ...@@ -13,6 +13,8 @@ MODSRC = m3d.F90 variables_3d.F90 advection_3d.F90 eqstate.F90 \
LIBSRC = start_macro.F90 hcc_check.F90 bdy_3d.F90 coordinates.F90 check_h.F90 bottom_friction_3d.F90 internal_pressure.F90 uu_momentum_3d.F90 vv_momentum_3d.F90 ww_momentum_3d.F90 uv_advect_3d.F90 uv_diffusion_3d.F90 stresses_3d.F90 ss_nn.F90 gotm.F90 slow_bottom_friction.F90 slow_advection.F90 slow_diffusion.F90 slow_terms.F90 stop_macro.F90 LIBSRC = start_macro.F90 hcc_check.F90 bdy_3d.F90 coordinates.F90 check_h.F90 bottom_friction_3d.F90 internal_pressure.F90 uu_momentum_3d.F90 vv_momentum_3d.F90 ww_momentum_3d.F90 uv_advect_3d.F90 uv_diffusion_3d.F90 stresses_3d.F90 ss_nn.F90 gotm.F90 slow_bottom_friction.F90 slow_advection.F90 slow_diffusion.F90 slow_terms.F90 stop_macro.F90
TEXSRC = m3d.F90 variables_3d.F90 coordinates.F90 check_h.F90 hcc_check.F90 uu_momentum_3d.F90 vv_momentum_3d.F90 ww_momentum_3d.F90 uv_advect_3d.F90 uv_diffusion_3d.F90 bottom_friction_3d.F90 internal_pressure.F90 ip_blumberg_mellor.F90 ip_blumberg_mellor_lin.F90 ip_z_interpol.F90 ip_song_wright.F90 ip_chu_fan.F90 advection_3d.F90 temperature.F90 salinity.F90 spm.F90 eqstate.F90 gotm.F90 ss_nn.F90 stresses_3d.F90 rivers.F90 bdy_3d.F90 slow_bottom_friction.F90 slow_advection.F90 slow_diffusion.F90 slow_terms.F90 start_macro.F90 stop_macro.F90
SRC = $(MODSRC) $(LIBSRC) SRC = $(MODSRC) $(LIBSRC)
MOD = \ MOD = \
...@@ -86,15 +88,15 @@ $(MOD): $(INCS) ...@@ -86,15 +88,15 @@ $(MOD): $(INCS)
$(OBJ): $(MOD) $(OBJ): $(MOD)
doc: $(SRC) doc: $(TEXSRC)
$(PROTEX) $(SRC) > $(DOCDIR)/3d.tex $(PROTEX) $(TEXSRC) > $(DOCDIR)/3d.tex
touch doc touch doc
clean: clean:
$(RM) $(LIB) $(MODDIR)/m3d.{m.mod} $(RM) $(LIB) $(MODDIR)/m3d.{m.mod}
realclean: clean realclean: clean
$(RM) *.o $(RM) *.o doc
distclean: realclean distclean: realclean
......
!$Id: bdy_3d.F90,v 1.7 2005-08-26 08:19:31 kbk Exp $ !$Id: bdy_3d.F90,v 1.8 2006-02-10 22:41:56 hb Exp $
#include "cppdefs.h" #include "cppdefs.h"
!----------------------------------------------------------------------- !-----------------------------------------------------------------------
!BOP !BOP
! !
! !MODULE: bdy_3d - boundary conditions 3D. ! !MODULE: bdy_3d - 3D boundary conditions \label{bdy-3d}
! !
! !INTERFACE: ! !INTERFACE:
module bdy_3d module bdy_3d
...@@ -32,6 +32,9 @@ ...@@ -32,6 +32,9 @@
! Original author(s): Karsten Bolding & Hans Burchard ! Original author(s): Karsten Bolding & Hans Burchard
! !
! $Log: bdy_3d.F90,v $ ! $Log: bdy_3d.F90,v $
! Revision 1.8 2006-02-10 22:41:56 hb
! Source code documentation extended
!
! Revision 1.7 2005-08-26 08:19:31 kbk ! Revision 1.7 2005-08-26 08:19:31 kbk
! zero gradient open boundaries for biological variables ! zero gradient open boundaries for biological variables
! !
...@@ -65,7 +68,8 @@ ...@@ -65,7 +68,8 @@
!----------------------------------------------------------------------- !-----------------------------------------------------------------------
!BOP !BOP
! !
! !IROUTINE: init_bdy_3d ! !IROUTINE: init_bdy_3d - initialising 3D boundary conditions
! \label{sec-init-bdy-3d}
! !
! !INTERFACE: ! !INTERFACE:
subroutine init_bdy_3d() subroutine init_bdy_3d()
...@@ -114,7 +118,8 @@ ...@@ -114,7 +118,8 @@
!----------------------------------------------------------------------- !-----------------------------------------------------------------------
!BOP !BOP
! !
! !IROUTINE: do_bdy_3d() ! !IROUTINE: do_bdy_3d() - updating 3D boundary conditions
! \label{sec-do-bdy-3d}
! !
! !INTERFACE: ! !INTERFACE:
subroutine do_bdy_3d(tag,field) subroutine do_bdy_3d(tag,field)
......
!$Id: bottom_friction_3d.F90,v 1.6 2006-02-04 11:47:26 hb Exp $ !$Id: bottom_friction_3d.F90,v 1.7 2006-02-10 22:41:56 hb Exp $
#include "cppdefs.h" #include "cppdefs.h"
!----------------------------------------------------------------------- !-----------------------------------------------------------------------
!BOP !BOP
! !
! !ROUTINE: bottom_friction_3d() - 3D-bottom friction ! !ROUTINE: bottom_friction_3d - bottom friction
! \label{sec-bottom-friction-3d} ! \label{sec-bottom-friction-3d}
! !
! !INTERFACE: ! !INTERFACE:
...@@ -11,6 +11,39 @@ ...@@ -11,6 +11,39 @@
! !
! !DESCRIPTION: ! !DESCRIPTION:
! !
! Based on the assumption that the velocity distribution in the bottom
! layer is logarithmic,
! the product of the drag coefficient with the
! absolute value of the current speed in the bottom layer,
!
! \begin{equation}
! r \sqrt{u_b^2+v_b^2}
! \end{equation}
!
! with the velocity components of the bottom layer, $u_b$ and $v_b$,
! and the drag coefficient
!
! \begin{equation}\label{r}
! r = \left(\frac{\kappa}{\ln \left(\frac{0.5h_1+z_0^b}{z_0^b}\right)}
! \right)^2,
! \end{equation}
!
! is calculated and
! provided as output parameters {\tt rru} (for U-points) and
! {\tt rrv} (for V-points). The layer height $h_1$ in (\ref{r}) is set to
! the thickness of the bottom layer in the respective U- or V-point.
!
! There are some experimental options for the interested user included
! here. It is possible to change the interpolation of $u$ to V-points
! and of $v$ to U-points from velocity-based interpolation (as done
! presently) to transport-based averaging (commented out). Furthermore,
! the user may activate some outcommented lines which allow the
! consideration of flow-depending bottom roughness length $z_0^b$
! according to (\ref{Defz0b}), see page \pageref{Defz0b}.
!
! For a derivation of (\ref{r}), see section \ref{SectionBedFric} on
! page \pageref{SectionBedFric}.
!
! !USES: ! !USES:
use parameters, only: kappa,avmmol use parameters, only: kappa,avmmol
use domain, only: iimin,iimax,jjmin,jjmax,kmax,au,av,min_depth use domain, only: iimin,iimax,jjmin,jjmax,kmax,au,av,min_depth
...@@ -28,6 +61,9 @@ ...@@ -28,6 +61,9 @@
! Original author(s): Hans Burchard & Karsten Bolding ! Original author(s): Hans Burchard & Karsten Bolding
! !
! $Log: bottom_friction_3d.F90,v $ ! $Log: bottom_friction_3d.F90,v $
! Revision 1.7 2006-02-10 22:41:56 hb
! Source code documentation extended
!
! Revision 1.6 2006-02-04 11:47:26 hb ! Revision 1.6 2006-02-04 11:47:26 hb
! Source code documentation extended ! Source code documentation extended
! !
......
!$Id: coordinates.F90,v 1.8 2005-10-06 09:54:01 hb Exp $ !$Id: coordinates.F90,v 1.9 2006-02-10 22:41:56 hb Exp $
#include "cppdefs.h" #include "cppdefs.h"
!----------------------------------------------------------------------- !-----------------------------------------------------------------------
!BOP !BOP
! !
! !ROUTINE: coordinates() - defines the vertical coordinate system. ! !ROUTINE: coordinates - defines the vertical coordinate
! \label{sec-coordinates}
! !
! !INTERFACE: ! !INTERFACE:
subroutine coordinates(cord_type,cord_relax,maxdepth) subroutine coordinates(cord_type,cord_relax,maxdepth)
! !
! !DESCRIPTION: ! !DESCRIPTION:
! !
! Here, the vertical layer distribution in T-, U- and V-points is updated
! during every macro time step. This is done for the old and the new
! layer thicknesses at every point. Calculation of the layer distribution
! in the U- and V-points is done indepently from the calculation in the
! T-points, since different methods for the calculation of the
! bathymetry values in the U- and V-points are possible, see routine
! {\tt uv\_depths} described on page \pageref{sec-uv-depth}.
!
! Here, three different methods for the vertical layer distribution
! are coded:
!
! \begin{enumerate}
! \item Classical $\sigma$ coordinates where layer interfaces for each
! layer index have a fixed relative position $\sigma_k$ in the water column,
! which may be even equidistant or non-equidistant, see equations
! (\ref{sigma}) and (\ref{formula_Antoine}).
! The surface and bottom zooming factors
! $d_u$ and $d_l$ are read in via the {\tt domain} namelist in {\tt getm.inp}
! as {\tt ddu} and {\tt ddl}.
! In the first call to coordinates, the relative interface positions
! {\tt dga} are calculated as a one-dimensional vector (in case of
! non-equidistant $\sigma$ coordinates), and those are then multiplied with
! the water depths in all T-, U- and V-points to get the layer thicknesses.
! \item Also $z$- (i.e.\ geopotential) coordinates are enabled in GETM
! in principle. However, they may not yet work and need further
! development. First of all, fixed $z$-levels are defined by means of
! zooming factors and the maximum water depth $H_{\max}$:
!
! \begin{equation}\label{formula_Antoine_zlevels}
! z_k = H_{\max}\left(\frac{\mbox{tanh}\left( (d_l+d_u)(1+\sigma_k)-d_l\right)
! +\mbox{tanh}(d_l)}{\mbox{tanh}(d_l)+\mbox{tanh}(d_u)}-1\right),
! \qquad k=0,\dots,N\qquad
! \end{equation}
!
! Then, layers are from the surface down filled into the T-point
! water column locally.
! When the last layer is shallower than {\tt hnmin} (hard coded as local
! variable), the two last layers are combined. The index of the lowest
! layer is then stored in the integer field {\tt kmin\_pmz}.
! layer thicknesses in U- and V-points are then taken as the minimum
! values of adjacent thicknesses in T-points, and bottom indices
! {\tt kumin\_pmz} and {\tt kvmin\_pmz} are taken as the maximum
! of adjacent {\tt kmin\_pmz} indices.
! \item The third and so far most powerful method are the genral
! vertical coordinates, discussed in section \ref{SectionGeneralCoordinates},
! see equations (\ref{sigma}) - (\ref{MLDTransform}), which is basically
! an interpolation between equidistant and non-equaidistant $\sigma$
! coordinates. During the first call, a three-dimensional field
! {\tt gga} containing the relative interface positions is calculated,
! which further down used together with the actual water depth in the
! T-, U- and V-points for calculating the updated old and new layer
! thicknesses.
!\end{enumerate}
!
! A fourth option will soon be the adaptive grids which have been
! conceptionally developed by \cite{BURCHARDea04}.
!
! !USES: ! !USES:
use halo_zones, only: update_3d_halo,wait_halo,H_TAG,HU_TAG,HV_TAG use halo_zones, only: update_3d_halo,wait_halo,H_TAG,HU_TAG,HV_TAG
use domain, only: iimin,iimax,jjmin,jjmax,kmax,H,HU,HV,az,au,av,min_depth use domain, only: iimin,iimax,jjmin,jjmax,kmax,H,HU,HV,az,au,av,min_depth
...@@ -32,6 +90,9 @@ ...@@ -32,6 +90,9 @@
! Original author(s): Hans Burchard & Karsten Bolding ! Original author(s): Hans Burchard & Karsten Bolding
! !
! $Log: coordinates.F90,v $ ! $Log: coordinates.F90,v $
! Revision 1.9 2006-02-10 22:41:56 hb
! Source code documentation extended
!
! Revision 1.8 2005-10-06 09:54:01 hb ! Revision 1.8 2005-10-06 09:54:01 hb
! added support for vertical slice model - via -DSLICE_MODEL ! added support for vertical slice model - via -DSLICE_MODEL
! !
......
!$Id: eqstate.F90,v 1.4 2004-01-02 13:54:24 kbk Exp $ !$Id: eqstate.F90,v 1.5 2006-02-10 22:41:56 hb Exp $
#include "cppdefs.h" #include "cppdefs.h"
!----------------------------------------------------------------------- !-----------------------------------------------------------------------
!BOP !BOP
! !
! !MODULE: eqstate ! !MODULE: eqstate \label{sec-eqstate}
! !
! !INTERFACE: ! !INTERFACE:
module eqstate module eqstate
...@@ -26,6 +26,9 @@ ...@@ -26,6 +26,9 @@
! Original author(s): Karsten Bolding & Hans Burchard ! Original author(s): Karsten Bolding & Hans Burchard
! !
! $Log: eqstate.F90,v $ ! $Log: eqstate.F90,v $
! Revision 1.5 2006-02-10 22:41:56 hb
! Source code documentation extended
!
! Revision 1.4 2004-01-02 13:54:24 kbk ! Revision 1.4 2004-01-02 13:54:24 kbk
! read equation of state info from namelist - Ruiz ! read equation of state info from namelist - Ruiz
! !
...@@ -134,7 +137,7 @@ ...@@ -134,7 +137,7 @@
!----------------------------------------------------------------------- !-----------------------------------------------------------------------
!BOP !BOP
! !
! !IROUTINE: do_eqstate() ! !IROUTINE: do_eqstate() - equation of state \label{sec-do-eqstate}
! !
! !INTERFACE: ! !INTERFACE:
subroutine do_eqstate() subroutine do_eqstate()
......
!$Id: gotm.F90,v 1.10 2006-01-09 10:35:40 lars Exp $ !$Id: gotm.F90,v 1.11 2006-02-10 22:41:56 hb Exp $
#include "cppdefs.h" #include "cppdefs.h"
!----------------------------------------------------------------------- !-----------------------------------------------------------------------
!BOP !BOP
! !
! !IROUTINE: gotm() - a wrapper to call GOTM. ! !IROUTINE: gotm() - a wrapper to call GOTM \label{sec-gotm}
! !
! !INTERFACE: ! !INTERFACE:
subroutine gotm() subroutine gotm()
...@@ -34,6 +34,9 @@ ...@@ -34,6 +34,9 @@
! Original author(s): Karsten Bolding & Hans Burchard ! Original author(s): Karsten Bolding & Hans Burchard
! !
! $Log: gotm.F90,v $ ! $Log: gotm.F90,v $
! Revision 1.11 2006-02-10 22:41:56 hb
! Source code documentation extended
!
! Revision 1.10 2006-01-09 10:35:40 lars ! Revision 1.10 2006-01-09 10:35:40 lars
! bug fix in call to do_turbulence() ! bug fix in call to do_turbulence()
! !
......
!$Id: hcc_check.F90,v 1.1 2004-05-04 09:23:51 kbk Exp $ !$Id: hcc_check.F90,v 1.2 2006-02-10 22:41:56 hb Exp $
#include "cppdefs.h" #include "cppdefs.h"
!----------------------------------------------------------------------- !-----------------------------------------------------------------------
!BOP !BOP
! !
! !ROUTINE: hcc_check() - hydrostatic consistency criteria ! !ROUTINE: hcc_check - hydrostatic consistency criteria
! !
! !INTERFACE: ! !INTERFACE:
subroutine hcc_check() subroutine hcc_check()
! !
! !DESCRIPTION: ! !DESCRIPTION:
! !
! This diagnostic routine calculates the hydrostatic consistency $h^c$
! in each T-point and each layer. $h^c$ is defined as:
!
! \begin{equation}\label{HCC}
! h^c_{i,j,k}=\max\left\{
! |\partial_xz_k| \frac{\Delta x}{\frac12(h_{i,j,k}+h_{i+1,j,k})},
! |\partial_yz_k| \frac{\Delta y}{\frac12(h_{i,j,k}+h_{i,j+1,k})}
! \right\}.
! \end{equation}
!
! For the numerical calculation it is used here that $\Delta x$ and
! $\Delta y$ can be cancelled out each.
! For $h^c \leq 1$, the grid box is hydrostatically consistent,
! else it is called hydrostatically inconsistent. In the latter case,
! numerical problems can be expected for terrain-following coordinates
! when stratification is strong.
!
! $h^c$ is stored in the 3d netcdf output file.
!
! !USES: ! !USES:
use domain, only: imin,imax,jmin,jmax,az,au,av,HU,HV use domain, only: imin,imax,jmin,jmax,az,au,av,HU,HV
use domain, only: iimin,iimax,jjmin,jjmax,kmax use domain, only: iimin,iimax,jjmin,jjmax,kmax
...@@ -20,6 +39,9 @@ ...@@ -20,6 +39,9 @@
! Original author(s): Karsten Bolding & Hans Burchard ! Original author(s): Karsten Bolding & Hans Burchard
! !
! $Log: hcc_check.F90,v $ ! $Log: hcc_check.F90,v $
! Revision 1.2 2006-02-10 22:41:56 hb
! Source code documentation extended
!
! Revision 1.1 2004-05-04 09:23:51 kbk ! Revision 1.1 2004-05-04 09:23:51 kbk
! hydrostatic consistency criteria stored in .3d.nc file ! hydrostatic consistency criteria stored in .3d.nc file
! !
...@@ -36,10 +58,10 @@ ...@@ -36,10 +58,10 @@
do j=jjmin,jjmax do j=jjmin,jjmax
do i=iimin,iimax do i=iimin,iimax
if (az(i,j) .ge. 1) then if (az(i,j) .ge. 1) then
du1=HU(i-1,j) du1=-HU(i-1,j)
du2=HU(i,j) du2=-HU(i,j)
dv1=HV(i,j-1) dv1=-HV(i,j-1)
dv2=HV(i,j) dv2=-HV(i,j)
do k=1,kmax do k=1,kmax
......
!$Id: internal_pressure.F90,v 1.13 2006-01-29 20:32:33 hb Exp $ !$Id: internal_pressure.F90,v 1.14 2006-02-10 22:41:56 hb Exp $
#include "cppdefs.h" #include "cppdefs.h"
!----------------------------------------------------------------------- !-----------------------------------------------------------------------
!BOP !BOP
...@@ -10,13 +10,51 @@ ...@@ -10,13 +10,51 @@
! !
! !DESCRIPTION: ! !DESCRIPTION:
! !
! \begin{verbatim} ! In GETM, various methods are provided for the calculation of the
! ip_method=1 ! old GETM version of Blumberg-Mellor, very stable ! internal pressure gradients terms in $x$- and $y$-direction.
! ip_method=2 ! new GETM version of 1, zero for linear profiles ! These terms which appear as layer-integrated terms in the
! ip_method=3 ! z-coordinates, linear interpolation ! equations for the layer-integrated momentum are for the
! ip_method=4 ! Song and Wright 1998 version ! eastward momentum $p_k$ (see equation (\ref{uEqvi})):
! ip_method=5 ! CHU and FAN 2003 version !
! \end{verbatim} ! \begin{equation}
! h_k\left(\frac12h_N(\partial^*_xb)_N
! +\sum_{j=k}^{N-1}\frac12(h_j+h_{j+1})(\partial^*_xb)_j
! \right)
! \end{equation}
!
! and for the northward layer-integrated momentum $q_k$
! (see equation (\ref{vEqvi})):
!
! \begin{equation}
! h_k\left(\frac12h_N(\partial^*_yb)_N
! +\sum_{j=k}^{N-1}\frac12(h_j+h_{j+1})(\partial^*_yb)_j
! \right)
! \end{equation}
!
! The major problem is how to calculate the horizontal (with respect
! to isogeopotentials) buoyancy gradients $\partial^*_xb$ and $\partial^*_yb$,
! which need to be defined at the interfaces positioned vertically
! between two velocity points.
!
! The methods for calculating the internal pressure gradient included in
! GETM are currently:
!
! \begin{enumerate}
! \item Method by \cite{MELLORea94}, see routine {\tt ip\_blumberg\_mellor}
! \item Modified \cite{MELLORea94} method, exact for linear density profiles
! with $z$-dependence only, see routine {\tt ip\_blumberg\_mellor\_lin}
! \item Calculation by mean of linear interpolation to $z$-levels, see routine
! {\tt ip\_z\_interpol}
! \item Method by \cite{SONG98}, see routine {\tt ip\_song\_wright}
! \item Method by \cite{CHUea03}, see routine {\tt ip\_chu\_fan}
! \end{enumerate}
!
! It is possible, by setting the compiler option {\tt SUBSTR\_INI\_PRESS},
! to substract the initial pressure gradient from all pressure
! gradients. This is only advisable for strong stratification
! without any initial internal pressure gradients. In this case
! any non-zero values of the resulting numerical initial pressure gradient
! are due to discretisation errors.
! !
! !USES: ! !USES:
use exceptions use exceptions
...@@ -55,6 +93,9 @@ ...@@ -55,6 +93,9 @@
! Original author(s): Hans Burchard & Karsten Bolding ! Original author(s): Hans Burchard & Karsten Bolding
! !
! $Log: internal_pressure.F90,v $ ! $Log: internal_pressure.F90,v $
! Revision 1.14 2006-02-10 22:41:56 hb
! Source code documentation extended
!
! Revision 1.13 2006-01-29 20:32:33 hb ! Revision 1.13 2006-01-29 20:32:33 hb
! Small LaTeX corrections to source code documentation ! Small LaTeX corrections to source code documentation
! !
...@@ -142,15 +183,18 @@ ...@@ -142,15 +183,18 @@
!----------------------------------------------------------------------- !-----------------------------------------------------------------------
!BOP !BOP
! !
! !IROUTINE: init_internal_pressure ! !IROUTINE: init_internal_pressure - initialising internal pressure gradient
! \label{sec-init-internal pressure}
! !
! !INTERFACE: ! !INTERFACE:
subroutine init_internal_pressure() subroutine init_internal_pressure()
IMPLICIT NONE IMPLICIT NONE
! !
! !DESCRIPTION: ! !DESCRIPTION:
! Reads the namelist and makes calls to the init functions of the !
! various model components. ! Here, some necessary memory is allocated (in case of the compiler option
! {\tt STATIC}), and information is written to the log-file of
! the simulation.
! !
! !REVISION HISTORY: ! !REVISION HISTORY:
! See the log for the module ! See the log for the module
...@@ -197,16 +241,25 @@ ...@@ -197,16 +241,25 @@
return return
end subroutine init_internal_pressure end subroutine init_internal_pressure
!EOC !EOC
!----------------------------------------------------------------------- !-----------------------------------------------------------------------
!BOP !BOP
! !
! !IROUTINE: do_internal_pressure() ! !IROUTINE: do_internal_pressure - internal pressure gradient
! \label{sec-do-internal-pressure}
! !
! !INTERFACE: ! !INTERFACE:
subroutine do_internal_pressure() subroutine do_internal_pressure()
! !
! !DESCRIPTION: ! !DESCRIPTION:
!
! Here, the chosen internal pressure gradient method is selected
! and (in case that the compiler option {\tt SUBSTR\_INI\_PRESS} is
! set), the initial pressure is calculated and subtracted from the
! updated internal pressure gradient.
!
! If GETM is executed as slice model (compiler option {\tt SLICE\_MODEL}
! is set, the internal pressure gradient for $j=2$ is copied to
! $j=3$.
! !
! !USES: ! !USES:
IMPLICIT NONE IMPLICIT NONE
...@@ -223,6 +276,7 @@ ...@@ -223,6 +276,7 @@
! !LOCAL VARIABLES: ! !LOCAL VARIABLES:
integer :: i,j,k integer :: i,j,k
logical, save :: first=.true. logical, save :: first=.true.
!EOP
!----------------------------------------------------------------------- !-----------------------------------------------------------------------
!BOC !BOC
#ifdef DEBUG #ifdef DEBUG
...@@ -251,14 +305,6 @@ ...@@ -251,14 +305,6 @@
stop 'do_internal_pressure()' stop 'do_internal_pressure()'
end select end select
#ifdef SLICE_MODEL
do i=iimin,iimax
do k=kmin(i,2),kmax
idpdx(i,3,k)=idpdx(i,2,k)
end do
end do
#endif
#ifdef SUBSTR_INI_PRESS #ifdef SUBSTR_INI_PRESS
if (first) then if (first) then
first = .false. first = .false.
...@@ -284,6 +330,14 @@ ...@@ -284,6 +330,14 @@
end if end if
#endif #endif