anchor.diff:                      endif
   
                       ! - - - - -   u/v grid calculations - - - - - 
  +                    ! contrary to the calculations above, these all take place on the parallelized terrain
                       z_u(:,i,:)   = (sum(dz_scl(1:(i-1))) + dz_scl(i)/2)   &
                                   + h1_u  *  SINH( (smooth_height/s1)**n -  ( (sum(dz_scl(1:(i-1)))+dz_scl(i)/2) /s1)**n ) / SINH((smooth_height/s1)**n)  &! large-scale terrain
                                   + h2_u  *  SINH( (smooth_height/s2)**n -  ( (sum(dz_scl(1:(i-1)))+dz_scl(i)/2) /s2)**n ) / SINH((smooth_height/s2)**n)   ! small terrain features
  
positive.code:                           - MAXVAL(h2)/s2 * COSH(smooth_height/s2)/SINH(smooth_height/s2)
  
              ! with the new (leuenberger et al 2010) Sleve formulation, the inveribiltiy criterion is as follows:
              allocate(gamma_n(this%kds : this%kde+1))
              i=kms                          
              gamma_n(i) =  1                                                     &    
  
negative.code:         h2_v =  h_v  - h1_v
  
          ! In case one wants to see how the terrain is split by smoothing, activate the block below and run in debug:
          ! if ((this_image()==1).and.(options%parameters%debug)) then
          !   call io_write("terrain_smooth_h1.nc", "h1", h1(:,:) )
          !   call io_write("terrain_smooth_h2.nc", "h2", h2(:,:) )
          !   call io_write("h1_u.nc", "h1_u", h1_u(:,:) )
          !   call io_write("h2_u.nc", "h2_u", h2_u(:,:) )
          ! endif
  
          if (this_image()==1) then
             print*, "       Max of full topography", MAXVAL(global_terrain )
  