Index: palm/trunk/SCRIPTS/mbuild
===================================================================
--- palm/trunk/SCRIPTS/mbuild	(revision 2364)
+++ palm/trunk/SCRIPTS/mbuild	(revision 2365)
@@ -27,4 +27,7 @@
 # -----------------
 # $Id$
+# Added lckea & lckeam. KIT/IMK-IFU Garmisch cluster. LRZ (SadiqHuq)
+# 
+# 2316 2017-07-20 07:53:42Z maronga
 # Removed qmake block as mrungui now runs as a python script
 # 
@@ -531,4 +534,5 @@
        # DETERMINE IP-ADDRES OF THE REMOTE-HOST
     case  $remote_host  in
+        (lcbwuni)        remote_address="129.13.82.89";;
         (lcbullhh)       remote_address=136.172.50.13;;
         (lccrayb)        remote_address=130.73.233.1;;
@@ -541,4 +545,6 @@
         (lckiaps)        remote_address=118.128.66.223;;
         (lckyut)         remote_address=133.5.4.37;;
+        (lclrz)          remote_address=129.187.20.240;;
+        (lckea*)         remote_address=172.27.80.109;;
         (lctit)          remote_address=10.1.6.170;;
         (lcxe6)          remote_address=129.177.20.113;;
@@ -547,4 +553,5 @@
         (ibmkisti)       remote_address=150.183.146.24;;
         (ibmku)          remote_address=133.5.4.129;;
+        (ibmmuc)         remote_address=129.187.11.197;;
         (ibms)           remote_address=150.183.5.101;;
         (nech)           remote_address=136.172.44.192;;
@@ -783,5 +790,10 @@
        elif [[ $(echo $remote_host | cut -c1-3) = ibm  &&  $(echo $string | cut -c1-3) = ibm ]]
        then
-          cpp_options="${cpp_options},-D__ibm"
+          if [[ $remote_host = ibmmuc ]]
+          then
+             cpp_options="${cpp_options}"
+          else
+             cpp_options="${cpp_options},-D__ibm"
+          fi
        elif [[ $(echo $remote_host | cut -c1-3) = nec  &&  $(echo $string | cut -c1-3) = nec ]]
        then
@@ -790,5 +802,10 @@
           if [[ $(echo $remote_host | cut -c1-3) = ibm ]]
           then
-             cpp_options="${cpp_options},-D__$string"
+             if [[ $remote_host = ibmmuc ]]
+             then
+                cpp_options="${cpp_options} -D__parallel"
+             else
+                cpp_options="${cpp_options},-D__$string"
+             fi
           else
              cpp_options="$cpp_options -D__$string "
Index: palm/trunk/SCRIPTS/mrun
===================================================================
--- palm/trunk/SCRIPTS/mrun	(revision 2364)
+++ palm/trunk/SCRIPTS/mrun	(revision 2365)
@@ -27,4 +27,8 @@
 # -----------------
 # $Id$
+# Added lckea & lckeam. KIT/IMK-IFU Garmisch cluster. LRZ (SadiqHuq)
+# Vertical grid nesting: set vnested_mode. -N Procs for fine and coarse grid.
+# 
+# 2303 2017-07-04 12:26:18Z raasch
 # bugfix: setting default value for write_binary
 # 
@@ -334,4 +338,6 @@
  module_calls=""
  mrun_script_name=$mc
+ vnested_dist=""
+ vnested_mode="vnested_twi"
  netcdf_inc=""
  netcdf_lib=""
@@ -341,4 +347,6 @@
  numprocs_atmos=0
  numprocs_ocean=0
+ numprocs_crse=0
+ numprocs_fine=0
  OOPT=""
  openmp=false
@@ -373,4 +381,5 @@
  run_coupled_model=false
  run_mode=""
+ run_vnested_model=false
  dashes="  ----------------------------------------------------------------------------"
  silent=false
@@ -475,5 +484,5 @@
     # READ SHELLSCRIPT-OPTIONS AND REBUILD THE MRUN-COMMAND STRING (MC),
     # WHICH WILL BE USED TO START RESTART-JOBS
- while  getopts  :a:bBc:Cd:Fg:G:h:H:i:kK:m:M:n:o:O:p:P:q:r:R:s:St:T:u:U:vw:xX:yY:zZ option
+ while  getopts  :a:bBc:Cd:Fg:G:h:H:i:kK:m:M:n:N:o:O:p:P:q:r:R:s:St:T:u:U:vw:xX:yY:zZ option
  do
    case  $option  in
@@ -495,4 +504,5 @@
        (M)   makefile=$OPTARG; mc="$mc -M$OPTARG";;
        (n)   node_usage=$OPTARG; mc="$mc -n$OPTARG";; 
+       (N)   run_vnested_model=true; vnested_dist=$OPTARG; mc="$mc -N'$OPTARG'";;
        (o)   output_list=$OPTARG; mc="$mc -o'$OPTARG'";;
        (O)   use_openmp=true; threads_per_task=$OPTARG; mc="$mc -O$OPTARG";;
@@ -550,4 +560,6 @@
     printf "\n        -M    Makefile name                            Makefile"
     printf "\n        -n    node usage (shared/not_shared)           depending on -h"
+    printf "\n        -N    Vertical grid nesting. Number of"
+    printf "\n              PE for Coarse and Fine grid"
     printf "\n        -o    OUTPUT control list                      \"\" "
     printf "\n        -O    threads per openMP task                  ---"
@@ -713,4 +725,34 @@
  fi
 
+      # NESTING (-N) selected and vnesting_mode specified
+ if [[ $run_vnested_model = true ]] 
+ then
+
+    if  [[ -n $vnested_dist ]]
+    then
+
+       numprocs_crse=`echo $vnested_dist | cut -d" " -s -f1`
+       numprocs_fine=`echo $vnested_dist | cut -d" " -s -f2`
+
+       if (( $numprocs_crse + $numprocs_fine != $numprocs ))
+       then
+
+          printf "\n  +++ number of processors does not fit to specification by \"-N\"."
+          printf "\n      PEs (total)     : $numprocs"
+          printf "\n      PEs (Coarse)    : $numprocs_crse"
+          printf "\n      PEs (Fine)      : $numprocs_fine"
+          locat=vnesting; exit
+
+       fi
+
+    else
+
+          printf "\n  +++ "Specify PE for fine and coarse grid: -N "nCGPE nFGPE \"-N\"."
+          locat=vnesting; exit
+
+    fi
+    vnested_dist=`echo "$numprocs_crse $numprocs_fine"`
+ fi
+
     # SAVE VALUES OF MRUN-OPTIONS SICHERN IN ORDER TO OVERWRITE
     # THOSE VALUES GIVEN IN THE CONFIGURATION-FILE
@@ -1195,4 +1237,6 @@
         (ibmh)       queue=cluster;;
         (ibmkisti)   queue=class.32plus;;
+        (ibmmuc*)    queue=test;;
+        (lcbwuni)    queue=develop;;
         (lcbullhh)   queue=compute;;
         (lccrayb)    queue=mpp1q;;
@@ -1203,4 +1247,5 @@
         (lckyuh)     queue=fx-single;;
         (lckyut)     queue=cx-single;;
+        (lclrz)      queue=mpp2;;
         (lctit)      queue=S;;
         (unics)      queue=unics;;
@@ -1764,5 +1809,5 @@
 
        # SET PREPROCESSOR-DIRECTIVES TO SELECT OPERATING SYSTEM SPECIFIC CODE
-    if [[ $(echo $localhost | cut -c1-3) = ibm ]]
+    if [[ $(echo $localhost | cut -c1-3) = ibm && $localhost != ibmmuc* ]]
     then
        cpp_options="${cpp_options},-D__ibm=__ibm"
@@ -1778,5 +1823,5 @@
 
        # SET DIRECTIVES GIVEN BY OPTION -K (E.G. parallel)
-    if [[ $(echo $localhost | cut -c1-3) = ibm ]]
+    if [[ $(echo $localhost | cut -c1-3) = ibm && $localhost != ibmmuc ]]
     then
        [[ -n $cond1 ]]  &&  cpp_options="${cpp_options},-D__$cond1=__$cond1"
@@ -1893,6 +1938,11 @@
     fi
  fi
- TEMPDIR=$tmp_user_catalog/${usern}.$kennung
-
+
+if [[ $localhost = ibmmuc* ]]
+then
+   TEMPDIR=$tmp_user_catalog/${USER}.$kennung
+else
+   TEMPDIR=$tmp_user_catalog/${usern}.$kennung
+fi
 
     # DETERMINE THE NAME OF THE DIRECTORY WHICH IS USED TO TEMPORARILY STORE DATA FOR RESTART RUNS
@@ -2792,4 +2842,41 @@
                    then
                       /opt/optibm/HPM_2_4_1/bin/hpmcount  a.out
+                   elif [[  $localhost = ibmmuc* ]]
+                   then
+                      ulimit -c unlimited        # only for debgingg
+                      echo $MP_NODES        > ~/job_queue/hostfile.$kennung
+                      echo $MP_PROCS       >> ~/job_queue/hostfile.$kennung
+                      cat $LOADL_HOSTFILE  >> ~/job_queue/hostfile.$kennung
+                      export MP_NODES=$nodes
+                      export MP_PROCS=$numprocs
+                      # export MPI_SINGLE_THREAD=no   # LRZ NetCDF
+                      # export MP_TASKS_PER_NODE=$tasks_per_node
+                      echo "Resource Info: " 
+                      echo "numprocs:         " $numprocs " MP_PROCS " $MP_PROCS
+                      echo "nodes:            " $nodes    " MP_NODES " $MP_NODES
+                      echo "tasks_per_node:   " $tasks_per_node
+                      echo "threads_per_task: " $threads_per_task
+                      export OMP_NUM_THREADS=1
+                      source /lrz/sys/share/modules/init/bash
+                      module li
+                      echo "runfile_atmos"
+                      
+                      if [[ $run_vnested_model = true ]]
+                      then
+                         
+                         printf "\n      Nested run ($numprocs_crse Coarse, $numprocs_fine Fine)"
+                         printf "\n      using $nested_mode nesting"
+                         printf "\n\n"
+                         
+                         echo "$vnested_mode $numprocs_crse $numprocs_fine"  >  runfile_atmos
+                         
+                         poe   ./a.out  <  runfile_atmos
+                         # mpiexec -n $numprocs  ./a.out  <  runfile_atmos
+                      else
+                         echo "precursor_atmos"  >  runfile_atmos
+                         poe   ./a.out -proc $numprocs -nodes  $nodes <  runfile_atmos  
+                         # mpiexec -n $numprocs  ./a.out  <  runfile_atmos  
+                      
+                      fi
                    else
                       if [[ $run_coupled_model = false ]]
@@ -2899,5 +2986,5 @@
                    printf "\n      threads per task: $threads_per_task  stacksize: unlimited"
                 fi
-                if [[ $run_coupled_model = false ]]
+                if [[ $run_coupled_model = false && $run_vnested_model = false ]]
                 then
                    if [[ "$ocean_file_appendix" = true ]]
@@ -2948,5 +3035,10 @@
                    then
                       mpirun_rsh -hostfile $PBS_NODEFILE -np `cat $PBS_NODEFILE | wc -l` a.out  < runfile_atmos
-
+                   elif [[ $host = lclrz || $host = lcbwuni  ]]
+                   then
+                      mpiexec  -n $ii  a.out  <  runfile_atmos  $ROPTeS
+                   elif [[ $host = lckea* ]]
+                   then
+                      srun  -n $ii  a.out  <  runfile_atmos  $ROPTeS
                    elif [[ $host = lckiaps ]]
                    then
@@ -2959,6 +3051,6 @@
                    fi
 
-                else
-
+                elif [[ $run_coupled_model = true ]]
+                then
                        # COUPLED RUN
                    (( iia = $numprocs_atmos / $threads_per_task ))
@@ -2992,4 +3084,23 @@
                    wait
 
+
+                elif [[ $run_vnested_model = true ]]
+                then
+                    printf "\n      Vertical Nested run ($numprocs_crse Coarse, $numprocs_fine Fine)"
+                    printf "\n      using $vnested_mode vnesting"
+                    printf "\n\n"
+
+                     echo "$vnested_mode $numprocs_crse $numprocs_fine"  >  runfile_atmos
+
+                     if [[  $host = lcbwuni  ||  $host = lclrz*  ]]
+                     then
+                         mpiexec  -n $ii  a.out  <  runfile_atmos  $ROPTeS
+                     elif [[ $host = lckea*  ]]
+                     then
+                         srun  -n $ii  a.out  <  runfile_atmos  $ROPTeS 
+                     else
+                        mpirun   -np $numprocs  ./a.out  $ROPTS  < runfile_atmos
+                     fi
+                     wait
                 fi
 
@@ -3657,4 +3768,5 @@
     [[ "$ocean_file_appendix" = true ]]  &&  mrun_com=${mrun_com}" -y"
     [[ $run_coupled_model = true ]]  &&  mrun_com=${mrun_com}" -Y \"$coupled_dist\""
+    [[ $run_vnested_model  = true ]]  &&  mrun_com=${mrun_com}" -N \"$vnested_dist\""
     [[ "$check_namelist_files" = false ]]  &&  mrun_com=${mrun_com}" -z"
     [[ "$combine_plot_fields" = false ]]  &&  mrun_com=${mrun_com}" -Z"
Index: palm/trunk/SCRIPTS/subjob
===================================================================
--- palm/trunk/SCRIPTS/subjob	(revision 2364)
+++ palm/trunk/SCRIPTS/subjob	(revision 2365)
@@ -28,4 +28,7 @@
 # -----------------
 # $Id$
+# Added lckea & lckeam. KIT/IMK-IFU Garmisch cluster. LRZ (SadiqHuq)
+# 
+# 2295 2017-06-27 14:25:52Z raasch
 # adjustments for using lcgeohu (cirrus @ HUB)
 # 
@@ -252,9 +255,15 @@
      (inferno)               local_address=130.75.105.5;   local_host=lcmuk;;
      (irifi)                 local_address=130.75.105.104; local_host=lcmuk;;
+   # (i*)                    local_address=129.187.11.197; local_host=ibmmuc;;
      (jaboticaba)            local_address=150.163.25.181; local_host=lcbr;;
+     (kea*)                  local_address=172.27.80.109;  local_host=lckeal;;
      (sno)                   local_address=130.75.105.113; local_host=lcmuk;;
      (kuma)                  local_address=130.75.105.115; local_host=lcmuk;;
      (levanto)               local_address=130.75.105.45;  local_host=lcmuk;;
      (login*)                local_address=118.128.66.201; local_host=lckiaps;;
+   # (login*)                local_address=129.187.11.197; local_host=ibmmuc;;
+     (lm*)                   local_address=129.187.11.197; local_host=ibmmuc;;
+     (lx*)                   local_address=129.187.20.240; local_host=lclrz;;
+     (mpp2*)                 local_address=129.187.20.105; local_host=lclrz;;
      (maestro)               local_address=130.75.105.2;   local_host=lcmuk;;
      (meller)                local_address=134.106.74.155; local_host=lcfor;;
@@ -278,4 +287,5 @@
      (tc*)                   local_address="ocean";        local_host=lcocean;; 
      (t2a*)                  local_address=10.1.6.165;     local_host=lctit;;
+     (uc1n*)                 local_address=129.13.82.89;   local_host=lcbwuni;;
      (urban*)                local_address=147.46.30.151   local_host=lcsb;;
      (vinessa)               local_address=130.75.105.112; local_host=lcmuk;;
@@ -409,4 +419,6 @@
         (ibmku)   queue=s4; remote_address=133.5.4.129; submcom=/usr/local/bin/llsubmit;;
         (ibms)    queue=p_normal; remote_address=150.183.5.101; submcom=/usr/lpp/LoadL/full/bin/llsubmit;;
+        (ibmmuc)  remote_address=129.187.11.197; submcom=/usr/bin/llsubmit;;
+        (lcbwuni) queue=develop; remote_address=129.13.82.89; submcom=/opt/moab/bin/msub;;
         (lcbullhh)    queue=compute; remote_address=136.172.50.13; submcom=/usr/bin/sbatch;;
         (lccrayb) queue=mpp1testq; remote_address=130.73.233.1; submcom="/opt/moab/default/bin/msub -E";;
@@ -417,8 +429,10 @@
         (lckyoto) remote_address=133.3.51.11; submcom=/thin/local/bin/qsub;;
         (lck)     remote_address=165.132.26.61; submcom=/usr/torque/bin/qsub;;
+        (lckeal)  queue=ivy; remote_address=172.27.80.109; submcom=/usr/bin/sbatch;;
         (lckiaps) remote_address=118.128.66.201; submcom=/opt/pbs/default/bin/qsub;;
         (lckordi) remote_address=210.219.61.8; submcom=/usr/torque/bin/qsub;;
         (lckyuh)  remote_address=133.5.4.33; submcom=/usr/bin/pjsub;;
         (lckyut)  remote_address=133.5.4.37; submcom=/usr/bin/pjsub;;
+        (lclrz)   remote_address=129.187.20.240; submcom=/usr/bin/sbatch;;
         (lcocean) remote_address="ocean"; submcom=qsub;;
         (lcsb)    remote_address=147.46.30.151; submcom=/usr/torque/bin/qsub;;
@@ -459,8 +473,16 @@
                      (*)                                     error=true;;
                  esac;;
+        (ibmmuc*) case  $ndq  in
+                     (test|micro|general|large|fat|fattest|special|tmp1|tmp2) error=false;;
+                     (*)                                     error=true;;
+                 esac;;
         (lcbullhh) case  $ndq  in
                      (compute|compute2|shared)  error=false;;
                      (*)                                     error=true;;
                  esac;;
+        (lcbwuni) case  $ndq  in
+                     (develop|singlenode|multinode|verylong|fat) error=false;;
+                     (*)                                     error=true;;
+                 esac;;
         (lccrayb) case  $ndq  in
                      (dataq|mpp1q|mpp1testq|mpp2q|mpp2testq|smp1q|smp1testq|specialm1q)   error=false;;
@@ -493,4 +515,8 @@
         (lckyut) case  $ndq  in
                      (cx-dbg|cx-single|cx-small|cx-middle|cx-large)  error=false;;
+                     (*)                                     error=true;;
+                 esac;;
+        (lclrz) case  $ndq  in
+                     (mpp1|mpp2|iuv2|myri)                   error=false;;
                      (*)                                     error=true;;
                  esac;;
@@ -744,5 +770,32 @@
     then
 
-       if [[ $remote_host != ibmkisti ]]
+       if [[ $remote_host == ibmmuc* ]]
+       then
+       cat > $job_to_send << %%END%%
+
+#!/bin/bash
+# @ job_type = parallel
+# @ job_name = $job_name
+# @ output = $remote_dayfile
+# @ error = $remote_dayfile
+# @ wall_clock_limit = $timestring
+$class
+$mcm_affinity_options
+$task_affinity
+$notify_user
+# @ network.MPI = sn_all,not_shared,us
+# @ notification = always
+# @ energy_policy_tag = table_kit_ifu
+# @ minimize_time_to_solution = yes
+# @ node = $nodes
+# @ total_tasks = $numprocs
+# @ node_topology = island
+# @ island_count = 1,2
+# @ environment = LD_LIBRARY_PATH=/lrz/sys/libraries/netcdf/4.2.1.1_impi4/lib:/lrz/sys/libraries/hdf5/1.8.15/ibmmpi/lib:/lrz/sys/libraries/fftw/3.3.3/avx/lib/
+# @ queue
+
+%%END%%
+
+       elif [[ $remote_host != ibmkisti ]]
        then
 
@@ -1248,4 +1301,80 @@
        submcom="$submcom -W group_list=$group_number -N $job_name -l walltime=$timestring -l select=$nodes:ncpus=$processes_per_node:mpiprocs=$tasks_per_node:mem=${Memory}gb -l place=scatter -o $remote_dayfile -j oe -et 1 -q $queue "
     fi
+
+ elif [[ $remote_host = lclrz ]]
+ then
+       cat > $job_to_send << %%END%%
+#!/bin/bash 
+#SBATCH -J $job_name
+#SBATCH -t $timestring
+#SBATCH -N $nodes
+#SBATCH --ntasks-per-node=$processes_per_node
+#SBATCH --get-user-env
+#SBATCH -o $remote_dayfile
+#SBATCH -e $remote_dayfile
+#SBATCH --mail-user=${email_notification}
+#SBATCH --clusters=$queue
+
+$init_cmds
+
+$module_calls
+
+%%END%%
+
+ elif [[ $remote_host = lckea* ]]
+ then
+       keal_tasks_per_core=1
+
+       if [[ $queue = haswell || $queue = ivy* ]]
+       then
+          if (( tasks_per_node > 20 ))
+          then
+             keal_tasks_per_core=2
+	  fi
+       fi
+
+       cat > $job_to_send << %%END%%
+#!/bin/bash 
+#SBATCH -J $job_name
+#SBATCH -t $timestring
+#SBATCH -N $nodes
+#SBATCH --ntasks-per-node=$processes_per_node
+#SBATCH --ntasks-per-core=$keal_tasks_per_core
+#SBATCH --mem-per-cpu=${memory}mb
+#SBATCH --get-user-env
+#SBATCH -o $remote_dayfile
+#SBATCH -e $remote_dayfile
+#SBATCH --mail-user=${email_notification}
+#SBATCH --mail-type=ALL
+#SBATCH --partition=$queue 
+
+export MV2_ENABLE_AFFINITY=0
+
+$init_cmds
+
+$module_calls
+
+%%END%%
+
+ elif [[ $remote_host=lcbwuni ]]
+ then
+    if [[ $email_notification = none ]]
+    then
+       email_directive=""
+    else
+       email_directive="#PBS -M $email_notification"
+    fi
+       cat > $job_to_send << %%END%%
+#!/bin/ksh
+#PBS -N $job_name
+#PBS -l walltime=$timestring
+#PBS -l nodes=${nodes}:ppn=$processes_per_node
+#PBS -l pmem=${memory}mb
+#PBS -m abe
+#PBS -o $remote_dayfile
+#PBS -j oe
+#PBS -q $queue
+$email_directive
+%%END%%
 
  else
@@ -1615,4 +1744,10 @@
              eval  $submcom  $job_on_remhost
           fi
+       elif [[ $local_host = lclrz || $local_host = lckea* ]]
+       then
+          eval  $submcom  $job_on_remhost
+       elif [[ $local_host = lcbwuni ]]
+       then
+          msub  -q $queue $job_on_remhost
        else
           qsub  $job_on_remhost
Index: palm/trunk/SOURCE/Makefile
===================================================================
--- palm/trunk/SOURCE/Makefile	(revision 2364)
+++ palm/trunk/SOURCE/Makefile	(revision 2365)
@@ -25,4 +25,7 @@
 # -----------------
 # $Id$
+# Added dependencies for vertical_nesting_mod
+# 
+# 2339 2017-08-07 13:55:26Z gronemeier
 # corrected timestamp in header
 # 
@@ -409,5 +412,5 @@
 	user_lpm_init.f90 user_lpm_set_attributes.f90 user_module.f90 \
 	user_parin.f90 user_read_restart_data.f90 \
-	user_spectra.f90 user_statistics.f90 virtual_flight_mod.f90 \
+	user_spectra.f90 user_statistics.f90 vertical_nesting_mod.f90 virtual_flight_mod.f90 \
 	wind_turbine_model_mod.f90 write_3d_binary.f90 write_var_list.f90
 
@@ -425,18 +428,4 @@
 # The following line is needed for palm_simple_install, don't remove it!
 #to_be_replaced_by_include
-
-#BOUNDS="-Rbc"  # Array bounds checking. Compromises performance seriously.
-
-F90 = ftn
-#COPT = -DMPI_REAL=MPI_DOUBLE_PRECISION -DMPI_2REAL=MPI_2DOUBLE_PRECISION -D__lc -D__parallel -D__fftw \
-#             -D__netcdf -D__netcdf4 -D__netcdf4_parallel -D__nopointer
-COPT = -DMPI_REAL=MPI_DOUBLE_PRECISION -DMPI_2REAL=MPI_2DOUBLE_PRECISION -D__lc -D__parallel -D__fftw \
-             -D__netcdf -D__netcdf4 -D__netcdf4_parallel
-
-#F90FLAGS = -O2 -e Fm -G2 -rm $(BOUNDS) -hnoomp
-F90FLAGS = -O2 -e F -G2 -rm $(BOUNDS) -hnoomp
-LDFLAGS = $(F90FLAGS) -dynamic
-
-# End of include
 
 .SUFFIXES:
@@ -469,5 +458,5 @@
 average_3d_data.o: modules.o cpulog_mod.o mod_kinds.o exchange_horiz_2d.o land_surface_model_mod.o \
 	radiation_model_mod.o urban_surface_mod.o
-boundary_conds.o: modules.o mod_kinds.o pmc_interface_mod.o surface_mod.o
+boundary_conds.o: modules.o mod_kinds.o pmc_interface_mod.o surface_mod.o vertical_nesting_mod.o
 buoyancy.o: modules.o mod_kinds.o
 calc_mean_profile.o: modules.o mod_kinds.o
@@ -480,5 +469,5 @@
 	microphysics_mod.o model_1d_mod.o netcdf_interface_mod.o plant_canopy_model_mod.o pmc_interface_mod.o \
 	radiation_model_mod.o spectra_mod.o subsidence_mod.o synthetic_turbulence_generator_mod.o \
-	wind_turbine_model_mod.o urban_surface_mod.o
+	wind_turbine_model_mod.o urban_surface_mod.o vertical_nesting_mod.o
 close_file.o: modules.o mod_kinds.o netcdf_interface_mod.o
 compute_vpt.o: modules.o mod_kinds.o
@@ -508,5 +497,5 @@
 diffusion_v.o: modules.o mod_kinds.o surface_mod.o
 diffusion_w.o: modules.o mod_kinds.o surface_mod.o
-diffusivities.o: modules.o mod_kinds.o surface_mod.o
+diffusivities.o: modules.o mod_kinds.o surface_mod.o vertical_nesting_mod.o
 disturb_field.o: modules.o cpulog_mod.o mod_kinds.o random_function_mod.o \
                  random_generator_parallel_mod.o
@@ -530,10 +519,11 @@
 init_advec.o: modules.o mod_kinds.o
 init_cloud_physics.o: modules.o mod_kinds.o
-init_coupling.o: modules.o mod_kinds.o
+init_coupling.o: modules.o mod_kinds.o vertical_nesting_mod.o
 init_dvrp.o: modules.o mod_kinds.o
-init_grid.o: modules.o mod_kinds.o advec_ws.o netcdf_interface_mod.o surface_mod.o
+init_grid.o: modules.o mod_kinds.o advec_ws.o netcdf_interface_mod.o surface_mod.o vertical_nesting_mod.o
 init_masks.o: modules.o mod_kinds.o netcdf_interface_mod.o
 init_ocean.o: modules.o eqn_state_seawater.o mod_kinds.o
-init_pegrid.o: modules.o mod_kinds.o spectra_mod.o synthetic_turbulence_generator_mod.o
+init_pegrid.o: modules.o mod_kinds.o spectra_mod.o synthetic_turbulence_generator_mod.o \
+   vertical_nesting_mod.o
 init_pt_anomaly.o: modules.o mod_kinds.o
 init_rankine.o: modules.o mod_kinds.o
@@ -637,5 +627,5 @@
 surface_layer_fluxes_mod.o: modules.o mod_kinds.o exchange_horiz_2d.o land_surface_model_mod.o \
 	urban_surface_mod.o surface_mod.o
-surface_mod.o: modules.o mod_kinds.o init_pegrid.o model_1d_mod.o 
+surface_mod.o: modules.o mod_kinds.o model_1d_mod.o 
 swap_timelevel.o: modules.o cpulog_mod.o mod_kinds.o land_surface_model_mod.o \
    pmc_interface_mod.o urban_surface_mod.o
@@ -648,5 +638,5 @@
 	spectra_mod.o user_actions.o microphysics_mod.o synthetic_turbulence_generator_mod.o \
         surface_layer_fluxes_mod.o surface_mod.o \
-	urban_surface_mod.o virtual_flight_mod.o wind_turbine_model_mod.o
+	urban_surface_mod.o vertical_nesting_mod.o virtual_flight_mod.o wind_turbine_model_mod.o
 time_integration_spinup.o: cpulog_mod.o modules.o land_surface_model_mod.o mod_kinds.o radiation_model_mod.o surface_layer_fluxes_mod.o urban_surface_mod.o 
 time_to_string.o: mod_kinds.o
@@ -689,4 +679,5 @@
 user_spectra.o: modules.o mod_kinds.o spectra_mod.o user_module.o
 user_statistics.o: modules.o mod_kinds.o netcdf_interface_mod.o user_module.o
+vertical_nesting_mod.o: modules.o mod_kinds.o surface_mod.o 
 virtual_flight_mod.o: modules.o cpulog_mod.o mod_kinds.o netcdf_interface_mod.o user_init_flight.o user_flight.o
 wind_turbine_model_mod.o: modules.o cpulog_mod.o mod_kinds.o
Index: palm/trunk/SOURCE/boundary_conds.f90
===================================================================
--- palm/trunk/SOURCE/boundary_conds.f90	(revision 2364)
+++ palm/trunk/SOURCE/boundary_conds.f90	(revision 2365)
@@ -25,4 +25,8 @@
 ! -----------------
 ! $Id$
+! Vertical grid nesting implemented: exclude setting vertical velocity to zero 
+! on fine grid (SadiqHuq)
+! 
+! 2320 2017-07-21 12:47:43Z suehring
 ! Remove unused control parameter large_scale_forcing from only-list
 ! 
@@ -166,5 +170,5 @@
     USE control_parameters,                                                    &
         ONLY:  bc_pt_t_val, bc_q_t_val, bc_s_t_val, constant_diffusion,        &
-               cloud_physics, dt_3d, humidity,                                 &
+               cloud_physics, coupling_mode, dt_3d, humidity,                  &
                ibc_pt_b, ibc_pt_t, ibc_q_b, ibc_q_t, ibc_s_b, ibc_s_t,         &
                ibc_sa_t, ibc_uv_b, ibc_uv_t, inflow_l, inflow_n, inflow_r,     &
@@ -237,6 +241,9 @@
     ENDIF
 
-    IF ( .NOT. nest_domain )  THEN
-       w_p(nzt:nzt+1,:,:) = 0.0_wp  ! nzt is not a prognostic level (but cf. pres)
+!
+!-- Vertical nesting: Vertical velocity not zero at the top of the fine grid
+    IF (  .NOT.  nest_domain  .AND.                                            &
+                 TRIM(coupling_mode) /= 'vnested_fine' )  THEN
+       w_p(nzt:nzt+1,:,:) = 0.0_wp  !< nzt is not a prognostic level (but cf. pres)
     ENDIF
 
Index: palm/trunk/SOURCE/check_parameters.f90
===================================================================
--- palm/trunk/SOURCE/check_parameters.f90	(revision 2364)
+++ palm/trunk/SOURCE/check_parameters.f90	(revision 2365)
@@ -25,4 +25,8 @@
 ! -----------------
 ! $Id$
+! Vertical grid nesting implemented: Check coupling mode. Generate file header
+! (SadiqHuq)
+! 
+! 2354 2017-08-17 10:49:36Z schwenkel
 ! Bugfix correlated to lsm_check_data_output_pr.
 ! If-statement for following checks is essential, otherwise units for lsm output 
@@ -561,4 +565,6 @@
     USE wind_turbine_model_mod,                                                &
         ONLY:  wtm_check_parameters, wind_turbine
+    USE vertical_nesting_mod,                                                  &
+        ONLY:  vnested, vnest_check_parameters
 
 
@@ -603,4 +609,6 @@
 !> @todo Check if any queries for other coupling modes (e.g. precursor_ocean) are missing
     IF ( coupling_mode /= 'uncoupled'            .AND.  &
+         coupling_mode /= 'vnested_crse'         .AND.  &
+         coupling_mode /= 'vnested_fine'         .AND.  &
          coupling_mode /= 'atmosphere_to_ocean'  .AND.  &
          coupling_mode /= 'ocean_to_atmosphere' )  THEN
@@ -611,5 +619,6 @@
 !
 !-- Check dt_coupling, restart_time, dt_restart, end_time, dx, dy, nx and ny
-    IF ( coupling_mode /= 'uncoupled')  THEN
+    IF ( coupling_mode /= 'uncoupled'  .AND.                                   &
+         coupling_mode(1:8) /= 'vnested_' )  THEN
 
        IF ( dt_coupling == 9999999.9_wp )  THEN
@@ -809,4 +818,8 @@
     IF ( coupling_mode == 'uncoupled' )  THEN
        coupling_string = ''
+    ELSEIF ( coupling_mode == 'vnested_crse' )  THEN
+       coupling_string = ' nested (coarse)'
+    ELSEIF ( coupling_mode == 'vnested_fine' )  THEN
+       coupling_string = ' nested (fine)'
     ELSEIF ( coupling_mode == 'atmosphere_to_ocean' )  THEN
        coupling_string = ' coupled (atmosphere)'
@@ -3858,4 +3871,7 @@
     ENDIF
 
+!-- Vertical nesting: check fine and coarse grid compatibility for data exchange
+    IF ( vnested )  CALL vnest_check_parameters
+
     CALL location_message( 'finished', .TRUE. )
 !
Index: palm/trunk/SOURCE/diffusivities.f90
===================================================================
--- palm/trunk/SOURCE/diffusivities.f90	(revision 2364)
+++ palm/trunk/SOURCE/diffusivities.f90	(revision 2365)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Vertical grid nesting implemented (SadiqHuq)
+! 
+! 2233 2017-05-30 18:08:54Z suehring
 !
 ! 2232 2017-05-30 17:47:52Z suehring
@@ -108,4 +111,7 @@
         ONLY :  bc_h
 
+    USE vertical_nesting_mod,                                                  &
+        ONLY:  vnest_boundary_conds_khkm, vnest_init
+
     IMPLICIT NONE
 
@@ -266,4 +272,7 @@
        kh(:,nyn+1,:) = kh(:,nyn,:)
     ENDIF
+!
+!-- Vertical nesting: set fine grid eddy viscosity top boundary condition
+    IF ( vnest_init )  CALL vnest_boundary_conds_khkm
 
  END SUBROUTINE diffusivities
Index: palm/trunk/SOURCE/init_coupling.f90
===================================================================
--- palm/trunk/SOURCE/init_coupling.f90	(revision 2364)
+++ palm/trunk/SOURCE/init_coupling.f90	(revision 2365)
@@ -25,4 +25,7 @@
 ! ------------------
 ! $Id$
+! Vertical nesting implemented (SadiqHuq)
+! 
+! 2298 2017-06-29 09:28:18Z raasch
 ! MPI2 coupling removed
 ! 
@@ -67,4 +70,6 @@
     USE pegrid
 
+    USE vertical_nesting_mod
+
     IMPLICIT NONE
 
@@ -96,4 +101,6 @@
        IF ( TRIM( coupling_mode ) == 'coupled_run' )  THEN
           i = 1
+       ELSEIF ( TRIM( coupling_mode ) == 'vnested_twi' )  THEN
+          i = 9
        ELSE
           i = 0
@@ -111,5 +118,5 @@
     i = bc_data(0)
 
-    IF ( i == 0 ) THEN
+    IF ( i == 0 )  THEN
        coupling_mode = 'uncoupled'
 !
@@ -121,4 +128,32 @@
           CLOSE ( 90 )
        ENDIF
+    ELSEIF ( i == 9 )  THEN
+
+!
+!--    Set a flag to identify runs with vertical nesting
+       vnested = .TRUE.
+       
+       comm_inter = MPI_COMM_WORLD
+       
+!
+!--    Split the total available PE's into two groups
+!--    numprocs for Coarse and Fine Grid are specified via mrun argument -N 
+       IF ( myid < bc_data(1) )  THEN
+          inter_color     = 0
+          numprocs        = bc_data(1)
+          coupling_mode   = 'vnested_crse'
+       ELSE
+          inter_color     = 1
+          numprocs        = bc_data(2)
+          coupling_mode   = 'vnested_fine'
+       ENDIF
+       
+       CALL MPI_COMM_SPLIT( MPI_COMM_WORLD, inter_color, 0, comm_palm, ierr )
+       comm2d = comm_palm
+       
+       OPEN( 90, FILE='VNESTING_PORT_OPENED', FORM='FORMATTED' )
+       WRITE ( 90, '(''TRUE'')' )
+       CLOSE ( 90 )
+      
     ELSE
        comm_inter = MPI_COMM_WORLD
@@ -154,3 +189,9 @@
     ENDIF
 
+    IF (  TRIM( coupling_mode ) == 'vnested_fine' )  THEN
+!
+!-- Set file extension for vertical nesting
+       coupling_char = '_N'
+    ENDIF
+
  END SUBROUTINE init_coupling
Index: palm/trunk/SOURCE/init_grid.f90
===================================================================
--- palm/trunk/SOURCE/init_grid.f90	(revision 2364)
+++ palm/trunk/SOURCE/init_grid.f90	(revision 2365)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Vertical nesting implemented (SadiqHuq)
+! 
+! 2319 2017-07-20 17:33:17Z suehring
 ! Remove print statements
 ! 
@@ -250,5 +253,6 @@
                canyon_height, canyon_wall_left, canyon_wall_south,             &
                canyon_width_x, canyon_width_y, constant_flux_layer,            &
-               coupling_char, dp_level_ind_b, dz, dz_max, dz_stretch_factor,   &
+               coupling_char, coupling_mode,                                   &
+               dp_level_ind_b, dz, dz_max, dz_stretch_factor,                  &
                dz_stretch_level, dz_stretch_level_index, grid_level, ibc_uv_b, &
                io_blocks, io_group, inflow_l, inflow_n, inflow_r, inflow_s,    &
@@ -285,4 +289,7 @@
     USE surface_mod,                                                           &
         ONLY:  get_topography_top_index, init_bc
+
+    USE vertical_nesting_mod,                                                  &
+        ONLY:  vnested, vnest_init_grid
 
     IMPLICIT NONE
@@ -2000,4 +2007,8 @@
     CALL exchange_horiz( l_wall, nbgp )     
 
+!
+!-- Vertical nesting: communicate vertical grid level arrays between fine and
+!-- coarse grid
+    IF ( vnested )  CALL vnest_init_grid
 
  END SUBROUTINE init_grid
Index: palm/trunk/SOURCE/init_pegrid.f90
===================================================================
--- palm/trunk/SOURCE/init_pegrid.f90	(revision 2364)
+++ palm/trunk/SOURCE/init_pegrid.f90	(revision 2365)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Vertical nesting implemented (SadiqHuq)
+! 
+! 2300 2017-06-29 13:31:14Z raasch
 ! host-specific settings removed
 ! 
@@ -228,4 +231,7 @@
         ONLY:  nxl_y, nxl_yd, nxl_z, nxr_y, nxr_yd, nxr_z, nyn_x, nyn_z, nys_x,&
                nys_z, nzb_x, nzb_y, nzb_yd, nzt_x, nzt_yd, nzt_y
+
+    USE vertical_nesting_mod,                                                  &
+        ONLY:  vnested, vnest_init_pegrid_domain, vnest_init_pegrid_rank
 
     IMPLICIT NONE
@@ -335,4 +341,9 @@
 
 !
+!-- Vertical nesting: store four lists that identify partner ranks to exchange 
+!-- data
+    IF ( vnested )  CALL vnest_init_pegrid_rank
+
+!
 !-- Determine sub-topologies for transpositions
 !-- Transposition from z to x:
@@ -642,5 +653,5 @@
     CALL MPI_TYPE_COMMIT( type_xy, ierr )
 
-    IF ( TRIM( coupling_mode ) /= 'uncoupled' )  THEN
+    IF ( TRIM( coupling_mode ) /= 'uncoupled' .AND. .NOT. vnested )  THEN
    
 !
@@ -743,4 +754,8 @@
     ENDIF
 
+!
+!-- Store partner grid point co-ordinates as lists.
+!-- Create custom MPI vector datatypes for contiguous data transfer
+    IF ( vnested )  CALL vnest_init_pegrid_domain
 
 #endif
Index: palm/trunk/SOURCE/local_stop.f90
===================================================================
--- palm/trunk/SOURCE/local_stop.f90	(revision 2364)
+++ palm/trunk/SOURCE/local_stop.f90	(revision 2365)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Vertical nesting implemented (SadiqHuq)
+! 
+! 2101 2017-01-05 16:42:31Z suehring
 !
 ! 2000 2016-08-20 18:09:15Z knoop
@@ -75,4 +78,7 @@
         ONLY:  nested_run
 
+    USE vertical_nesting_mod,                                                  &
+        ONLY:  vnested
+
 #if defined( __parallel )
     IF ( coupling_mode == 'uncoupled' )  THEN
@@ -92,4 +98,11 @@
           ENDIF
        ENDIF
+    ELSEIF ( coupling_mode(1:8) == 'vnested_' )  THEN
+
+       PRINT*, '+++ local_stop:'
+       PRINT*, '     model "', TRIM( coupling_mode ), '" terminated'
+!
+!--    Abort both coarse and fine grid
+       CALL MPI_ABORT( MPI_COMM_WORLD, 9999, ierr )
     ELSE
 
Index: palm/trunk/SOURCE/parin.f90
===================================================================
--- palm/trunk/SOURCE/parin.f90	(revision 2364)
+++ palm/trunk/SOURCE/parin.f90	(revision 2365)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Vertical grid nesting: add vnest_start_time to d3par (SadiqHuq)
+! 
+! 2339 2017-08-07 13:55:26Z gronemeier
 ! corrected timestamp in header
 ! 
@@ -354,4 +357,6 @@
         ONLY:  wtm_parin
 
+    USE vertical_nesting_mod,                                                  &
+        ONLY:  vnest_start_time
 
     IMPLICIT NONE
@@ -445,5 +450,5 @@
              skip_time_do2d_xy, skip_time_do2d_xz, skip_time_do2d_yz,          &
              skip_time_do3d, skip_time_domask, synchronous_exchange,           &
-             termination_time_needed, z_max_do2d
+             termination_time_needed, vnest_start_time, z_max_do2d
 
 
Index: palm/trunk/SOURCE/read_var_list.f90
===================================================================
--- palm/trunk/SOURCE/read_var_list.f90	(revision 2364)
+++ palm/trunk/SOURCE/read_var_list.f90	(revision 2365)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Vertical grid nesting implemented (SadiqHuq)
+! 
+! 2339 2017-08-07 13:55:26Z gronemeier
 ! corrected timestamp in header
 ! 
@@ -251,4 +254,7 @@
     USE synthetic_turbulence_generator_mod,                                    &
         ONLY:  stg_read_restart_data
+
+    USE vertical_nesting_mod,                                                  &
+        ONLY:  vnest_init
 
     IMPLICIT NONE
@@ -785,5 +791,7 @@
              READ ( 13 )  vg_vertical_gradient_level_ind
           CASE ( 'virtual_flight' )
-             READ ( 13 ) virtual_flight
+             READ ( 13 )  virtual_flight
+          CASE ( 'vnest_init' )
+             READ ( 13 )  vnest_init
           CASE ( 'volume_flow_area' )
              READ ( 13 )  volume_flow_area
Index: palm/trunk/SOURCE/time_integration.f90
===================================================================
--- palm/trunk/SOURCE/time_integration.f90	(revision 2364)
+++ palm/trunk/SOURCE/time_integration.f90	(revision 2365)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Vertical grid nesting implemented (SadiqHuq)
+! 
+! 2320 2017-07-21 12:47:43Z suehring
 ! Set bottom boundary conditions after nesting interpolation and anterpolation
 ! 
@@ -399,4 +402,10 @@
         ONLY:  wind_turbine, wtm_forces
 
+    USE vertical_nesting_mod,                                                  &
+        ONLY:  vnested, vnest_anterpolate, vnest_anterpolate_e,                &
+               vnest_boundary_conds, vnest_boundary_conds_khkm,                & 
+               vnest_deallocate, vnest_init, vnest_init_fine,                  &
+               vnest_start_time
+
     IMPLICIT NONE
 
@@ -429,5 +438,5 @@
 !-- Data exchange between coupled models in case that a call has been omitted 
 !-- at the end of the previous run of a job chain.
-    IF ( coupling_mode /= 'uncoupled'  .AND.  run_coupled )  THEN
+    IF ( coupling_mode /= 'uncoupled'  .AND.  run_coupled .AND. .NOT. vnested)  THEN
 !
 !--    In case of model termination initiated by the local model the coupler 
@@ -458,4 +467,14 @@
        CALL cpu_log( log_point_s(10), 'timesteps', 'start' )
 !
+!--    Vertical nesting: initialize fine grid
+       IF ( vnested ) THEN
+          IF ( .NOT. vnest_init  .AND.  simulated_time >= vnest_start_time )  THEN
+             CALL cpu_log( log_point(80), 'vnest_init', 'start' )
+             CALL vnest_init_fine
+             vnest_init = .TRUE.
+             CALL cpu_log( log_point(80), 'vnest_init', 'stop' )
+          ENDIF
+       ENDIF
+!
 !--    Determine ug, vg and w_subs in dependence on data from external file 
 !--    LSF_DATA
@@ -624,4 +643,12 @@
 !--       Swap the time levels in preparation for the next time step.
           CALL swap_timelevel
+
+!
+!--       Vertical nesting: Interpolate fine grid data to the coarse grid
+          IF ( vnest_init ) THEN
+             CALL cpu_log( log_point(81), 'vnest_anterpolate', 'start' )
+             CALL vnest_anterpolate
+             CALL cpu_log( log_point(81), 'vnest_anterpolate', 'stop' )
+          ENDIF
 
           IF ( nested_run )  THEN
@@ -736,5 +763,28 @@
           IF ( intermediate_timestep_count == 1  .OR. &
                 call_psolver_at_all_substeps )  THEN
-             CALL pres
+
+             IF (  vnest_init ) THEN
+!
+!--             Compute pressure in the CG, interpolate top boundary conditions
+!--             to the FG and then compute pressure in the FG
+                IF ( coupling_mode == 'vnested_crse' )  CALL pres
+
+                CALL cpu_log( log_point(82), 'vnest_bc', 'start' )
+                CALL vnest_boundary_conds
+                CALL cpu_log( log_point(82), 'vnest_bc', 'stop' )
+ 
+                IF ( coupling_mode == 'vnested_fine' )  CALL pres
+
+!--             Anterpolate TKE, satisfy Germano Identity
+                CALL cpu_log( log_point(83), 'vnest_anter_e', 'start' )
+                CALL vnest_anterpolate_e
+                CALL cpu_log( log_point(83), 'vnest_anter_e', 'stop' )
+
+             ELSE
+
+                CALL pres
+
+             ENDIF
+
           ENDIF
 
@@ -910,5 +960,6 @@
 !
 !--    Data exchange between coupled models
-       IF ( coupling_mode /= 'uncoupled'  .AND.  run_coupled )  THEN
+       IF ( coupling_mode /= 'uncoupled'  .AND.  run_coupled                   &
+                                          .AND. .NOT. vnested )  THEN
           time_coupling = time_coupling + dt_3d
 
@@ -1154,4 +1205,7 @@
     ENDDO   ! time loop
 
+!-- Vertical nesting: Deallocate variables initialized for vertical nesting    
+    IF ( vnest_init )  CALL vnest_deallocate
+
     IF ( myid == 0 )  CALL finish_progress_bar
 
Index: palm/trunk/SOURCE/timestep.f90
===================================================================
--- palm/trunk/SOURCE/timestep.f90	(revision 2364)
+++ palm/trunk/SOURCE/timestep.f90	(revision 2365)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Vertical grid nesting: Sync fine and coarse grid timestep (SadiqHuq)
+! 
+! 2258 2017-06-08 07:55:13Z suehring
 ! Bugfix, add pre-preprocessor directives to enable non-parrallel mode
 !
@@ -137,4 +140,7 @@
         ONLY:  flow_statistics_called, hom, u_max, u_max_ijk, v_max, v_max_ijk,&
                w_max, w_max_ijk
+
+    USE vertical_nesting_mod,                                                  &
+        ONLY:  vnested, vnest_timestep_sync
 
     IMPLICIT NONE
@@ -376,4 +382,8 @@
     ENDIF
 
+!
+!-- Vertical nesting: coarse and fine grid timestep has to be identical    
+    IF ( vnested )  CALL vnest_timestep_sync
+
     CALL cpu_log( log_point(12), 'calculate_timestep', 'stop' )
 
Index: palm/trunk/SOURCE/vertical_nesting_mod.f90
===================================================================
--- palm/trunk/SOURCE/vertical_nesting_mod.f90	(revision 2365)
+++ palm/trunk/SOURCE/vertical_nesting_mod.f90	(revision 2365)
@@ -0,0 +1,3826 @@
+!> @file vertical_nesting_mod.f90
+!------------------------------------------------------------------------------!
+! This file is part of PALM.
+!
+! PALM is free software: you can redistribute it and/or modify it under the
+! terms of the GNU General Public License as published by the Free Software
+! Foundation, either version 3 of the License, or (at your option) any later
+! version.
+!
+! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
+! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
+! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
+!
+! You should have received a copy of the GNU General Public License along with
+! PALM. If not, see <http://www.gnu.org/licenses/>.
+!
+! Copyright 1997-2017 Leibniz Universitaet Hannover
+!------------------------------------------------------------------------------!
+!
+! Current revisions:
+! -----------------
+! 
+! 
+! Former revisions:
+! -----------------
+! $Id$
+! Initial revision (SadiqHuq)
+! 
+! 
+! 
+! 
+! Description:
+! ------------
+! Module for vertical nesting.
+!
+! Definition of parameters and variables for vertical nesting
+!------------------------------------------------------------------------------!
+ MODULE vertical_nesting_mod
+
+    USE kinds
+
+    IMPLICIT NONE
+
+    INTEGER(iwp),DIMENSION(3,2)              :: bdims=0,bdims_rem=0
+    INTEGER(iwp)                             :: cg_nprocs, fg_nprocs
+    INTEGER(iwp),DIMENSION(:,:,:),ALLOCATABLE:: c2f_dims_cg, f2c_dims_cg
+    INTEGER(iwp),DIMENSION(:),ALLOCATABLE    :: c2f_dims_fg, f2c_dims_fg
+    INTEGER(iwp)                             :: TYPE_VNEST_BC, TYPE_VNEST_ANTER
+
+    INTEGER(iwp),DIMENSION(:,:),ALLOCATABLE  :: f_rnk_lst, c_rnk_lst
+    INTEGER(iwp),DIMENSION(3)                :: cfratio
+
+    INTEGER(iwp)                             :: nxc, nxf, nyc, nyf, nzc, nzf
+    INTEGER(iwp)                             :: ngp_c, ngp_f
+
+    INTEGER(iwp)                             :: target_idex, n_cell_c, n_cell_f
+    INTEGER(iwp),DIMENSION(2)                :: pdims_partner
+    INTEGER(iwp),DIMENSION(2)                :: offset,map_coord
+
+    REAL(wp)                                :: dxc, dyc, dxf, dyf, dzc, dzf,dtc,dtf
+
+    REAL(wp),DIMENSION(:,:,:), ALLOCATABLE  :: work3d 
+    REAL(wp),DIMENSION(:,:),   ALLOCATABLE  :: work2dshf
+    REAL(wp),DIMENSION(:,:),   ALLOCATABLE  :: work2dusws
+    REAL(wp),DIMENSION(:,:),   ALLOCATABLE  :: work2dvsws
+    REAL(wp),DIMENSION(:,:),   ALLOCATABLE  :: work2dts
+    REAL(wp),DIMENSION(:,:),   ALLOCATABLE  :: work2dus
+    REAL(wp),DIMENSION(:,:),   ALLOCATABLE  :: work2dz0
+
+
+    REAL(wp), DIMENSION(:),    ALLOCATABLE  :: zuc, zuf, zwc, zwf
+    REAL(wp), DIMENSION(:,:,:), POINTER      :: interpol3d,anterpol3d 
+!   REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  interpol3d 
+
+    LOGICAL ::  vnest_init = .FALSE., vnested = .FALSE., &
+	    vnest_twi = .FALSE., vnest_couple_rk3 = .FALSE.
+
+    ! PARIN
+    REAL(wp) ::  vnest_start_time = 9999999.9
+
+    SAVE
+
+!-- Public functions
+    PUBLIC vnest_init_fine, vnest_boundary_conds, vnest_anterpolate,          &
+	   vnest_boundary_conds_khkm, vnest_anterpolate_e,                    &
+	   vnest_init_pegrid_rank, vnest_init_pegrid_domain, vnest_init_grid, &
+	   vnest_timestep_sync, vnest_deallocate
+
+!-- Public constants and variables
+    PUBLIC vnested, vnest_init, vnest_twi, vnest_couple_rk3,                  &
+           vnest_start_time
+
+    PRIVATE bdims, bdims_rem,                                                 &
+            work3d, work2dshf, work2dusws, work2dvsws,                        &
+            work2dts, work2dus, work2dz0,                                     &
+            dxc, dyc, dxf, dyf, dzc, dzf, dtc, dtf,                           &
+            zuc, zuf, zwc, zwf, interpol3d, anterpol3d,                       &
+            cg_nprocs, fg_nprocs,                                             &
+            c2f_dims_cg, c2f_dims_fg, f2c_dims_cg, f2c_dims_fg,               &
+            f_rnk_lst, c_rnk_lst, cfratio, pdims_partner,                     &
+            nxc, nxf, nyc, nyf, nzc, nzf,                                     &
+            ngp_c, ngp_f, target_idex, n_cell_c, n_cell_f,                    &
+	    offset, map_coord, TYPE_VNEST_BC, TYPE_VNEST_ANTER                
+
+    INTERFACE vnest_anterpolate
+       MODULE PROCEDURE vnest_anterpolate
+    END INTERFACE vnest_anterpolate
+
+    INTERFACE vnest_anterpolate_e
+       MODULE PROCEDURE vnest_anterpolate_e
+    END INTERFACE vnest_anterpolate_e
+
+    INTERFACE vnest_boundary_conds
+       MODULE PROCEDURE vnest_boundary_conds
+    END INTERFACE vnest_boundary_conds
+
+    INTERFACE vnest_boundary_conds_khkm
+       MODULE PROCEDURE vnest_boundary_conds_khkm
+    END INTERFACE vnest_boundary_conds_khkm
+
+    INTERFACE vnest_check_parameters
+       MODULE PROCEDURE vnest_check_parameters
+    END INTERFACE vnest_check_parameters
+
+    INTERFACE vnest_deallocate
+       MODULE PROCEDURE vnest_deallocate
+    END INTERFACE vnest_deallocate
+
+    INTERFACE vnest_init_fine
+       MODULE PROCEDURE vnest_init_fine
+    END INTERFACE vnest_init_fine
+
+    INTERFACE vnest_init_grid
+       MODULE PROCEDURE vnest_init_grid
+    END INTERFACE vnest_init_grid
+
+    INTERFACE vnest_init_pegrid_domain
+       MODULE PROCEDURE vnest_init_pegrid_domain
+    END INTERFACE vnest_init_pegrid_domain
+
+    INTERFACE vnest_init_pegrid_rank
+       MODULE PROCEDURE vnest_init_pegrid_rank
+    END INTERFACE vnest_init_pegrid_rank
+
+    INTERFACE vnest_timestep_sync
+       MODULE PROCEDURE vnest_timestep_sync
+    END INTERFACE vnest_timestep_sync
+
+ CONTAINS
+    
+    
+    
+    SUBROUTINE vnest_init_fine
+   
+        !--------------------------------------------------------------------------------!
+        ! Description:
+        ! ------------
+        ! At the specified vnest_start_time initialize the Fine Grid based on the coarse
+        ! grid values
+        !------------------------------------------------------------------------------!
+   
+   
+        USE arrays_3d
+        USE control_parameters
+        USE grid_variables
+        USE indices
+        USE interfaces
+        USE pegrid
+        USE surface_mod,                                                        &
+            ONLY :  surf_def_h, surf_def_v
+        
+   
+        IMPLICIT NONE
+   
+        REAL(wp)    ::  time_since_reference_point_rem
+   
+        INTEGER(iwp)   :: i, j, k,im,jn,ko
+        INTEGER(iwp)                             ::  if, jf, kf
+   
+#if defined( __parallel )
+   
+        if (myid ==0 )print *, ' TIME TO INIT FINE from COARSE', simulated_time
+   
+        !
+        !-- In case of model termination initiated by the remote model
+        !-- (terminate_coupled_remote > 0), initiate termination of the local model.
+        !-- The rest of the coupler must then be skipped because it would cause an MPI
+        !-- intercomminucation hang.
+        !-- If necessary, the coupler will be called at the beginning of the next
+        !-- restart run.
+   
+        IF ( myid == 0) THEN
+            CALL MPI_SENDRECV( terminate_coupled,        1, MPI_INTEGER, &
+                target_id, 0,                             &
+                terminate_coupled_remote, 1, MPI_INTEGER, &
+                target_id, 0,                             &
+                comm_inter, status, ierr )
+        ENDIF
+        CALL MPI_BCAST( terminate_coupled_remote, 1, MPI_INTEGER, 0, comm2d, &
+            ierr )
+   
+        IF ( terminate_coupled_remote > 0 )  THEN
+            WRITE( message_string, * ) 'remote model "',                         &
+                TRIM( coupling_mode_remote ),             &
+                '" terminated',                           &
+                '&with terminate_coupled_remote = ',      &
+                terminate_coupled_remote,                 &
+                '&local model  "', TRIM( coupling_mode ), &
+                '" has',                                  &
+                '&terminate_coupled = ',                  &
+                terminate_coupled
+            CALL message( 'vnest_init_fine', 'PA0310', 1, 2, 0, 6, 0 )
+            RETURN
+        ENDIF
+   
+   
+        !
+        !-- Exchange the current simulated time between the models,
+        !-- currently just for total_2ding
+        IF ( myid == 0 ) THEN
+   
+            CALL MPI_SEND( time_since_reference_point, 1, MPI_REAL, target_id, &
+                11, comm_inter, ierr )
+            CALL MPI_RECV( time_since_reference_point_rem, 1, MPI_REAL,        &
+                target_id, 11, comm_inter, status, ierr )
+   
+        ENDIF
+   
+        CALL MPI_BCAST( time_since_reference_point_rem, 1, MPI_REAL, 0, comm2d, &
+            ierr )
+   
+   
+        IF ( coupling_mode == 'vnested_crse' )  THEN
+!-- Send data to fine grid for initialization
+   
+            offset(1) = ( pdims_partner(1) / pdims(1) ) * pcoord(1)
+            offset(2) = ( pdims_partner(2) / pdims(2) ) * pcoord(2)
+   
+            do j = 0,   ( pdims_partner(2) / pdims(2) ) - 1
+                do i = 0,   ( pdims_partner(1) / pdims(1) ) - 1
+                    map_coord(1) = i+offset(1)
+                    map_coord(2) = j+offset(2)
+   
+                    target_idex = f_rnk_lst(map_coord(1),map_coord(2)) + numprocs
+   
+                    CALL MPI_RECV( bdims_rem,   6, MPI_INTEGER, target_idex, 10, &
+                        comm_inter,status, ierr )
+   
+                    bdims (1,1) = bdims_rem (1,1) / cfratio(1)
+                    bdims (1,2) = bdims_rem (1,2) / cfratio(1)
+                    bdims (2,1) = bdims_rem (2,1) / cfratio(2)
+                    bdims (2,2) = bdims_rem (2,2) / cfratio(2)
+                    bdims (3,1) = bdims_rem (3,1)
+                    bdims (3,2) = bdims_rem (3,2) / cfratio(3)
+   
+   
+                    CALL MPI_SEND( bdims, 6, MPI_INTEGER, target_idex, 9, &
+                        comm_inter, ierr )
+   
+   
+                    n_cell_c = (bdims(1,2)-bdims(1,1)+3) * &
+                        (bdims(2,2)-bdims(2,1)+3) * &
+                        (bdims(3,2)-bdims(3,1)+3)
+   
+                    CALL MPI_SEND( u( bdims(3,1):bdims(3,2)+2, &
+                        bdims(2,1)-1:bdims(2,2)+1, &
+                        bdims(1,1)-1:bdims(1,2)+1),&
+                        n_cell_c, MPI_REAL, target_idex,  &
+                        101, 	comm_inter, ierr)
+   
+                    CALL MPI_SEND( v( bdims(3,1):bdims(3,2)+2, &
+                        bdims(2,1)-1:bdims(2,2)+1, &
+                        bdims(1,1)-1:bdims(1,2)+1),&
+                        n_cell_c, MPI_REAL, target_idex,  &
+                        102, 	comm_inter, ierr)
+   
+                    CALL MPI_SEND( w( bdims(3,1):bdims(3,2)+2, &
+                        bdims(2,1)-1:bdims(2,2)+1, &
+                        bdims(1,1)-1:bdims(1,2)+1),&
+                        n_cell_c, MPI_REAL, target_idex,  &
+                        103, 	comm_inter, ierr)
+   
+                    CALL MPI_SEND( pt(bdims(3,1):bdims(3,2)+2, &
+                        bdims(2,1)-1:bdims(2,2)+1, &
+                        bdims(1,1)-1:bdims(1,2)+1),&
+                        n_cell_c, MPI_REAL, target_idex,  &
+                        105, 	comm_inter, ierr)
+   
+            IF ( humidity )  THEN
+                    CALL MPI_SEND( q(bdims(3,1):bdims(3,2)+2, &
+                        bdims(2,1)-1:bdims(2,2)+1, &
+                        bdims(1,1)-1:bdims(1,2)+1),&
+                        n_cell_c, MPI_REAL, target_idex,  &
+                        116, 	comm_inter, ierr)
+            ENDIF
+ 
+                     CALL MPI_SEND( e( bdims(3,1):bdims(3,2)+2, &
+                        bdims(2,1)-1:bdims(2,2)+1, &
+                        bdims(1,1)-1:bdims(1,2)+1),&
+                        n_cell_c, MPI_REAL, target_idex,  &
+                        104, 	comm_inter, ierr)
+   
+                   CALL MPI_SEND(kh( bdims(3,1):bdims(3,2)+2, &
+                        bdims(2,1)-1:bdims(2,2)+1, &
+                        bdims(1,1)-1:bdims(1,2)+1),&
+                        n_cell_c, MPI_REAL, target_idex,  &
+                        106, 	comm_inter, ierr)
+   
+                    CALL MPI_SEND(km( bdims(3,1):bdims(3,2)+2, &
+                        bdims(2,1)-1:bdims(2,2)+1, &
+                        bdims(1,1)-1:bdims(1,2)+1),&
+                        n_cell_c, MPI_REAL, target_idex,  &
+                        107, 	comm_inter, ierr)
+   
+!-- Send Surface fluxes
+            IF ( use_surface_fluxes )  THEN
+   
+                   n_cell_c = (bdims(1,2)-bdims(1,1)+3) * &
+                        (bdims(2,2)-bdims(2,1)+3)
+   
+!-- WARNING
+!-- shf,z0 not interpolated
+!-- line commented in interpolate_to_fine_flux
+!-- FG needs to read it's own data file   
+!MERGE-WIP                  CALL MPI_SEND(surf_def_h(0)%shf ( bdims(2,1)-1:bdims(2,2)+1, &
+!MERGE-WIP                      bdims(1,1)-1:bdims(1,2)+1),&
+!MERGE-WIP                      n_cell_c, MPI_REAL, target_idex,  &
+!MERGE-WIP                      109, 	comm_inter, ierr   )
+!MERGE-WIP 
+!MERGE-WIP                  CALL MPI_SEND(surf_def_h(0)%usws( bdims(2,1)-1:bdims(2,2)+1, &
+!MERGE-WIP                      bdims(1,1)-1:bdims(1,2)+1),&
+!MERGE-WIP                      n_cell_c, MPI_REAL, target_idex,  &
+!MERGE-WIP                      110, 	comm_inter, ierr   )
+!MERGE-WIP 
+!MERGE-WIP                  CALL MPI_SEND(surf_def_h(0)%vsws( bdims(2,1)-1:bdims(2,2)+1, &
+!MERGE-WIP                      bdims(1,1)-1:bdims(1,2)+1),&
+!MERGE-WIP                      n_cell_c, MPI_REAL, target_idex,  &
+!MERGE-WIP                      111, 	comm_inter, ierr   )
+!MERGE-WIP 
+!MERGE                   CALL MPI_SEND(ts  ( bdims(2,1)-1:bdims(2,2)+1, &
+!MERGE                       bdims(1,1)-1:bdims(1,2)+1),&
+!MERGE                       n_cell_c, MPI_REAL, target_idex,  &
+!MERGE                       112, 	comm_inter, ierr   )
+!MERGE  
+!MERGE                   CALL MPI_SEND(us  ( bdims(2,1)-1:bdims(2,2)+1, &
+!MERGE                       bdims(1,1)-1:bdims(1,2)+1),&
+!MERGE                       n_cell_c, MPI_REAL, target_idex,  &
+!MERGE                       113, 	comm_inter, ierr   )
+!MERGE  
+!MERGE                   CALL MPI_SEND(z0  ( bdims(2,1)-1:bdims(2,2)+1, &
+!MERGE                       bdims(1,1)-1:bdims(1,2)+1),&
+!MERGE                       n_cell_c, MPI_REAL, target_idex,  &
+!MERGE                       114, 	comm_inter, ierr   )
+            ENDIF
+   
+   
+   
+   
+                end do
+            end do
+   
+        ELSEIF ( coupling_mode == 'vnested_fine' )  THEN
+!-- Receive data from coarse grid for initialization
+   
+            offset(1) = pcoord(1) / ( pdims(1)/pdims_partner(1) )
+            offset(2) = pcoord(2) / ( pdims(2)/pdims_partner(2) )
+            map_coord(1) = offset(1)
+            map_coord(2) = offset(2)
+            target_idex = c_rnk_lst(map_coord(1),map_coord(2))
+   
+            bdims (1,1) = nxl
+            bdims (1,2) = nxr
+            bdims (2,1) = nys
+            bdims (2,2) = nyn
+            bdims (3,1) = nzb
+            bdims (3,2) = nzt
+   
+            CALL MPI_SEND( bdims,       6, MPI_INTEGER, target_idex, 10, &
+                comm_inter, ierr )
+   
+            CALL MPI_RECV( bdims_rem,   6, MPI_INTEGER, target_idex, 9, &
+                comm_inter,status, ierr )
+   
+            n_cell_c = (bdims_rem(1,2)-bdims_rem(1,1)+3) * &
+                (bdims_rem(2,2)-bdims_rem(2,1)+3) * &
+                (bdims_rem(3,2)-bdims_rem(3,1)+3)
+   
+            ALLOCATE( work3d ( bdims_rem(3,1)  :bdims_rem(3,2)+2, &
+                bdims_rem(2,1)-1:bdims_rem(2,2)+1, &
+                bdims_rem(1,1)-1:bdims_rem(1,2)+1))
+   
+   
+            CALL MPI_RECV( work3d,n_cell_c, MPI_REAL, target_idex, 101, &
+                comm_inter,status, ierr )
+            interpol3d => u
+            call interpolate_to_fine_u ( 101 )
+   
+            CALL MPI_RECV( work3d,n_cell_c, MPI_REAL, target_idex, 102, &
+                comm_inter,status, ierr )
+            interpol3d => v
+            call interpolate_to_fine_v ( 102 )
+   
+            CALL MPI_RECV( work3d,n_cell_c, MPI_REAL, target_idex, 103, &
+                comm_inter,status, ierr )
+            interpol3d => w
+            call interpolate_to_fine_w ( 103 )
+   
+            CALL MPI_RECV( work3d,n_cell_c, MPI_REAL, target_idex, 105, &
+                comm_inter,status, ierr )
+            interpol3d => pt
+            call interpolate_to_fine_s ( 105 )
+   
+            IF ( humidity )  THEN
+            CALL MPI_RECV( work3d,n_cell_c, MPI_REAL, target_idex, 116, &
+                comm_inter,status, ierr )
+            interpol3d => q
+            call interpolate_to_fine_s ( 116 )
+            ENDIF
+
+            CALL MPI_RECV( work3d,n_cell_c, MPI_REAL, target_idex, 104, &
+                comm_inter,status, ierr )
+            interpol3d => e
+            call interpolate_to_fine_s ( 104 )
+   
+!-- kh,km no target attribute, use of pointer not possible
+            CALL MPI_RECV( work3d,n_cell_c, MPI_REAL, target_idex, 106, &
+                comm_inter,status, ierr )
+            call interpolate_to_fine_kh ( 106 )
+   
+            CALL MPI_RECV( work3d,n_cell_c, MPI_REAL, target_idex, 107, &
+                comm_inter,status, ierr )
+            call interpolate_to_fine_km ( 107 )
+   
+            DEALLOCATE(   work3d       )
+            NULLIFY   (   interpol3d   )
+   
+!-- Recv Surface Fluxes	    
+            IF ( use_surface_fluxes )  THEN
+            n_cell_c = (bdims_rem(1,2)-bdims_rem(1,1)+3) * &
+                (bdims_rem(2,2)-bdims_rem(2,1)+3)
+   
+            ALLOCATE( work2dshf  ( bdims_rem(2,1)-1:bdims_rem(2,2)+1, &
+                bdims_rem(1,1)-1:bdims_rem(1,2)+1) )
+            ALLOCATE( work2dusws ( bdims_rem(2,1)-1:bdims_rem(2,2)+1, &
+                bdims_rem(1,1)-1:bdims_rem(1,2)+1) )
+            ALLOCATE( work2dvsws ( bdims_rem(2,1)-1:bdims_rem(2,2)+1, &
+                bdims_rem(1,1)-1:bdims_rem(1,2)+1) )
+            ALLOCATE( work2dts   ( bdims_rem(2,1)-1:bdims_rem(2,2)+1, &
+                bdims_rem(1,1)-1:bdims_rem(1,2)+1) )
+            ALLOCATE( work2dus   ( bdims_rem(2,1)-1:bdims_rem(2,2)+1, &
+                bdims_rem(1,1)-1:bdims_rem(1,2)+1) )
+            ALLOCATE( work2dz0   ( bdims_rem(2,1)-1:bdims_rem(2,2)+1, &
+                bdims_rem(1,1)-1:bdims_rem(1,2)+1) )
+   
+!MERGE-WIP          CALL MPI_RECV( work2dshf ,n_cell_c, MPI_REAL, target_idex, 109, &
+!MERGE-WIP              comm_inter,status, ierr )
+!MERGE-WIP 
+!MERGE-WIP          CALL MPI_RECV( work2dusws,n_cell_c, MPI_REAL, target_idex, 110, &
+!MERGE-WIP              comm_inter,status, ierr )
+!MERGE-WIP 
+!MERGE-WIP          CALL MPI_RECV( work2dvsws,n_cell_c, MPI_REAL, target_idex, 111, &
+!MERGE-WIP              comm_inter,status, ierr )
+!MERGE-WIP 
+!MERGE           CALL MPI_RECV( work2dts  ,n_cell_c, MPI_REAL, target_idex, 112, &
+!MERGE               comm_inter,status, ierr )
+!MERGE  
+!MERGE           CALL MPI_RECV( work2dus  ,n_cell_c, MPI_REAL, target_idex, 113, &
+!MERGE               comm_inter,status, ierr )
+!MERGE  
+!MERGE           CALL MPI_RECV( work2dz0  ,n_cell_c, MPI_REAL, target_idex, 114, &
+!MERGE               comm_inter,status, ierr )
+!MERGE  
+!MERGE-WIP           CALL interpolate_to_fine_flux ( 108 )
+   
+            DEALLOCATE(   work2dshf    )
+            DEALLOCATE(   work2dusws   )
+            DEALLOCATE(   work2dvsws   )
+            DEALLOCATE(   work2dts     )
+            DEALLOCATE(   work2dus     )
+            DEALLOCATE(   work2dz0     )
+          ENDIF
+   
+          IF ( .NOT. constant_diffusion ) THEN
+             DO kf =  bdims(3,1)+1,bdims(3,2)+1
+                 DO jf =  bdims(2,1),bdims(2,2)
+                     DO if =  bdims(1,1),bdims(1,2)
+   
+                         IF ( e(kf,jf,if) < 0.0 ) THEN
+                              e(kf,jf,if) = 1E-15_wp
+                         END IF
+   
+                     END DO
+                 END DO
+             END DO
+          ENDIF
+   
+           w(nzt+1,:,:) = w(nzt,:,:)
+          
+           CALL exchange_horiz( u, nbgp )
+           CALL exchange_horiz( v, nbgp )
+           CALL exchange_horiz( w, nbgp )
+           CALL exchange_horiz( pt, nbgp )
+           IF ( .NOT. constant_diffusion )  CALL exchange_horiz( e, nbgp )
+           IF ( humidity )  CALL exchange_horiz( q, nbgp )
+   
+           !
+           !--         Velocity boundary conditions at the bottom boundary
+           IF ( ibc_uv_b == 0 ) THEN
+               u(nzb,:,:) = 0.0_wp
+               v(nzb,:,:) = 0.0_wp
+           ELSE
+               u(nzb,:,:) = u(nzb+1,:,:)
+               v(nzb,:,:) = v(nzb+1,:,:)
+           END IF
+           
+           
+           w(nzb,:,:) = 0.0_wp
+
+	   !
+           !-- Temperature boundary conditions at the bottom boundary
+           IF ( ibc_pt_b /= 0 ) THEN
+               pt(nzb,:,:) = pt(nzb+1,:,:)
+           END IF
+           
+           !
+           !-- Bottom boundary condition for the turbulent kinetic energy
+           !-- Generally a Neumann condition with de/dz=0 is assumed
+           IF ( .NOT. constant_diffusion ) THEN
+               e(nzb,:,:) = e(nzb+1,:,:)
+           END IF
+           
+           !
+           !-- Bottom boundary condition for turbulent diffusion coefficients
+           km(nzb,:,:) = km(nzb+1,:,:)
+           kh(nzb,:,:) = kh(nzb+1,:,:)
+           
+           !diffusivities required
+           IF ( .NOT. humidity ) THEN
+               CALL diffusivities( pt, pt_reference )
+           ELSE
+               CALL diffusivities( vpt, pt_reference )
+           ENDIF
+   
+
+!
+!-- Reset Fine Grid top Boundary Condition
+!-- At the top of the FG, the scalars always follow Dirichlet condition
+
+            ibc_pt_t = 0
+
+!-- Initialize old time levels
+            pt_p = pt; u_p = u; v_p = v; w_p = w
+            IF ( .NOT. constant_diffusion ) e_p = e
+            IF ( humidity ) THEN
+               ibc_q_t = 0
+               q_p     = q
+            ENDIF
+   
+        ENDIF
+   
+
+        if (myid==0) print *, '** Fine Initalized ** simulated_time:', simulated_time
+#endif
+    CONTAINS
+   
+       SUBROUTINE interpolate_to_fine_w( tag )
+      
+#if defined( __parallel )
+      
+           USE arrays_3d
+           USE control_parameters
+           USE grid_variables
+           USE indices
+           USE pegrid
+           
+      
+           IMPLICIT NONE
+      
+           INTEGER(iwp), intent(in) ::  tag
+           INTEGER(iwp)                             ::  i, j, k
+           INTEGER(iwp)                             ::  if, jf, kf
+           INTEGER(iwp)                             ::  bottomx, topx
+           INTEGER(iwp)                             ::  bottomy, topy
+           INTEGER(iwp)                             ::  bottomz, topz
+           REAL(wp)                                ::  eps, alpha, eminus, edot, eplus
+           REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  wprs, wprf
+      
+           INTEGER(iwp) :: nzbottom, nztop
+      
+      
+           nzbottom = bdims_rem (3,1)
+           nztop  = bdims_rem (3,2)
+      
+           ALLOCATE( wprf(nzbottom:nztop, bdims_rem(2,1)-1: bdims_rem(2,2)+1,nxl:nxr) )
+           ALLOCATE( wprs(nzbottom:nztop,nys:nyn,nxl:nxr) )
+      
+      
+           !
+           !-- Initialisation of the velocity component w
+           !
+           !-- Interpolation in x-direction
+           DO k = nzbottom, nztop
+               DO j = bdims_rem(2,1)-1, bdims_rem(2,2)+1
+                   DO i = bdims_rem(1,1),bdims_rem(1,2)
+      
+                       bottomx = (nxf+1)/(nxc+1) * i
+                       topx    = (nxf+1)/(nxc+1) * (i+1) - 1
+      
+                       DO if = bottomx, topx
+      
+                           eps    = ( if * dxf + 0.5 * dxf - i * dxc - 0.5 * dxc ) / dxc
+                           alpha  = ( ( dxf / dxc )**2.0 - 1.0 ) / 24.0
+                           eminus = eps * ( eps - 1.0 ) / 2.0 + alpha
+                           edot  = ( 1.0 - eps**2.0 ) - 2.0 * alpha
+                           eplus  = eps * ( eps + 1.0 ) / 2.0 + alpha
+                      
+                           wprf(k,j,if) = eminus * work3d(k,j,i-1) &
+                               + edot  * work3d(k,j,i)   &
+                               + eplus  * work3d(k,j,i+1)
+                       END DO
+      
+                   END DO
+               END DO
+           END DO
+      
+           !
+           !-- Interpolation in y-direction (quadratic, Clark and Farley)
+           DO k = nzbottom, nztop
+               DO j = bdims_rem(2,1), bdims_rem(2,2)
+                
+                   bottomy = (nyf+1)/(nyc+1) * j
+                   topy    = (nyf+1)/(nyc+1) * (j+1) - 1
+      
+                   DO if = nxl, nxr
+                       DO jf = bottomy, topy
+                      
+                           eps    = ( jf * dyf + 0.5 * dyf - j * dyc - 0.5 * dyc ) / dyc
+                           alpha  = ( ( dyf / dyc )**2.0 - 1.0 ) / 24.0
+                           eminus = eps * ( eps - 1.0 ) / 2.0 + alpha
+                           edot  = ( 1.0 - eps**2.0 ) - 2.0 * alpha
+                           eplus  = eps * ( eps + 1.0 ) / 2.0 + alpha
+      
+                           wprs(k,jf,if) = eminus * wprf(k,j-1,if) &
+                               + edot  * wprf(k,j,if)   &
+                               + eplus  * wprf(k,j+1,if)
+      
+                       END DO
+                   END DO
+      
+               END DO
+           END DO
+      
+           !
+           !-- Interpolation in z-direction (linear)
+      
+           DO k = nzbottom, nztop-1
+      
+               bottomz = (dzc/dzf) * k
+               topz    = (dzc/dzf) * (k+1) - 1
+      
+               DO jf = nys, nyn
+                   DO if = nxl, nxr
+                       DO kf = bottomz, topz
+      
+                           w(kf,jf,if) = wprs(k,jf,if) + ( zwf(kf) - zwc(k) ) &
+                               * ( wprs(k+1,jf,if) - wprs(k,jf,if) ) / dzc
+      
+                       END DO
+                   END DO
+               END DO
+      
+           END DO
+      
+           DO jf = nys, nyn
+               DO if = nxl, nxr
+      
+                   w(nzt,jf,if) = wprs(nztop,jf,if)
+       
+               END DO
+           END DO
+           !
+           !    w(nzb:nzt+1,nys:nyn,nxl:nxr) = 0
+      
+           DEALLOCATE( wprf, wprs )
+      
+#endif
+       END SUBROUTINE interpolate_to_fine_w
+      
+       SUBROUTINE interpolate_to_fine_u( tag )
+      
+#if defined( __parallel )
+      
+           USE arrays_3d
+           USE control_parameters
+           USE grid_variables
+           USE indices
+           USE pegrid
+           
+      
+           IMPLICIT NONE
+      
+           INTEGER(iwp), intent(in) ::  tag
+           INTEGER(iwp)                             ::  i, j, k
+           INTEGER(iwp)                             ::  if, jf, kf
+           INTEGER(iwp)                             ::  bottomx, topx
+           INTEGER(iwp)                             ::  bottomy, topy
+           INTEGER(iwp)                             ::  bottomz, topz
+           REAL(wp)                                ::  eps, alpha, eminus, edot, eplus
+           REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  uprs, uprf
+      
+           INTEGER(iwp) :: nzbottom, nztop
+      
+      
+           nzbottom = bdims_rem (3,1)
+           nztop  = bdims_rem (3,2)
+      
+           ALLOCATE( uprf(nzbottom:nztop+2,nys:nyn,bdims_rem(1,1)-1:bdims_rem(1,2)+1) )
+           ALLOCATE( uprs(nzb+1:nzt+1,nys:nyn,bdims_rem(1,1)-1:bdims_rem(1,2)+1) )
+      
+           !
+           !-- Initialisation of the velocity component uf
+      
+           !
+           !-- Interpolation in y-direction (quadratic, Clark and Farley)
+      
+           DO k = nzbottom, nztop+2
+               DO j = bdims_rem(2,1), bdims_rem(2,2)
+      
+                   bottomy = (nyf+1)/(nyc+1) * j
+                   topy    = (nyf+1)/(nyc+1) * (j+1) - 1
+      
+                   DO i = bdims_rem(1,1)-1, bdims_rem(1,2)+1
+                       DO jf = bottomy, topy
+      
+                           eps    = ( jf * dyf + 0.5 * dyf - j * dyc - 0.5 * dyc ) / dyc
+                           alpha  = ( ( dyf / dyc )**2.0 - 1.0 ) / 24.0
+                           eminus = eps * ( eps - 1.0 ) / 2.0 + alpha
+                           edot  = ( 1.0 - eps**2.0 ) - 2.0 * alpha
+                           eplus  = eps * ( eps + 1.0 ) / 2.0 + alpha
+      
+                           uprf(k,jf,i) = eminus * work3d(k,j-1,i) &
+                               + edot  * work3d(k,j,i)   &
+                               + eplus  * work3d(k,j+1,i)
+      
+                       END DO
+                   END DO
+      
+               END DO
+           END DO
+      
+           !
+           !-- Interpolation in z-direction (quadratic, Clark and Farley)
+      
+           DO k = nzbottom+1, nztop
+      
+               bottomz = (dzc/dzf) * (k-1) + 1
+               topz    = (dzc/dzf) * k
+      
+               DO jf = nys, nyn
+                   DO i = bdims_rem(1,1)-1, bdims_rem(1,2)+1
+                       DO kf = bottomz, topz
+                      
+                           eps    = ( zuf(kf) - zuc(k) ) / dzc
+                           alpha  = ( ( dzf / dzc )**2.0 - 1.0 ) / 24.0
+                           eminus = eps * ( eps - 1.0 ) / 2.0 + alpha
+                           edot  = ( 1.0 - eps**2.0 ) - 2.0 * alpha
+                           eplus  = eps * ( eps + 1.0 ) / 2.0 + alpha
+      
+                           uprs(kf,jf,i) = eminus * uprf(k-1,jf,i) &
+                               + edot  * uprf(k,jf,i)   &
+                               + eplus  * uprf(k+1,jf,i)
+      
+                       END DO
+                   END DO
+               END DO
+      
+           END DO
+      
+           DO jf = nys, nyn
+               DO i = bdims_rem(1,1)-1, bdims_rem(1,2)+1
+      
+                   eps    = ( zuf(nzt+1) - zuc(nztop+1) ) / dzc
+                   alpha  = ( ( dzf / dzc )**2.0 - 1.0 ) / 24.0
+                   eminus = eps * ( eps - 1.0 ) / 2.0 + alpha
+                   edot  = ( 1.0 - eps**2.0 ) - 2.0 * alpha
+                   eplus  = eps * ( eps + 1.0 ) / 2.0 + alpha
+      
+                   uprs(nzt+1,jf,i)  = eminus * uprf(nztop,jf,i)   &
+                       + edot  * uprf(nztop+1,jf,i) &
+                       + eplus  * uprf(nztop+2,jf,i)
+      
+               END DO
+           END DO
+      
+           !
+           !-- Interpolation in x-direction (linear)
+      
+           DO kf = nzb+1, nzt+1
+               DO jf = nys, nyn
+                   DO i = bdims_rem(1,1), bdims_rem(1,2)
+      
+                       bottomx = (nxf+1)/(nxc+1) * i
+                       topx    = (nxf+1)/(nxc+1) * (i+1) - 1
+      
+                       DO if = bottomx, topx
+                           u(kf,jf,if)  = uprs(kf,jf,i) + ( if * dxf - i * dxc ) &
+                               * ( uprs(kf,jf,i+1) - uprs(kf,jf,i) ) / dxc
+                       END DO
+      
+                   END DO
+               END DO
+           END DO
+           !
+           !-- Determination of uf at the bottom boundary
+      
+      
+      
+           DEALLOCATE( uprf, uprs )
+      
+#endif
+       END SUBROUTINE interpolate_to_fine_u
+      
+      
+       SUBROUTINE interpolate_to_fine_v( tag )
+      
+#if defined( __parallel )
+      
+           USE arrays_3d
+           USE control_parameters
+           USE grid_variables
+           USE indices
+           USE pegrid
+           
+      
+           IMPLICIT NONE
+          
+           INTEGER(iwp), intent(in) ::  tag
+           INTEGER(iwp)                             ::  i, j, k
+           INTEGER(iwp)                             ::  if, jf, kf
+           INTEGER(iwp)                             ::  bottomx, topx
+           INTEGER(iwp)                             ::  bottomy, topy
+           INTEGER(iwp)                             ::  bottomz, topz
+           REAL(wp)                                ::  eps, alpha, eminus, edot, eplus
+           REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  vprs, vprf
+      
+           INTEGER(iwp) :: nzbottom, nztop
+      
+      
+           nzbottom = bdims_rem (3,1)
+           nztop  = bdims_rem (3,2)
+      
+           ALLOCATE( vprf(nzbottom:nztop+2,bdims_rem(2,1)-1:bdims_rem(2,2)+1,nxl:nxr) )
+           ALLOCATE( vprs(nzb+1:nzt+1, bdims_rem(2,1)-1:bdims_rem(2,2)+1,nxl:nxr) )
+           !
+           !-- Initialisation of the velocity component vf
+      
+           !
+           !-- Interpolation in x-direction (quadratic, Clark and Farley)
+      
+           DO k = nzbottom, nztop+2
+               DO j = bdims_rem(2,1)-1, bdims_rem(2,2)+1
+                   DO i = bdims_rem(1,1), bdims_rem(1,2)
+      
+                       bottomx = (nxf+1)/(nxc+1) * i
+                       topx    = (nxf+1)/(nxc+1) * (i+1) - 1
+      
+                       DO if = bottomx, topx
+      
+                           eps    = ( if * dxf + 0.5 * dxf - i * dxc - 0.5 * dxc ) / dxc
+                           alpha  = ( ( dxf / dxc )**2.0 - 1.0 ) / 24.0
+                           eminus = eps * ( eps - 1.0 ) / 2.0 + alpha
+                           edot  = ( 1.0 - eps**2.0 ) - 2.0 * alpha
+                           eplus  = eps * ( eps + 1.0 ) / 2.0 + alpha
+      
+                           vprf(k,j,if) = eminus * work3d(k,j,i-1) &
+                               + edot  * work3d(k,j,i)   &
+                               + eplus  * work3d(k,j,i+1)
+      
+                       END DO
+      
+                   END DO
+               END DO
+           END DO
+      
+           !
+           !-- Interpolation in z-direction (quadratic, Clark and Farley)
+      
+           DO k = nzbottom+1, nztop
+             
+               bottomz = (dzc/dzf) * (k-1) + 1
+               topz    = (dzc/dzf) * k
+      
+               DO j = bdims_rem(2,1)-1, bdims_rem(2,2)+1
+                   DO if = nxl, nxr
+                       DO kf = bottomz, topz
+      
+                           eps    = ( zuf(kf) - zuc(k) ) / dzc
+                           alpha  = ( ( dzf / dzc )**2.0 - 1.0 ) / 24.0
+                           eminus = eps * ( eps - 1.0 ) / 2.0 + alpha
+                           edot  = ( 1.0 - eps**2.0 ) - 2.0 * alpha
+                           eplus  = eps * ( eps + 1.0 ) / 2.0 + alpha
+      
+                           vprs(kf,j,if) = eminus * vprf(k-1,j,if) &
+                               + edot  * vprf(k,j,if)   &
+                               + eplus  * vprf(k+1,j,if)
+      
+                       END DO
+                   END DO
+               END DO
+      
+           END DO
+             
+           DO j = bdims_rem(2,1)-1, bdims_rem(2,2)+1
+               DO if = nxl, nxr
+      
+                   eps    = ( zuf(nzt+1) - zuc(nztop+1) ) / dzc
+                   alpha  = ( ( dzf / dzc )**2.0 - 1.0 ) / 24.0
+                   eminus = eps * ( eps - 1.0 ) / 2.0 + alpha
+                   edot  = ( 1.0 - eps**2.0 ) - 2.0 * alpha
+                   eplus  = eps * ( eps + 1.0 ) / 2.0 + alpha
+      
+                   vprs(nzt+1,j,if)  = eminus * vprf(nztop,j,if)   &
+                       + edot  * vprf(nztop+1,j,if) &
+                       + eplus  * vprf(nztop+2,j,if)
+      
+               END DO
+           END DO
+      
+           !
+           !-- Interpolation in y-direction (linear)
+      
+           DO kf = nzb+1, nzt+1
+               DO j = bdims_rem(2,1), bdims_rem(2,2)
+      
+                   bottomy = (nyf+1)/(nyc+1) * j
+                   topy    = (nyf+1)/(nyc+1) * (j+1) - 1
+      
+                   DO if = nxl, nxr
+                       DO jf = bottomy, topy
+                           v (kf,jf,if) = vprs(kf,j,if) + ( jf * dyf - j * dyc ) &
+                               * ( vprs(kf,j+1,if) - vprs(kf,j,if) ) / dyc
+                       END DO
+                   END DO
+        
+               END DO
+           END DO
+      
+           !
+           !-- Determination of vf at the bottom boundary
+      
+      
+           DEALLOCATE( vprf, vprs )
+      
+#endif
+       END SUBROUTINE interpolate_to_fine_v
+      
+      
+       SUBROUTINE interpolate_to_fine_s( tag )
+      
+#if defined( __parallel )
+      
+           USE arrays_3d
+           USE control_parameters
+           USE grid_variables
+           USE indices
+           USE pegrid
+           
+      
+           IMPLICIT NONE
+      
+          
+           INTEGER(iwp), intent(in) ::  tag
+           INTEGER(iwp)                             ::  i, j, k
+           INTEGER(iwp)                             ::  if, jf, kf
+           INTEGER(iwp)                             ::  bottomx, topx
+           INTEGER(iwp)                             ::  bottomy, topy
+           INTEGER(iwp)                             ::  bottomz, topz
+           REAL(wp)                                ::  eps, alpha, eminus, edot, eplus
+           REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  ptprs, ptprf
+      
+           INTEGER(iwp) :: nzbottom, nztop
+      
+      
+           nzbottom = bdims_rem (3,1)
+           nztop  = bdims_rem (3,2)
+      
+           ALLOCATE( ptprf(nzbottom:nztop+2,bdims_rem(2,1)-1:bdims_rem(2,2)+1,nxl:nxr) )
+           ALLOCATE( ptprs(nzbottom:nztop+2,nys:nyn,nxl:nxr) )
+      
+           !
+           !-- Initialisation of scalar variables
+      
+           !
+           !-- Interpolation in x-direction (quadratic, Clark and Farley)
+      
+           DO k = nzbottom, nztop+2
+               DO j = bdims_rem(2,1)-1, bdims_rem(2,2)+1
+                   DO i = bdims_rem(1,1), bdims_rem(1,2)
+      
+                       bottomx = (nxf+1)/(nxc+1) * i
+                       topx    = (nxf+1)/(nxc+1) * (i+1) - 1
+      
+                       DO if = bottomx, topx
+      
+                           eps    = ( if * dxf + 0.5 * dxf - i * dxc - 0.5 * dxc ) / dxc
+                           alpha  = ( ( dxf / dxc )**2.0 - 1.0 ) / 24.0
+                           eminus = eps * ( eps - 1.0 ) / 2.0 + alpha
+                           edot  = ( 1.0 - eps**2.0 ) - 2.0 * alpha
+                           eplus  = eps * ( eps + 1.0 ) / 2.0 + alpha
+      
+                           ptprf(k,j,if) = eminus * work3d(k,j,i-1) &
+                               + edot  * work3d(k,j,i)   &
+                               + eplus  * work3d(k,j,i+1)
+                       END DO
+      
+                   END DO
+               END DO
+           END DO
+      
+           !
+           !-- Interpolation in y-direction (quadratic, Clark and Farley)
+      
+           DO k = nzbottom, nztop+2
+               DO j = bdims_rem(2,1), bdims_rem(2,2)
+        
+                   bottomy = (nyf+1)/(nyc+1) * j
+                   topy    = (nyf+1)/(nyc+1) * (j+1) - 1
+      
+                   DO if = nxl, nxr
+                       DO jf = bottomy, topy
+      
+                           eps    = ( jf * dyf + 0.5 * dyf - j * dyc - 0.5 * dyc ) / dyc
+                           alpha  = ( ( dyf / dyc )**2.0 - 1.0 ) / 24.0
+                           eminus = eps * ( eps - 1.0 ) / 2.0 + alpha
+                           edot  = ( 1.0 - eps**2.0 ) - 2.0 * alpha
+                           eplus  = eps * ( eps + 1.0 ) / 2.0 + alpha
+      
+                           ptprs(k,jf,if) = eminus * ptprf(k,j-1,if) &
+                               + edot  * ptprf(k,j,if)   &
+                               + eplus  * ptprf(k,j+1,if)
+      
+                       END DO
+                   END DO
+      
+               END DO
+           END DO
+      
+           !
+           !-- Interpolation in z-direction (quadratic, Clark and Farley)
+      
+           DO k = nzbottom+1, nztop
+            
+               bottomz = (dzc/dzf) * (k-1) + 1
+               topz    = (dzc/dzf) * k
+      
+               DO jf = nys, nyn
+                   DO if = nxl, nxr
+                       DO kf = bottomz, topz
+                       
+                           eps    = ( zuf(kf) - zuc(k) ) / dzc
+                           alpha  = ( ( dzf / dzc )**2.0 - 1.0 ) / 24.0
+                           eminus = eps * ( eps - 1.0 ) / 2.0 + alpha
+                           edot  = ( 1.0 - eps**2.0 ) - 2.0 * alpha
+                           eplus  = eps * ( eps + 1.0 ) / 2.0 + alpha
+       
+                           interpol3d(kf,jf,if) = eminus * ptprs(k-1,jf,if) &
+                               + edot  * ptprs(k,jf,if)   &
+                               + eplus  * ptprs(k+1,jf,if)
+      
+                       END DO
+                   END DO
+               END DO
+      
+           END DO
+               
+           DO jf = nys, nyn
+               DO if = nxl, nxr
+      
+                   eps    = ( zuf(nzt+1) - zuc(nztop+1) ) / dzc
+                   alpha  = ( ( dzf / dzc )**2.0 - 1.0 ) / 24.0
+                   eminus = eps * ( eps - 1.0 ) / 2.0 + alpha
+                   edot  = ( 1.0 - eps**2.0 ) - 2.0 * alpha
+                   eplus  = eps * ( eps + 1.0 ) / 2.0 + alpha
+      
+                   interpol3d(nzt+1,jf,if) = eminus * ptprs(nztop,jf,if)   &
+                       + edot  * ptprs(nztop+1,jf,if) &
+                       + eplus  * ptprs(nztop+2,jf,if)
+      
+               END DO
+           END DO
+      
+      
+           DEALLOCATE( ptprf, ptprs ) 
+      
+#endif
+       END SUBROUTINE interpolate_to_fine_s
+      
+      
+       SUBROUTINE interpolate_to_fine_kh( tag )
+      
+#if defined( __parallel )
+      
+           USE arrays_3d
+           USE control_parameters
+           USE grid_variables
+           USE indices
+           USE pegrid
+           
+      
+           IMPLICIT NONE
+      
+          
+           INTEGER(iwp), intent(in) ::  tag
+           INTEGER(iwp)                             ::  i, j, k
+           INTEGER(iwp)                             ::  if, jf, kf
+           INTEGER(iwp)                             ::  bottomx, topx
+           INTEGER(iwp)                             ::  bottomy, topy
+           INTEGER(iwp)                             ::  bottomz, topz
+           REAL(wp)                                ::  eps, alpha, eminus, edot, eplus
+           REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  uprs, uprf
+           REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  vprs, vprf
+           REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  wprs, wprf
+           REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  ptprs, ptprf
+           REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  eprs, eprf
+           REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  kmprs, kmprf
+           REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  khprs, khprf
+           REAL(wp), DIMENSION(:,:), ALLOCATABLE   ::  shfpr, uswspr, vswspr
+           REAL(wp), DIMENSION(:,:), ALLOCATABLE   ::  tspr, uspr, z0pr
+      
+           INTEGER(iwp) :: nzbottom, nztop
+      
+      
+           nzbottom = bdims_rem (3,1)
+           nztop  = bdims_rem (3,2)
+           !   nztop  = blk_dim_rem (3,2)+1
+      
+      
+           ALLOCATE( ptprf(nzbottom:nztop+2,bdims_rem(2,1)-1:bdims_rem(2,2)+1,nxl:nxr) )
+           ALLOCATE( ptprs(nzbottom:nztop+2,nys:nyn,nxl:nxr) )
+      
+      
+           !
+           !-- Initialisation of scalar variables
+      
+           !
+           !-- Interpolation in x-direction (quadratic, Clark and Farley)
+      
+           DO k = nzbottom, nztop+2
+               DO j = bdims_rem(2,1)-1, bdims_rem(2,2)+1
+                   DO i = bdims_rem(1,1), bdims_rem(1,2)
+      
+                       bottomx = (nxf+1)/(nxc+1) * i
+                       topx    = (nxf+1)/(nxc+1) * (i+1) - 1
+      
+                       DO if = bottomx, topx
+      
+                           eps    = ( if * dxf + 0.5 * dxf - i * dxc - 0.5 * dxc ) / dxc
+                           alpha  = ( ( dxf / dxc )**2.0 - 1.0 ) / 24.0
+                           eminus = eps * ( eps - 1.0 ) / 2.0 + alpha
+                           edot  = ( 1.0 - eps**2.0 ) - 2.0 * alpha
+                           eplus  = eps * ( eps + 1.0 ) / 2.0 + alpha
+      
+                           ptprf(k,j,if) = eminus * work3d(k,j,i-1) &
+                               + edot  * work3d(k,j,i)   &
+                               + eplus  * work3d(k,j,i+1)
+                       END DO
+      
+                   END DO
+               END DO
+           END DO
+      
+           !
+           !-- Interpolation in y-direction (quadratic, Clark and Farley)
+      
+           DO k = nzbottom, nztop+2
+               DO j = bdims_rem(2,1), bdims_rem(2,2)
+        
+                   bottomy = (nyf+1)/(nyc+1) * j
+                   topy    = (nyf+1)/(nyc+1) * (j+1) - 1
+      
+                   DO if = nxl, nxr
+                       DO jf = bottomy, topy
+      
+                           eps    = ( jf * dyf + 0.5 * dyf - j * dyc - 0.5 * dyc ) / dyc
+                           alpha  = ( ( dyf / dyc )**2.0 - 1.0 ) / 24.0
+                           eminus = eps * ( eps - 1.0 ) / 2.0 + alpha
+                           edot  = ( 1.0 - eps**2.0 ) - 2.0 * alpha
+                           eplus  = eps * ( eps + 1.0 ) / 2.0 + alpha
+      
+                           ptprs(k,jf,if) = eminus * ptprf(k,j-1,if) &
+                               + edot  * ptprf(k,j,if)   &
+                               + eplus  * ptprf(k,j+1,if)
+      
+                       END DO
+                   END DO
+      
+               END DO
+           END DO
+      
+           !
+           !-- Interpolation in z-direction (quadratic, Clark and Farley)
+      
+           DO k = nzbottom+1, nztop
+            
+               bottomz = (dzc/dzf) * (k-1) + 1
+               topz    = (dzc/dzf) * k
+      
+               DO jf = nys, nyn
+                   DO if = nxl, nxr
+                       DO kf = bottomz, topz
+                       
+                           eps    = ( zuf(kf) - zuc(k) ) / dzc
+                           alpha  = ( ( dzf / dzc )**2.0 - 1.0 ) / 24.0
+                           eminus = eps * ( eps - 1.0 ) / 2.0 + alpha
+                           edot  = ( 1.0 - eps**2.0 ) - 2.0 * alpha
+                           eplus  = eps * ( eps + 1.0 ) / 2.0 + alpha
+       
+                           kh(kf,jf,if) = eminus * ptprs(k-1,jf,if) &
+                               + edot  * ptprs(k,jf,if)   &
+                               + eplus  * ptprs(k+1,jf,if)
+      
+                       END DO
+                   END DO
+               END DO
+      
+           END DO
+               
+           DO jf = nys, nyn
+               DO if = nxl, nxr
+      
+                   eps    = ( zuf(nzt+1) - zuc(nztop+1) ) / dzc
+                   alpha  = ( ( dzf / dzc )**2.0 - 1.0 ) / 24.0
+                   eminus = eps * ( eps - 1.0 ) / 2.0 + alpha
+                   edot  = ( 1.0 - eps**2.0 ) - 2.0 * alpha
+                   eplus  = eps * ( eps + 1.0 ) / 2.0 + alpha
+      
+                   kh(nzt+1,jf,if) = eminus * ptprs(nztop,jf,if)   &
+                       + edot  * ptprs(nztop+1,jf,if) &
+                       + eplus  * ptprs(nztop+2,jf,if)
+      
+               END DO
+           END DO
+      
+      
+           DEALLOCATE( ptprf, ptprs ) 
+      
+#endif
+       END SUBROUTINE interpolate_to_fine_kh
+      
+       SUBROUTINE interpolate_to_fine_km( tag )
+      
+#if defined( __parallel )
+      
+           USE arrays_3d
+           USE control_parameters
+           USE grid_variables
+           USE indices
+           USE pegrid
+           
+      
+           IMPLICIT NONE
+      
+          
+           INTEGER(iwp), intent(in) ::  tag
+           INTEGER(iwp)                             ::  i, j, k
+           INTEGER(iwp)                             ::  if, jf, kf
+           INTEGER(iwp)                             ::  bottomx, topx
+           INTEGER(iwp)                             ::  bottomy, topy
+           INTEGER(iwp)                             ::  bottomz, topz
+           REAL(wp)                                ::  eps, alpha, eminus, edot, eplus
+           REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  uprs, uprf
+           REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  vprs, vprf
+           REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  wprs, wprf
+           REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  ptprs, ptprf
+           REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  eprs, eprf
+           REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  kmprs, kmprf
+           REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  khprs, khprf
+           REAL(wp), DIMENSION(:,:), ALLOCATABLE   ::  shfpr, uswspr, vswspr
+           REAL(wp), DIMENSION(:,:), ALLOCATABLE   ::  tspr, uspr, z0pr
+      
+           INTEGER(iwp) :: nzbottom, nztop
+      
+      
+           nzbottom = bdims_rem (3,1)
+           nztop  = bdims_rem (3,2)
+           !   nztop  = blk_dim_rem (3,2)+1
+      
+      
+           ALLOCATE( ptprf(nzbottom:nztop+2,bdims_rem(2,1)-1:bdims_rem(2,2)+1,nxl:nxr) )
+           ALLOCATE( ptprs(nzbottom:nztop+2,nys:nyn,nxl:nxr) )
+      
+      
+           !
+           !-- Initialisation of scalar variables
+      
+           !
+           !-- Interpolation in x-direction (quadratic, Clark and Farley)
+      
+           DO k = nzbottom, nztop+2
+               DO j = bdims_rem(2,1)-1, bdims_rem(2,2)+1
+                   DO i = bdims_rem(1,1), bdims_rem(1,2)
+      
+                       bottomx = (nxf+1)/(nxc+1) * i
+                       topx    = (nxf+1)/(nxc+1) * (i+1) - 1
+      
+                       DO if = bottomx, topx
+      
+                           eps    = ( if * dxf + 0.5 * dxf - i * dxc - 0.5 * dxc ) / dxc
+                           alpha  = ( ( dxf / dxc )**2.0 - 1.0 ) / 24.0
+                           eminus = eps * ( eps - 1.0 ) / 2.0 + alpha
+                           edot  = ( 1.0 - eps**2.0 ) - 2.0 * alpha
+                           eplus  = eps * ( eps + 1.0 ) / 2.0 + alpha
+      
+                           ptprf(k,j,if) = eminus * work3d(k,j,i-1) &
+                               + edot  * work3d(k,j,i)   &
+                               + eplus  * work3d(k,j,i+1)
+                       END DO
+      
+                   END DO
+               END DO
+           END DO
+      
+           !
+           !-- Interpolation in y-direction (quadratic, Clark and Farley)
+      
+           DO k = nzbottom, nztop+2
+               DO j = bdims_rem(2,1), bdims_rem(2,2)
+        
+                   bottomy = (nyf+1)/(nyc+1) * j
+                   topy    = (nyf+1)/(nyc+1) * (j+1) - 1
+      
+                   DO if = nxl, nxr
+                       DO jf = bottomy, topy
+      
+                           eps    = ( jf * dyf + 0.5 * dyf - j * dyc - 0.5 * dyc ) / dyc
+                           alpha  = ( ( dyf / dyc )**2.0 - 1.0 ) / 24.0
+                           eminus = eps * ( eps - 1.0 ) / 2.0 + alpha
+                           edot  = ( 1.0 - eps**2.0 ) - 2.0 * alpha
+                           eplus  = eps * ( eps + 1.0 ) / 2.0 + alpha
+      
+                           ptprs(k,jf,if) = eminus * ptprf(k,j-1,if) &
+                               + edot  * ptprf(k,j,if)   &
+                               + eplus  * ptprf(k,j+1,if)
+      
+                       END DO
+                   END DO
+      
+               END DO
+           END DO
+      
+           !
+           !-- Interpolation in z-direction (quadratic, Clark and Farley)
+      
+           DO k = nzbottom+1, nztop
+            
+               bottomz = (dzc/dzf) * (k-1) + 1
+               topz    = (dzc/dzf) * k
+      
+               DO jf = nys, nyn
+                   DO if = nxl, nxr
+                       DO kf = bottomz, topz
+                       
+                           eps    = ( zuf(kf) - zuc(k) ) / dzc
+                           alpha  = ( ( dzf / dzc )**2.0 - 1.0 ) / 24.0
+                           eminus = eps * ( eps - 1.0 ) / 2.0 + alpha
+                           edot  = ( 1.0 - eps**2.0 ) - 2.0 * alpha
+                           eplus  = eps * ( eps + 1.0 ) / 2.0 + alpha
+       
+                           km(kf,jf,if) = eminus * ptprs(k-1,jf,if) &
+                               + edot  * ptprs(k,jf,if)   &
+                               + eplus  * ptprs(k+1,jf,if)
+      
+                       END DO
+                   END DO
+               END DO
+      
+           END DO
+               
+           DO jf = nys, nyn
+               DO if = nxl, nxr
+      
+                   eps    = ( zuf(nzt+1) - zuc(nztop+1) ) / dzc
+                   alpha  = ( ( dzf / dzc )**2.0 - 1.0 ) / 24.0
+                   eminus = eps * ( eps - 1.0 ) / 2.0 + alpha
+                   edot  = ( 1.0 - eps**2.0 ) - 2.0 * alpha
+                   eplus  = eps * ( eps + 1.0 ) / 2.0 + alpha
+      
+                   km(nzt+1,jf,if)  = eminus * ptprs(nztop,jf,if)   &
+                       + edot  * ptprs(nztop+1,jf,if) &
+                       + eplus  * ptprs(nztop+2,jf,if)
+      
+               END DO
+           END DO
+      
+      
+           DEALLOCATE( ptprf, ptprs ) 
+      
+#endif
+       END SUBROUTINE interpolate_to_fine_km
+      
+      
+      
+      
+       SUBROUTINE interpolate_to_fine_flux( tag )
+      
+#if defined( __parallel )
+      
+           USE arrays_3d
+           USE control_parameters
+           USE grid_variables
+           USE indices
+           USE pegrid
+           
+      
+           IMPLICIT NONE
+      
+          
+           INTEGER(iwp), intent(in) ::  tag
+           INTEGER(iwp)                             ::  i, j, k
+           INTEGER(iwp)                             ::  if, jf, kf
+           INTEGER(iwp)                             ::  bottomx, topx
+           INTEGER(iwp)                             ::  bottomy, topy
+           INTEGER(iwp)                             ::  bottomz, topz
+           REAL(wp)                                ::  eps, alpha, eminus, edot, eplus
+           REAL(wp), DIMENSION(:,:), ALLOCATABLE   ::  shfpr, uswspr, vswspr
+           REAL(wp), DIMENSION(:,:), ALLOCATABLE   ::  tspr, uspr, z0pr
+      
+      
+           INTEGER(iwp) :: nzbottom, nztop
+      
+           ALLOCATE( shfpr (bdims_rem(2,1)-1:bdims_rem(2,2)+1,nxl:nxr) )
+           ALLOCATE( uswspr(bdims_rem(2,1)-1:bdims_rem(2,2)+1,nxl:nxr) )
+           ALLOCATE( vswspr(bdims_rem(2,1)-1:bdims_rem(2,2)+1,nxl:nxr) )
+           ALLOCATE( tspr  (bdims_rem(2,1)-1:bdims_rem(2,2)+1,nxl:nxr) )
+           ALLOCATE( uspr  (bdims_rem(2,1)-1:bdims_rem(2,2)+1,nxl:nxr) )
+           ALLOCATE( z0pr  (bdims_rem(2,1)-1:bdims_rem(2,2)+1,nxl:nxr) )
+      
+           !
+           !-- Initialisation of scalar variables (2D)
+      
+           !
+           !-- Interpolation in x-direction (quadratic, Clark and Farley)
+      
+           DO j = bdims_rem(2,1)-1, bdims_rem(2,2)+1
+               DO i = bdims_rem(1,1), bdims_rem(1,2)
+              
+                   bottomx = (nxf+1)/(nxc+1) * i
+                   topx    = (nxf+1)/(nxc+1) * (i+1) - 1
+      
+                   DO if = bottomx, topx
+      
+                       eps    = ( if * dxf + 0.5 * dxf - i * dxc - 0.5 * dxc ) / dxc
+                       alpha  = ( ( dxf / dxc )**2.0 - 1.0 ) / 24.0
+                       eminus = eps * ( eps - 1.0 ) / 2.0 + alpha
+                       edot  = ( 1.0 - eps**2.0 ) - 2.0 * alpha
+                       eplus  = eps * ( eps + 1.0 ) / 2.0 + alpha
+      
+                       shfpr(j,if)  = eminus * work2dshf(j,i-1) &
+                           + edot  * work2dshf(j,i)   &
+                           + eplus  * work2dshf(j,i+1)
+      
+                       uswspr(j,if) = eminus * work2dusws(j,i-1) &
+                           + edot  * work2dusws(j,i)   &
+                           + eplus  * work2dusws(j,i+1)
+      
+                       vswspr(j,if) = eminus * work2dvsws(j,i-1) &
+                           + edot  * work2dvsws(j,i)   &
+                           + eplus  * work2dvsws(j,i+1)
+      
+                       tspr(j,if)   = eminus * work2dts(j,i-1) &
+                           + edot  * work2dts(j,i)   &
+                           + eplus  * work2dts(j,i+1)
+      
+                       uspr(j,if)   = eminus * work2dus(j,i-1) &
+                           + edot  * work2dus(j,i)   &
+                           + eplus  * work2dus(j,i+1)
+      
+                       z0pr(j,if)   = eminus * work2dz0(j,i-1) &
+                           + edot  * work2dz0(j,i)   &
+                           + eplus  * work2dz0(j,i+1)
+      
+                   END DO
+      
+               END DO
+           END DO
+      
+           !
+           !-- Interpolation in y-direction (quadratic, Clark and Farley)
+      
+           DO j = bdims_rem(2,1), bdims_rem(2,2)
+             
+               bottomy = (nyf+1)/(nyc+1) * j
+               topy    = (nyf+1)/(nyc+1) * (j+1) - 1
+      
+               DO if = nxl, nxr
+                   DO jf = bottomy, topy
+      
+                       eps    = ( jf * dyf + 0.5 * dyf - j * dyc - 0.5 * dyc ) / dyc
+                       alpha  = ( ( dyf / dyc )**2.0 - 1.0 ) / 24.0
+                       eminus = eps * ( eps - 1.0 ) / 2.0 + alpha
+                       edot  = ( 1.0 - eps**2.0 ) - 2.0 * alpha
+                       eplus  = eps * ( eps + 1.0 ) / 2.0 + alpha
+      
+!-- WARNING
+!-- shf,z0 not interpolated
+!-- line commented in interpolate_to_fine_flux
+!-- FG needs to read it's own data file      
+!MERGE-WIP!                    surf_def_h(0)%shf(jf,if)  = eminus * shfpr(j-1,if) &
+!MERGE-WIP!                        + edot  * shfpr(j,if)   &
+!MERGE-WIP!                        + eplus  * shfpr(j+1,if)
+!MERGE-WIP    
+!MERGE-WIP                     surf_def_h(0)%usws(jf,if) = eminus * uswspr(j-1,if) &
+!MERGE-WIP                         + edot  * uswspr(j,if)   &
+!MERGE-WIP                         + eplus  * uswspr(j+1,if)
+!MERGE-WIP    
+!MERGE-WIP                     surf_def_h(0)%vsws(jf,if) = eminus * vswspr(j-1,if) &
+!MERGE-WIP                         + edot  * vswspr(j,if)   &
+!MERGE-WIP                         + eplus  * vswspr(j+1,if)
+!MERGE-WIP    
+!MERGE                      ts(jf,if)   = eminus * tspr(j-1,if) &
+!MERGE                          + edot  * tspr(j,if)   &
+!MERGE                          + eplus  * tspr(j+1,if)
+!MERGE     
+!MERGE                      us(jf,if)   = eminus * uspr(j-1,if) &
+!MERGE                          + edot  * uspr(j,if)   &
+!MERGE                          + eplus  * uspr(j+1,if)
+!MERGE     
+!MERGE!                     z0(jf,if)   = eminus * z0pr(j-1,if) &
+!MERGE!                         + edot  * z0pr(j,if)   &
+!MERGE!                         + eplus  * z0pr(j+1,if)
+      
+                   END DO
+               END DO
+      
+           END DO
+      
+      
+           DEALLOCATE( shfpr, uswspr, vswspr )
+           DEALLOCATE( tspr, uspr, z0pr )
+      
+      
+#endif
+       END SUBROUTINE interpolate_to_fine_flux
+   
+   
+    END SUBROUTINE vnest_init_fine
+   
+    SUBROUTINE vnest_boundary_conds
+        !------------------------------------------------------------------------------!
+        ! Description:
+        ! ------------
+        ! Boundary conditions for the prognostic quantities.
+        ! One additional bottom boundary condition is applied for the TKE (=(u*)**2)
+        ! in prandtl_fluxes. The cyclic lateral boundary conditions are implicitly
+        ! handled in routine exchange_horiz. Pressure boundary conditions are
+        ! explicitly set in routines pres, poisfft, poismg and sor.
+        !------------------------------------------------------------------------------!
+    
+        USE arrays_3d
+        USE control_parameters
+        USE grid_variables
+        USE indices
+        USE pegrid
+        
+    
+        IMPLICIT NONE
+    
+        INTEGER(iwp) ::  i, j, k
+        INTEGER(iwp) ::  if, jf
+    
+        REAL(wp)    ::  c_max, denom
+    
+#if defined( __parallel )
+    
+        !
+        !-- vnest: top boundary conditions
+    
+        IF ( coupling_mode == 'vnested_crse' )  THEN
+        !-- Send data to fine grid for TOP BC
+    
+            offset(1) = ( pdims_partner(1) / pdims(1) ) * pcoord(1)
+            offset(2) = ( pdims_partner(2) / pdims(2) ) * pcoord(2)
+    
+            do j = 0,   ( pdims_partner(2) / pdims(2) ) - 1
+                do i = 0,   ( pdims_partner(1) / pdims(1) ) - 1
+                    map_coord(1) = i+offset(1)
+                    map_coord(2) = j+offset(2)
+    
+                    target_idex = f_rnk_lst(map_coord(1),map_coord(2)) + numprocs
+      
+                    bdims (1,1) =  c2f_dims_cg (0,map_coord(1),map_coord(2))
+                    bdims (1,2) =  c2f_dims_cg (1,map_coord(1),map_coord(2))
+                    bdims (2,1) =  c2f_dims_cg (2,map_coord(1),map_coord(2))
+                    bdims (2,2) =  c2f_dims_cg (3,map_coord(1),map_coord(2))
+                    bdims (3,1) =  c2f_dims_cg (4,map_coord(1),map_coord(2))
+                    bdims (3,2) =  c2f_dims_cg (5,map_coord(1),map_coord(2))
+  
+                    n_cell_c =      ( (bdims(1,2)-bdims(1,1)) + 3 ) * &
+                        ( (bdims(2,2)-bdims(2,1)) + 3 ) * &
+                        ( (bdims(3,2)-bdims(3,1)) + 1 )
+    
+                    CALL MPI_SEND(u (bdims(3,1), bdims(2,1)-1, bdims(1,1)-1), &
+                        1, TYPE_VNEST_BC, target_idex,    &
+                        201,    comm_inter, ierr)
+    
+                    CALL MPI_SEND(v(bdims(3,1), bdims(2,1)-1, bdims(1,1)-1),&
+                        1, TYPE_VNEST_BC, target_idex,    &
+                        202,    comm_inter, ierr)
+
+                   CALL MPI_SEND(w(bdims(3,1), bdims(2,1)-1, bdims(1,1)-1),&
+                        1, TYPE_VNEST_BC, target_idex,    &
+                        203,    comm_inter, ierr)
+    
+                    CALL MPI_SEND(pt(bdims(3,1), bdims(2,1)-1, bdims(1,1)-1),&
+                        1, TYPE_VNEST_BC, target_idex,    &
+                        205,    comm_inter, ierr)
+    
+                    IF ( humidity )  THEN 
+                    CALL MPI_SEND(q(bdims(3,1), bdims(2,1)-1, bdims(1,1)-1),&
+                        1, TYPE_VNEST_BC, target_idex,    &
+                        209,    comm_inter, ierr)
+                    ENDIF 
+ 
+                end do
+            end do
+    
+        ELSEIF ( coupling_mode == 'vnested_fine' )  THEN
+        !-- Receive data from coarse grid for TOP BC
+    
+            offset(1) = pcoord(1) / ( pdims(1)/pdims_partner(1) )
+            offset(2) = pcoord(2) / ( pdims(2)/pdims_partner(2) )
+            map_coord(1) = offset(1)
+            map_coord(2) = offset(2)
+            target_idex = c_rnk_lst(map_coord(1),map_coord(2))
+    
+            bdims_rem (1,1) = c2f_dims_fg(0)
+            bdims_rem (1,2) = c2f_dims_fg(1)
+            bdims_rem (2,1) = c2f_dims_fg(2)
+            bdims_rem (2,2) = c2f_dims_fg(3)
+            bdims_rem (3,1) = c2f_dims_fg(4)
+            bdims_rem (3,2) = c2f_dims_fg(5)
+ 
+            n_cell_c =                                    &
+                ( (bdims_rem(1,2)-bdims_rem(1,1)) + 3 ) * &
+                ( (bdims_rem(2,2)-bdims_rem(2,1)) + 3 ) * &
+                ( (bdims_rem(3,2)-bdims_rem(3,1)) + 1 )
+    
+            ALLOCATE( work3d  (                    &
+                bdims_rem(3,1)  :bdims_rem(3,2)  , &
+                bdims_rem(2,1)-1:bdims_rem(2,2)+1, &
+                bdims_rem(1,1)-1:bdims_rem(1,2)+1))
+    
+    
+            CALL MPI_RECV( work3d ,n_cell_c, MPI_REAL, target_idex, 201, &
+                comm_inter,status, ierr )
+            interpol3d => u
+            call vnest_set_topbc_u
+    
+            CALL MPI_RECV( work3d ,n_cell_c, MPI_REAL, target_idex, 202, &
+                comm_inter,status, ierr )
+            interpol3d => v
+            call vnest_set_topbc_v
+    
+            CALL MPI_RECV( work3d ,n_cell_c, MPI_REAL, target_idex, 203, &
+                comm_inter,status, ierr )
+            interpol3d => w
+            call vnest_set_topbc_w
+    
+            CALL MPI_RECV( work3d,n_cell_c, MPI_REAL, target_idex, 205, &
+                comm_inter,status, ierr )
+            interpol3d => pt
+            call vnest_set_topbc_s
+    
+            IF ( humidity )  THEN
+    		CALL MPI_RECV( work3d,n_cell_c, MPI_REAL, target_idex, 209, &
+                        comm_inter,status, ierr )
+                interpol3d => q
+                call vnest_set_topbc_s
+    
+                CALL exchange_horiz_2d(q (nzt+1,:,:) )
+            ENDIF
+    
+!-- TKE Neumann BC for FG top
+            DO jf = nys, nyn
+                DO if = nxl, nxr
+                   e(nzt+1,jf,if) = e(nzt,jf,if)
+                END DO
+            END DO
+
+!
+!-- w level nzt+1 does not impact results. Only to avoid jumps while
+!-- plotting profiles
+            w(nzt+1,:,:) = w(nzt,:,:)
+
+            CALL exchange_horiz_2d(u (nzt+1,:,:) )
+            CALL exchange_horiz_2d(v (nzt+1,:,:) )
+            CALL exchange_horiz_2d(pt(nzt+1,:,:) )
+            CALL exchange_horiz_2d(e (nzt+1,:,:) )
+            CALL exchange_horiz_2d(w (nzt+1,:,:) )
+            CALL exchange_horiz_2d(w (nzt  ,:,:) )
+    
+            NULLIFY    (   interpol3d  )
+            DEALLOCATE (   work3d      )
+    
+        ENDIF
+    
+    
+#endif
+    CONTAINS
+    
+       SUBROUTINE vnest_set_topbc_w
+       
+#if defined( __parallel )
+       
+           USE arrays_3d
+           USE control_parameters
+           USE grid_variables
+           USE indices
+           USE pegrid
+           
+       
+           IMPLICIT NONE
+       
+           INTEGER(iwp)                             ::  i, j, k
+           INTEGER(iwp)                             ::  if, jf
+           INTEGER(iwp)                             ::  bottomx, topx
+           INTEGER(iwp)                             ::  bottomy, topy
+           REAL(wp)                                ::  eps, alpha, eminus, edot, eplus
+           REAL(wp), DIMENSION(:,:), ALLOCATABLE   ::  wprf
+         
+           
+           ALLOCATE( wprf(bdims_rem(2,1)-1:bdims_rem(2,2)+1,nxl:nxr) )
+       
+           !
+           !-- Determination of a boundary condition for the vertical velocity component w:
+           !-- In this case only interpolation in x- and y- direction is necessary, as the
+           !-- boundary w-node of the fine grid coincides with a w-node in the coarse grid.
+           !-- For both interpolations the scheme of Clark and Farley is used.
+       
+           !
+           !-- Interpolation in x-direction
+       
+           DO j = bdims_rem(2,1)-1, bdims_rem(2,2)+1
+       
+               DO i = bdims_rem(1,1), bdims_rem(1,2)
+       
+                   bottomx = (nxf+1)/(nxc+1) * i
+                   topx    = (nxf+1)/(nxc+1) * (i+1) - 1
+       
+                   DO if = bottomx, topx
+       
+                       eps = (if * dxf + 0.5 * dxf - i * dxc - 0.5 * dxc) / dxc
+                       alpha = ( (dxf/dxc)**2.0 - 1.0) / 24.0
+                       eminus = eps * ( eps - 1.0 ) / 2.0 + alpha
+                       edot  = ( 1.0 - eps**2.0 ) - 2.0 * alpha
+                       eplus  = eps * ( eps + 1.0 ) / 2.0 + alpha
+                       wprf(j,if) = eminus * work3d(bdims_rem(3,1),j,i-1) &
+                           + edot  * work3d(bdims_rem(3,1),j,i)   &
+                           + eplus  * work3d(bdims_rem(3,1),j,i+1)
+       
+                   END DO
+       
+               END DO
+       
+           END DO
+              
+           !
+           !-- Interpolation in y-direction
+       
+           DO j = bdims_rem(2,1), bdims_rem(2,2)
+       
+               bottomy = (nyf+1)/(nyc+1) * j
+               topy    = (nyf+1)/(nyc+1) * (j+1) - 1
+       
+               DO if = nxl, nxr
+       
+                   DO jf = bottomy, topy
+       
+                       eps = (jf * dyf + 0.5 * dyf - j * dyc - 0.5 * dyc) / dyc
+       
+                       alpha = ( (dyf/dyc)**2.0 - 1.0) / 24.0
+       
+                       eminus = eps * ( eps - 1.0 ) / 2.0 + alpha
+       
+                       edot = ( 1.0 - eps**2.0 ) - 2.0 * alpha
+       
+                       eplus = eps * ( eps + 1.0 ) / 2.0 + alpha
+       
+                       w(nzt,jf,if) = eminus * wprf(j-1,if) &
+                           + edot  * wprf(j,if)   &
+                           + eplus  * wprf(j+1,if)
+       
+                   END DO
+       
+               END DO
+       
+           END DO
+
+           DEALLOCATE( wprf )
+#endif 
+       
+       END SUBROUTINE vnest_set_topbc_w
+       
+       
+       SUBROUTINE vnest_set_topbc_u
+       
+#if defined( __parallel )
+       
+           USE arrays_3d
+           USE control_parameters
+           USE grid_variables
+           USE indices
+           USE pegrid
+           
+       
+           IMPLICIT NONE
+       
+           INTEGER(iwp)                             ::  i, j, k
+           INTEGER(iwp)                             ::  if, jf
+           INTEGER(iwp)                             ::  bottomx, topx
+           INTEGER(iwp)                             ::  bottomy, topy
+           REAL(wp)                                ::  eps, alpha, eminus, edot, eplus
+           REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  uprf
+           REAL(wp), DIMENSION(:,:), ALLOCATABLE   ::  uprs
+         
+           
+       
+           ALLOCATE( uprf(bdims_rem(3,1):bdims_rem(3,2),nys:nyn,bdims_rem(1,1)-1:bdims_rem(1,2)+1) )
+           ALLOCATE( uprs(nys:nyn,bdims_rem(1,1)-1:bdims_rem(1,2)+1) )
+       
+       
+           !
+           !-- Interpolation in y-direction
+       
+           DO k = bdims_rem(3,1), bdims_rem(3,2)
+               DO j = bdims_rem(2,1), bdims_rem(2,2)
+       
+                   bottomy = (nyf+1)/(nyc+1) * j
+                   topy    = (nyf+1)/(nyc+1) * (j+1) - 1
+       
+                   DO i = bdims_rem(1,1)-1, bdims_rem(1,2)+1
+                       DO jf = bottomy, topy
+       
+                           eps = (jf * dyf + 0.5 * dyf - j * dyc - 0.5 * dyc) / dyc
+                           alpha = ( (dyf/dyc)**2.0 - 1.0) / 24.0
+                           eminus = eps * ( eps - 1.0 ) / 2.0 + alpha
+                           edot  = ( 1.0 - eps**2.0 ) - 2.0 * alpha
+                           eplus  = eps * ( eps + 1.0 ) / 2.0 + alpha
+       
+                           uprf(k,jf,i) = eminus * work3d(k,j-1,i) &
+                               + edot  * work3d(k,j,i)   &
+                               + eplus  * work3d(k,j+1,i)
+                       END DO
+                   END DO
+       
+               END DO
+           END DO
+       
+           !
+           !-- Interpolation in z-direction
+       
+           DO jf = nys, nyn
+               DO i = bdims_rem(1,1)-1, bdims_rem(1,2)+1
+                   eps = ( zuf(nzt+1) - zuc(bdims_rem(3,1)+1) ) / dzc
+                   alpha = ( (dzf/dzc)**2.0 - 1.0) / 24.0
+                   eminus = eps * ( eps - 1.0 ) / 2.0 + alpha
+                   edot  = ( 1.0 - eps**2.0 ) - 2.0 * alpha
+                   eplus  = eps * ( eps + 1.0 ) / 2.0 + alpha
+                   uprs(jf,i) = eminus * uprf(bdims_rem(3,1),jf,i)   &
+                       + edot  * uprf(bdims_rem(3,1)+1,jf,i) &
+                       + eplus  * uprf(bdims_rem(3,1)+2,jf,i)
+               END DO
+           END DO
+       
+           !
+           !-- Interpolation in x-direction
+       
+           DO jf = nys, nyn
+               DO i = bdims_rem(1,1), bdims_rem(1,2)
+       
+                   bottomx = (nxf+1)/(nxc+1) * i
+                   topx    = (nxf+1)/(nxc+1) * (i+1) - 1
+       
+                   DO if = bottomx, topx
+                       u(nzt+1,jf,if) = uprs(jf,i) + ( if * dxf - i * dxc ) * ( uprs(jf,i+1) - uprs(jf,i) ) / dxc
+                   END DO
+       
+               END DO
+           END DO           
+       
+       
+       
+           DEALLOCATE ( uprf, uprs )
+#endif 
+       
+       END SUBROUTINE vnest_set_topbc_u
+       
+       
+       SUBROUTINE vnest_set_topbc_v
+       
+#if defined( __parallel )
+       
+           USE arrays_3d
+           USE control_parameters
+           USE grid_variables
+           USE indices
+           USE pegrid
+           
+       
+           IMPLICIT NONE
+       
+           INTEGER(iwp)                             ::  i, j, k
+           INTEGER(iwp)                             ::  if, jf
+           INTEGER(iwp)                             ::  bottomx, topx
+           INTEGER(iwp)                             ::  bottomy, topy
+           REAL(wp)                                ::  eps, alpha, eminus, edot, eplus
+           REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  vprf
+           REAL(wp), DIMENSION(:,:), ALLOCATABLE   ::  vprs
+         
+           
+       
+           ALLOCATE( vprf(bdims_rem(3,1):bdims_rem(3,2),bdims_rem(2,1)-1:bdims_rem(2,2)+1,nxl:nxr) )
+           ALLOCATE( vprs(bdims_rem(2,1)-1:bdims_rem(2,2)+1,nxl:nxr) )
+           !
+           !-- Determination of a boundary condition for the horizontal velocity component v:
+           !-- Interpolation in x- and z-direction is carried out by using the scheme,
+           !-- which was derived by Clark and Farley (1984). In y-direction a
+           !-- linear interpolation is carried out.
+       
+           !
+           !-- Interpolation in x-direction
+       
+           DO k = bdims_rem(3,1), bdims_rem(3,2)
+               DO j = bdims_rem(2,1)-1, bdims_rem(2,2)+1
+                   DO i = bdims_rem(1,1), bdims_rem(1,2)
+       
+                       bottomx = (nxf+1)/(nxc+1) * i
+                       topx    = (nxf+1)/(nxc+1) * (i+1) - 1
+       
+                       DO if = bottomx, topx
+       
+                           eps = (if * dxf + 0.5 * dxf - i * dxc - 0.5 * dxc) / dxc
+                           alpha = ( (dxf/dxc)**2.0 - 1.0) / 24.0
+                           eminus = eps * ( eps - 1.0 ) / 2.0 + alpha
+                           edot  = ( 1.0 - eps**2.0 ) - 2.0 * alpha
+                           eplus  = eps * ( eps + 1.0 ) / 2.0 + alpha
+                           vprf(k,j,if) = eminus * work3d(k,j,i-1) &
+                               + edot  * work3d(k,j,i)   &
+                               + eplus  * work3d(k,j,i+1)
+                       END DO
+       
+                   END DO
+               END DO
+           END DO
+       
+           !
+           !-- Interpolation in z-direction
+       
+           DO j = bdims_rem(2,1)-1, bdims_rem(2,2)+1
+               DO if = nxl, nxr
+       
+                   eps = ( zuf(nzt+1) - zuc(bdims_rem(3,1)+1) ) / dzc
+                   alpha = ( (dzf/dzc)**2.0 - 1.0) / 24.0
+                   eminus = eps * ( eps - 1.0 ) / 2.0 + alpha
+                   edot  = ( 1.0 - eps**2.0 ) - 2.0 * alpha
+                   eplus  = eps * ( eps + 1.0 ) / 2.0 + alpha
+                   vprs(j,if) = eminus * vprf(bdims_rem(3,1),j,if)   &
+                       + edot  * vprf(bdims_rem(3,1)+1,j,if) &
+                       + eplus  * vprf(bdims_rem(3,1)+2,j,if)
+       
+               END DO
+           END DO
+       
+           !
+           !-- Interpolation in y-direction
+       
+           DO j = bdims_rem(2,1), bdims_rem(2,2)
+               DO if = nxl, nxr
+       
+                   bottomy = (nyf+1)/(nyc+1) * j
+                   topy    = (nyf+1)/(nyc+1) * (j+1) - 1
+       
+                   DO jf = bottomy, topy
+       
+                       v(nzt+1,jf,if) = vprs(j,if) + ( jf * dyf - j * dyc ) * ( vprs(j+1,if) - vprs(j,if) ) / dyc
+       
+                   END DO
+               END DO
+           END DO      
+       
+       
+           DEALLOCATE ( vprf, vprs)
+       
+       
+#endif 
+       
+       END SUBROUTINE vnest_set_topbc_v
+       
+       
+       SUBROUTINE vnest_set_topbc_s
+       
+#if defined( __parallel )
+       
+           USE arrays_3d
+           USE control_parameters
+           USE grid_variables
+           USE indices
+           USE pegrid
+           
+       
+           IMPLICIT NONE
+       
+           INTEGER(iwp)                             ::  i, j, k
+           INTEGER(iwp)                             ::  if, jf
+           INTEGER(iwp)                             ::  bottomx, topx
+           INTEGER(iwp)                             ::  bottomy, topy
+           REAL(wp)                                ::  eps, alpha, eminus, edot, eplus
+           REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  ptprf, ptprs
+         
+           
+       
+           ALLOCATE( ptprf(bdims_rem(3,1):bdims_rem(3,2),bdims_rem(2,1)-1:bdims_rem(2,2)+1,nxl:nxr) )
+           ALLOCATE( ptprs(bdims_rem(3,1):bdims_rem(3,2),nys:nyn,nxl:nxr) )
+       
+           !
+           !-- Determination of a boundary condition for the potential temperature pt:
+           !-- The scheme derived by Clark and Farley can be used in all three dimensions.
+       
+           !
+           !-- Interpolation in x-direction
+       
+           DO k = bdims_rem(3,1), bdims_rem(3,2)
+       
+               DO j = bdims_rem(2,1)-1, bdims_rem(2,2)+1
+       
+                   DO i = bdims_rem(1,1), bdims_rem(1,2)
+       
+                       bottomx = (nxf+1)/(nxc+1) * i
+                       topx    = (nxf+1)/(nxc+1) *(i+1) - 1
+       
+                       DO if = bottomx, topx
+       
+                           eps = (if * dxf + 0.5 * dxf - i * dxc - 0.5 * dxc) / dxc
+       
+                           alpha = ( (dxf/dxc)**2.0 - 1.0) / 24.0
+       
+                           eminus = eps * (eps - 1.0 ) / 2.0 + alpha
+       
+                           edot = ( 1.0 - eps**2.0 ) - 2.0 * alpha
+       
+                           eplus = eps * ( eps + 1.0 ) / 2.0 + alpha
+       
+                           ptprf(k,j,if) = eminus * work3d(k,j,i-1) &
+                               + edot  * work3d(k,j,i)   &
+                               + eplus  * work3d(k,j,i+1)
+                       END DO
+       
+                   END DO
+       
+               END DO
+       
+           END DO
+       
+           !
+           !-- Interpolation in y-direction
+       
+           DO k = bdims_rem(3,1), bdims_rem(3,2)
+       
+               DO j = bdims_rem(2,1), bdims_rem(2,2)
+       
+                   bottomy = (nyf+1)/(nyc+1) * j
+                   topy    = (nyf+1)/(nyc+1) * (j+1) - 1
+       
+                   DO if = nxl, nxr
+       
+                       DO jf = bottomy, topy
+       
+                           eps = (jf * dyf + 0.5 * dyf - j * dyc - 0.5 * dyc) / dyc
+       
+                           alpha = ( (dyf/dyc)**2.0 - 1.0) / 24.0
+                   
+                           eminus = eps * (eps - 1.0) / 2.0 + alpha
+       
+                           edot  = ( 1.0 - eps**2.0 ) - 2.0 * alpha
+       
+                           eplus  = eps * ( eps + 1.0 ) / 2.0 + alpha
+       
+                           ptprs(k,jf,if) = eminus * ptprf(k,j-1,if) &
+                               + edot  * ptprf(k,j,if)   &
+                               + eplus  * ptprf(k,j+1,if)
+                       END DO
+       
+                   END DO
+       
+               END DO
+       
+           END DO
+       
+           !
+           !-- Interpolation in z-direction
+       
+           DO jf = nys, nyn
+               DO if = nxl, nxr
+       
+                   eps = ( zuf(nzt+1) - zuc(bdims_rem(3,1)+1) ) / dzc
+       
+                   alpha = ( (dzf/dzc)**2.0 - 1.0) / 24.0
+       
+                   eminus = eps * ( eps - 1.0 ) / 2.0 + alpha
+       
+                   edot  = ( 1.0 - eps**2.0 ) - 2.0 * alpha
+       
+                   eplus  = eps * ( eps + 1.0 ) / 2.0 + alpha
+       
+                   interpol3d (nzt+1,jf,if) = eminus * ptprs(bdims_rem(3,1),jf,if)   &
+                       + edot  * ptprs(bdims_rem(3,1)+1,jf,if) &
+                       + eplus  * ptprs(bdims_rem(3,1)+2,jf,if)
+       
+               END DO
+           END DO          
+           
+           DEALLOCATE ( ptprf, ptprs )
+       
+       
+#endif 
+       
+       END SUBROUTINE vnest_set_topbc_s
+    END SUBROUTINE vnest_boundary_conds
+    
+ 
+    SUBROUTINE vnest_boundary_conds_khkm
+    
+        !--------------------------------------------------------------------------------!
+        ! Description:
+        ! ------------
+        ! Boundary conditions for the prognostic quantities.
+        ! One additional bottom boundary condition is applied for the TKE (=(u*)**2)
+        ! in prandtl_fluxes. The cyclic lateral boundary conditions are implicitly
+        ! handled in routine exchange_horiz. Pressure boundary conditions are
+        ! explicitly set in routines pres, poisfft, poismg and sor.
+        !------------------------------------------------------------------------------!
+    
+        USE arrays_3d
+        USE control_parameters
+        USE grid_variables
+        USE indices
+        USE pegrid
+        
+    
+        IMPLICIT NONE
+    
+        INTEGER(iwp) ::  i, j, k
+        INTEGER(iwp) ::  if, jf
+    
+        REAL(wp)    ::  c_max, denom
+    
+#if defined( __parallel )
+    
+        IF ( coupling_mode == 'vnested_crse' )  THEN
+        ! Send data to fine grid for TOP BC
+    
+            offset(1) = ( pdims_partner(1) / pdims(1) ) * pcoord(1)
+            offset(2) = ( pdims_partner(2) / pdims(2) ) * pcoord(2)
+    
+            do j = 0,   ( pdims_partner(2) / pdims(2) ) - 1
+                do i = 0,   ( pdims_partner(1) / pdims(1) ) - 1
+                    map_coord(1) = i+offset(1)
+                    map_coord(2) = j+offset(2)
+    
+                    target_idex = f_rnk_lst(map_coord(1),map_coord(2)) + numprocs
+    
+                    CALL MPI_RECV( bdims_rem,   6, MPI_INTEGER, target_idex, 10, &
+                        comm_inter,status, ierr )
+    
+                    bdims (1,1) = bdims_rem (1,1) / cfratio(1)
+                    bdims (1,2) = bdims_rem (1,2) / cfratio(1)
+                    bdims (2,1) = bdims_rem (2,1) / cfratio(2)
+                    bdims (2,2) = bdims_rem (2,2) / cfratio(2)
+                    bdims (3,1) = bdims_rem (3,2) / cfratio(3)
+                    bdims (3,2) = bdims     (3,1) + 2
+    
+                    CALL MPI_SEND( bdims,       6, MPI_INTEGER, target_idex, 9, &
+                        comm_inter, ierr )
+    
+    
+                    n_cell_c =      ( (bdims(1,2)-bdims(1,1)) + 3 ) * &
+                        ( (bdims(2,2)-bdims(2,1)) + 3 ) * &
+                        ( (bdims(3,2)-bdims(3,1)) + 1 )
+    
+                    CALL MPI_SEND(kh(bdims(3,1)  :bdims(3,2)  , &
+                        bdims(2,1)-1:bdims(2,2)+1, &
+                        bdims(1,1)-1:bdims(1,2)+1),&
+                        n_cell_c, MPI_REAL, target_idex,    &
+                        207, 	comm_inter, ierr)
+    
+                    CALL MPI_SEND(km(bdims(3,1)  :bdims(3,2)  , &
+                        bdims(2,1)-1:bdims(2,2)+1, &
+                        bdims(1,1)-1:bdims(1,2)+1),&
+                        n_cell_c, MPI_REAL, target_idex,    &
+                        208, 	comm_inter, ierr)
+    
+    
+    
+                end do
+            end do
+    
+        ELSEIF ( coupling_mode == 'vnested_fine' )  THEN
+                      ! Receive data from coarse grid for TOP BC
+    
+            offset(1) = pcoord(1) / ( pdims(1)/pdims_partner(1) )
+            offset(2) = pcoord(2) / ( pdims(2)/pdims_partner(2) )
+            map_coord(1) = offset(1)
+            map_coord(2) = offset(2)
+            target_idex = c_rnk_lst(map_coord(1),map_coord(2))
+    
+            bdims (1,1) = nxl
+            bdims (1,2) = nxr
+            bdims (2,1) = nys
+            bdims (2,2) = nyn
+            bdims (3,1) = nzb
+            bdims (3,2) = nzt
+    
+            CALL MPI_SEND( bdims,       6, MPI_INTEGER, target_idex, 10, &
+                comm_inter, ierr )
+    
+            CALL MPI_RECV( bdims_rem,   6, MPI_INTEGER, target_idex, 9, &
+                comm_inter,status, ierr )
+    
+            n_cell_c =       ( (bdims_rem(1,2)-bdims_rem(1,1)) + 3 ) * &
+                ( (bdims_rem(2,2)-bdims_rem(2,1)) + 3 ) * &
+                ( (bdims_rem(3,2)-bdims_rem(3,1)) + 1 )
+    
+            ALLOCATE( work3d  ( bdims_rem(3,1)  :bdims_rem(3,2)  , &
+                bdims_rem(2,1)-1:bdims_rem(2,2)+1, &
+                bdims_rem(1,1)-1:bdims_rem(1,2)+1))
+    
+    
+            CALL MPI_RECV( work3d,n_cell_c, MPI_REAL, target_idex, 207, &
+                comm_inter,status, ierr )
+    
+        ! Neumann BC for FG kh
+        DO jf = nys, nyn
+            DO if = nxl, nxr
+               kh(nzt+1,jf,if) = kh(nzt,jf,if)
+            END DO
+        END DO
+    
+        CALL MPI_RECV( work3d,n_cell_c, MPI_REAL, target_idex, 208, &
+             comm_inter,status, ierr )
+    
+        ! Neumann BC for FG kh
+        DO jf = nys, nyn
+            DO if = nxl, nxr
+               km(nzt+1,jf,if) = km(nzt,jf,if)
+            END DO
+        END DO
+    
+    
+        !
+        !-- The following evaluation can only be performed, if the fine grid is situated below the inversion
+        !!    DO jf = nys-1, nyn+1
+        !!       DO if = nxl-1, nxr+1
+        !!
+        !!          km(nzt+1,jf,if) = 0.1 * l_grid(nzt+1) * SQRT( e(nzt+1,jf,if) )
+        !!          kh(nzt+1,jf,if) = 3.0 * km(nzt+1,jf,if)
+        !!
+        !!       END DO
+        !!    END DO
+    
+        CALL exchange_horiz_2d(km(nzt+1,:,:) )
+        CALL exchange_horiz_2d(kh(nzt+1,:,:) )
+    
+        DEALLOCATE (   work3d      )
+    
+        ENDIF
+    
+#endif 
+    
+    CONTAINS
+ 
+       SUBROUTINE vnest_set_topbc_kh
+       
+#if defined( __parallel )
+       
+           USE arrays_3d
+           USE control_parameters
+           USE grid_variables
+           USE indices
+           USE pegrid
+           
+       
+           IMPLICIT NONE
+       
+           INTEGER(iwp)                             ::  i, j, k
+           INTEGER(iwp)                             ::  if, jf
+           INTEGER(iwp)                             ::  bottomx, topx
+           INTEGER(iwp)                             ::  bottomy, topy
+           REAL(wp)                                ::  eps, alpha, eminus, edot, eplus
+           REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  ptprf, ptprs
+         
+           
+       
+           ALLOCATE( ptprf(bdims_rem(3,1):bdims_rem(3,2),bdims_rem(2,1)-1:bdims_rem(2,2)+1,nxl:nxr) )
+           ALLOCATE( ptprs(bdims_rem(3,1):bdims_rem(3,2),nys:nyn,nxl:nxr) )
+       
+           !
+           !-- Determination of a boundary condition for the potential temperature pt:
+           !-- The scheme derived by Clark and Farley can be used in all three dimensions.
+       
+           !
+           !-- Interpolation in x-direction
+       
+           DO k = bdims_rem(3,1), bdims_rem(3,2)
+       
+               DO j = bdims_rem(2,1)-1, bdims_rem(2,2)+1
+       
+                   DO i = bdims_rem(1,1), bdims_rem(1,2)
+       
+                       bottomx = (nxf+1)/(nxc+1) * i
+                       topx    = (nxf+1)/(nxc+1) *(i+1) - 1
+       
+                       DO if = bottomx, topx
+       
+                           eps = (if * dxf + 0.5 * dxf - i * dxc - 0.5 * dxc) / dxc
+       
+                           alpha = ( (dxf/dxc)**2.0 - 1.0) / 24.0
+       
+                           eminus = eps * (eps - 1.0 ) / 2.0 + alpha
+       
+                           edot = ( 1.0 - eps**2.0 ) - 2.0 * alpha
+       
+                           eplus = eps * ( eps + 1.0 ) / 2.0 + alpha
+       
+                           ptprf(k,j,if) = eminus * work3d(k,j,i-1) &
+                               + edot  * work3d(k,j,i)   &
+                               + eplus  * work3d(k,j,i+1)
+                       END DO
+       
+                   END DO
+       
+               END DO
+       
+           END DO
+       
+           !
+           !-- Interpolation in y-direction
+       
+           DO k = bdims_rem(3,1), bdims_rem(3,2)
+       
+               DO j = bdims_rem(2,1), bdims_rem(2,2)
+       
+                   bottomy = (nyf+1)/(nyc+1) * j
+                   topy    = (nyf+1)/(nyc+1) * (j+1) - 1
+       
+                   DO if = nxl, nxr
+       
+                       DO jf = bottomy, topy
+       
+                           eps = (jf * dyf + 0.5 * dyf - j * dyc - 0.5 * dyc) / dyc
+       
+                           alpha = ( (dyf/dyc)**2.0 - 1.0) / 24.0
+                   
+                           eminus = eps * (eps - 1.0) / 2.0 + alpha
+       
+                           edot  = ( 1.0 - eps**2.0 ) - 2.0 * alpha
+       
+                           eplus  = eps * ( eps + 1.0 ) / 2.0 + alpha
+       
+                           ptprs(k,jf,if) = eminus * ptprf(k,j-1,if) &
+                               + edot  * ptprf(k,j,if)   &
+                               + eplus  * ptprf(k,j+1,if)
+                       END DO
+       
+                   END DO
+       
+               END DO
+       
+           END DO
+       
+           !
+           !-- Interpolation in z-direction
+       
+           DO jf = nys, nyn
+               DO if = nxl, nxr
+       
+                   eps = ( zuf(nzt+1) - zuc(bdims_rem(3,1)+1) ) / dzc
+       
+                   alpha = ( (dzf/dzc)**2.0 - 1.0) / 24.0
+       
+                   eminus = eps * ( eps - 1.0 ) / 2.0 + alpha
+       
+                   edot  = ( 1.0 - eps**2.0 ) - 2.0 * alpha
+       
+                   eplus  = eps * ( eps + 1.0 ) / 2.0 + alpha
+       
+                   kh (nzt+1,jf,if) = eminus * ptprs(bdims_rem(3,1),jf,if)   &
+                       + edot  * ptprs(bdims_rem(3,1)+1,jf,if) &
+                       + eplus  * ptprs(bdims_rem(3,1)+2,jf,if)
+       
+               END DO
+           END DO          
+           
+           DEALLOCATE ( ptprf, ptprs )
+       
+       
+#endif 
+       
+       END SUBROUTINE vnest_set_topbc_kh
+       
+       SUBROUTINE vnest_set_topbc_km
+       
+#if defined( __parallel )
+       
+           USE arrays_3d
+           USE control_parameters
+           USE grid_variables
+           USE indices
+           USE pegrid
+           
+       
+           IMPLICIT NONE
+       
+           INTEGER(iwp)                             ::  i, j, k
+           INTEGER(iwp)                             ::  if, jf
+           INTEGER(iwp)                             ::  bottomx, topx
+           INTEGER(iwp)                             ::  bottomy, topy
+           REAL(wp)                                ::  eps, alpha, eminus, edot, eplus
+           REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  ptprf, ptprs
+         
+           
+       
+           ALLOCATE( ptprf(bdims_rem(3,1):bdims_rem(3,2),bdims_rem(2,1)-1:bdims_rem(2,2)+1,nxl:nxr) )
+           ALLOCATE( ptprs(bdims_rem(3,1):bdims_rem(3,2),nys:nyn,nxl:nxr) )
+       
+           !
+           !-- Determination of a boundary condition for the potential temperature pt:
+           !-- The scheme derived by Clark and Farley can be used in all three dimensions.
+       
+           !
+           !-- Interpolation in x-direction
+       
+           DO k = bdims_rem(3,1), bdims_rem(3,2)
+       
+               DO j = bdims_rem(2,1)-1, bdims_rem(2,2)+1
+       
+                   DO i = bdims_rem(1,1), bdims_rem(1,2)
+       
+                       bottomx = (nxf+1)/(nxc+1) * i
+                       topx    = (nxf+1)/(nxc+1) *(i+1) - 1
+       
+                       DO if = bottomx, topx
+       
+                           eps = (if * dxf + 0.5 * dxf - i * dxc - 0.5 * dxc) / dxc
+       
+                           alpha = ( (dxf/dxc)**2.0 - 1.0) / 24.0
+       
+                           eminus = eps * (eps - 1.0 ) / 2.0 + alpha
+       
+                           edot = ( 1.0 - eps**2.0 ) - 2.0 * alpha
+       
+                           eplus = eps * ( eps + 1.0 ) / 2.0 + alpha
+       
+                           ptprf(k,j,if) = eminus * work3d(k,j,i-1) &
+                               + edot  * work3d(k,j,i)   &
+                               + eplus  * work3d(k,j,i+1)
+                       END DO
+       
+                   END DO
+       
+               END DO
+       
+           END DO
+       
+           !
+           !-- Interpolation in y-direction
+       
+           DO k = bdims_rem(3,1), bdims_rem(3,2)
+       
+               DO j = bdims_rem(2,1), bdims_rem(2,2)
+       
+                   bottomy = (nyf+1)/(nyc+1) * j
+                   topy    = (nyf+1)/(nyc+1) * (j+1) - 1
+       
+                   DO if = nxl, nxr
+       
+                       DO jf = bottomy, topy
+       
+                           eps = (jf * dyf + 0.5 * dyf - j * dyc - 0.5 * dyc) / dyc
+       
+                           alpha = ( (dyf/dyc)**2.0 - 1.0) / 24.0
+                   
+                           eminus = eps * (eps - 1.0) / 2.0 + alpha
+       
+                           edot  = ( 1.0 - eps**2.0 ) - 2.0 * alpha
+       
+                           eplus  = eps * ( eps + 1.0 ) / 2.0 + alpha
+       
+                           ptprs(k,jf,if) = eminus * ptprf(k,j-1,if) &
+                               + edot  * ptprf(k,j,if)   &
+                               + eplus  * ptprf(k,j+1,if)
+                       END DO
+       
+                   END DO
+       
+               END DO
+       
+           END DO
+       
+           !
+           !-- Interpolation in z-direction
+       
+           DO jf = nys, nyn
+               DO if = nxl, nxr
+       
+                   eps = ( zuf(nzt+1) - zuc(bdims_rem(3,1)+1) ) / dzc
+       
+                   alpha = ( (dzf/dzc)**2.0 - 1.0) / 24.0
+       
+                   eminus = eps * ( eps - 1.0 ) / 2.0 + alpha
+       
+                   edot  = ( 1.0 - eps**2.0 ) - 2.0 * alpha
+       
+                   eplus  = eps * ( eps + 1.0 ) / 2.0 + alpha
+       
+                   km (nzt+1,jf,if) = eminus * ptprs(bdims_rem(3,1),jf,if)   &
+                       + edot  * ptprs(bdims_rem(3,1)+1,jf,if) &
+                       + eplus  * ptprs(bdims_rem(3,1)+2,jf,if)
+       
+               END DO
+           END DO          
+           
+           DEALLOCATE ( ptprf, ptprs )
+       
+       
+#endif 
+       
+       END SUBROUTINE vnest_set_topbc_km
+    
+ 
+    END SUBROUTINE vnest_boundary_conds_khkm
+   
+   
+   
+    SUBROUTINE vnest_anterpolate
+   
+        !--------------------------------------------------------------------------------!
+        ! Description:
+        ! ------------
+        ! Anterpolate data from fine grid to coarse grid.
+        !------------------------------------------------------------------------------!
+   
+        USE arrays_3d
+        USE control_parameters
+        USE grid_variables
+        USE indices
+        USE interfaces
+        USE pegrid
+        USE surface_mod,                                                           &
+            ONLY :  bc_h
+        
+   
+        IMPLICIT NONE
+   
+        REAL(wp)    ::  time_since_reference_point_rem
+        INTEGER(iwp)   :: i, j, k, im, jn, ko
+
+!---    INTEGER(iwp) ::  j  !< grid index y direction
+!--     INTEGER(iwp) ::  k  !< grid index z direction
+        INTEGER(iwp) ::  kb !< variable to set respective boundary value, depends on facing. 
+        INTEGER(iwp) ::  l  !< running index boundary type, for up- and downward-facing walls
+        INTEGER(iwp) ::  m  !< running index surface elements
+   
+#if defined( __parallel )
+   
+   
+        !
+        !-- In case of model termination initiated by the remote model
+        !-- (terminate_coupled_remote > 0), initiate termination of the local model.
+        !-- The rest of the coupler must then be skipped because it would cause an MPI
+        !-- intercomminucation hang.
+        !-- If necessary, the coupler will be called at the beginning of the next
+        !-- restart run.
+   
+        IF ( myid == 0) THEN
+            CALL MPI_SENDRECV( terminate_coupled,        1, MPI_INTEGER, &
+                target_id, 0,                             &
+                terminate_coupled_remote, 1, MPI_INTEGER, &
+                target_id, 0,                             &
+                comm_inter, status, ierr )
+        ENDIF
+        CALL MPI_BCAST( terminate_coupled_remote, 1, MPI_INTEGER, 0, comm2d, &
+            ierr )
+   
+        IF ( terminate_coupled_remote > 0 )  THEN
+            WRITE( message_string, * ) 'remote model "',                         &
+                TRIM( coupling_mode_remote ),             &
+                '" terminated',                           &
+                '&with terminate_coupled_remote = ',      &
+                terminate_coupled_remote,                 &
+                '&local model  "', TRIM( coupling_mode ), &
+                '" has',                                  &
+                '&terminate_coupled = ',                  &
+                terminate_coupled
+            CALL message( 'vnest_anterpolate', 'PA0310', 1, 2, 0, 6, 0 )
+            RETURN
+        ENDIF
+   
+   
+        !
+        !-- Exchange the current simulated time between the models
+   
+        IF ( myid == 0 ) THEN
+   
+            CALL MPI_SEND( time_since_reference_point, 1, MPI_REAL, target_id, &
+                11, comm_inter, ierr )
+            CALL MPI_RECV( time_since_reference_point_rem, 1, MPI_REAL,        &
+                target_id, 11, comm_inter, status, ierr )
+   
+        ENDIF
+   
+        CALL MPI_BCAST( time_since_reference_point_rem, 1, MPI_REAL, 0, comm2d, &
+            ierr )
+   
+        IF ( coupling_mode == 'vnested_crse' )  THEN
+                      ! Receive data from fine grid for anterpolation
+   
+            offset(1) = ( pdims_partner(1) / pdims(1) ) * pcoord(1)
+            offset(2) = ( pdims_partner(2) / pdims(2) ) * pcoord(2)
+   
+            do j = 0,   ( pdims_partner(2) / pdims(2) ) - 1
+                do i = 0,   ( pdims_partner(1) / pdims(1) ) - 1
+                    map_coord(1) = i+offset(1)
+                    map_coord(2) = j+offset(2)
+   
+                    target_idex = f_rnk_lst(map_coord(1),map_coord(2)) + numprocs
+   
+                    CALL MPI_RECV( bdims_rem,   6, MPI_INTEGER, target_idex, 10, &
+                        comm_inter,status, ierr )
+   
+                    bdims (1,1) = bdims_rem (1,1) / cfratio(1)
+                    bdims (1,2) = bdims_rem (1,2) / cfratio(1)
+                    bdims (2,1) = bdims_rem (2,1) / cfratio(2)
+                    bdims (2,2) = bdims_rem (2,2) / cfratio(2)
+                    bdims (3,1) = bdims_rem (3,1)
+                    bdims (3,2) = bdims_rem (3,2) / cfratio(3)
+   
+                    CALL MPI_SEND( bdims, 6, MPI_INTEGER, target_idex, 9, &
+                        comm_inter, ierr )
+   
+                    n_cell_c =                      &
+                        (bdims(1,2)-bdims(1,1)+1) * &
+                        (bdims(2,2)-bdims(2,1)+1) * &
+                        (bdims(3,2)-bdims(3,1)+0)
+   
+                    CALL MPI_RECV( u(            &
+                        bdims(3,1)+1:bdims(3,2), &
+                        bdims(2,1)  :bdims(2,2), &
+                        bdims(1,1)  :bdims(1,2)),&
+                        n_cell_c, MPI_REAL, target_idex, 101, &
+                        comm_inter,status, ierr )
+   
+                    CALL MPI_RECV( v(            &
+                        bdims(3,1)+1:bdims(3,2), &
+                        bdims(2,1)  :bdims(2,2), &
+                        bdims(1,1)  :bdims(1,2)),&
+                        n_cell_c, MPI_REAL, target_idex, 102, &
+                        comm_inter,status, ierr )
+   
+                    CALL MPI_RECV(pt(            &
+                        bdims(3,1)+1:bdims(3,2), &
+                        bdims(2,1)  :bdims(2,2), &
+                        bdims(1,1)  :bdims(1,2)),&
+                        n_cell_c, MPI_REAL, target_idex, 105, &
+                        comm_inter,status, ierr )
+   
+   	        IF ( humidity ) THEN
+                    CALL MPI_RECV(q(            &
+                        bdims(3,1)+1:bdims(3,2), &
+                        bdims(2,1)  :bdims(2,2), &
+                        bdims(1,1)  :bdims(1,2)),&
+                        n_cell_c, MPI_REAL, target_idex, 106, &
+                        comm_inter,status, ierr )
+   		ENDIF
+   
+                    CALL MPI_RECV( w(              &
+                        bdims(3,1)  :bdims(3,2)-1, &
+                        bdims(2,1)  :bdims(2,2),   &
+                        bdims(1,1)  :bdims(1,2)),  &
+                        n_cell_c, MPI_REAL, target_idex, 103, &
+                        comm_inter,status, ierr )
+   
+                end do
+            end do
+   
+   
+   
+            !
+            !-- Boundary conditions for the velocity components u and v
+   
+   
+            IF ( ibc_uv_b == 0 ) THEN
+                u(nzb,:,:) = 0.0_wp 
+                v(nzb,:,:) = 0.0_wp 
+            ELSE
+                u(nzb,:,:) = u(nzb+1,:,:)
+                v(nzb,:,:) = v(nzb+1,:,:)
+            END IF
+            !
+            !-- Boundary conditions for the velocity components w
+
+            w(nzb,:,:) = 0.0_wp
+   
+!
+!-- Temperature at bottom boundary.
+!-- Neumann, zero-gradient
+            IF ( ibc_pt_b == 1 )  THEN 
+               DO  l = 0, 1 
+!
+!--       Set kb, for upward-facing surfaces value at topography top (k-1) is set,
+!--       for downward-facing surfaces at topography bottom (k+1). 
+                    kb = MERGE( -1, 1, l == 0 )
+                    DO  m = 1, bc_h(l)%ns
+                       i = bc_h(l)%i(m)     
+                       j = bc_h(l)%j(m)
+                       k = bc_h(l)%k(m)
+                       pt(k+kb,j,i) = pt(k,j,i)
+                    ENDDO
+                 ENDDO
+            ENDIF
+   
+   
+            CALL exchange_horiz( u, nbgp )
+            CALL exchange_horiz( v, nbgp )
+            CALL exchange_horiz( w, nbgp )
+            CALL exchange_horiz( pt, nbgp )
+   
+   
+        ELSEIF ( coupling_mode == 'vnested_fine' )  THEN
+                      ! Send data to coarse grid for anterpolation
+   
+            offset(1) = pcoord(1) / ( pdims(1)/pdims_partner(1) )
+            offset(2) = pcoord(2) / ( pdims(2)/pdims_partner(2) )
+            map_coord(1) = offset(1)
+            map_coord(2) = offset(2)
+            target_idex = c_rnk_lst(map_coord(1),map_coord(2))
+   
+!-- Limit anterpolation level to nzt - z nesting ratio (a pseudo-buffer layer)
+            bdims (1,1) = nxl
+            bdims (1,2) = nxr
+            bdims (2,1) = nys
+            bdims (2,2) = nyn
+            bdims (3,1) = nzb
+            bdims (3,2) = nzt-cfratio(3) 
+   
+            CALL MPI_SEND( bdims, 6, MPI_INTEGER, target_idex, 10, &
+                comm_inter, ierr )
+   
+            CALL MPI_RECV( bdims_rem,   6, MPI_INTEGER, target_idex, 9, &
+                comm_inter,status, ierr )
+   
+   
+            ALLOCATE( work3d (                   &
+                bdims_rem(3,1)+1:bdims_rem(3,2), &
+                bdims_rem(2,1)  :bdims_rem(2,2), &
+                bdims_rem(1,1)  :bdims_rem(1,2)))
+   
+   
+            anterpol3d => u
+   	   
+            CALL anterpolate_to_crse_u ( 101 )
+            CALL MPI_SEND( work3d, 1, TYPE_VNEST_ANTER, target_idex,  &
+                101, 	comm_inter, ierr)
+   
+            anterpol3d => v
+   
+            CALL anterpolate_to_crse_v ( 102 )
+            CALL MPI_SEND( work3d, 1, TYPE_VNEST_ANTER, target_idex,  &
+                102, 	comm_inter, ierr)
+   
+            anterpol3d => pt
+   
+            CALL anterpolate_to_crse_s ( 105 )
+            CALL MPI_SEND( work3d, 1, TYPE_VNEST_ANTER, target_idex,  &
+                105, 	comm_inter, ierr)
+   
+   
+          IF ( humidity ) THEN
+   
+            anterpol3d => q
+   
+            CALL anterpolate_to_crse_s ( 106 )
+            CALL MPI_SEND( work3d, 1, TYPE_VNEST_ANTER, target_idex,  &
+                106, 	comm_inter, ierr)
+          ENDIF
+   
+   
+            DEALLOCATE(   work3d                    )
+            ALLOCATE(   work3d ( bdims_rem(3,1)  :bdims_rem(3,2)-1, &
+                                 bdims_rem(2,1)  :bdims_rem(2,2), &
+                                 bdims_rem(1,1)  :bdims_rem(1,2)))
+            anterpol3d => w
+            CALL anterpolate_to_crse_w ( 103 )
+            CALL MPI_SEND( work3d, 1, TYPE_VNEST_ANTER, target_idex,  &
+                103, 	comm_inter, ierr)
+   
+            NULLIFY   (   anterpol3d                )
+            DEALLOCATE(   work3d                    )
+   
+        ENDIF
+   
+   
+#endif
+
+    CONTAINS
+       SUBROUTINE anterpolate_to_crse_u( tag )
+      
+#if defined( __parallel )
+      
+           USE arrays_3d
+           USE control_parameters
+           USE grid_variables
+           USE indices
+           USE pegrid
+           
+      
+           IMPLICIT NONE
+      
+           INTEGER(iwp)               ::  i, j, k
+           INTEGER(iwp)               ::  if, jf, kf
+           INTEGER(iwp)               ::  bottomx, topx
+           INTEGER(iwp)               ::  bottomy, topy
+           INTEGER(iwp)               ::  bottomz, topz
+           REAL(wp)                  ::  aweight
+           INTEGER(iwp), intent(in)   ::  tag
+      
+           !
+           !-- Anterpolation of the velocity components u
+           !-- only values in yz-planes that coincide in the fine and
+           !-- the coarse grid are considered
+      
+           DO k = bdims_rem(3,1)+1, bdims_rem(3,2)
+      
+               bottomz = (dzc/dzf) * (k-1) + 1
+               topz    = (dzc/dzf) * k
+      
+               DO j = bdims_rem(2,1),bdims_rem(2,2)
+       
+                   bottomy = (nyf+1) / (nyc+1) * j
+                   topy    = (nyf+1) / (nyc+1) * (j+1) - 1
+      
+                   DO i = bdims_rem(1,1),bdims_rem(1,2)
+      
+                       if = (nxf+1) / (nxc+1) * i
+      
+                       aweight   = 0.0
+      
+                       DO kf = bottomz, topz
+                           DO jf = bottomy, topy
+      
+                               aweight   = aweight + anterpol3d(kf,jf,if) *        &
+                                   (dzf/dzc) * (dyf/dyc)
+      
+                           END DO
+                       END DO
+      
+                       work3d(k,j,i)   = aweight
+      
+                   END DO
+      
+               END DO
+      
+           END DO
+      
+      
+#endif
+      
+       END SUBROUTINE anterpolate_to_crse_u
+      
+      
+       SUBROUTINE anterpolate_to_crse_v( tag )
+      
+#if defined( __parallel )
+      
+           USE arrays_3d
+           USE control_parameters
+           USE grid_variables
+           USE indices
+           USE pegrid
+           
+      
+           IMPLICIT NONE
+      
+           INTEGER(iwp)               ::  i, j, k
+           INTEGER(iwp)               ::  if, jf, kf
+           INTEGER(iwp)               ::  bottomx, topx
+           INTEGER(iwp)               ::  bottomy, topy
+           INTEGER(iwp)               ::  bottomz, topz
+           REAL(wp)                  ::  aweight
+           INTEGER(iwp), intent(in) ::  tag
+           !
+           !-- Anterpolation of the velocity components v
+           !-- only values in xz-planes that coincide in the fine and
+           !-- the coarse grid are considered
+      
+           DO k = bdims_rem(3,1)+1, bdims_rem(3,2)
+      
+               bottomz = (dzc/dzf) * (k-1) + 1
+               topz    = (dzc/dzf) * k
+      
+               DO j = bdims_rem(2,1), bdims_rem(2,2)
+      
+                   jf = (nyf+1) / (nyc+1) * j
+      
+                   DO i = bdims_rem(1,1), bdims_rem(1,2)
+      
+                       bottomx = (nxf+1) / (nxc+1) * i
+                       topx    = (nxf+1) / (nxc+1) * (i+1) - 1
+      
+                       aweight   = 0.0
+      
+                       DO kf = bottomz, topz
+                           DO if = bottomx, topx
+      
+                               aweight   = aweight + anterpol3d(kf,jf,if) *        &
+                                   (dzf/dzc) * (dxf/dxc)
+      
+      
+                           END DO
+                       END DO
+      
+                       work3d(k,j,i)   = aweight
+      
+                   END DO
+               END DO
+           END DO
+      
+      
+#endif
+      
+       END SUBROUTINE anterpolate_to_crse_v
+      
+      
+       SUBROUTINE anterpolate_to_crse_w( tag )
+      
+#if defined( __parallel )
+      
+           USE arrays_3d
+           USE control_parameters
+           USE grid_variables
+           USE indices
+           USE pegrid
+           
+      
+           IMPLICIT NONE
+      
+           INTEGER(iwp)               ::  i, j, k
+           INTEGER(iwp)               ::  if, jf, kf
+           INTEGER(iwp)               ::  bottomx, topx
+           INTEGER(iwp)               ::  bottomy, topy
+           INTEGER(iwp)               ::  bottomz, topz
+           REAL(wp)                  ::  aweight
+           INTEGER(iwp), intent(in) ::  tag
+           !
+           !-- Anterpolation of the velocity components w
+           !-- only values in xy-planes that coincide in the fine and
+           !-- the coarse grid are considered
+      
+           DO k = bdims_rem(3,1), bdims_rem(3,2)-1
+      
+               kf = cfratio(3) * k
+      
+               DO j = bdims_rem(2,1), bdims_rem(2,2)
+      
+                   bottomy = (nyf+1) / (nyc+1) * j
+                   topy    = (nyf+1) / (nyc+1) * (j+1) - 1
+      
+                   DO i = bdims_rem(1,1), bdims_rem(1,2)
+      
+                       bottomx = (nxf+1) / (nxc+1) * i
+                       topx    = (nxf+1) / (nxc+1) * (i+1) - 1
+      
+                       aweight   = 0.0
+      
+                       DO jf = bottomy, topy
+                           DO if = bottomx, topx
+      
+                               aweight   = aweight + anterpol3d (kf,jf,if) *        &
+                                   (dxf/dxc) * (dyf/dyc)
+      
+                           END DO
+                       END DO
+      
+                       work3d(k,j,i)   = aweight
+      
+                   END DO
+      
+               END DO
+      
+           END DO
+      
+#endif
+      
+       END SUBROUTINE anterpolate_to_crse_w
+      
+      
+       SUBROUTINE anterpolate_to_crse_s( tag )
+      
+#if defined( __parallel )
+      
+           USE arrays_3d
+           USE control_parameters
+           USE grid_variables
+           USE indices
+           USE pegrid
+           
+      
+           IMPLICIT NONE
+      
+           INTEGER(iwp)               ::  i, j, k
+           INTEGER(iwp)               ::  if, jf, kf
+           INTEGER(iwp)               ::  bottomx, topx
+           INTEGER(iwp)               ::  bottomy, topy
+           INTEGER(iwp)               ::  bottomz, topz
+           REAL(wp)                  ::  aweight
+           INTEGER(iwp), intent(in)   ::  tag
+      
+           !
+           !-- Anterpolation of the potential temperature pt
+           !-- all fine grid values are considered
+      
+           DO k = bdims_rem(3,1)+1, bdims_rem(3,2)
+      
+               bottomz = (dzc/dzf) * (k-1) + 1
+               topz    = (dzc/dzf) * k
+      
+               DO j = bdims_rem(2,1), bdims_rem(2,2)
+      
+                   bottomy = (nyf+1) / (nyc+1) * j
+                   topy    = (nyf+1) / (nyc+1) * (j+1) - 1
+      
+                   DO i = bdims_rem(1,1), bdims_rem(1,2)
+      
+                       bottomx = (nxf+1) / (nxc+1) * i
+                       topx    = (nxf+1) / (nxc+1) * (i+1) - 1
+      
+                       aweight   = 0.0
+          
+                       DO kf = bottomz, topz
+                           DO jf = bottomy, topy
+                               DO if = bottomx, topx
+      
+                                   aweight   = aweight  + anterpol3d(kf,jf,if) *                  &
+                                       (dzf/dzc) * (dyf/dyc) * (dxf/dxc)
+      
+                               END DO
+                           END DO
+                       END DO
+      
+                       work3d(k,j,i)   = aweight
+      
+                   END DO
+      
+               END DO
+      
+           END DO
+      
+#endif
+      
+       END SUBROUTINE anterpolate_to_crse_s
+    END SUBROUTINE vnest_anterpolate
+   
+   
+   
+    SUBROUTINE vnest_anterpolate_e
+   
+        !--------------------------------------------------------------------------------!
+        ! Description:
+        ! ------------
+        ! Anterpolate TKE from fine grid to coarse grid.
+        !------------------------------------------------------------------------------!
+   
+        USE arrays_3d
+        USE control_parameters
+        USE grid_variables
+        USE indices
+        USE interfaces
+        USE pegrid
+        
+   
+        IMPLICIT NONE
+   
+        REAL(wp)    ::  time_since_reference_point_rem
+        INTEGER(iwp)   :: i, j, k, im, jn, ko
+   
+#if defined( __parallel )
+   
+        !
+        !-- In case of model termination initiated by the remote model
+        !-- (terminate_coupled_remote > 0), initiate termination of the local model.
+        !-- The rest of the coupler must then be skipped because it would cause an MPI
+        !-- intercomminucation hang.
+        !-- If necessary, the coupler will be called at the beginning of the next
+        !-- restart run.
+   
+        IF ( myid == 0) THEN
+            CALL MPI_SENDRECV( terminate_coupled,        1, MPI_INTEGER, &
+                target_id, 0,                             &
+                terminate_coupled_remote, 1, MPI_INTEGER, &
+                target_id, 0,                             &
+                comm_inter, status, ierr )
+        ENDIF
+        CALL MPI_BCAST( terminate_coupled_remote, 1, MPI_INTEGER, 0, comm2d, &
+            ierr )
+   
+        IF ( terminate_coupled_remote > 0 )  THEN
+            WRITE( message_string, * ) 'remote model "',                         &
+                TRIM( coupling_mode_remote ),             &
+                '" terminated',                           &
+                '&with terminate_coupled_remote = ',      &
+                terminate_coupled_remote,                 &
+                '&local model  "', TRIM( coupling_mode ), &
+                '" has',                                  &
+                '&terminate_coupled = ',                  &
+                terminate_coupled
+            CALL message( 'vnest_anterpolate_e', 'PA0310', 1, 2, 0, 6, 0 )
+            RETURN
+        ENDIF
+   
+   
+        !
+        !-- Exchange the current simulated time between the models
+        IF ( myid == 0 ) THEN
+   
+            CALL MPI_SEND( time_since_reference_point, 1, MPI_REAL, target_id, &
+                11, comm_inter, ierr )
+            CALL MPI_RECV( time_since_reference_point_rem, 1, MPI_REAL,        &
+                target_id, 11, comm_inter, status, ierr )
+   
+        ENDIF
+   
+        CALL MPI_BCAST( time_since_reference_point_rem, 1, MPI_REAL, 0, comm2d, &
+            ierr )
+   
+        IF ( coupling_mode == 'vnested_crse' )  THEN
+                      ! Receive data from fine grid for anterpolation
+   
+            offset(1) = ( pdims_partner(1) / pdims(1) ) * pcoord(1)
+            offset(2) = ( pdims_partner(2) / pdims(2) ) * pcoord(2)
+   
+            do j = 0,   ( pdims_partner(2) / pdims(2) ) - 1
+                do i = 0,   ( pdims_partner(1) / pdims(1) ) - 1
+                    map_coord(1) = i+offset(1)
+                    map_coord(2) = j+offset(2)
+   
+                    target_idex = f_rnk_lst(map_coord(1),map_coord(2)) + numprocs
+  
+                    bdims (1,1) = f2c_dims_cg (0,map_coord(1),map_coord(2)) 
+                    bdims (1,2) = f2c_dims_cg (1,map_coord(1),map_coord(2)) 
+                    bdims (2,1) = f2c_dims_cg (2,map_coord(1),map_coord(2)) 
+                    bdims (2,2) = f2c_dims_cg (3,map_coord(1),map_coord(2)) 
+                    bdims (3,1) = f2c_dims_cg (4,map_coord(1),map_coord(2))
+                    bdims (3,2) = f2c_dims_cg (5,map_coord(1),map_coord(2)) 
+
+
+                     n_cell_c = (bdims(1,2)-bdims(1,1)+1) * &
+                         (bdims(2,2)-bdims(2,1)+1) * &
+                         (bdims(3,2)-bdims(3,1)+0)
+
+
+                    CALL MPI_RECV( e( bdims(3,1)+1:bdims(3,2), &
+                        bdims(2,1)  :bdims(2,2), &
+                        bdims(1,1)  :bdims(1,2)),&
+                        n_cell_c, MPI_REAL, target_idex, 104, &
+                        comm_inter,status, ierr )
+                end do
+            end do
+   
+   
+            !
+            !-- Boundary conditions
+   
+            IF ( .NOT. constant_diffusion ) THEN
+                e(nzb,:,:)   = e(nzb+1,:,:)
+            END IF
+   
+            IF ( .NOT. constant_diffusion )  CALL exchange_horiz( e, nbgp )
+   
+        ELSEIF ( coupling_mode == 'vnested_fine' )  THEN
+                      ! Send data to coarse grid for anterpolation
+   
+            offset(1) = pcoord(1) / ( pdims(1)/pdims_partner(1) )
+            offset(2) = pcoord(2) / ( pdims(2)/pdims_partner(2) )
+            map_coord(1) = offset(1)
+            map_coord(2) = offset(2)
+            target_idex = c_rnk_lst(map_coord(1),map_coord(2))
+
+            bdims_rem (1,1) =  f2c_dims_fg (0)
+            bdims_rem (1,2) =  f2c_dims_fg (1)
+            bdims_rem (2,1) =  f2c_dims_fg (2)
+            bdims_rem (2,2) =  f2c_dims_fg (3)
+            bdims_rem (3,1) =  f2c_dims_fg (4)
+            bdims_rem (3,2) =  f2c_dims_fg (5)
+
+            ALLOCATE( work3d (                   &
+                bdims_rem(3,1)+1:bdims_rem(3,2), &
+                bdims_rem(2,1)  :bdims_rem(2,2), &
+                bdims_rem(1,1)  :bdims_rem(1,2)))
+   
+            anterpol3d => e
+   
+            CALL anterpolate_to_crse_e ( 104 )
+
+            CALL MPI_SEND( work3d, 1, TYPE_VNEST_ANTER, target_idex,  &
+                           104, 	comm_inter, ierr)
+
+            NULLIFY   (   anterpol3d                )
+            DEALLOCATE(   work3d                    )
+        ENDIF
+   
+#endif
+   
+    CONTAINS 
+   
+   
+   
+   
+   
+       SUBROUTINE anterpolate_to_crse_e( tag )
+      
+#if defined( __parallel )
+      
+           USE arrays_3d
+           USE control_parameters
+           USE grid_variables
+           USE indices
+           USE pegrid
+           
+      
+           IMPLICIT NONE
+      
+           INTEGER(iwp)               ::  i, j, k
+           INTEGER(iwp)               ::  if, jf, kf
+           INTEGER(iwp)               ::  bottomx, topx
+           INTEGER(iwp)               ::  bottomy, topy
+           INTEGER(iwp)               ::  bottomz, topz
+           REAL(wp)                  ::  aweight_a, aweight_b, aweight_c, aweight_d, aweight_e
+           REAL(wp)                  ::  energ
+           INTEGER(iwp), intent(in)   ::  tag
+      
+          
+           DO k = bdims_rem(3,1)+1, bdims_rem(3,2)
+      
+               bottomz = (dzc/dzf) * (k-1) + 1
+               topz    = (dzc/dzf) * k
+             
+               DO j = bdims_rem(2,1), bdims_rem(2,2)
+      
+                   bottomy = (nyf+1) / (nyc+1) * j
+                   topy    = (nyf+1) / (nyc+1) * (j+1) - 1
+      
+                   DO i = bdims_rem(1,1), bdims_rem(1,2)
+      
+                       bottomx = (nxf+1) / (nxc+1) * i
+                       topx    = (nxf+1) / (nxc+1) * (i+1) - 1
+      
+                       aweight_a   = 0.0
+                       aweight_b   = 0.0
+                       aweight_c   = 0.0
+                       aweight_d   = 0.0
+                       aweight_e   = 0.0
+      
+                       DO kf = bottomz, topz
+                           DO jf = bottomy, topy
+                               DO if = bottomx, topx
+      
+                                   aweight_a = aweight_a + anterpol3d(kf,jf,if)  *                     &
+                                       (dzf/dzc) * (dyf/dyc) * (dxf/dxc)
+                              
+      
+                                   energ = ( 0.5 * ( u(kf,jf,if)   + u(kf,jf,if+1) ) )**2.0 +          &
+                                       ( 0.5 * ( v(kf,jf,if)   + v(kf,jf+1,if) ) )**2.0 +              &
+                                       ( 0.5 * ( w(kf-1,jf,if) + w(kf,jf,if) ) )**2.0
+      
+                                   aweight_b   =   aweight_b + energ         *                         &
+                                       (dzf/dzc) * (dyf/dyc) * (dxf/dxc)
+       
+                                   aweight_c   =   aweight_c + 0.5 * ( u(kf,jf,if) + u(kf,jf,if+1) ) * &
+                                       (dzf/dzc) * (dyf/dyc) * (dxf/dxc)
+      
+                                   aweight_d   =   aweight_d + 0.5 * ( v(kf,jf,if) + v(kf,jf+1,if) ) * &
+                                       (dzf/dzc) * (dyf/dyc) * (dxf/dxc)
+       
+                                   aweight_e   =   aweight_e + 0.5 * ( w(kf-1,jf,if) + w(kf,jf,if) ) * &
+                                       (dzf/dzc) * (dyf/dyc) * (dxf/dxc)
+      
+      
+                               END DO
+                           END DO
+                       END DO
+      
+                       work3d(k,j,i)  = aweight_a + 0.5 * ( aweight_b -  &
+                           aweight_c**2.0 -                              &
+                           aweight_d**2.0 -                              &
+                           aweight_e**2.0                 )
+      
+                   END DO
+      
+               END DO
+      
+           END DO
+      
+       
+#endif
+      
+       END SUBROUTINE anterpolate_to_crse_e
+    END SUBROUTINE vnest_anterpolate_e
+
+    SUBROUTINE vnest_init_pegrid_rank
+! Domain decomposition and exchange of grid variables between coarse and fine
+! Given processor coordinates as index f_rnk_lst(pcoord(1), pcoord(2))
+!   returns the rank. A single coarse block will have to send data to multiple
+!   fine blocks. In the coarse grid the pcoords of the remote block is first found and then using
+!   f_rnk_lst the target_idex is identified.
+! blk_dim stores the index limits of a given block. blk_dim_remote is received
+! from the asscoiated nest partner.
+! cf_ratio(1:3) is the ratio between fine and coarse grid: nxc/nxf, nyc/nyf and
+! ceiling(dxc/dxf) 
+
+
+       USE control_parameters,                                                    &
+           ONLY:  coupling_mode, coupling_mode_remote, coupling_topology, dz
+
+       USE grid_variables,                                                        &
+           ONLY:  dx, dy
+           
+       USE indices,                                                               &
+	   ONLY: nbgp, nx, ny, nz
+
+       USE kinds
+         
+       USE pegrid
+ 
+
+       IMPLICIT NONE
+
+       INTEGER(iwp)                    :: dest_rnk
+       INTEGER(iwp)               ::  i                        !<
+
+       IF (myid == 0) THEN
+           IF ( coupling_mode == 'vnested_crse') THEN
+               CALL MPI_SEND( pdims, 2, MPI_INTEGER, numprocs, 33, comm_inter, &
+                   ierr )
+               CALL MPI_RECV( pdims_partner, 2, MPI_INTEGER, numprocs, 66,      &
+                   comm_inter, status, ierr )
+           ELSEIF ( coupling_mode == 'vnested_fine') THEN
+               CALL MPI_RECV( pdims_partner, 2, MPI_INTEGER, 0, 33,      &
+                   comm_inter, status, ierr )
+               CALL MPI_SEND( pdims, 2, MPI_INTEGER, 0, 66, comm_inter, &
+                   ierr )
+           ENDIF
+       ENDIF
+
+
+       IF ( coupling_mode == 'vnested_crse') THEN
+           CALL MPI_BCAST( pdims_partner, 2, MPI_INTEGER, 0, comm2d, ierr )
+           ALLOCATE( c_rnk_lst( 0:(pdims(1)-1)    ,0:(pdims(2)-1) )   )
+           ALLOCATE( f_rnk_lst( 0:(pdims_partner(1)-1)  ,0:(pdims_partner(2)-1) ) )
+           do i=0,numprocs-1
+               CALL MPI_CART_COORDS( comm2d, i, ndim, pcoord, ierr )
+               call MPI_Cart_rank(comm2d, pcoord, dest_rnk, ierr)
+               c_rnk_lst(pcoord(1),pcoord(2)) = dest_rnk
+           end do
+       ELSEIF ( coupling_mode == 'vnested_fine') THEN
+           CALL MPI_BCAST( pdims_partner, 2, MPI_INTEGER, 0, comm2d, ierr )
+           ALLOCATE( c_rnk_lst( 0:(pdims_partner(1)-1)  ,0:(pdims_partner(2)-1) ) )
+           ALLOCATE( f_rnk_lst( 0:(pdims(1)-1)    ,0:(pdims(2)-1) )   )
+
+           do i=0,numprocs-1
+               CALL MPI_CART_COORDS( comm2d, i, ndim, pcoord, ierr )
+               call MPI_Cart_rank(comm2d, pcoord, dest_rnk, ierr)
+               f_rnk_lst(pcoord(1),pcoord(2)) = dest_rnk
+           enddo
+       ENDIF
+
+
+       IF ( coupling_mode == 'vnested_crse') THEN
+           if (myid == 0) then
+               CALL MPI_SEND( c_rnk_lst, pdims(1)*pdims(2),     MPI_INTEGER, numprocs, 0, comm_inter, ierr )
+               CALL MPI_RECV( f_rnk_lst, pdims_partner(1)*pdims_partner(2), MPI_INTEGER, numprocs, 4, comm_inter,status, ierr )
+           end if
+           CALL MPI_BCAST( f_rnk_lst, pdims_partner(1)*pdims_partner(2), MPI_INTEGER, 0, comm2d, ierr )
+       ELSEIF ( coupling_mode == 'vnested_fine') THEN
+           if (myid == 0) then
+               CALL MPI_RECV( c_rnk_lst, pdims_partner(1)*pdims_partner(2),  MPI_INTEGER, 0, 0, comm_inter,status, ierr )
+               CALL MPI_SEND( f_rnk_lst, pdims(1)*pdims(2),      MPI_INTEGER, 0, 4, comm_inter, ierr )
+           end if
+           CALL MPI_BCAST( c_rnk_lst, pdims_partner(1)*pdims_partner(2), MPI_INTEGER, 0, comm2d, ierr )
+       ENDIF
+
+       !-- Reason for MPI error unknown; solved if three lines duplicated
+       CALL MPI_CART_COORDS( comm2d, myid, ndim, pcoord, ierr )
+       CALL MPI_CART_SHIFT( comm2d, 0, 1, pleft, pright, ierr )
+       CALL MPI_CART_SHIFT( comm2d, 1, 1, psouth, pnorth, ierr )
+
+
+ 
+    END SUBROUTINE vnest_init_pegrid_rank
+
+
+    SUBROUTINE vnest_init_pegrid_domain
+
+       USE control_parameters,                                                    &
+           ONLY:  coupling_mode, coupling_mode_remote, coupling_topology, dz
+
+       USE grid_variables,                                                        &
+           ONLY:  dx, dy
+           
+       USE indices,                                                               &
+	   ONLY: nbgp, nx, ny, nz, nxl, nxr, nys, nyn, nzb, nzt, &
+	         nxlg, nxrg, nysg, nyng
+
+       USE kinds
+         
+       USE pegrid
+       
+       IMPLICIT NONE
+
+       INTEGER(iwp)              :: dest_rnk
+       INTEGER(iwp) ::  i, j                               !<
+       INTEGER(iwp) ::  tempx, tempy                               !<
+       INTEGER(iwp) :: TYPE_INT_YZ,  SIZEOFREAL
+       INTEGER(iwp) :: MTV_X,MTV_Y,MTV_Z,MTV_RX,MTV_RY,MTV_RZ
+   
+       !
+       !--    Pass the number of grid points of the coarse model to
+       !--    the nested model and vice versa
+       IF ( coupling_mode ==  'vnested_crse' )  THEN
+
+           nxc = nx
+           nyc = ny
+           nzc = nz
+           dxc = dx
+           dyc = dy
+           dzc = dz
+           cg_nprocs = numprocs
+
+           IF ( myid == 0 )  THEN
+
+               CALL MPI_SEND( nxc, 1, MPI_INTEGER  , numprocs, 1,  comm_inter,  &
+                   ierr )
+               CALL MPI_SEND( nyc, 1, MPI_INTEGER  , numprocs, 2,  comm_inter,  &
+                   ierr )
+               CALL MPI_SEND( nzc, 1, MPI_INTEGER  , numprocs, 3,  comm_inter,  &
+                   ierr )
+               CALL MPI_SEND( dxc, 1, MPI_REAL     , numprocs, 4,  comm_inter,  &
+                   ierr )
+               CALL MPI_SEND( dyc, 1, MPI_REAL     , numprocs, 5,  comm_inter,  &
+                   ierr )
+               CALL MPI_SEND( dzc, 1, MPI_REAL     , numprocs, 6,  comm_inter,  &
+                   ierr )
+               CALL MPI_SEND( pdims, 2, MPI_INTEGER, numprocs, 7,  comm_inter,  &
+                   ierr )
+               CALL MPI_SEND( cg_nprocs, 1, MPI_INTEGER, numprocs, 8,  comm_inter,  &
+                   ierr )
+               CALL MPI_RECV( nxf, 1, MPI_INTEGER,   numprocs, 21, comm_inter,  &
+                   status, ierr )
+               CALL MPI_RECV( nyf, 1, MPI_INTEGER,   numprocs, 22, comm_inter,  &
+                   status, ierr )
+               CALL MPI_RECV( nzf, 1, MPI_INTEGER,   numprocs, 23, comm_inter,  &
+                   status, ierr )
+               CALL MPI_RECV( dxf, 1, MPI_REAL,      numprocs, 24, comm_inter,  &
+                   status, ierr )
+               CALL MPI_RECV( dyf, 1, MPI_REAL,      numprocs, 25, comm_inter,  &
+                   status, ierr )
+               CALL MPI_RECV( dzf, 1, MPI_REAL,      numprocs, 26, comm_inter,  &
+                   status, ierr )
+               CALL MPI_RECV( pdims_partner, 2, MPI_INTEGER,                     &
+                   numprocs, 27, comm_inter, status, ierr )
+               CALL MPI_RECV( fg_nprocs, 1, MPI_INTEGER,                     &
+                   numprocs, 28, comm_inter, status, ierr )
+           ENDIF
+
+           CALL MPI_BCAST( nxf, 1,   MPI_INTEGER, 0, comm2d, ierr )
+           CALL MPI_BCAST( nyf, 1,   MPI_INTEGER, 0, comm2d, ierr )
+           CALL MPI_BCAST( nzf, 1,   MPI_INTEGER, 0, comm2d, ierr )
+           CALL MPI_BCAST( dxf, 1,   MPI_REAL,    0, comm2d, ierr )
+           CALL MPI_BCAST( dyf, 1,   MPI_REAL,    0, comm2d, ierr )
+           CALL MPI_BCAST( dzf, 1,   MPI_REAL,    0, comm2d, ierr )
+           CALL MPI_BCAST( pdims_partner, 2, MPI_INTEGER, 0, comm2d, ierr )
+           CALL MPI_BCAST( fg_nprocs,  1, MPI_INTEGER, 0, comm2d, ierr )
+
+       ELSEIF ( coupling_mode ==  'vnested_fine' )  THEN
+
+           nxf = nx
+           nyf = ny
+           nzf = nz
+           dxf = dx
+           dyf = dy
+           dzf = dz
+           fg_nprocs = numprocs
+
+           IF ( myid == 0 ) THEN
+
+               CALL MPI_RECV( nxc, 1, MPI_INTEGER, 0, 1, comm_inter, status, &
+                   ierr )
+               CALL MPI_RECV( nyc, 1, MPI_INTEGER, 0, 2, comm_inter, status, &
+                   ierr )
+               CALL MPI_RECV( nzc, 1, MPI_INTEGER, 0, 3, comm_inter, status, &
+                   ierr )
+               CALL MPI_RECV( dxc, 1, MPI_REAL,    0, 4, comm_inter, status, &
+                   ierr )
+               CALL MPI_RECV( dyc, 1, MPI_REAL,    0, 5, comm_inter, status, &
+                   ierr )
+               CALL MPI_RECV( dzc, 1, MPI_REAL,    0, 6, comm_inter, status, &
+                   ierr )
+               CALL MPI_RECV( pdims_partner, 2, MPI_INTEGER, 0, 7, comm_inter, &
+                   status, ierr )
+               CALL MPI_RECV( cg_nprocs,  1, MPI_INTEGER, 0, 8, comm_inter, &
+                   status, ierr )
+               CALL MPI_SEND( nxf, 1, MPI_INTEGER, 0, 21,  comm_inter, ierr )
+               CALL MPI_SEND( nyf, 1, MPI_INTEGER, 0, 22,  comm_inter, ierr )
+               CALL MPI_SEND( nzf, 1, MPI_INTEGER, 0, 23, comm_inter, ierr )
+               CALL MPI_SEND( dxf, 1, MPI_REAL,    0, 24, comm_inter, ierr )
+               CALL MPI_SEND( dyf, 1, MPI_REAL,    0, 25, comm_inter, ierr )
+               CALL MPI_SEND( dzf, 1, MPI_REAL,    0, 26, comm_inter, ierr )
+               CALL MPI_SEND( pdims,2,MPI_INTEGER, 0, 27, comm_inter, ierr )
+               CALL MPI_SEND( fg_nprocs,1,MPI_INTEGER, 0, 28, comm_inter, ierr )
+           ENDIF
+
+           CALL MPI_BCAST( nxc, 1, MPI_INTEGER, 0, comm2d, ierr)
+           CALL MPI_BCAST( nyc, 1, MPI_INTEGER, 0, comm2d, ierr)
+           CALL MPI_BCAST( nzc, 1, MPI_INTEGER, 0, comm2d, ierr)
+           CALL MPI_BCAST( dxc, 1, MPI_REAL,    0, comm2d, ierr)
+           CALL MPI_BCAST( dyc, 1, MPI_REAL,    0, comm2d, ierr)
+           CALL MPI_BCAST( dzc, 1, MPI_REAL,    0, comm2d, ierr)
+           CALL MPI_BCAST( pdims_partner, 2, MPI_INTEGER, 0, comm2d, ierr)
+           CALL MPI_BCAST( cg_nprocs,  1, MPI_INTEGER, 0, comm2d, ierr)
+
+       ENDIF
+  
+       ngp_c = ( nxc+1 + 2 * nbgp ) * ( nyc+1 + 2 * nbgp )
+       ngp_f = ( nxf+1 + 2 * nbgp ) * ( nyf+1 + 2 * nbgp )
+
+       IF ( coupling_mode(1:8) == 'vnested_')  coupling_topology = 1
+
+
+!-- Nesting Ratio: For each coarse grid cell how many fine grid cells exist
+       cfratio(1) = INT ( (nxf+1) / (nxc+1) )
+       cfratio(2) = INT ( (nyf+1) / (nyc+1) )
+       cfratio(3) = CEILING (  dzc    /  dzf  )
+
+!-- target_id is used only for exhange of information like simulated_time
+!-- which are then MPI_BCAST to other processors in the group
+          IF ( myid == 0 )  THEN
+   
+              IF ( TRIM( coupling_mode ) == 'vnested_crse' )  THEN
+                  target_id = numprocs
+              ELSE IF ( TRIM( coupling_mode ) == 'vnested_fine' )  THEN
+                  target_id = 0
+              ENDIF
+   
+          ENDIF
+
+!-- Store partner grid dimenstions and create MPI derived types
+
+         IF ( coupling_mode == 'vnested_crse' )  THEN
+
+            offset(1) = ( pdims_partner(1) / pdims(1) ) * pcoord(1)
+            offset(2) = ( pdims_partner(2) / pdims(2) ) * pcoord(2)
+      
+	    tempx =  ( pdims_partner(1) / pdims(1) ) - 1
+	    tempy =  ( pdims_partner(2) / pdims(2) ) - 1
+            ALLOCATE( c2f_dims_cg (0:5,offset(1):tempx+offset(1),offset(2):tempy+offset(2) ) )
+            ALLOCATE( f2c_dims_cg (0:5,offset(1):tempx+offset(1),offset(2):tempy+offset(2) ) )
+
+            do j = 0,   ( pdims_partner(2) / pdims(2) ) - 1
+                do i = 0,   ( pdims_partner(1) / pdims(1) ) - 1
+                   map_coord(1) = i+offset(1)
+                   map_coord(2) = j+offset(2)
+      
+                   target_idex = f_rnk_lst(map_coord(1),map_coord(2)) + numprocs
+
+                   CALL MPI_RECV( bdims_rem,   6, MPI_INTEGER, target_idex, 10, &
+                       comm_inter,status, ierr )
+
+!-- Store the CG dimensions that correspond to the FG partner; needed for FG top BC
+!-- One CG can have multiple FG partners. The 3D array is mapped by partner proc co-ord
+                   c2f_dims_cg (0,map_coord(1),map_coord(2)) = bdims_rem (1,1) / cfratio(1)
+                   c2f_dims_cg (1,map_coord(1),map_coord(2)) = bdims_rem (1,2) / cfratio(1)
+                   c2f_dims_cg (2,map_coord(1),map_coord(2)) = bdims_rem (2,1) / cfratio(2)
+                   c2f_dims_cg (3,map_coord(1),map_coord(2)) = bdims_rem (2,2) / cfratio(2)
+                   c2f_dims_cg (4,map_coord(1),map_coord(2)) = bdims_rem (3,2) / cfratio(3)
+                   c2f_dims_cg (5,map_coord(1),map_coord(2)) =(bdims_rem (3,2) / cfratio(3)) + 2
+ 
+!-- Store the CG dimensions that  correspond to the FG partner; needed for anterpolation
+                   f2c_dims_cg (0,map_coord(1),map_coord(2)) = bdims_rem (1,1) / cfratio(1)
+                   f2c_dims_cg (1,map_coord(1),map_coord(2)) = bdims_rem (1,2) / cfratio(1)
+                   f2c_dims_cg (2,map_coord(1),map_coord(2)) = bdims_rem (2,1) / cfratio(2)
+                   f2c_dims_cg (3,map_coord(1),map_coord(2)) = bdims_rem (2,2) / cfratio(2)
+                   f2c_dims_cg (4,map_coord(1),map_coord(2)) = bdims_rem (3,1)
+                   f2c_dims_cg (5,map_coord(1),map_coord(2)) =(bdims_rem (3,2)-cfratio(3))/ cfratio(3)
+   
+                   CALL MPI_SEND( c2f_dims_cg (:,map_coord(1),map_coord(2)), 6, &
+			           MPI_INTEGER, target_idex, 100, comm_inter, ierr )
+      
+                   CALL MPI_SEND( f2c_dims_cg (:,map_coord(1),map_coord(2)), 6, &
+			           MPI_INTEGER, target_idex, 101, comm_inter, ierr )
+      
+                end do
+            end do
+     
+!-- A derived data type to pack 3 Z-levels of CG to set FG top BC
+            MTV_X = ( nxr - nxl + 1 ) + 2*nbgp
+            MTV_Y = ( nyn - nys + 1 ) + 2*nbgp
+            MTV_Z = nzt+1 - nzb +1
+            
+            MTV_RX = ( c2f_dims_cg (1,offset(1),offset(2)) - c2f_dims_cg (0,offset(1),offset(2)) ) +1+2
+            MTV_RY = ( c2f_dims_cg (3,offset(1),offset(2)) - c2f_dims_cg (2,offset(1),offset(2)) ) +1+2
+            MTV_RZ = ( c2f_dims_cg (5,offset(1),offset(2)) - c2f_dims_cg (4,offset(1),offset(2)) ) +1
+            
+            CALL MPI_TYPE_EXTENT(MPI_REAL, SIZEOFREAL, IERR)
+
+            CALL MPI_TYPE_VECTOR ( MTV_RY, MTV_RZ, MTV_Z, MPI_REAL, TYPE_INT_YZ, IERR)
+            CALL MPI_TYPE_HVECTOR( MTV_RX, 1, MTV_Z*MTV_Y*SIZEOFREAL, &
+	                                      TYPE_INT_YZ, TYPE_VNEST_BC, IERR)
+            CALL MPI_TYPE_FREE(TYPE_INT_YZ, IERR)
+            CALL MPI_TYPE_COMMIT(TYPE_VNEST_BC, IERR)			   
+
+
+         ELSEIF ( coupling_mode == 'vnested_fine' )  THEN
+
+           ALLOCATE( c2f_dims_fg (0:5) )
+           ALLOCATE( f2c_dims_fg (0:5) )
+    
+           offset(1) = pcoord(1) / ( pdims(1)/pdims_partner(1) )
+           offset(2) = pcoord(2) / ( pdims(2)/pdims_partner(2) )
+           map_coord(1) = offset(1)
+           map_coord(2) = offset(2)
+           target_idex = c_rnk_lst(map_coord(1),map_coord(2))
+    
+           bdims (1,1) = nxl
+           bdims (1,2) = nxr
+           bdims (2,1) = nys
+           bdims (2,2) = nyn
+           bdims (3,1) = nzb
+           bdims (3,2) = nzt
+
+           CALL MPI_SEND( bdims,    6, MPI_INTEGER, target_idex, 10, &
+               comm_inter, ierr )
+           
+!-- Store the CG dimensions that correspond to the FG partner; needed for FG top BC
+!-- One FG can have only one CG partner
+           CALL MPI_RECV( c2f_dims_fg,   6, MPI_INTEGER, target_idex, 100, &
+               comm_inter,status, ierr )
+           
+           CALL MPI_RECV( f2c_dims_fg,   6, MPI_INTEGER, target_idex, 101, &
+               comm_inter,status, ierr )
+ 
+!-- Store the CG dimensions that  correspond to the FG partner; needed for anterpolation
+
+           n_cell_c = (f2c_dims_fg(1)-f2c_dims_fg(0)+1) * &
+                      (f2c_dims_fg(3)-f2c_dims_fg(2)+1) * &
+                      (f2c_dims_fg(5)-f2c_dims_fg(4)+0)
+
+           CALL MPI_TYPE_CONTIGUOUS(n_cell_c, MPI_REAL, TYPE_VNEST_ANTER, IERR)
+           CALL MPI_TYPE_COMMIT(TYPE_VNEST_ANTER, ierr) 
+
+        ENDIF
+   
+       END SUBROUTINE vnest_init_pegrid_domain
+
+
+       SUBROUTINE vnest_init_grid
+ 
+          USE arrays_3d,                                                             &
+              ONLY:  zu, zw
+              
+          USE control_parameters,                                                    &
+              ONLY:  coupling_mode
+             
+          USE indices,                                                               &
+              ONLY:  nzt
+         
+          USE kinds
+          
+          USE pegrid
+
+	  IMPLICIT NONE
+      
+          !-- Allocate and Exchange zuc and zuf, zwc and zwf
+          IF ( coupling_mode(1:8)  == 'vnested_' )  THEN
+          
+             ALLOCATE( zuc(0:nzc+1), zuf(0:nzf+1) )
+             ALLOCATE( zwc(0:nzc+1), zwf(0:nzf+1) )
+          
+             IF ( coupling_mode ==  'vnested_crse' )  THEN
+          
+                   zuc = zu
+                   zwc = zw
+                IF ( myid == 0 )  THEN
+          
+                   CALL MPI_SEND( zuc, nzt+2, MPI_REAL, numprocs, 41, comm_inter,  &
+                                  ierr )
+                   CALL MPI_RECV( zuf, nzf+2, MPI_REAL, numprocs, 42, comm_inter,  &
+                                  status, ierr )
+          
+                   CALL MPI_SEND( zwc, nzt+2, MPI_REAL, numprocs, 43, comm_inter,  &
+                                  ierr )
+                   CALL MPI_RECV( zwf, nzf+2, MPI_REAL, numprocs, 44, comm_inter,  &
+                                  status, ierr )
+          
+                ENDIF
+          
+                CALL MPI_BCAST( zuf,nzf+2,MPI_REAL,    0, comm2d, ierr ) 
+                CALL MPI_BCAST( zwf,nzf+2,MPI_REAL,    0, comm2d, ierr ) 
+          
+             ELSEIF ( coupling_mode ==  'vnested_fine' )  THEN
+          
+                   zuf = zu
+                   zwf = zw
+                IF ( myid == 0 )  THEN
+          
+                   CALL MPI_RECV( zuc,nzc+2, MPI_REAL, 0, 41, comm_inter, status, &
+                                  ierr )
+                   CALL MPI_SEND( zuf,nzt+2, MPI_REAL, 0, 42, comm_inter, ierr )
+          
+                   CALL MPI_RECV( zwc,nzc+2, MPI_REAL, 0, 43, comm_inter, status, &
+                                  ierr )
+                   CALL MPI_SEND( zwf,nzt+2, MPI_REAL, 0, 44, comm_inter, ierr )
+                ENDIF
+          
+                CALL MPI_BCAST( zuc,nzc+2,MPI_REAL,  0, comm2d, ierr ) 
+                CALL MPI_BCAST( zwc,nzc+2,MPI_REAL,  0, comm2d, ierr ) 
+          
+             ENDIF
+          ENDIF
+
+       END SUBROUTINE vnest_init_grid
+
+
+       SUBROUTINE vnest_check_parameters
+
+          USE arrays_3d,                                                             &
+              ONLY:  zu, zw
+
+          USE control_parameters,                                                    &
+              ONLY:  coupling_mode
+             
+          USE indices,                                                               &
+              ONLY:  nzt
+         
+          USE kinds
+          
+          USE pegrid
+
+	  IMPLICIT NONE
+
+
+	  IF (myid==0) PRINT*, '*** vnest: check parameters not implemented yet ***'
+      
+
+       END SUBROUTINE vnest_check_parameters
+
+
+       SUBROUTINE vnest_timestep_sync
+
+         USE control_parameters,                                                    &
+             ONLY:  coupling_mode, dt_3d, dt_coupling
+      
+         USE interfaces
+      
+         USE kinds
+
+         USE pegrid
+
+         IMPLICIT NONE
+
+         IF ( coupling_mode == 'vnested_crse')  THEN
+         	dtc = dt_3d
+         	if (myid == 0) then
+                   CALL MPI_SEND( dt_3d, 1, MPI_REAL, target_id,                     &
+                        31, comm_inter, ierr )
+                   CALL MPI_RECV( dtf, 1, MPI_REAL,                                  &
+                        target_id, 32, comm_inter, status, ierr )
+
+         	endif
+         	CALL MPI_BCAST( dtf, 1,   MPI_REAL,    0, comm2d, ierr ) 
+         ELSE
+         	dtf = dt_3d
+         	if (myid == 0) then
+                   CALL MPI_RECV( dtc, 1, MPI_REAL,                                  &
+                        target_id, 31, comm_inter, status, ierr )
+                   CALL MPI_SEND( dt_3d, 1, MPI_REAL, target_id,                     &
+                        32, comm_inter, ierr )
+
+         	endif
+         	CALL MPI_BCAST( dtc, 1,   MPI_REAL,    0, comm2d, ierr ) 
+         
+         ENDIF
+!-- Identical timestep for coarse and fine grids
+          dt_3d = MIN( dtc, dtf )
+!-- Nest coupling at every timestep
+          dt_coupling = dt_3d
+
+       END SUBROUTINE vnest_timestep_sync
+       
+       SUBROUTINE vnest_deallocate
+          USE control_parameters,                                                    &
+              ONLY:  coupling_mode
+             
+          IMPLICIT NONE
+          
+          IF ( ALLOCATED(c_rnk_lst) ) DEALLOCATE (c_rnk_lst)
+          IF ( ALLOCATED(f_rnk_lst) ) DEALLOCATE (f_rnk_lst)
+          
+          IF ( coupling_mode == 'vnested_crse')  THEN
+             IF ( ALLOCATED (c2f_dims_cg) ) DEALLOCATE (c2f_dims_cg)
+             IF ( ALLOCATED (f2c_dims_cg) ) DEALLOCATE (f2c_dims_cg)
+          ELSEIF( coupling_mode == 'vnested_fine' )  THEN
+             IF ( ALLOCATED (c2f_dims_fg) ) DEALLOCATE (c2f_dims_fg)
+             IF ( ALLOCATED (f2c_dims_fg) ) DEALLOCATE (f2c_dims_fg)
+          ENDIF
+       
+       END SUBROUTINE vnest_deallocate
+
+ END MODULE vertical_nesting_mod
Index: palm/trunk/SOURCE/write_var_list.f90
===================================================================
--- palm/trunk/SOURCE/write_var_list.f90	(revision 2364)
+++ palm/trunk/SOURCE/write_var_list.f90	(revision 2365)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Vertical nesting implemented (SadiqHuq)
+! 
+! 2339 2017-08-07 13:55:26Z gronemeier
 ! corrected timestamp in header
 ! 
@@ -224,4 +227,6 @@
                v_max_ijk, w_max, w_max_ijk
 
+    USE vertical_nesting_mod,                                                  &
+        ONLY:  vnest_init
 
     IMPLICIT NONE
@@ -695,4 +700,6 @@
     WRITE ( 14 )  'virtual_flight                '
     WRITE ( 14 )  virtual_flight
+    WRITE ( 14 )  'vnest_init                    '
+    WRITE ( 14 )  vnest_init
     WRITE ( 14 )  'volume_flow_area              '
     WRITE ( 14 )  volume_flow_area
Index: palm/trunk/UTIL/combine_plot_fields.f90
===================================================================
--- palm/trunk/UTIL/combine_plot_fields.f90	(revision 2364)
+++ palm/trunk/UTIL/combine_plot_fields.f90	(revision 2365)
@@ -21,8 +21,11 @@
 ! -----------------
 ! 
-!
+! 
 ! Former revisions:
 ! -----------------
 ! $Id$
+! Vertical grid nesting implemented (SadiqHuq)
+!
+! 1809 2016-04-05 20:13:28Z raasch
 !
 ! 1808 2016-04-05 19:44:00Z raasch
@@ -120,4 +123,5 @@
     LOGICAL  ::  avs_output, compressed, found, iso2d_output, netcdf_output,   &
                  netcdf_parallel, netcdf_0, netcdf_1
+    LOGICAL  ::  vnest
 
     REAL(wp) ::  cpu_start_time, cpu_end_time, dx, simulated_time
@@ -126,4 +130,6 @@
     REAL(spk), DIMENSION(:,:,:), ALLOCATABLE ::  pf3d, pf3d_tmp
 
+
+
     PRINT*, ''
     PRINT*, ''
@@ -133,7 +139,11 @@
 !-- Find out if a coupled run has been carried out
     INQUIRE( FILE='COUPLING_PORT_OPENED', EXIST=found )
+    INQUIRE( FILE='VNESTING_PORT_OPENED', EXIST=vnest )
     IF ( found )  THEN
        models = 2
        PRINT*, '    coupled run'
+    ELSEIF ( vnest )  THEN
+       models = 2
+       PRINT*, '    Vertically nested grid coupling'
     ELSE
        models = 1
@@ -160,10 +170,21 @@
           PRINT*, '*** combine_plot_fields ***'
           IF ( model == 2 )  THEN
-             model_string = '_O'
-             PRINT*, '    now combining ocean data'
-             PRINT*, '    ========================'
+             IF ( vnest )  THEN
+                model_string = '_N'
+                PRINT*, '    now combining FINE data'
+                PRINT*, '    ========================'
+             ELSE
+                model_string = '_O'
+                PRINT*, '    now combining ocean data'
+                PRINT*, '    ========================'
+             ENDIF
           ELSE
-             PRINT*, '    now combining atmosphere data'
-             PRINT*, '    ============================='
+             IF ( vnest )  THEN
+                PRINT*, '    now combining COARSE data'
+                PRINT*, '    ============================='
+             ELSE
+                PRINT*, '    now combining atmosphere data'
+                PRINT*, '    ============================='
+             ENDIF
           ENDIF
        ENDIF
