coordinates.F90 17.4 KB
Newer Older
kbk's avatar
kbk committed
1
!$Id: coordinates.F90,v 1.11 2006-03-01 15:54:08 kbk Exp $
gotm's avatar
gotm committed
2 3 4 5
#include "cppdefs.h"
!-----------------------------------------------------------------------
!BOP
!
hb's avatar
hb committed
6 7
! !ROUTINE:  coordinates - defines the vertical coordinate
! \label{sec-coordinates}
gotm's avatar
gotm committed
8 9
!
! !INTERFACE:
kbk's avatar
kbk committed
10
   subroutine coordinates(cord_type,cord_relax,maxdepth)
gotm's avatar
gotm committed
11 12 13
!
! !DESCRIPTION:
!
hb's avatar
hb committed
14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70
! 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}.
!
gotm's avatar
gotm committed
71
! !USES:
kbk's avatar
kbk committed
72
   use halo_zones, only: update_3d_halo,wait_halo,H_TAG,HU_TAG,HV_TAG
kbk's avatar
kbk committed
73 74
   use domain, only: iimin,iimax,jjmin,jjmax,kmax,H,HU,HV,az,au,av,min_depth
   use domain, only: ga,ddu,ddl,d_gamma,gamma_surf
gotm's avatar
gotm committed
75
   use variables_3d, only: dt,kmin,kumin,kvmin,ho,hn,huo,hun,hvo,hvn
kbk's avatar
kbk committed
76
   use variables_3d, only: kmin_pmz,kumin_pmz,kvmin_pmz
gotm's avatar
gotm committed
77 78 79 80
   use variables_3d, only: sseo,ssen,ssuo,ssun,ssvo,ssvn
   IMPLICIT NONE
!
! !INPUT PARAMETERS:
kbk's avatar
kbk committed
81 82
   integer, intent(in)                 :: cord_type
   REALTYPE, intent(in)                :: cord_relax
kbk's avatar
kbk committed
83
   REALTYPE, intent(in)                :: maxdepth
gotm's avatar
gotm committed
84 85 86 87 88
!
! !INPUT/OUTPUT PARAMETERS:
!
! !OUTPUT PARAMETERS:
!
kbk's avatar
kbk committed
89 90 91
! !REVISION HISTORY:
!  Original author(s): Hans Burchard & Karsten Bolding
!
gotm's avatar
gotm committed
92
! !LOCAL VARIABLES:
kbk's avatar
kbk committed
93 94 95 96 97 98
   integer         :: i,j,k,rc,kk
   REALTYPE        :: tmp,kmaxm1,alpha
   REALTYPE        :: HH,zz,r,hnmin=0.01
   logical, save   :: first=.true.,equiv_sigma=.false.
   logical         :: kminset
   REALTYPE, save, dimension(:),     allocatable  :: dga,be,sig,zlev
kbk's avatar
kbk committed
99
   REALTYPE, save, dimension(:,:,:), allocatable  :: gga
gotm's avatar
gotm committed
100 101 102 103 104 105 106 107 108 109 110
!EOP
!-----------------------------------------------------------------------
!BOC
#ifdef DEBUG
   integer, save :: Ncall = 0
   Ncall = Ncall+1
   write(debug,*) 'coordinates() # ',Ncall
#endif

   if (first) then
      first = .false.
kbk's avatar
kbk committed
111
      if (.not. allocated(ga)) allocate(ga(0:kmax),stat=rc)
gotm's avatar
gotm committed
112 113 114 115 116 117 118 119 120 121 122 123 124
      if (rc /= 0) stop 'coordinates: Error allocating (ga)'
      select case (cord_type)
         case (1) ! sigma coordinates
            if (ddu .le. _ZERO_ .and. ddl .le. _ZERO_) then
               equiv_sigma=.true.
               ga(0) = -_ONE_
               do k=1,kmax
                  ga(k) = ga(k-1) + _ONE_/kmax
               end do
               ga(kmax) = _ZERO_
            else
               ! Non-equidistant sigma coordinates
               ! This zooming routine is from Antoine Garapon, ICCH, DK
kbk's avatar
kbk committed
125 126
               if (ddu .lt. _ZERO_) ddu=_ZERO_
               if (ddl .lt. _ZERO_) ddl=_ZERO_
gotm's avatar
gotm committed
127 128 129 130 131 132 133 134 135 136 137 138 139
               allocate(dga(0:kmax),stat=rc)
               if (rc /= 0) STOP 'coordinates: Error allocating (dga)'
               ga(0)= -_ONE_
               dga(0)= _ZERO_
               do k=1,kmax
                  ga(k)=tanh((ddl+ddu)*k/float(kmax)-ddl)+tanh(ddl)
                  ga(k)=ga(k)/(tanh(ddl)+tanh(ddu)) - _ONE_
                  dga(k)=ga(k)-ga(k-1)
               end do
            end if
            kmin=1
            kumin=1
            kvmin=1
kbk's avatar
kbk committed
140 141 142
            kmin_pmz=1
            kumin_pmz=1
            kvmin_pmz=1
gotm's avatar
gotm committed
143
         case (2) ! z-level
kbk's avatar
kbk committed
144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211
            allocate(zlev(0:kmax),stat=rc)
            if (rc /= 0) stop 'coordinates: Error allocating (zlev)'
            allocate(gga(I3DFIELD),stat=rc)  ! dimensionless gamma-coordinate
            if (rc /= 0) stop 'coordinates: Error allocating memory (gga)'
            ! Calculate profile of zlevels
            zlev(0)=-maxdepth
            do k=1,kmax
               zlev(k)=tanh((ddl+ddu)*k/float(kmax)-ddl)+tanh(ddl)
               zlev(k)=zlev(k)/(tanh(ddl)+tanh(ddu)) - _ONE_
               zlev(k)=zlev(k)*maxdepth
            end do
            kmin_pmz=1
            do j=jjmin,jjmax
               do i=iimin,iimax
                 HH=max(sseo(i,j)+H(i,j),min_depth)
                 kminset=.true.
                 if (az(i,j).eq.0) then
                    do k=1,kmax
                       hn(i,j,k)=HH/float(kmax)
                       gga(i,j,k)=hn(i,j,k)/HH
                    end do
                    kmin_pmz(i,j)=1
                  else
                     k=kmax
111                     if (zlev(k-1)-(k-1)*hnmin.ge.-H(i,j)) then
                          hn(i,j,k)=zlev(k)-zlev(k-1)
                    else
                       hn(i,j,k)=max(zlev(k)-(-H(i,j)+(k-1)*hnmin),hnmin)
                       if (kminset.and.hn(i,j,k).eq.hnmin) then
                          kmin_pmz(i,j)=k+1
                          kminset=.false.
                       end if
                    end if
                    gga(i,j,k)=hn(i,j,k)/HH
                    k=k-1
                    if (k.ge.1) goto 111
                 end if
               end do
            end do
            do i=iimin-1,iimax
               do j=jjmin,jjmax
                  do k=1,kmax
                     huo(i,j,k)=min(hn(i,j,k),hn(i+1,j,k))
                     hun(i,j,k)=huo(i,j,k)
                  end do
                 kumin_pmz(i,j)=max(kmin_pmz(i,j),kmin_pmz(i+1,j))
               end do
            end do
            do i=iimin,iimax
               do j=jjmin-1,jjmax
                  do k=1,kmax
                     hvo(i,j,k)=min(hn(i,j,k),hn(i+1,j,k))
                     hvn(i,j,k)=hvo(i,j,k)
                  end do
                 kvmin_pmz(i,j)=max(kmin_pmz(i,j),kmin_pmz(i,j+1))
               end do
            end do
            do k=1,kmax
               do j=jjmin-1,jjmax
                  do i=iimin,iimax
                     hvo(i,j,k)=min(hn(i,j,k),hn(i+1,j,k))
                     hvn(i,j,k)=hvo(i,j,k)
                  end do
               end do
            end do
            kmin=1
            kumin=1
            kvmin=1
gotm's avatar
gotm committed
212 213 214 215 216 217 218 219 220 221 222 223
         case (3) ! general vertical coordinates
            do k=0,kmax
               ga(k) = k
            end do
            allocate(sig(0:kmax),stat=rc)    ! dimensionless sigma-coordinate
            if (rc /= 0) STOP 'coordinates: Error allocating (sig)'
            allocate(be(0:kmax),stat=rc)     ! dimensionless beta-coordinate
            if (rc /= 0) STOP 'coordinates: Error allocating (be)'
            allocate(gga(I3DFIELD),stat=rc)  ! dimensionless gamma-coordinate
            if (rc /= 0) stop 'coordinates: Error allocating memory (gga)'
            be(0)=  -_ONE_
            sig(0)= -_ONE_
kbk's avatar
kbk committed
224 225
            if (ddu .lt. _ZERO_) ddu=_ZERO_
            if (ddl .lt. _ZERO_) ddl=_ZERO_
gotm's avatar
gotm committed
226 227 228 229 230
            do k=1,kmax
               be(k)=tanh((ddl+ddu)*k/float(kmax)-ddl)+tanh(ddl)
               be(k)=be(k)/(tanh(ddl)+tanh(ddu))-_ONE_
               sig(k)=k/float(kmax)-_ONE_
            end do
kbk's avatar
kbk committed
231
            if (gamma_surf) then
gotm's avatar
gotm committed
232 233 234 235
               kk=kmax
            else
               kk=1
            end if
236 237
            do j=jjmin-HALO,jjmax+HALO
               do i=iimin-HALO,iimax+HALO
kbk's avatar
kbk committed
238
                  HH=max(sseo(i,j)+H(i,j),min_depth)
kbk's avatar
kbk committed
239 240 241
                  alpha=min(&
                           ((be(kk)-be(kk-1))-D_gamma/HH&
                            *(sig(kk)-sig(kk-1)))&
gotm's avatar
gotm committed
242 243 244
                            /((be(kk)-be(kk-1))-(sig(kk)-sig(kk-1))),_ONE_)
                  gga(i,j,0)=-_ONE_
                  do k=1,kmax
245
                     gga(i,j,k)=alpha*sig(k)+(1.-alpha)*be(k)
gotm's avatar
gotm committed
246 247 248 249 250 251 252 253 254 255 256 257 258 259 260
                     if (gga(i,j,k) .lt. gga(i,j,k-1)) then
                        STDERR kk,(be(kk)-be(kk-1)),(sig(kk)-sig(kk-1))
                        STDERR D_gamma,HH
                        STDERR alpha
                        STDERR k-1,gga(i,j,k-1),be(k-1),sig(k-1)
                        STDERR k,gga(i,j,k),be(k),sig(k)
                        stop 'coordinates'
                     end if
                  end do
               end do
            end do

            kmin=1
            kumin=1
            kvmin=1
kbk's avatar
kbk committed
261 262 263
            kmin_pmz=1
            kumin_pmz=1
            kvmin_pmz=1
gotm's avatar
gotm committed
264 265

! Here, the initial layer distribution is calculated.
266 267 268 269 270 271
            do k=1,kmax
               do j=jjmin-HALO,jjmax+HALO
                  do i=iimin-HALO,iimax+HALO
                     HH=max(sseo(i,j)+H(i,j),min_depth)
                     hn(i,j,k)=HH*(gga(i,j,k)-gga(i,j,k-1))
                  end do
gotm's avatar
gotm committed
272 273
               end do
            end do
274 275 276 277 278 279 280 281 282

            do k=1,kmax
               do j=jjmin-HALO,jjmax+HALO
                  do i=iimin-HALO,iimax+HALO-1
                     HH=max(ssuo(i,j)+HU(i,j),min_depth)
                     huo(i,j,k)=HH*0.5*            &
                      (gga(i,j,k)-gga(i,j,k-1)+gga(i+1,j,k)-gga(i+1,j,k-1))
                     hun(i,j,k)=huo(i,j,k)
                  end do
gotm's avatar
gotm committed
283 284 285
               end do
            end do

286 287 288 289 290 291 292 293
            do k=1,kmax
               do j=jjmin-HALO,jjmax+HALO-1
                  do i=iimin-HALO,iimax+HALO
                     HH=max(ssvo(i,j)+HV(i,j),min_depth)
                     hvo(i,j,k)=HH*0.5*            &
                      (gga(i,j,k)-gga(i,j,k-1)+gga(i,j+1,k)-gga(i,j+1,k-1))
                     hvn(i,j,k)=hvo(i,j,k)
                  end do
gotm's avatar
gotm committed
294 295
               end do
            end do
kbk's avatar
kbk committed
296

gotm's avatar
gotm committed
297
         case default
kbk's avatar
kbk committed
298

gotm's avatar
gotm committed
299 300 301 302 303 304 305 306
      end select

   end if ! first

   select case (cord_type)
      case (1) ! sigma coordinates
         if (equiv_sigma) then
            kmaxm1= _ONE_/float(kmax)
307 308
            do j=jjmin-HALO,jjmax+HALO
               do i=iimin-HALO,iimax+HALO
gotm's avatar
gotm committed
309 310 311 312 313
                  ho(i,j,:)=(sseo(i,j)+H(i,j))*kmaxm1
                  hn(i,j,:)=(ssen(i,j)+H(i,j))*kmaxm1
               end do
            end do

314 315
            do j=jjmin-HALO,jjmax+HALO
               do i=iimin-HALO,iimax+HALO-1
gotm's avatar
gotm committed
316 317 318 319 320
                  huo(i,j,:)=(ssuo(i,j)+HU(i,j))*kmaxm1
                  hun(i,j,:)=(ssun(i,j)+HU(i,j))*kmaxm1
               end do
            end do

321 322
            do j=jjmin-HALO,jjmax+HALO-1
               do i=iimin-HALO,iimax+HALO
gotm's avatar
gotm committed
323 324 325 326 327 328
                  hvo(i,j,:)=(ssvo(i,j)+HV(i,j))*kmaxm1
                  hvn(i,j,:)=(ssvn(i,j)+HV(i,j))*kmaxm1
               end do
            end do

         else ! non-equivdistant
329 330 331

            do j=jjmin-HALO,jjmax+HALO
               do i=iimin-HALO,iimax+HALO
gotm's avatar
gotm committed
332 333 334 335 336
                  ho(i,j,1:kmax)=(sseo(i,j)+H(i,j))*dga(1:kmax)
                  hn(i,j,1:kmax)=(ssen(i,j)+H(i,j))*dga(1:kmax)
               end do
            end do

337 338
            do j=jjmin-HALO,jjmax+HALO
               do i=iimin-HALO,iimax+HALO-1
gotm's avatar
gotm committed
339 340 341 342 343
                  huo(i,j,1:kmax)=(ssuo(i,j)+HU(i,j))*dga(1:kmax)
                  hun(i,j,1:kmax)=(ssun(i,j)+HU(i,j))*dga(1:kmax)
               end do
            end do

344 345
            do j=jjmin-HALO,jjmax+HALO-1
               do i=iimin-HALO,iimax+HALO
gotm's avatar
gotm committed
346 347 348 349 350
                  hvo(i,j,1:kmax)=(ssvo(i,j)+HV(i,j))*dga(1:kmax)
                  hvn(i,j,1:kmax)=(ssvn(i,j)+HV(i,j))*dga(1:kmax)
               end do
            end do
         end if
kbk's avatar
kbk committed
351

gotm's avatar
gotm committed
352
      case (2) ! z-level
kbk's avatar
kbk committed
353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399
!        Maybe use mask and not loop over k, e.g. ho(i,j,:)=hn(i,j,:)
!kbk#define USE_ARRAY_SYNTAX
         do j=jjmin-HALO,jjmax+HALO
            do i=iimin-HALO,iimax+HALO
               HH=ssen(i,j)+H(i,j)
#ifdef USE_ARRAY_SYNTAX
               ho(i,j,:)=hn(i,j,:)
               hn(i,j,:)=gga(i,j,:)*HH
#else
               do k=1,kmax
                  ho(i,j,k)=hn(i,j,k)
                  hn(i,j,k)=gga(i,j,k)*HH
               end do
#endif
            end do
         end do
         do j=jjmin-HALO,jjmax+HALO
            do i=iimin-HALO,iimax+HALO-1
               HH=ssun(i,j)+HU(i,j)
#ifdef USE_ARRAY_SYNTAX
               huo(i,j,:)=hun(i,j,:)
               hun(i,j,:)=0.5*(gga(i,j,:)+gga(i+1,j,:))*HH
#else
               do k=1,kmax
                  huo(i,j,k)=hun(i,j,k)
                  hun(i,j,k)=0.5*(gga(i,j,k)+gga(i+1,j,k))*HH
               end do
#endif
            end do
         end do
         do j=jjmin-HALO,jjmax+HALO-1
            do i=iimin-HALO,iimax+HALO
               HH=ssvn(i,j)+HV(i,j)
#ifdef USE_ARRAY_SYNTAX
               hvo(i,j,:)=hvn(i,j,:)
               hvn(i,j,:)=0.5*(gga(i,j,:)+gga(i,j+1,:))*HH
#else
               do k=1,kmax
                  hvo(i,j,k)=hvn(i,j,k)
                  hvn(i,j,k)=0.5*(gga(i,j,k)+gga(i,j+1,k))*HH
               end do
#endif

#undef USE_ARRAY_SYNTAX
            end do
         end do

gotm's avatar
gotm committed
400 401 402 403 404 405
      case (3) ! general vertical coordinates

! The general vertical coordinates can be relaxed towards the new layer
! thicknesses by the following relaxation time scale r. This should
! later be generalised also for sigma coordinates.

406 407
         do j=jjmin-HALO,jjmax+HALO
            do i=iimin-HALO,iimax+HALO
kbk's avatar
kbk committed
408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424
               r=cord_relax/dt*H(i,j)/maxdepth
               HH=ssen(i,j)+H(i,j)
               if (HH .lt. D_gamma) then
                  do k=1,kmax
                     ho(i,j,k)=hn(i,j,k)
                     hn(i,j,k)=HH/float(kmax)
                  end do
               else
                  zz=-H(i,j)
                  do k=1,kmax-1
                     ho(i,j,k)=hn(i,j,k)
                     hn(i,j,k)=(ho(i,j,k)*r+HH*(gga(i,j,k)-gga(i,j,k-1)))/(r+1.)
                     zz=zz+hn(i,j,k)
                  end do
                  ho(i,j,kmax)=hn(i,j,kmax)
                  hn(i,j,kmax)=ssen(i,j)-zz
               end if
gotm's avatar
gotm committed
425 426 427
            end do
         end do

428 429 430
         do j=jjmin-HALO,jjmax+HALO
            do i=iimin-HALO,iimax+HALO-1
!KBK               if (au(i,j) .gt. 0) then
kbk's avatar
kbk committed
431 432 433
                  r=cord_relax/dt*HU(i,j)/maxdepth
                  zz=-HU(i,j)
                  HH=ssun(i,j)+HU(i,j)
gotm's avatar
gotm committed
434
                  do k=1,kmax-1
kbk's avatar
kbk committed
435
                     huo(i,j,k)=hun(i,j,k)
gotm's avatar
gotm committed
436
                     hun(i,j,k)=(huo(i,j,k)*r+HH*0.5*(gga(i,j,k)-gga(i,j,k-1) &
kbk's avatar
kbk committed
437
                               +gga(i+1,j,k)-gga(i+1,j,k-1)))/(r+1.)
kbk's avatar
kbk committed
438
                     zz=zz+hun(i,j,k)
gotm's avatar
gotm committed
439
                  end do
kbk's avatar
kbk committed
440
                  huo(i,j,kmax)=hun(i,j,kmax)
gotm's avatar
gotm committed
441
                  hun(i,j,kmax)=ssun(i,j)-zz
442
!KBK               end if
gotm's avatar
gotm committed
443 444 445
            end do
         end do

446 447 448
         do j=jjmin-HALO,jjmax+HALO-1
            do i=iimin-HALO,iimax+HALO
!KBK               if (av(i,j).gt.0) then
kbk's avatar
kbk committed
449 450 451
                  r=cord_relax/dt*HV(i,j)/maxdepth
                  zz=-HV(i,j)
                  HH=ssvn(i,j)+HV(i,j)
gotm's avatar
gotm committed
452
                  do k=1,kmax-1
kbk's avatar
kbk committed
453
                     hvo(i,j,k)=hvn(i,j,k)
gotm's avatar
gotm committed
454
                     hvn(i,j,k)=(hvo(i,j,k)*r+HH*0.5*(gga(i,j,k)-gga(i,j,k-1) &
kbk's avatar
kbk committed
455
                               +gga(i,j+1,k)-gga(i,j+1,k-1)))/(r+1.)
gotm's avatar
gotm committed
456 457
                     zz=zz+hvn(i,j,k)
                  end do
kbk's avatar
kbk committed
458
                  hvo(i,j,kmax)=hvn(i,j,kmax)
gotm's avatar
gotm committed
459
                  hvn(i,j,kmax)=ssvn(i,j)-zz
460
!KBK               end if
gotm's avatar
gotm committed
461 462 463 464 465 466
            end do
         end do

      case default

   end select
kbk's avatar
kbk committed
467

468 469 470 471 472 473 474 475 476 477 478
#ifdef SLICE_MODEL
   do i=iimin,iimax
      do k=kvmin(i,2),kmax
         hvo(i,1,k)=hvo(i,2,k)
         hvo(i,3,k)=hvo(i,2,k)
         hvn(i,1,k)=hvn(i,2,k)
         hvn(i,3,k)=hvn(i,2,k)
      end do
   end do
#endif

gotm's avatar
gotm committed
479
#ifdef DEBUG
480
   write(debug,*) 'Leaving coordinates()'
gotm's avatar
gotm committed
481 482 483 484 485 486 487 488 489
   write(debug,*)
#endif
   return
   end subroutine coordinates
!EOC

!-----------------------------------------------------------------------
! Copyright (C) 2001 - Hans Burchard and Karsten Bolding               !
!-----------------------------------------------------------------------