Index: /palm/trunk/DOC/app/chapter_4.2.html
===================================================================
--- /palm/trunk/DOC/app/chapter_4.2.html	(revision 101)
+++ /palm/trunk/DOC/app/chapter_4.2.html	(revision 102)
@@ -3,5 +3,5 @@
 <meta content="text/html; charset=ISO-8859-1" http-equiv="content-type"><title>PALM chapter 4.2</title></head>
 <body><h3 style="line-height: 100%;"><a name="Kapitel4.2"></a>4.2 <a href="#Laufparameter">Runtime
-parameters</a> and <a href="#Paketparameter">package
+parameters</a>, <a href="#particle_parameters">particle parameters</a>, and <a href="#Paketparameter">package
 parameters</a></h3>
 <h3 style="margin-bottom: 0cm; line-height: 100%;"><a name="Laufparameter"></a>Runtime parameters:</h3>
@@ -2093,9 +2093,7 @@
 is <a href="http://www.muk.uni-hannover.de/institut/software/iso2d_beschreibung.html#YRIGHT">yright</a>).</p>
 </td> </tr> </tbody></table><br>
-<br><h3 style="line-height: 100%;"><a name="Paketparameter"></a>Package
+<br><h3 style="line-height: 100%;"><a name="particle_parameters"></a>Particle
 parameters: </h3>
-Package (<span style="font-weight: bold;">mrun</span>
-option -p): <span style="font-weight: bold;"><a name="particles_package"></a>particles</span>&nbsp;&nbsp;&nbsp;&nbsp;
-NAMELIST group name: <span style="font-weight: bold;">particles_par<br>
+<span style="font-weight: bold;"></span><span style="font-weight: bold;"></span>NAMELIST group name: <span style="font-weight: bold;">particles_par<br>
 </span><table style="text-align: left; width: 100%;" border="1" cellpadding="2" cellspacing="2"> <tbody><tr>
 <td style="vertical-align: top;"><font size="4"><b>Parameter
@@ -2654,5 +2652,6 @@
 and performance of PALM may decrease.</p> </td> </tr>
 </tbody></table><span style="font-weight: bold;"><br>
-<br></span><span style="font-weight: bold;">Package
+</span><br><h3 style="line-height: 100%;"><a name="Paketparameter"></a>Package
+parameters: </h3><br><span style="font-weight: bold;">Package
 (<span style="font-weight: bold;">mrun</span> option
 -p): <span style="font-weight: bold;"><a name="dvrp_graphics"></a>dvrp_graphics</span>
Index: /palm/trunk/SCRIPTS/mrun
===================================================================
--- /palm/trunk/SCRIPTS/mrun	(revision 101)
+++ /palm/trunk/SCRIPTS/mrun	(revision 102)
@@ -120,4 +120,8 @@
      #                     removed (incl. option -f)
      # 19/06/07 - Siggi  - time limit for cdata jobs increased
+     # 25/07/07 - Siggi  - two sets of executables can be started on lcmuk,
+     #                     if new option -Y (coupling) is given,
+     #                     output of executables to aout_output removed,
+     #                     messages are immediately written to stdout instead
 
 
@@ -207,4 +211,5 @@
  remotecall=false
  remote_username=""
+ run_coupled_model=false
  run_mode=""
  scirocco=false
@@ -226,5 +231,5 @@
 
  typeset -i  iec=0 iic=0 iin=0 ioc=0 iout=0 memory=0 stagein_anz=0 stageout_anz=0
- typeset -i  cputime i ii icycle inode ival jobges jobsek maxcycle minuten nodes pes sekunden tp1
+ typeset -i  cputime i ii iii icycle inode ival jobges jobsek maxcycle minuten nodes pes sekunden tp1
 
  typeset  -R30 calltime
@@ -344,5 +349,5 @@
     # SHELLSCRIPT-OPTIONEN EINLESEN UND KOMMANDO NEU ZUSAMMENSETZEN, FALLS ES
     # FUER FOLGEJOBS BENOETIGT WIRD
- while  getopts  :a:AbBc:Cd:D:FG:h:H:i:IkK:m:M:n:o:Op:P:q:r:R:s:St:T:u:U:vxX: option
+ while  getopts  :a:AbBc:Cd:D:FG:h:H:i:IkK:m:M:n:o:Op:P:q:r:R:s:St:T:u:U:vxX:Y option
  do
    case  $option  in
@@ -382,4 +387,5 @@
        (x)   do_trace=true;set -x; mc="$mc -x";;
        (X)   numprocs=$OPTARG; mc="$mc -X$OPTARG";;
+       (Y)   run_coupled_model=true; mc="$mc -Y";;
        (\?)  printf "\n  +++ unknown option $OPTARG \n"
              printf "\n  --> type \"$0 ?\" for available options \n"
@@ -436,4 +442,5 @@
     printf "\n        -x    tracing of mrun for debug purposes       ---"
     printf "\n        -X    # of processors (on parallel machines)   1"
+    printf "\n        -Y    run coupled model                        ---"
     printf "\n "
     printf "\n      Possible values of positional parameter <modus>:"
@@ -2618,12 +2625,12 @@
                 if [[ "$cond1" = hpmcount  ||  "$cond2" = hpmcount ]]
                 then
-                   /opt/optibm/HPM_2_4_1/bin/hpmcount  a.out  $ROPTS  >  aout_output  2>&1
+                   /opt/optibm/HPM_2_4_1/bin/hpmcount  a.out  $ROPTS
                 else
-                   poe ./a.out  $ROPTS  >  aout_output  2>&1
+                   poe ./a.out  $ROPTS
                 fi
              else
                 if [[ $localhost = ibmb  ||  $localhost = ibmh  ||  $localhost = ibms ]]
                 then
-                   poe  a.out  -procs $numprocs  -nodes 1  -rmpool 0  $ROPTS  >  aout_output  2>&1
+                   poe  a.out  -procs $numprocs  -nodes 1  -rmpool 0  $ROPTS
                 elif [[ $localhost = ibmy ]]
                 then
@@ -2640,5 +2647,5 @@
                    fi
                    export MP_HOSTFILE=hostfile
-                   ./a.out  -procs $tasks_per_node  $ROPTS  >  aout_output  2>&1
+                   ./a.out  -procs $tasks_per_node  $ROPTS
                 else
                    if [[ "$host_file" = "" ]]
@@ -2651,5 +2658,5 @@
                    fi
                    export MP_HOSTFILE=$host_file
-                   poe  a.out  -procs $numprocs  -tasks_per_node $numprocs  $ROPTS  >  aout_output  2>&1
+                   poe  a.out  -procs $numprocs  -tasks_per_node $numprocs  $ROPTS
                 fi
              fi
@@ -2663,5 +2670,5 @@
                 #  export F_TRACE=YES|FMT1|FMT2  # output of ftrace informations to job protocol
                 echo "*** execution on single node with mpirun"
-                mpirun  -np $numprocs  ./a.out  $ROPTS  >  aout_output  2>&1
+                mpirun  -np $numprocs  ./a.out  $ROPTS
              else
                 (( i = 0 ))
@@ -2699,5 +2706,5 @@
                 echo "export use_seperate_pe_for_dvrp_output"  >>  mpi_exec_shell
                 echo " "                                 >>  mpi_exec_shell
-                echo "exec  ./a.out > aout_output  2>&1" >>  mpi_exec_shell
+                echo "exec  ./a.out"                     >>  mpi_exec_shell
 
                 chmod u+x  mpi_exec_shell
@@ -2731,14 +2738,27 @@
                 printf "\n      threads per task: $threads_per_task"
              fi
-             printf "\n\n"
-             if [[ $scirocco = true ]]
-             then
-                mpiexec  -machinefile hostfile  -n $ii  a.out  $ROPTS   >  aout_output  2>&1
+             if [[ $run_coupled_model = false ]]
+             then
+                printf "\n\n"
+                mpiexec  -machinefile hostfile  -n $ii  a.out  $ROPTS
              else
-                mpirun  -machinefile hostfile  -np $ii  a.out  $ROPTS   >  aout_output  2>&1
-             fi
+                ((  iii = ii / 2 ))
+                printf "\n      coupled run ($iii atmosphere, $iii ocean)"
+                printf "\n\n"
+                mpiexec  -machinefile hostfile  -n $iii  -env coupling_mode atmosphere_to_ocean  a.out  $ROPTS  &
+                mpiexec  -machinefile hostfile  -n $iii  -env coupling_mode ocean_to_atmosphere  a.out  $ROPTS  &
+                wait
+             fi
+
+#             if [[ $scirocco = true ]]
+#             then
+#                mpiexec  -machinefile hostfile  -n $ii  a.out  $ROPTS
+#             else
+#                mpirun  -machinefile hostfile  -np $ii  a.out  $ROPTS
+#             fi
+
           elif [[ $host = decalpha ]]
           then
-             dmpirun  -np $numprocs  a.out  $ROPTS  >  aout_output  2>&1
+             dmpirun  -np $numprocs  a.out  $ROPTS
           elif [[ $host = lctit ]]
           then
@@ -2760,12 +2780,12 @@
              if [[ "$QUEUE" = interactive ]]
              then
-                mpirun  -np $numprocs  a.out  $ROPTS  >  aout_output  2>&1
+                mpirun  -np $numprocs  a.out  $ROPTS
              else
-                n1ge  -fore  -g 1S060156  -mpi $numprocs  -mem 4  -N palm  -q $queue  a.out  $ROPTS  >  aout_output  2>&1
+                n1ge  -fore  -g 1S060156  -mpi $numprocs  -mem 4  -N palm  -q $queue  a.out  $ROPTS
              fi
 
              rm  $runfile
           else
-             mpprun  -n $numprocs  a.out  $ROPTS  >  aout_output  2>&1
+             mpprun  -n $numprocs  a.out  $ROPTS
           fi
           [[ $? != 0 ]]  &&  execution_error=true
@@ -2778,5 +2798,5 @@
           fi
        else
-          a.out  $ROPTS  >  aout_output  2>&1
+          a.out  $ROPTS
        fi
     fi
@@ -2785,10 +2805,10 @@
 
           # ABBRUCH BEI LAUFZEITFEHLER
-       [[ ! ( "$cond1" = debug  ||  "$cond2" = debug ) ]]  &&  cat  aout_output
+#       [[ ! ( "$cond1" = debug  ||  "$cond2" = debug ) ]]  &&  cat  aout_output*
        printf "\n  +++ runtime error occured"
        locat=execution
        exit
     else
-       [[ ! ( "$cond1" = debug  ||  "$cond2" = debug ) ]]  &&  cat  aout_output
+#       [[ ! ( "$cond1" = debug  ||  "$cond2" = debug ) ]]  &&  cat  aout_output*
        printf "\n$striche\n  *** execution finished \n"
     fi
@@ -3544,4 +3564,5 @@
     [[ $delete_temporary_catalog = false ]]  &&  mrun_com=${mrun_com}" -B"
     [[ $node_usage != default ]]  &&  mrun_com=${mrun_com}" -n $node_usage"
+    [[ $run_coupled_model = true ]]  &&  mrun_com=${mrun_com}" -Y"
     if [[ $do_remote = true ]]
     then
Index: /palm/trunk/SCRIPTS/subjob
===================================================================
--- /palm/trunk/SCRIPTS/subjob	(revision 101)
+++ /palm/trunk/SCRIPTS/subjob	(revision 102)
@@ -82,4 +82,5 @@
      # 01/03/07 - Siggi - adjustments for RIAM machines gate and NEC-SX8 (n-sx)
      # 12/04/07 - Siggi - option -f (filetransfer protocol) removed, scp only
+     # 27/07/07 - Siggi - autan admitted
 
 
@@ -144,4 +145,5 @@
  case  $local_host  in
      (atmos)                 local_addres=172.20.25.35;   local_host=lcide;;
+     (autan)                 local_addres=130.75.105.57;  local_host=lcmuk;;
      (bora)                  local_addres=130.75.105.103; local_host=lcmuk;;
      (breg*-en0|berni*-en0)  local_addres=130.73.230.10;  local_host=ibmb;;
Index: /palm/trunk/SOURCE/CURRENT_MODIFICATIONS
===================================================================
--- /palm/trunk/SOURCE/CURRENT_MODIFICATIONS	(revision 101)
+++ /palm/trunk/SOURCE/CURRENT_MODIFICATIONS	(revision 102)
@@ -1,4 +1,20 @@
 New:
 ---
+Allows runs for a coupled atmosphere-ocean LES,
+coupling frequency is controlled by new d3par-parameter dt_coupling,
+the coupling mode (atmosphere_to_ocean or ocean_to_atmosphere) for the
+respective processes is read from environment variable coupling_mode,
+which is set by the mpiexec-command,
+communication between the two models is done using the intercommunicator
+comm_inter,
+local files opened by the ocean model get the additional suffic "_O"
+
+A momentum flux can be set as top boundary condition using the new
+inipar parameter top_momentumflux_u|v.
+
+boundary_conds, check_open, check_parameters, diffusion_u, diffusion_v, flow_statistics, header, init_pegrid, init_3d_model, modules, palm, parin, prognostic_equations, read_var_list, read_3d_binary, swap_timelevel, time_integration, write_var_list, write_3d_binary
+
+New:
+surface_coupler
 
 
Index: /palm/trunk/SOURCE/Makefile
===================================================================
--- /palm/trunk/SOURCE/Makefile	(revision 101)
+++ /palm/trunk/SOURCE/Makefile	(revision 102)
@@ -4,9 +4,12 @@
 # Actual revisions:
 # -----------------
-# +eqn_state_seawater, init_ocean
+# +surface_coupler
 #
 # Former revisions:
 # -----------------
 # $Id$
+#
+# 96 2007-06-04 08:07:41Z raasch
+# +eqn_state_seawater, init_ocean
 #
 # 82 2007-04-16 15:40:52Z raasch
@@ -61,6 +64,6 @@
 	random_gauss.f90 read_3d_binary.f90 read_var_list.f90 run_control.f90 \
 	singleton.f90 sor.f90 spline_x.f90 spline_y.f90 spline_z.f90 \
-	sum_up_3d_data.f90 swap_timelevel.f90 temperton_fft.f90 \
-	time_integration.f90 time_to_string.f90 timestep.f90 \
+	sum_up_3d_data.f90 surface_coupler.f90 swap_timelevel.f90 \
+	temperton_fft.f90 time_integration.f90 time_to_string.f90 timestep.f90 \
 	timestep_scheme_steering.f90 transpose.f90 user_interface.f90 \
 	wall_fluxes.f90 write_3d_binary.f90 write_compressed.f90 \
@@ -90,7 +93,7 @@
 	production_e.o prognostic_equations.o random_function.o random_gauss.o \
 	read_3d_binary.o read_var_list.o run_control.o singleton.o sor.o \
-	spline_x.o spline_y.o spline_z.o sum_up_3d_data.o swap_timelevel.o \
-	temperton_fft.o time_integration.o time_to_string.o timestep.o \
-	timestep_scheme_steering.o transpose.o user_interface.o \
+	spline_x.o spline_y.o spline_z.o sum_up_3d_data.o surface_coupler.o \
+	swap_timelevel.o temperton_fft.o time_integration.o time_to_string.o \
+        timestep.o timestep_scheme_steering.o transpose.o user_interface.o \
 	wall_fluxes.o write_3d_binary.o write_compressed.o write_var_list.o
 
@@ -221,4 +224,5 @@
 spline_z.o: modules.o
 sum_up_3d_data.o: modules.o
+surface_coupler.o: modules.o
 swap_timelevel.o: modules.o
 temperton_fft.o: modules.o
Index: /palm/trunk/SOURCE/boundary_conds.f90
===================================================================
--- /palm/trunk/SOURCE/boundary_conds.f90	(revision 101)
+++ /palm/trunk/SOURCE/boundary_conds.f90	(revision 102)
@@ -4,5 +4,5 @@
 ! Actual revisions:
 ! -----------------
-! 
+! Boundary conditions for temperature adjusted for coupled runs
 !
 ! Former revisions:
@@ -93,5 +93,7 @@
 
 !
-!--    Temperature at bottom boundary
+!--    Temperature at bottom boundary.
+!--    In case of coupled runs (ibc_pt_b = 2) the temperature is given by
+!--    the sea surface temperature of the coupled ocean model.
        IF ( ibc_pt_b == 0 )  THEN
           DO  i = nxl-1, nxr+1
@@ -100,5 +102,5 @@
              ENDDO
           ENDDO
-       ELSE
+       ELSEIF ( ibc_pt_b == 1 )  THEN
           DO  i = nxl-1, nxr+1
              DO  j = nys-1, nyn+1
Index: /palm/trunk/SOURCE/check_open.f90
===================================================================
--- /palm/trunk/SOURCE/check_open.f90	(revision 101)
+++ /palm/trunk/SOURCE/check_open.f90	(revision 102)
@@ -4,5 +4,6 @@
 ! Actual revisions:
 ! -----------------
-! 
+! Different filenames are used in case of a coupled simulation,
+! coupling_char added to all relevant filenames
 !
 ! Former revisions:
@@ -146,14 +147,15 @@
        CASE ( 11 )
 
-          OPEN ( 11, FILE='PARIN', FORM='FORMATTED', STATUS='OLD' )
+          OPEN ( 11, FILE='PARIN'//coupling_char, FORM='FORMATTED', &
+                     STATUS='OLD' )
 
        CASE ( 13 )
 
           IF ( myid_char == '' )  THEN
-             OPEN ( 13, FILE='BININ'//myid_char, FORM='UNFORMATTED', &
-                        STATUS='OLD' )
-          ELSE
-             OPEN ( 13, FILE='BININ/'//myid_char, FORM='UNFORMATTED', &
-                        STATUS='OLD' )
+             OPEN ( 13, FILE='BININ'//coupling_char//myid_char, &
+                        FORM='UNFORMATTED', STATUS='OLD' )
+          ELSE
+             OPEN ( 13, FILE='BININ'//coupling_char//'/'//myid_char, &
+                        FORM='UNFORMATTED', STATUS='OLD' )
           ENDIF
 
@@ -161,9 +163,9 @@
 
           IF ( myid_char == '' )  THEN
-             OPEN ( 14, FILE='BINOUT'//myid_char, FORM='UNFORMATTED', &
-                        POSITION='APPEND' )
+             OPEN ( 14, FILE='BINOUT'//coupling_char//myid_char, &
+                        FORM='UNFORMATTED', POSITION='APPEND' )
           ELSE
              IF ( myid == 0  .AND. .NOT. openfile(file_id)%opened_before )  THEN
-                CALL local_system( 'mkdir  BINOUT' )
+                CALL local_system( 'mkdir  BINOUT' // coupling_char )
              ENDIF
 #if defined( __parallel )
@@ -173,36 +175,36 @@
              CALL MPI_BARRIER( comm2d, ierr )
 #endif
-             OPEN ( 14, FILE='BINOUT/'//myid_char_14, FORM='UNFORMATTED', &
-                        POSITION='APPEND' )
+             OPEN ( 14, FILE='BINOUT'//coupling_char//'/'//myid_char_14, &
+                        FORM='UNFORMATTED', POSITION='APPEND' )
           ENDIF
 
        CASE ( 15 )
 
-          OPEN ( 15, FILE='RUN_CONTROL', FORM='FORMATTED' )
+          OPEN ( 15, FILE='RUN_CONTROL'//coupling_char, FORM='FORMATTED' )
 
        CASE ( 16 )
 
-          OPEN ( 16, FILE='LIST_PROFIL', FORM='FORMATTED' )
+          OPEN ( 16, FILE='LIST_PROFIL'//coupling_char, FORM='FORMATTED' )
 
        CASE ( 17 )
 
-          OPEN ( 17, FILE='LIST_PROFIL_1D', FORM='FORMATTED' )
+          OPEN ( 17, FILE='LIST_PROFIL_1D'//coupling_char, FORM='FORMATTED' )
 
        CASE ( 18 )
 
-          OPEN ( 18, FILE='CPU_MEASURES', FORM='FORMATTED' )
+          OPEN ( 18, FILE='CPU_MEASURES'//coupling_char, FORM='FORMATTED' )
 
        CASE ( 19 )
 
-          OPEN ( 19, FILE='HEADER', FORM='FORMATTED' )
+          OPEN ( 19, FILE='HEADER'//coupling_char, FORM='FORMATTED' )
 
        CASE ( 20 )
 
           IF ( myid == 0  .AND. .NOT. openfile(file_id)%opened_before )  THEN
-             CALL local_system( 'mkdir  DATA_LOG' )
+             CALL local_system( 'mkdir  DATA_LOG' // coupling_char )
           ENDIF
           IF ( myid_char == '' )  THEN
-             OPEN ( 20, FILE='DATA_LOG/_0000', FORM='UNFORMATTED', &
-                        POSITION='APPEND' )
+             OPEN ( 20, FILE='DATA_LOG'//coupling_char//'/_0000', &
+                        FORM='UNFORMATTED', POSITION='APPEND' )
           ELSE
 #if defined( __parallel )
@@ -212,6 +214,6 @@
              CALL MPI_BARRIER( comm2d, ierr )
 #endif
-             OPEN ( 20, FILE='DATA_LOG/'//myid_char, FORM='UNFORMATTED', &
-                        POSITION='APPEND' )
+             OPEN ( 20, FILE='DATA_LOG'//coupling_char//'/'//myid_char, &
+                        FORM='UNFORMATTED', POSITION='APPEND' )
           ENDIF
 
@@ -219,9 +221,9 @@
 
           IF ( data_output_2d_on_each_pe )  THEN
-             OPEN ( 21, FILE='PLOT2D_XY'//myid_char, FORM='UNFORMATTED', &
-                        POSITION='APPEND' )
-          ELSE
-             OPEN ( 21, FILE='PLOT2D_XY', FORM='UNFORMATTED', &
-                        POSITION='APPEND' )
+             OPEN ( 21, FILE='PLOT2D_XY'//coupling_char//myid_char, &
+                        FORM='UNFORMATTED', POSITION='APPEND' )
+          ELSE
+             OPEN ( 21, FILE='PLOT2D_XY'//coupling_char, &
+                        FORM='UNFORMATTED', POSITION='APPEND' )
           ENDIF
 
@@ -250,6 +252,6 @@
 !--          Create output file for local parameters
              IF ( iso2d_output )  THEN
-                OPEN ( 27, FILE='PLOT2D_XY_LOCAL', FORM='FORMATTED', &
-                           DELIM='APOSTROPHE' )
+                OPEN ( 27, FILE='PLOT2D_XY_LOCAL'//coupling_char, &
+                           FORM='FORMATTED', DELIM='APOSTROPHE' )
                 openfile(27)%opened = .TRUE.
              ENDIF
@@ -260,8 +262,8 @@
 
           IF ( data_output_2d_on_each_pe )  THEN
-             OPEN ( 22, FILE='PLOT2D_XZ'//myid_char, FORM='UNFORMATTED', &
-                        POSITION='APPEND' )
-          ELSE
-             OPEN ( 22, FILE='PLOT2D_XZ', FORM='UNFORMATTED', &
+             OPEN ( 22, FILE='PLOT2D_XZ'//coupling_char//myid_char, &
+                        FORM='UNFORMATTED', POSITION='APPEND' )
+          ELSE
+             OPEN ( 22, FILE='PLOT2D_XZ'//coupling_char, FORM='UNFORMATTED', &
                         POSITION='APPEND' )
           ENDIF
@@ -289,6 +291,6 @@
 !
 !--          Create output file for local parameters
-             OPEN ( 28, FILE='PLOT2D_XZ_LOCAL', FORM='FORMATTED', &
-                        DELIM='APOSTROPHE' )
+             OPEN ( 28, FILE='PLOT2D_XZ_LOCAL'//coupling_char, &
+                        FORM='FORMATTED', DELIM='APOSTROPHE' )
              openfile(28)%opened = .TRUE.
 
@@ -298,8 +300,8 @@
 
           IF ( data_output_2d_on_each_pe )  THEN
-             OPEN ( 23, FILE='PLOT2D_YZ'//myid_char, FORM='UNFORMATTED', &
-                        POSITION='APPEND' )
-          ELSE
-             OPEN ( 23, FILE='PLOT2D_YZ', FORM='UNFORMATTED', &
+             OPEN ( 23, FILE='PLOT2D_YZ'//coupling_char//myid_char, &
+                        FORM='UNFORMATTED', POSITION='APPEND' )
+          ELSE
+             OPEN ( 23, FILE='PLOT2D_YZ'//coupling_char, FORM='UNFORMATTED', &
                         POSITION='APPEND' )
           ENDIF
@@ -327,6 +329,6 @@
 !
 !--          Create output file for local parameters
-             OPEN ( 29, FILE='PLOT2D_YZ_LOCAL', FORM='FORMATTED', &
-                        DELIM='APOSTROPHE' )
+             OPEN ( 29, FILE='PLOT2D_YZ_LOCAL'//coupling_char, &
+                        FORM='FORMATTED', DELIM='APOSTROPHE' )
              openfile(29)%opened = .TRUE.
 
@@ -335,5 +337,6 @@
        CASE ( 30 )
 
-          OPEN ( 30, FILE='PLOT3D_DATA'//myid_char, FORM='UNFORMATTED' )
+          OPEN ( 30, FILE='PLOT3D_DATA'//coupling_char//myid_char, &
+                     FORM='UNFORMATTED' )
 !
 !--       Write coordinate file for AVS
@@ -508,5 +511,6 @@
              WRITE ( suffix, '(''_'',I1)' )  file_id - 40
           ENDIF
-          OPEN ( file_id, FILE='PLOT1D_DATA'//TRIM( suffix ), FORM='FORMATTED' )
+          OPEN ( file_id, FILE='PLOT1D_DATA'//coupling_char//TRIM( suffix ), &
+                          FORM='FORMATTED' )
 !
 !--       Write contents comments at the top of the file
@@ -521,6 +525,6 @@
              WRITE ( suffix, '(''_'',I1)' )  file_id - 50
           ENDIF
-          OPEN ( file_id, FILE='PLOTTS_DATA'//TRIM( suffix ), FORM='FORMATTED',&
-                          RECL=496 )
+          OPEN ( file_id, FILE='PLOTTS_DATA'//coupling_char//TRIM( suffix ), &
+                          FORM='FORMATTED', RECL=496 )
 !
 !--       Write PROFIL parameter file for output of time series
@@ -539,6 +543,6 @@
 !
 !--       Write RAHMEN parameter
-          OPEN ( 90, FILE='PLOTTS_PAR'//TRIM( suffix ), FORM='FORMATTED', &
-                     DELIM='APOSTROPHE' )
+          OPEN ( 90, FILE='PLOTTS_PAR'//coupling_char//TRIM( suffix ), &
+                     FORM='FORMATTED', DELIM='APOSTROPHE' )
           WRITE ( 90, RAHMEN )
 !
@@ -575,9 +579,9 @@
 
           IF ( myid_char == '' )  THEN
-             OPEN ( 80, FILE='PARTICLE_INFOS'//myid_char, FORM='FORMATTED', &
-                        POSITION='APPEND' )
+             OPEN ( 80, FILE='PARTICLE_INFOS'//coupling_char//myid_char, &
+                        FORM='FORMATTED', POSITION='APPEND' )
           ELSE
              IF ( myid == 0  .AND.  .NOT. openfile(80)%opened_before )  THEN
-                CALL local_system( 'mkdir  PARTICLE_INFOS' )
+                CALL local_system( 'mkdir  PARTICLE_INFOS' // coupling_char )
              ENDIF
 #if defined( __parallel )
@@ -592,6 +596,6 @@
              ENDIF
 #endif
-             OPEN ( 80, FILE='PARTICLE_INFOS/'//myid_char, FORM='FORMATTED', &
-                        POSITION='APPEND' )
+             OPEN ( 80, FILE='PARTICLE_INFOS'//coupling_char//'/'//myid_char, &
+                        FORM='FORMATTED', POSITION='APPEND' )
           ENDIF
 
@@ -602,20 +606,20 @@
        CASE ( 81 )
 
-             OPEN ( 81, FILE='PLOTSP_X_PAR', FORM='FORMATTED', &
+             OPEN ( 81, FILE='PLOTSP_X_PAR'//coupling_char, FORM='FORMATTED', &
                         DELIM='APOSTROPHE', RECL=1500, POSITION='APPEND' )
 
        CASE ( 82 )
 
-             OPEN ( 82, FILE='PLOTSP_X_DATA', FORM='FORMATTED', &
+             OPEN ( 82, FILE='PLOTSP_X_DATA'//coupling_char, FORM='FORMATTED', &
                         POSITION = 'APPEND' )
 
        CASE ( 83 )
 
-             OPEN ( 83, FILE='PLOTSP_Y_PAR', FORM='FORMATTED', &
+             OPEN ( 83, FILE='PLOTSP_Y_PAR'//coupling_char, FORM='FORMATTED', &
                         DELIM='APOSTROPHE', RECL=1500, POSITION='APPEND' )
 
        CASE ( 84 )
 
-             OPEN ( 84, FILE='PLOTSP_Y_DATA', FORM='FORMATTED', &
+             OPEN ( 84, FILE='PLOTSP_Y_DATA'//coupling_char, FORM='FORMATTED', &
                         POSITION='APPEND' )
 
@@ -623,9 +627,9 @@
 
           IF ( myid_char == '' )  THEN
-             OPEN ( 85, FILE='PARTICLE_DATA'//myid_char, FORM='UNFORMATTED', &
-                        POSITION='APPEND' )
+             OPEN ( 85, FILE='PARTICLE_DATA'//coupling_char//myid_char, &
+                        FORM='UNFORMATTED', POSITION='APPEND' )
           ELSE
              IF ( myid == 0  .AND.  .NOT. openfile(85)%opened_before )  THEN
-                CALL local_system( 'mkdir  PARTICLE_DATA' )
+                CALL local_system( 'mkdir  PARTICLE_DATA' // coupling_char )
              ENDIF
 #if defined( __parallel )
@@ -635,6 +639,6 @@
              CALL MPI_BARRIER( comm2d, ierr )
 #endif
-             OPEN ( 85, FILE='PARTICLE_DATA/'//myid_char, FORM='UNFORMATTED', &
-                        POSITION='APPEND' )
+             OPEN ( 85, FILE='PARTICLE_DATA'//coupling_char//'/'//myid_char, &
+                        FORM='UNFORMATTED', POSITION='APPEND' )
           ENDIF
 
@@ -656,8 +660,8 @@
 !--       Set filename depending on unit number
           IF ( file_id == 101 )  THEN
-             filename = 'DATA_2D_XY_NETCDF'
+             filename = 'DATA_2D_XY_NETCDF' // coupling_char
              av = 0
           ELSE
-             filename = 'DATA_2D_XY_AV_NETCDF'
+             filename = 'DATA_2D_XY_AV_NETCDF' // coupling_char
              av = 1
           ENDIF
@@ -716,8 +720,8 @@
 !--       Set filename depending on unit number
           IF ( file_id == 102 )  THEN
-             filename = 'DATA_2D_XZ_NETCDF'
+             filename = 'DATA_2D_XZ_NETCDF' // coupling_char
              av = 0
           ELSE
-             filename = 'DATA_2D_XZ_AV_NETCDF'
+             filename = 'DATA_2D_XZ_AV_NETCDF' // coupling_char
              av = 1
           ENDIF
@@ -776,8 +780,8 @@
 !--       Set filename depending on unit number
           IF ( file_id == 103 )  THEN
-             filename = 'DATA_2D_YZ_NETCDF'
+             filename = 'DATA_2D_YZ_NETCDF' // coupling_char
              av = 0
           ELSE
-             filename = 'DATA_2D_YZ_AV_NETCDF'
+             filename = 'DATA_2D_YZ_AV_NETCDF' // coupling_char
              av = 1
           ENDIF
@@ -834,12 +838,16 @@
        CASE ( 104 )
 !
+!--       Set filename
+          filename = 'DATA_1D_PR_NETCDF' // coupling_char
+
+!
 !--       Inquire, if there is a NetCDF file from a previuos run. This should
 !--       be opened for extension, if its variables match the actual run.
-          INQUIRE( FILE='DATA_1D_PR_NETCDF', EXIST=netcdf_extend )
+          INQUIRE( FILE=filename, EXIST=netcdf_extend )
 
           IF ( netcdf_extend )  THEN
 !
 !--          Open an existing NetCDF file for output 
-             nc_stat = NF90_OPEN( 'DATA_1D_PR_NETCDF', NF90_WRITE, id_set_pr )
+             nc_stat = NF90_OPEN( filename, NF90_WRITE, id_set_pr )
              IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 29 )
 !
@@ -854,5 +862,5 @@
                 nc_stat = NF90_CLOSE( id_set_pr )
                 IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 30 )
-                CALL local_system( 'rm DATA_1D_PR_NETCDF' )
+                CALL local_system( 'rm ' // TRIM( filename ) )
              ENDIF
 
@@ -864,5 +872,5 @@
              IF ( netcdf_64bit )  THEN
 #if defined( __netcdf_64bit )
-                nc_stat = NF90_CREATE( 'DATA_1D_PR_NETCDF',                    &
+                nc_stat = NF90_CREATE( filename,                               &
                                        OR( NF90_NOCLOBBER, NF90_64BIT_OFFSET ),&
                                        id_set_pr )
@@ -870,10 +878,8 @@
                 IF ( myid == 0 )  PRINT*, '+++ WARNING: NetCDF: no 64-bit ', &
                                                'offset allowed on this machine'
-                nc_stat = NF90_CREATE( 'DATA_1D_PR_NETCDF', NF90_NOCLOBBER, &
-                                       id_set_pr )
+                nc_stat = NF90_CREATE( filename, NF90_NOCLOBBER, id_set_pr )
 #endif
              ELSE
-                nc_stat = NF90_CREATE( 'DATA_1D_PR_NETCDF', NF90_NOCLOBBER, &
-                                       id_set_pr )
+                nc_stat = NF90_CREATE( filename, NF90_NOCLOBBER, id_set_pr )
              ENDIF
              IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 31 )
@@ -886,12 +892,16 @@
        CASE ( 105 )
 !
+!--       Set filename
+          filename = 'DATA_1D_TS_NETCDF' // coupling_char
+
+!
 !--       Inquire, if there is a NetCDF file from a previuos run. This should
 !--       be opened for extension, if its variables match the actual run.
-          INQUIRE( FILE='DATA_1D_TS_NETCDF', EXIST=netcdf_extend )
+          INQUIRE( FILE=filename, EXIST=netcdf_extend )
 
           IF ( netcdf_extend )  THEN
 !
 !--          Open an existing NetCDF file for output 
-             nc_stat = NF90_OPEN( 'DATA_1D_TS_NETCDF', NF90_WRITE, id_set_ts )
+             nc_stat = NF90_OPEN( filename, NF90_WRITE, id_set_ts )
              IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 32 )
 !
@@ -906,5 +916,5 @@
                 nc_stat = NF90_CLOSE( id_set_ts )
                 IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 33 )
-                CALL local_system( 'rm DATA_1D_TS_NETCDF' )
+                CALL local_system( 'rm ' // TRIM( filename ) )
              ENDIF
 
@@ -916,5 +926,5 @@
              IF ( netcdf_64bit )  THEN
 #if defined( __netcdf_64bit )
-                nc_stat = NF90_CREATE( 'DATA_1D_TS_NETCDF',                    &
+                nc_stat = NF90_CREATE( filename,                               &
                                        OR( NF90_NOCLOBBER, NF90_64BIT_OFFSET ),&
                                        id_set_ts )
@@ -922,10 +932,8 @@
                 IF ( myid == 0 )  PRINT*, '+++ WARNING: NetCDF: no 64-bit ', &
                                                'offset allowed on this machine'
-                nc_stat = NF90_CREATE( 'DATA_1D_TS_NETCDF', NF90_NOCLOBBER, &
-                                       id_set_ts )
+                nc_stat = NF90_CREATE( filename, NF90_NOCLOBBER, id_set_ts )
 #endif
              ELSE
-                nc_stat = NF90_CREATE( 'DATA_1D_TS_NETCDF', NF90_NOCLOBBER, &
-                                       id_set_ts )
+                nc_stat = NF90_CREATE( filename, NF90_NOCLOBBER, id_set_ts )
              ENDIF
              IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 34 )
@@ -941,8 +949,8 @@
 !--       Set filename depending on unit number
           IF ( file_id == 106 )  THEN
-             filename = 'DATA_3D_NETCDF'
+             filename = 'DATA_3D_NETCDF' // coupling_char
              av = 0
           ELSE
-             filename = 'DATA_3D_AV_NETCDF'
+             filename = 'DATA_3D_AV_NETCDF' // coupling_char
              av = 1
           ENDIF
@@ -1000,12 +1008,16 @@
        CASE ( 107 )
 !
+!--       Set filename
+          filename = 'DATA_1D_SP_NETCDF' // coupling_char
+
+!
 !--       Inquire, if there is a NetCDF file from a previuos run. This should
 !--       be opened for extension, if its variables match the actual run.
-          INQUIRE( FILE='DATA_1D_SP_NETCDF', EXIST=netcdf_extend )
+          INQUIRE( FILE=filename, EXIST=netcdf_extend )
 
           IF ( netcdf_extend )  THEN
 !
 !--          Open an existing NetCDF file for output 
-             nc_stat = NF90_OPEN( 'DATA_1D_SP_NETCDF', NF90_WRITE, id_set_sp )
+             nc_stat = NF90_OPEN( filename, NF90_WRITE, id_set_sp )
              IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 38 )
 !
@@ -1020,5 +1032,5 @@
                 nc_stat = NF90_CLOSE( id_set_sp )
                 IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 39 )
-                CALL local_system( 'rm DATA_1D_SP_NETCDF' )
+                CALL local_system( 'rm ' // TRIM( filename ) )
              ENDIF
 
@@ -1030,5 +1042,5 @@
              IF ( netcdf_64bit )  THEN
 #if defined( __netcdf_64bit )
-                nc_stat = NF90_CREATE( 'DATA_1D_SP_NETCDF',                    &
+                nc_stat = NF90_CREATE( filename,                               &
                                        OR( NF90_NOCLOBBER, NF90_64BIT_OFFSET ),&
                                        id_set_sp )
@@ -1036,10 +1048,8 @@
                 IF ( myid == 0 )  PRINT*, '+++ WARNING: NetCDF: no 64-bit ', &
                                                'offset allowed on this machine'
-                nc_stat = NF90_CREATE( 'DATA_1D_SP_NETCDF', NF90_NOCLOBBER, &
-                                       id_set_sp )
+                nc_stat = NF90_CREATE( filename, NF90_NOCLOBBER, id_set_sp )
 #endif
              ELSE
-                nc_stat = NF90_CREATE( 'DATA_1D_SP_NETCDF', NF90_NOCLOBBER, &
-                                       id_set_sp )
+                nc_stat = NF90_CREATE( filename, NF90_NOCLOBBER, id_set_sp )
              ENDIF
              IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 40 )
@@ -1054,7 +1064,7 @@
 
           IF ( myid_char == '' )  THEN
-             filename = 'DATA_PRT_NETCDF'
-          ELSE
-             filename = 'DATA_PRT_NETCDF/' // myid_char
+             filename = 'DATA_PRT_NETCDF' // coupling_char
+          ELSE
+             filename = 'DATA_PRT_NETCDF' // coupling_char // '/' // myid_char
           ENDIF
 !
@@ -1091,5 +1101,6 @@
                 IF ( myid == 0  .AND. .NOT. openfile(file_id)%opened_before ) &
                 THEN    ! needs modification in case of non-extendable sets
-                   CALL local_system( 'mkdir  DATA_PRT_NETCDF/' )
+                   CALL local_system( 'mkdir  DATA_PRT_NETCDF' // &
+                                       coupling_char // '/' )
                 ENDIF
 #if defined( __parallel )
@@ -1126,12 +1137,16 @@
        CASE ( 109 )
 !
+!--       Set filename
+          filename = 'DATA_1D_PTS_NETCDF' // coupling_char
+
+!
 !--       Inquire, if there is a NetCDF file from a previuos run. This should
 !--       be opened for extension, if its variables match the actual run.
-          INQUIRE( FILE='DATA_1D_PTS_NETCDF', EXIST=netcdf_extend )
+          INQUIRE( FILE=filename, EXIST=netcdf_extend )
 
           IF ( netcdf_extend )  THEN
 !
 !--          Open an existing NetCDF file for output 
-             nc_stat = NF90_OPEN( 'DATA_1D_PTS_NETCDF', NF90_WRITE, id_set_pts )
+             nc_stat = NF90_OPEN( filename, NF90_WRITE, id_set_pts )
              IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 393 )
 !
@@ -1146,5 +1161,5 @@
                 nc_stat = NF90_CLOSE( id_set_pts )
                 IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 394 )
-                CALL local_system( 'rm DATA_1D_PTS_NETCDF' )
+                CALL local_system( 'rm ' // TRIM( filename ) )
              ENDIF
 
@@ -1156,5 +1171,5 @@
              IF ( netcdf_64bit )  THEN
 #if defined( __netcdf_64bit )
-                nc_stat = NF90_CREATE( 'DATA_1D_PTS_NETCDF',                   &
+                nc_stat = NF90_CREATE( filename,                               &
                                        OR( NF90_NOCLOBBER, NF90_64BIT_OFFSET ),&
                                        id_set_pts )
@@ -1162,10 +1177,8 @@
                 IF ( myid == 0 )  PRINT*, '+++ WARNING: NetCDF: no 64-bit ', &
                                                'offset allowed on this machine'
-                nc_stat = NF90_CREATE( 'DATA_1D_PTS_NETCDF', NF90_NOCLOBBER, &
-                                       id_set_pts )
+                nc_stat = NF90_CREATE( filename, NF90_NOCLOBBER, id_set_pts )
 #endif
              ELSE
-                nc_stat = NF90_CREATE( 'DATA_1D_PTS_NETCDF', NF90_NOCLOBBER, &
-                                       id_set_pts )
+                nc_stat = NF90_CREATE( filename, NF90_NOCLOBBER, id_set_pts )
              ENDIF
              IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 395 )
Index: /palm/trunk/SOURCE/check_parameters.f90
===================================================================
--- /palm/trunk/SOURCE/check_parameters.f90	(revision 101)
+++ /palm/trunk/SOURCE/check_parameters.f90	(revision 102)
@@ -4,5 +4,6 @@
 ! Actual revisions:
 ! -----------------
-! 
+! Check coupling_mode and set default (obligatory) values (like boundary
+! conditions for temperature and fluxes) in case of coupled runs
 !
 ! Former revisions:
@@ -69,4 +70,5 @@
     CHARACTER (LEN=8)   ::  date
     CHARACTER (LEN=10)  ::  time
+    CHARACTER (LEN=40)  ::  coupling_string
     CHARACTER (LEN=100) ::  action
 
@@ -85,4 +87,16 @@
           PRINT*, '    is set before running PALM'
        ENDIF
+    ENDIF
+
+!
+!-- Check the coupling mode
+    IF ( coupling_mode /= 'uncoupled'            .AND.  &
+         coupling_mode /= 'atmosphere_to_ocean'  .AND.  &
+         coupling_mode /= 'ocean_to_atmosphere' )  THEN
+       IF ( myid == 0 )  THEN
+          PRINT*, '+++ check_parameters:'
+          PRINT*, '    illegal coupling mode: ', TRIM( coupling_mode )
+       ENDIF
+       CALL local_stop
     ENDIF
 
@@ -93,9 +107,17 @@
     run_date = date(7:8)//'-'//date(5:6)//'-'//date(3:4)
     run_time = time(1:2)//':'//time(3:4)//':'//time(5:6)
-
-    WRITE ( run_description_header, '(A,2X,A,2X,A,A,A,I2.2,2X,A,A,2X,A,1X,A)' )&
-              TRIM( version ), TRIM( revision ), &
-              'run: ', TRIM( run_identifier ), '.', &
-              runnr, 'host: ', TRIM( host ), run_date, run_time
+    IF ( coupling_mode == 'uncoupled' )  THEN
+       coupling_string = ''
+    ELSEIF ( coupling_mode == 'atmosphere_to_ocean' )  THEN
+       coupling_string = ' coupled (atmosphere)'
+    ELSEIF ( coupling_mode == 'ocean_to_atmosphere' )  THEN
+       coupling_string = ' coupled (ocean)'
+    ENDIF       
+
+    WRITE ( run_description_header,                                        &
+                             '(A,2X,A,2X,A,A,A,I2.2,A,2X,A,A,2X,A,1X,A)' ) &
+              TRIM( version ), TRIM( revision ), 'run: ',                  &
+              TRIM( run_identifier ), '.', runnr, TRIM( coupling_string ), &
+              'host: ', TRIM( host ), run_date, run_time
 
 !
@@ -1029,15 +1051,20 @@
 !
 !-- Boundary conditions for potential temperature
-    IF ( bc_pt_b == 'dirichlet' )  THEN
-       ibc_pt_b = 0
-    ELSEIF ( bc_pt_b == 'neumann' )  THEN
-       ibc_pt_b = 1
+    IF ( coupling_mode == 'atmosphere_to_ocean' )  THEN
+       ibc_pt_b = 2
     ELSE
-       IF ( myid == 0 )  THEN
-          PRINT*, '+++ check_parameters:'
-          PRINT*, '    unknown boundary condition: bc_pt_b = ', bc_pt_b
-       ENDIF
-       CALL local_stop
-    ENDIF
+       IF ( bc_pt_b == 'dirichlet' )  THEN
+          ibc_pt_b = 0
+       ELSEIF ( bc_pt_b == 'neumann' )  THEN
+          ibc_pt_b = 1
+       ELSE
+          IF ( myid == 0 )  THEN
+             PRINT*, '+++ check_parameters:'
+             PRINT*, '    unknown boundary condition: bc_pt_b = ', bc_pt_b
+          ENDIF
+          CALL local_stop
+       ENDIF
+    ENDIF
+
     IF ( bc_pt_t == 'dirichlet' )  THEN
        ibc_pt_t = 0
Index: /palm/trunk/SOURCE/diffusion_u.f90
===================================================================
--- /palm/trunk/SOURCE/diffusion_u.f90	(revision 101)
+++ /palm/trunk/SOURCE/diffusion_u.f90	(revision 102)
@@ -4,5 +4,5 @@
 ! Actual revisions:
 ! -----------------
-! 
+! Momentumflux at top (uswst) included as boundary condition
 ! 
 ! Former revisions:
@@ -54,5 +54,6 @@
 ! Call for all grid points
 !------------------------------------------------------------------------------!
-    SUBROUTINE diffusion_u( ddzu, ddzw, km, km_damp_y, tend, u, usws, v, w )
+    SUBROUTINE diffusion_u( ddzu, ddzw, km, km_damp_y, tend, u, usws, uswst, &
+                            v, w )
 
        USE control_parameters
@@ -66,5 +67,5 @@
        REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt+1), km_damp_y(nys-1:nyn+1)
        REAL    ::  tend(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1)
-       REAL, DIMENSION(:,:),   POINTER ::  usws
+       REAL, DIMENSION(:,:),   POINTER ::  usws, uswst
        REAL, DIMENSION(:,:,:), POINTER ::  km, u, v, w
        REAL, DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::  usvs
@@ -155,5 +156,5 @@
 !--          Compute vertical diffusion. In case of simulating a Prandtl layer,
 !--          index k starts at nzb_u_inner+2.
-             DO  k = nzb_diff_u(j,i), nzt
+             DO  k = nzb_diff_u(j,i), nzt_diff
 !
 !--             Interpolate eddy diffusivities on staggered gridpoints
@@ -195,6 +196,26 @@
                       & + ( kmzp * ( w(k,j,i)   - w(k,j,i-1)   ) * ddx       &
                       &   ) * ddzw(k)                                        &
-                      & + ( kmzp * ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1)   &
+                      & + ( kmzp * ( u(k+1,j,i) - u(k,j,i)     ) * ddzu(k+1) &
                       &   + usws(j,i)                                        &
+                      &   ) * ddzw(k)
+             ENDIF
+
+!
+!--          Vertical diffusion at the first gridpoint below the top boundary,
+!--          if the momentum flux at the top is prescribed by the user
+             IF ( use_top_fluxes )  THEN
+                k = nzt
+!
+!--             Interpolate eddy diffusivities on staggered gridpoints
+                kmzp = 0.25 * &
+                       ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
+                kmzm = 0.25 * &
+                       ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
+
+                tend(k,j,i) = tend(k,j,i)                                    &
+                      & - ( kmzm * ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx       &
+                      &   ) * ddzw(k)                                        &
+                      & + ( -uswst(j,i)                                      &
+                      &   - kmzm * ( u(k,j,i)   - u(k-1,j,i)   ) * ddzu(k)   &
                       &   ) * ddzw(k)
              ENDIF
@@ -210,5 +231,5 @@
 !------------------------------------------------------------------------------!
     SUBROUTINE diffusion_u_ij( i, j, ddzu, ddzw, km, km_damp_y, tend, u, usws, &
-                               v, w )
+                               uswst, v, w )
 
        USE control_parameters
@@ -223,5 +244,5 @@
        REAL    ::  tend(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1)
        REAL, DIMENSION(nzb:nzt+1)      ::  usvs
-       REAL, DIMENSION(:,:),   POINTER ::  usws
+       REAL, DIMENSION(:,:),   POINTER ::  usws, uswst
        REAL, DIMENSION(:,:,:), POINTER ::  km, u, v, w
 
@@ -303,5 +324,5 @@
 !--    Compute vertical diffusion. In case of simulating a Prandtl layer,
 !--    index k starts at nzb_u_inner+2.
-       DO  k = nzb_diff_u(j,i), nzt
+       DO  k = nzb_diff_u(j,i), nzt_diff
 !
 !--       Interpolate eddy diffusivities on staggered gridpoints
@@ -339,9 +360,29 @@
                       & + ( kmzp * ( w(k,j,i)   - w(k,j,i-1)   ) * ddx       &
                       &   ) * ddzw(k)                                        &
-                      & + ( kmzp * ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1)   &
+                      & + ( kmzp * ( u(k+1,j,i) - u(k,j,i)     ) * ddzu(k+1) &
                       &   + usws(j,i)                                        &
                       &   ) * ddzw(k)
        ENDIF
 
+!
+!--    Vertical diffusion at the first gridpoint below the top boundary,
+!--    if the momentum flux at the top is prescribed by the user
+       IF ( use_top_fluxes )  THEN
+          k = nzt
+!
+!--       Interpolate eddy diffusivities on staggered gridpoints
+          kmzp = 0.25 * &
+                 ( km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) )
+          kmzm = 0.25 * &
+                 ( km(k,j,i)+km(k-1,j,i)+km(k,j,i-1)+km(k-1,j,i-1) )
+
+          tend(k,j,i) = tend(k,j,i)                                          &
+                      & - ( kmzm * ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx       &
+                      &   ) * ddzw(k)                                        &
+                      & + ( -uswst(j,i)                                      &
+                      &   - kmzm * ( u(k,j,i)   - u(k-1,j,i)   ) * ddzu(k)   &
+                      &   ) * ddzw(k)
+       ENDIF
+
     END SUBROUTINE diffusion_u_ij
 
Index: /palm/trunk/SOURCE/diffusion_v.f90
===================================================================
--- /palm/trunk/SOURCE/diffusion_v.f90	(revision 101)
+++ /palm/trunk/SOURCE/diffusion_v.f90	(revision 102)
@@ -4,5 +4,5 @@
 ! Actual revisions:
 ! -----------------
-! 
+! Momentumflux at top (vswst) included as boundary condition
 !
 ! Former revisions:
@@ -52,5 +52,6 @@
 ! Call for all grid points
 !------------------------------------------------------------------------------!
-    SUBROUTINE diffusion_v( ddzu, ddzw, km, km_damp_x, tend, u, v, vsws, w )
+    SUBROUTINE diffusion_v( ddzu, ddzw, km, km_damp_x, tend, u, v, vsws, &
+                            vswst, w )
 
        USE control_parameters
@@ -64,5 +65,5 @@
        REAL    ::  ddzu(1:nzt+1), ddzw(1:nzt+1), km_damp_x(nxl-1:nxr+1)
        REAL    ::  tend(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1)
-       REAL, DIMENSION(:,:),   POINTER ::  vsws
+       REAL, DIMENSION(:,:),   POINTER ::  vsws, vswst
        REAL, DIMENSION(:,:,:), POINTER ::  km, u, v, w
        REAL, DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::  vsus
@@ -153,5 +154,5 @@
 !--          Compute vertical diffusion. In case of simulating a Prandtl
 !--          layer, index k starts at nzb_v_inner+2.
-             DO  k = nzb_diff_v(j,i), nzt
+             DO  k = nzb_diff_v(j,i), nzt_diff
 !
 !--             Interpolate eddy diffusivities on staggered gridpoints
@@ -193,6 +194,26 @@
                       & + ( kmzp * ( w(k,j,i) - w(k,j-1,i)     ) * ddy       &
                       &   ) * ddzw(k)                                        &
-                      & + ( kmzp * ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1)     &
+                      & + ( kmzp * ( v(k+1,j,i) - v(k,j,i)     ) * ddzu(k+1) &
                       &   + vsws(j,i)                                        &
+                      &   ) * ddzw(k)
+             ENDIF
+
+!
+!--          Vertical diffusion at the first gridpoint below the top boundary,
+!--          if the momentum flux at the top is prescribed by the user
+             IF ( use_top_fluxes )  THEN
+                k = nzt
+!
+!--             Interpolate eddy diffusivities on staggered gridpoints
+                kmzp = 0.25 * &
+                       ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
+                kmzm = 0.25 * &
+                       ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
+
+                tend(k,j,i) = tend(k,j,i)                                    &
+                      & - ( kmzm *  ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy      &
+                      &   ) * ddzw(k)                                        &
+                      & + ( -vswst(j,i)                                      &
+                      &   - kmzm * ( v(k,j,i)   - v(k-1,j,i)    ) * ddzu(k)  &
                       &   ) * ddzw(k)
              ENDIF
@@ -208,5 +229,5 @@
 !------------------------------------------------------------------------------!
     SUBROUTINE diffusion_v_ij( i, j, ddzu, ddzw, km, km_damp_x, tend, u, v, &
-                               vsws, w )
+                               vsws, vswst, w )
 
        USE control_parameters
@@ -221,5 +242,5 @@
        REAL    ::  tend(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1)
        REAL, DIMENSION(nzb:nzt+1)      ::  vsus
-       REAL, DIMENSION(:,:),   POINTER ::  vsws
+       REAL, DIMENSION(:,:),   POINTER ::  vsws, vswst
        REAL, DIMENSION(:,:,:), POINTER ::  km, u, v, w
 
@@ -302,5 +323,5 @@
 !--    Compute vertical diffusion. In case of simulating a Prandtl layer,
 !--    index k starts at nzb_v_inner+2.
-       DO  k = nzb_diff_v(j,i), nzt
+       DO  k = nzb_diff_v(j,i), nzt_diff
 !
 !--       Interpolate eddy diffusivities on staggered gridpoints
@@ -338,9 +359,29 @@
                       & + ( kmzp * ( w(k,j,i) - w(k,j-1,i)     ) * ddy       &
                       &   ) * ddzw(k)                                        &
-                      & + ( kmzp * ( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1)     &
+                      & + ( kmzp * ( v(k+1,j,i) - v(k,j,i)     ) * ddzu(k+1) &
                       &   + vsws(j,i)                                        &
                       &   ) * ddzw(k)
        ENDIF
 
+!
+!--    Vertical diffusion at the first gridpoint below the top boundary,
+!--    if the momentum flux at the top is prescribed by the user
+       IF ( use_top_fluxes )  THEN
+          k = nzt
+!
+!--       Interpolate eddy diffusivities on staggered gridpoints
+          kmzp = 0.25 * &
+                 ( km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) )
+          kmzm = 0.25 * &
+                 ( km(k,j,i)+km(k-1,j,i)+km(k,j-1,i)+km(k-1,j-1,i) )
+
+          tend(k,j,i) = tend(k,j,i)                                          &
+                      & - ( kmzm *  ( w(k-1,j,i) - w(k-1,j-1,i) ) * ddy      &
+                      &   ) * ddzw(k)                                        &
+                      & + ( -vswst(j,i)                                      &
+                      &   - kmzm * ( v(k,j,i)   - v(k-1,j,i)    ) * ddzu(k)  &
+                      &   ) * ddzw(k)
+       ENDIF
+
     END SUBROUTINE diffusion_v_ij
 
Index: /palm/trunk/SOURCE/flow_statistics.f90
===================================================================
--- /palm/trunk/SOURCE/flow_statistics.f90	(revision 101)
+++ /palm/trunk/SOURCE/flow_statistics.f90	(revision 102)
@@ -4,5 +4,5 @@
 ! Actual revisions:
 ! -----------------
-! 
+! Prescribed momentum fluxes at the top surface are used
 !
 ! Former revisions:
@@ -381,5 +381,5 @@
 !--          ----  strictly speaking the following k-loop would have to be
 !--                split up according to the staggered grid.
-             DO  k = nzb_diff_s_outer(j,i)-1, nzt
+             DO  k = nzb_diff_s_outer(j,i)-1, nzt_diff
 !
 !--             Momentum flux w"u"
@@ -398,7 +398,4 @@
                                  + ( w(k,j,i)   - w(k,j-1,i) ) * ddy           &
                                                                ) * rmask(j,i,sr)
-             ENDDO
-
-             DO  k = nzb_diff_s_outer(j,i)-1, nzt_diff
 !
 !--             Heat flux w"pt"
@@ -490,4 +487,8 @@
 !--          Subgridscale fluxes at the top surface
              IF ( use_top_fluxes )  THEN
+                sums_l(nzt,12,tn) = sums_l(nzt,12,tn) + &
+                                    uswst(j,i) * rmask(j,i,sr)    ! w"u"
+                sums_l(nzt,14,tn) = sums_l(nzt,14,tn) + &
+                                    vswst(j,i) * rmask(j,i,sr)    ! w"v"
                 sums_l(nzt,16,tn) = sums_l(nzt,16,tn) + &
                                     tswst(j,i)  * rmask(j,i,sr)   ! w"pt"
Index: /palm/trunk/SOURCE/header.f90
===================================================================
--- /palm/trunk/SOURCE/header.f90	(revision 101)
+++ /palm/trunk/SOURCE/header.f90	(revision 102)
@@ -4,5 +4,5 @@
 ! Actual revisions:
 ! -----------------
-! 
+! Output of informations for coupled model runs (boundary conditions etc.)
 !
 ! Former revisions:
@@ -117,6 +117,8 @@
     host_chr = host(1:10)
     ver_rev = TRIM( version ) // '  ' // TRIM( revision )
-    WRITE ( io, 100 )  ver_rev, TRIM( run_classification ), run_date, &
-                       run_identifier, run_time, runnr, ADJUSTR( host_chr )
+    WRITE ( io, 100 )  ver_rev, TRIM( run_classification )
+    IF ( coupling_mode /= 'uncoupled' )  WRITE ( io, 101 )  coupling_mode
+    WRITE ( io, 102 )  run_date, run_identifier, run_time, runnr, &
+                       ADJUSTR( host_chr )
 #if defined( __parallel )
     IF ( npex == -1  .AND.  pdims(2) /= 1 )  THEN
@@ -130,7 +132,7 @@
     ENDIF
     IF ( threads_per_task == 1 )  THEN
-       WRITE ( io, 101 )  numprocs, pdims(1), pdims(2), TRIM( char1 )
+       WRITE ( io, 103 )  numprocs, pdims(1), pdims(2), TRIM( char1 )
     ELSE
-       WRITE ( io, 102 )  numprocs*threads_per_task, numprocs, &
+       WRITE ( io, 104 )  numprocs*threads_per_task, numprocs, &
                           threads_per_task, pdims(1), pdims(2), TRIM( char1 )
     ENDIF
@@ -139,11 +141,11 @@
          npex == -1  .AND.  pdims(2) == 1 )                      &
     THEN
-       WRITE ( io, 104 )
+       WRITE ( io, 106 )
     ELSEIF ( pdims(2) == 1 )  THEN
-       WRITE ( io, 105 )  'x'
+       WRITE ( io, 107 )  'x'
     ELSEIF ( pdims(1) == 1 )  THEN
-       WRITE ( io, 105 )  'y'
-    ENDIF
-    IF ( use_seperate_pe_for_dvrp_output )  WRITE ( io, 103 )
+       WRITE ( io, 107 )  'y'
+    ENDIF
+    IF ( use_seperate_pe_for_dvrp_output )  WRITE ( io, 105 )
 #endif
     WRITE ( io, 99 )
@@ -378,6 +380,8 @@
     IF ( ibc_pt_b == 0 )  THEN
        runten = TRIM( runten ) // ' pt(0)   = pt_surface'
-    ELSE
+    ELSEIF ( ibc_pt_b == 1 )  THEN
        runten = TRIM( runten ) // ' pt(0)   = pt(1)'
+    ELSEIF ( ibc_pt_b == 2 )  THEN
+       runten = TRIM( runten ) // ' pt(0) = from coupled model'
     ENDIF
     IF ( ibc_pt_t == 0 )  THEN
@@ -457,6 +461,10 @@
     IF ( use_top_fluxes )  THEN
        WRITE ( io, 304 )
-       IF ( constant_top_heatflux )  THEN
-          WRITE ( io, 306 )  top_heatflux
+       IF ( coupling_mode == 'uncoupled' )  THEN
+          IF ( constant_top_heatflux )  THEN
+             WRITE ( io, 306 )  top_heatflux
+          ENDIF
+       ELSEIF ( coupling_mode == 'ocean_to_atmosphere' )  THEN
+          WRITE ( io, 316 )
        ENDIF
        IF ( ocean  .AND.  constant_top_salinityflux )  THEN
@@ -1193,17 +1201,19 @@
 100 FORMAT (/1X,'*************************',11X,42('-')/        &
             1X,'* ',A,' *',11X,A/                               &
-            1X,'*************************',11X,42('-')//        &
-            ' Date:            ',A8,11X,'Run:       ',A20/      &
+            1X,'*************************',11X,42('-'))
+101 FORMAT (37X,'coupled run: ',A/ &
+            37X,42('-'))
+102 FORMAT (/' Date:            ',A8,11X,'Run:       ',A20/      &
             ' Time:            ',A8,11X,'Run-No.:   ',I2.2/     &
             ' Run on host:   ',A10)
 #if defined( __parallel )
-101 FORMAT (' Number of PEs:',7X,I4,11X,'Processor grid (x,y): (',I3,',',I3, &
+103 FORMAT (' Number of PEs:',7X,I4,11X,'Processor grid (x,y): (',I3,',',I3, &
               ')',1X,A)
-102 FORMAT (' Number of PEs:',7X,I4,11X,'Tasks:',I4,'   threads per task:',I4/ &
+104 FORMAT (' Number of PEs:',7X,I4,11X,'Tasks:',I4,'   threads per task:',I4/ &
               37X,'Processor grid (x,y): (',I3,',',I3,')',1X,A)
-103 FORMAT (37X,'One additional PE is used to handle'/37X,'the dvrp output!')
-104 FORMAT (37X,'A 1d-decomposition along x is forced'/ &
+105 FORMAT (37X,'One additional PE is used to handle'/37X,'the dvrp output!')
+106 FORMAT (37X,'A 1d-decomposition along x is forced'/ &
             37X,'because the job is running on an SMP-cluster')
-105 FORMAT (37X,'A 1d-decomposition along ',A,' is used')
+107 FORMAT (37X,'A 1d-decomposition along ',A,' is used')
 #endif
 110 FORMAT (/' Numerical Schemes:'/ &
@@ -1332,4 +1342,6 @@
 314 FORMAT ('       Predefined scalar value at the surface')
 315 FORMAT ('       Humidity / scalar flux at top surface is 0.0')
+316 FORMAT ('       Sensible heatflux and momentum flux from coupled ', &
+                    'atmosphere model')
 317 FORMAT (//' Lateral boundaries:'/ &
             '       left/right:  ',A/    &
Index: /palm/trunk/SOURCE/init_3d_model.f90
===================================================================
--- /palm/trunk/SOURCE/init_3d_model.f90	(revision 101)
+++ /palm/trunk/SOURCE/init_3d_model.f90	(revision 102)
@@ -7,5 +7,5 @@
 ! Actual revisions:
 ! -----------------
-! 
+! Flux initialization in case of coupled runs, +momentum fluxes at top boundary
 !
 ! Former revisions:
@@ -102,5 +102,7 @@
               ts(nys-1:nyn+1,nxl-1:nxr+1), tswst_1(nys-1:nyn+1,nxl-1:nxr+1),  &
               us(nys-1:nyn+1,nxl-1:nxr+1), usws_1(nys-1:nyn+1,nxl-1:nxr+1),   &
-              vsws_1(nys-1:nyn+1,nxl-1:nxr+1), z0(nys-1:nyn+1,nxl-1:nxr+1) )
+              uswst_1(nys-1:nyn+1,nxl-1:nxr+1),                               &
+              vsws_1(nys-1:nyn+1,nxl-1:nxr+1),                                &
+              vswst_1(nys-1:nyn+1,nxl-1:nxr+1), z0(nys-1:nyn+1,nxl-1:nxr+1) )
 
     IF ( timestep_scheme(1:5) /= 'runge' )  THEN
@@ -111,4 +113,6 @@
                  tswst_2(nys-1:nyn+1,nxl-1:nxr+1), &
                  usws_2(nys-1:nyn+1,nxl-1:nxr+1),  &
+                 uswst_2(nys-1:nyn+1,nxl-1:nxr+1), &
+                 vswst_2(nys-1:nyn+1,nxl-1:nxr+1), &
                  vsws_2(nys-1:nyn+1,nxl-1:nxr+1) )
     ENDIF
@@ -260,5 +264,7 @@
        tswst_m => tswst_1;  tswst => tswst_2
        usws_m  => usws_1;   usws  => usws_2
+       uswst_m => uswst_1;  uswst => uswst_2
        vsws_m  => vsws_1;   vsws  => vsws_2
+       vswst_m => vswst_1;  vswst => vswst_2
        e_m  => e_1;   e  => e_2;   e_p  => e_3;   te_m  => e_3
        kh_m => kh_1;  kh => kh_2
@@ -287,5 +293,7 @@
        tswst => tswst_1
        usws  => usws_1
+       uswst => uswst_1
        vsws  => vsws_1
+       vswst => vswst_1
        e     => e_1;   e_p  => e_2;   te_m  => e_3;   e_m  => e_3
        kh    => kh_1
@@ -381,4 +389,6 @@
              vsws = 0.0
           ENDIF
+          uswst = top_momentumflux_u
+          vswst = top_momentumflux_v
 
 !
@@ -494,10 +504,12 @@
                            ! e = (u*/cm)**2 at k=nzb+1
           ENDIF
-          e    = 0.0
-          rif  = 0.0
-          ts   = 0.0
-          us   = 0.0
-          usws = 0.0
-          vsws = 0.0
+          e     = 0.0
+          rif   = 0.0
+          ts    = 0.0
+          us    = 0.0
+          usws  = 0.0
+          uswst = top_momentumflux_u
+          vsws  = 0.0
+          vswst = top_momentumflux_v
           IF ( humidity  .OR.  passive_scalar ) qs = 0.0
 
@@ -617,5 +629,12 @@
                 saswst = top_salinityflux
              ENDIF
-         ENDIF
+          ENDIF
+
+!
+!--       Initialization in case of a coupled model run
+          IF ( coupling_mode == 'ocean_to_atmosphere' )  THEN
+             tswst = 0.0
+             IF ( ASSOCIATED( tswst_m ) )  tswst_m = tswst
+          ENDIF
 
        ENDIF
Index: /palm/trunk/SOURCE/init_pegrid.f90
===================================================================
--- /palm/trunk/SOURCE/init_pegrid.f90	(revision 101)
+++ /palm/trunk/SOURCE/init_pegrid.f90	(revision 102)
@@ -4,5 +4,6 @@
 ! Actual revisions:
 ! -----------------
-! 
+! Intercommunicator (comm_inter) and derived data type (type_xy) for
+! coupled model runs created
 !
 ! Former revisions:
@@ -480,4 +481,60 @@
        CALL MPI_SEND( ibuf, 12, MPI_INTEGER, 0, myid, comm2d, ierr )       
     ENDIF
+#endif
+
+#if defined( __mpi2 )
+!
+!-- In case of coupled runs, get the port name on PE0 of the atmosphere model
+!-- and pass it to PE0 of the ocean model
+    IF ( myid == 0 )  THEN
+
+       IF ( coupling_mode == 'atmosphere_to_ocean' )  THEN
+
+          CALL MPI_OPEN_PORT( MPI_INFO_NULL, port_name, ierr )
+          CALL MPI_PUBLISH_NAME( 'palm_coupler', MPI_INFO_NULL, port_name, &
+                                 ierr )
+
+       ELSEIF ( coupling_mode == 'ocean_to_atmosphere' )  THEN
+
+          CALL MPI_LOOKUP_NAME( 'palm_coupler', MPI_INFO_NULL, port_name, ierr )
+
+       ENDIF
+
+    ENDIF
+
+!
+!-- In case of coupled runs, establish the connection between the atmosphere
+!-- and the ocean model and define the intercommunicator (comm_inter)
+    CALL MPI_BARRIER( comm2d, ierr )
+    IF ( coupling_mode == 'atmosphere_to_ocean' )  THEN
+
+       print*, '... before COMM_ACCEPT'
+       CALL MPI_COMM_ACCEPT( port_name, MPI_INFO_NULL, 0, MPI_COMM_WORLD, &
+                             comm_inter, ierr )
+       print*, '--- ierr = ', ierr
+       print*, '--- comm_inter atmosphere = ', comm_inter
+
+    ELSEIF ( coupling_mode == 'ocean_to_atmosphere' )  THEN
+
+       IF ( myid == 0 )  PRINT*, '*** read: ', port_name, '  ierr = ', ierr
+       print*, '... before COMM_CONNECT'
+       CALL MPI_COMM_CONNECT( port_name, MPI_INFO_NULL, 0, MPI_COMM_WORLD, &
+                              comm_inter, ierr )
+       print*, '--- ierr = ', ierr
+       print*, '--- comm_inter ocean      = ', comm_inter
+
+    ENDIF
+
+!
+!-- In case of coupled runs, create a new MPI derived datatype for the
+!-- exchange of surface (xy) data .
+!-- Gridpoint number for the exchange of ghost points (xy-plane)
+    ngp_xy  = ( nxr - nxl + 3 ) * ( nyn - nys + 3 )
+
+!
+!-- Define a new MPI derived datatype for the exchange of ghost points in
+!-- y-direction for 2D-arrays (line)
+    CALL MPI_TYPE_VECTOR( ngp_xy, 1, nzt-nzb+2, MPI_REAL, type_xy, ierr )
+    CALL MPI_TYPE_COMMIT( type_xy, ierr )
 #endif
 
Index: /palm/trunk/SOURCE/modules.f90
===================================================================
--- /palm/trunk/SOURCE/modules.f90	(revision 101)
+++ /palm/trunk/SOURCE/modules.f90	(revision 102)
@@ -5,5 +5,6 @@
 ! Actual revisions:
 ! -----------------
-! 
+! +comm_inter, coupling_char, coupling_mode, dt_coupling, ngp_xy, port_name,
+! time_coupling, type_xy, uswst*, vswst*
 !
 ! Former revisions:
@@ -107,9 +108,11 @@
     REAL, DIMENSION(:,:), ALLOCATABLE, TARGET ::                               &
           qsws_1, qsws_2, qswst_1, qswst_2, rif_1, rif_2, saswsb_1, saswst_1,  &
-          shf_1, shf_2, tswst_1, tswst_2, usws_1, usws_2, vsws_1, vsws_2
+          shf_1, shf_2, tswst_1, tswst_2, usws_1, usws_2, uswst_1, uswst_2,    &
+          vsws_1, vsws_2, vswst_1, vswst_2
 
     REAL, DIMENSION(:,:), POINTER ::                                           &
           qsws, qsws_m, qswst, qswst_m, rif, rif_m, saswsb, saswst, shf,       &
-          shf_m, tswst, tswst_m, usws, usws_m, vsws, vsws_m
+          shf_m, tswst, tswst_m, usws, uswst, usws_m, uswst_m, vsws, vswst,    &
+          vsws_m, vswst_m
 
     REAL, DIMENSION(:,:,:), ALLOCATABLE ::                                     &
@@ -229,4 +232,5 @@
 
     CHARACTER (LEN=1)   ::  cycle_mg = 'w', timestep_reason = ' '
+    CHARACTER (LEN=2)   ::  coupling_char = ''
     CHARACTER (LEN=5)   ::  write_binary = 'false'
     CHARACTER (LEN=6)   ::  grid_matching = 'match'
@@ -246,4 +250,5 @@
                             bc_sa_t = 'neumann', &
                             bc_uv_b = 'dirichlet', bc_uv_t = 'dirichlet', &
+                            coupling_mode = 'uncoupled', &
                             dissipation_1d = 'as_in_3d_model', &
                             fft_method = 'system-specific', &
@@ -348,5 +353,6 @@
              disturbance_level_t = -9999999.9, &
              dt = -1.0, dt_averaging_input = 0.0, &
-             dt_averaging_input_pr = 9999999.9, dt_data_output = 9999999.9, &
+             dt_averaging_input_pr = 9999999.9, dt_coupling = 9999999.9, &
+             dt_data_output = 9999999.9, &
              dt_data_output_av = 9999999.9, dt_disturb = 9999999.9, &
              dt_dopr = 9999999.9, dt_dopr_listing = 9999999.9, &
@@ -382,5 +388,6 @@
              surface_scalarflux = 0.0, surface_waterflux = 0.0, &
              s_surface = 0.0, s_surface_initial_change = 0.0, &
-             termination_time_needed = -1.0, time_disturb = 0.0, &
+             termination_time_needed = -1.0, time_coupling = 0.0, &
+             time_disturb = 0.0, &
              time_dopr = 0.0, time_dopr_av = 0.0, time_dopr_listing = 0.0, &
              time_dopts = 0.0, time_dosp = 0.0, time_dosp_av = 0.0, &
@@ -389,5 +396,6 @@
              time_do_sla = 0.0, time_dvrp = 0.0, time_prel = 0.0, &
              time_restart = 9999999.9, time_run_control = 0.0, &
-             top_heatflux = 9999999.9, top_salinityflux = 9999999.9, &
+             top_heatflux = 9999999.9, top_momentumflux_u = 0.0, &
+             top_momentumflux_v = 0.0, top_salinityflux = 9999999.9, &
              ug_surface = 0.0, u_gtrans = 0.0, &
              ups_limit_e = 0.0, ups_limit_pt = 0.0, ups_limit_u = 0.0, &
@@ -916,10 +924,15 @@
 
 #if defined( __parallel )
-    INTEGER ::  comm1dx, comm1dy, comm2d, comm_palm, ierr, myidx, myidy,       &
-                ndim = 2, ngp_y, pleft, pnorth, pright, psouth,                &
+#if defined( __mpi2 )
+    CHARACTER (LEN=MPI_MAX_PORT_NAME) ::  port_name
+#endif
+
+    INTEGER ::  comm1dx, comm1dy, comm2d, comm_inter, comm_palm, ierr, myidx,  &
+                myidy, ndim = 2, ngp_xy, ngp_y, pleft, pnorth, pright, psouth, &
                 sendrecvcount_xy, sendrecvcount_yz, sendrecvcount_zx,          &
                 sendrecvcount_zyd, sendrecvcount_yxd,                          &
-                type_x, type_x_int, ibuf(12), pcoord(2), pdims(2),             &
-                status(MPI_STATUS_SIZE)
+                type_x, type_x_int, type_xy
+
+    INTEGER ::  ibuf(12), pcoord(2), pdims(2), status(MPI_STATUS_SIZE)
 
     INTEGER, DIMENSION(:), ALLOCATABLE ::  ngp_yz, type_xz
Index: /palm/trunk/SOURCE/palm.f90
===================================================================
--- /palm/trunk/SOURCE/palm.f90	(revision 101)
+++ /palm/trunk/SOURCE/palm.f90	(revision 102)
@@ -4,5 +4,5 @@
 ! Actual revisions:
 ! -----------------
-! 
+! Get coupling mode from environment variable
 !
 ! Former revisions:
@@ -43,4 +43,5 @@
     USE arrays_3d
     USE constants
+    USE control_parameters
     USE cpulog
     USE dvrp_variables
@@ -53,5 +54,4 @@
     USE spectrum
     USE statistics
-    USE control_parameters
 
     IMPLICIT NONE
@@ -63,5 +63,5 @@
     INTEGER           ::  i, run_description_header_i(80)
 
-    version = 'PALM 3.3'
+    version = 'PALM 3.3a'
 
 #if defined( __parallel )
@@ -75,4 +75,13 @@
 #endif
 
+#if defined( __mpi2 )
+!
+!-- Get information about the coupling mode from the environment variable
+!-- which has been set by the mpiexec command
+    CALL local_getenv( 'coupling_mode', 13, coupling_mode, i )
+    IF ( i == 0 )  coupling_mode = 'uncoupled'
+    IF ( coupling_mode == 'ocean_to_atmosphere' )  coupling_char = '_O'
+#endif
+
 !
 !-- Initialize measuring of the CPU-time remaining to the run
@@ -104,5 +113,13 @@
 !
 !-- Open a file for debug output
-    OPEN( 9, FILE='DEBUG'//myid_char, FORM='FORMATTED' )
+    OPEN( 9, FILE='DEBUG'//TRIM( coupling_char )//myid_char, FORM='FORMATTED' )
+
+#if defined( __mpi2 )
+!
+!-- TEST OUTPUT (TO BE REMOVED)
+    WRITE(9,*) '*** coupling_mode = "', TRIM( coupling_mode ), '"'
+    CALL LOCAL_FLUSH( 9 )
+    print*, '*** PE', myid, '  ', TRIM( coupling_mode )
+#endif
 
 !
@@ -167,4 +184,16 @@
     CALL cpu_log( log_point(4), 'last actions', 'stop' )
 
+#if defined( __mpi2 )
+!
+!-- Test exchange via intercommunicator
+    IF ( coupling_mode == 'atmosphere_to_ocean' )  THEN
+       i = 12345 + myid
+       CALL MPI_SEND( i, 1, MPI_INTEGER, myid, 11, comm_inter, ierr )
+    ELSEIF ( coupling_mode == 'ocean_to_atmosphere' )  THEN
+       CALL MPI_RECV( i, 1, MPI_INTEGER, myid, 11, comm_inter, status, ierr )
+       PRINT*, '### myid: ', myid, '   received from atmosphere:  i = ', i
+    ENDIF
+#endif
+
 !
 !-- Take final CPU-time for CPU-time analysis
Index: /palm/trunk/SOURCE/parin.f90
===================================================================
--- /palm/trunk/SOURCE/parin.f90	(revision 101)
+++ /palm/trunk/SOURCE/parin.f90	(revision 102)
@@ -4,5 +4,5 @@
 ! Actual revisions:
 ! -----------------
-! 
+! +top_momentumflux_u|v in inipar, +dt_coupling in d3par
 !
 ! Former revisions:
@@ -91,5 +91,6 @@
              surface_waterflux, s_surface, s_surface_initial_change, &
              s_vertical_gradient, s_vertical_gradient_level, top_heatflux, &
-             top_salinityflux, timestep_scheme, topography, ug_surface, &
+             top_momentumflux_u, top_momentumflux_v, top_salinityflux, &
+             timestep_scheme, topography, ug_surface, &
              ug_vertical_gradient, ug_vertical_gradient_level, ups_limit_e, &
              ups_limit_pt, ups_limit_u, ups_limit_v, ups_limit_w, &
@@ -109,5 +110,5 @@
                        disturbance_level_t, do2d_at_begin, do3d_at_begin, &
                        do3d_compress, do3d_comp_prec, dt, dt_averaging_input, &
-                       dt_averaging_input_pr, dt_data_output, &
+                       dt_averaging_input_pr, dt_coupling, dt_data_output, &
                        dt_data_output_av, dt_disturb, dt_dopr, &
                        dt_dopr_listing, dt_dots, dt_do2d_xy, dt_do2d_xz, &
Index: /palm/trunk/SOURCE/prognostic_equations.f90
===================================================================
--- /palm/trunk/SOURCE/prognostic_equations.f90	(revision 101)
+++ /palm/trunk/SOURCE/prognostic_equations.f90	(revision 102)
@@ -4,5 +4,5 @@
 ! Actual revisions:
 ! -----------------
-! 
+! +uswst, vswst as arguments in calls of diffusion_u|v
 !
 ! Former revisions:
@@ -148,8 +148,8 @@
           IF ( tsc(2) == 2.0  .AND.  timestep_scheme(1:8) == 'leapfrog' )  THEN
              CALL diffusion_u( i, j, ddzu, ddzw, km_m, km_damp_y, tend, u_m,   &
-                               usws_m, v_m, w_m )
+                               usws_m, uswst_m, v_m, w_m )
           ELSE
              CALL diffusion_u( i, j, ddzu, ddzw, km, km_damp_y, tend, u, usws, &
-                               v, w )
+                               uswst, v, w )
           ENDIF
           CALL coriolis( i, j, 1 )
@@ -216,8 +216,8 @@
           IF ( tsc(2) == 2.0  .AND.  timestep_scheme(1:8) == 'leapfrog' )  THEN
              CALL diffusion_v( i, j, ddzu, ddzw, km_m, km_damp_x, tend, u_m, &
-                               v_m, vsws_m, w_m )
+                               v_m, vsws_m, vswst_m, w_m )
           ELSE
              CALL diffusion_v( i, j, ddzu, ddzw, km, km_damp_x, tend, u, v,  &
-                               vsws, w )
+                               vsws, vswst, w )
           ENDIF
           CALL coriolis( i, j, 2 )
@@ -817,8 +817,8 @@
              THEN
                 CALL diffusion_u( i, j, ddzu, ddzw, km_m, km_damp_y, tend,  &
-                                  u_m, usws_m, v_m, w_m )
+                                  u_m, usws_m, uswst_m, v_m, w_m )
              ELSE
                 CALL diffusion_u( i, j, ddzu, ddzw, km, km_damp_y, tend, u, &
-                                  usws, v, w )
+                                  usws, uswst, v, w )
              ENDIF
              CALL coriolis( i, j, 1 )
@@ -868,8 +868,8 @@
              THEN
                 CALL diffusion_v( i, j, ddzu, ddzw, km_m, km_damp_x, tend,     &
-                                  u_m, v_m, vsws_m, w_m )
+                                  u_m, v_m, vsws_m, vswst_m, w_m )
              ELSE
                 CALL diffusion_v( i, j, ddzu, ddzw, km, km_damp_x, tend, u, v, &
-                                  vsws, w )
+                                  vsws, vswst, w )
              ENDIF
              CALL coriolis( i, j, 2 )
@@ -1257,8 +1257,8 @@
     ENDIF
     IF ( tsc(2) == 2.0  .AND.  timestep_scheme(1:8) == 'leapfrog' )  THEN
-       CALL diffusion_u( ddzu, ddzw, km_m, km_damp_y, tend, u_m, usws_m, v_m, &
-                         w_m )
+       CALL diffusion_u( ddzu, ddzw, km_m, km_damp_y, tend, u_m, usws_m, &
+                         uswst_m, v_m, w_m )
     ELSE
-       CALL diffusion_u( ddzu, ddzw, km, km_damp_y, tend, u, usws, v, w )
+       CALL diffusion_u( ddzu, ddzw, km, km_damp_y, tend, u, usws, uswst, v, w )
     ENDIF
     CALL coriolis( 1 )
@@ -1330,7 +1330,7 @@
     IF ( tsc(2) == 2.0  .AND.  timestep_scheme(1:8) == 'leapfrog' )  THEN
        CALL diffusion_v( ddzu, ddzw, km_m, km_damp_x, tend, u_m, v_m, vsws_m, &
-                         w_m )
+                         vswst_m, w_m )
     ELSE
-       CALL diffusion_v( ddzu, ddzw, km, km_damp_x, tend, u, v, vsws, w )
+       CALL diffusion_v( ddzu, ddzw, km, km_damp_x, tend, u, v, vsws, vswst, w )
     ENDIF
     CALL coriolis( 2 )
Index: /palm/trunk/SOURCE/read_3d_binary.f90
===================================================================
--- /palm/trunk/SOURCE/read_3d_binary.f90	(revision 101)
+++ /palm/trunk/SOURCE/read_3d_binary.f90	(revision 102)
@@ -4,5 +4,5 @@
 ! Actual revisions:
 ! -----------------
-! 
+! +uswst, uswst_m, vswst, vswst_m
 !
 ! Former revisions:
@@ -333,6 +333,10 @@
           CASE ( 'usws' )
              READ ( 13 )  usws
+          CASE ( 'uswst' )
+             READ ( 13 )  uswst
           CASE ( 'usws_m' )
              READ ( 13 )  usws_m
+          CASE ( 'uswst_m' )
+             READ ( 13 )  uswst_m
           CASE ( 'us_av' )
              ALLOCATE( us_av(nys-1:nyn+1,nxl-1:nxr+1) )
@@ -366,6 +370,10 @@
           CASE ( 'vsws' )
              READ ( 13 )  vsws
+          CASE ( 'vswst' )
+             READ ( 13 )  vswst
           CASE ( 'vsws_m' )
              READ ( 13 )  vsws_m
+          CASE ( 'vswst_m' )
+             READ ( 13 )  vswst_m
           CASE ( 'w' )
              READ ( 13 )  w
Index: /palm/trunk/SOURCE/read_var_list.f90
===================================================================
--- /palm/trunk/SOURCE/read_var_list.f90	(revision 101)
+++ /palm/trunk/SOURCE/read_var_list.f90	(revision 102)
@@ -4,5 +4,5 @@
 ! Actual revisions:
 ! -----------------
-! 
+! +time_coupling, top_momentumflux_u|v
 !
 ! Former revisions:
@@ -358,4 +358,6 @@
           CASE ( 's_vertical_gradient_level' )
              READ ( 13 )  s_vertical_gradient_level
+          CASE ( 'time_coupling' )
+             READ ( 13 )  time_coupling
           CASE ( 'time_disturb' )
              READ ( 13 )  time_disturb
@@ -396,4 +398,8 @@
           CASE ( 'top_heatflux' )
              READ ( 13 )  top_heatflux
+          CASE ( 'top_momentumflux_u' )
+             READ ( 13 )  top_momentumflux_u
+          CASE ( 'top_momentumflux_v' )
+             READ ( 13 )  top_momentumflux_v
           CASE ( 'top_salinityflux' )
              READ ( 13 )  top_salinityflux
Index: /palm/trunk/SOURCE/surface_coupler.f90
===================================================================
--- /palm/trunk/SOURCE/surface_coupler.f90	(revision 102)
+++ /palm/trunk/SOURCE/surface_coupler.f90	(revision 102)
@@ -0,0 +1,142 @@
+ SUBROUTINE surface_coupler
+
+!------------------------------------------------------------------------------!
+! Actual revisions:
+! -----------------
+! 
+!
+! Former revisions:
+! ------------------
+! $Id$
+!
+! Initial revision
+!
+! Description:
+! ------------
+! Data exchange at the interface between coupled models
+!------------------------------------------------------------------------------!
+
+    USE arrays_3d
+    USE control_parameters
+    USE cpulog
+    USE grid_variables
+    USE indices
+    USE interfaces
+    USE pegrid
+
+    IMPLICIT NONE
+
+    INTEGER ::  k
+
+    REAL    ::  simulated_time_remote
+
+#if defined( __parallel )  &&  defined( __mpi2 )
+
+    CALL cpu_log( log_point(39), 'surface_coupler', 'start' )
+
+!
+!-- First exchange the current simulated time between the models,
+!-- currently just for testing
+    CALL MPI_SEND( simulated_time, 1, MPI_REAL, myid, 11, comm_inter, ierr )
+    CALL MPI_RECV( simulated_time_remote, 1, MPI_REAL, myid, 11, &
+                   comm_inter, status, ierr )
+    WRITE ( 9, * )  simulated_time, ' remote: ', simulated_time_remote
+    CALL local_flush( 9 )
+
+!
+!-- Exchange the interface data
+    IF ( coupling_mode == 'atmosphere_to_ocean' )  THEN
+
+!
+!--    Send heat flux at bottom surface to the ocean model
+       WRITE ( 9, * )  '*** send shf to ocean'
+       CALL local_flush( 9 )
+       CALL MPI_SEND( shf(nys-1,nxl-1), ngp_xy, MPI_REAL, myid, 12, &
+                      comm_inter, ierr )
+       WRITE ( 9, * )  '    ready'
+       CALL local_flush( 9 )
+
+!
+!--    Receive temperature at the bottom surface from the ocean model
+       WRITE ( 9, * )  '*** receive pt from ocean'
+       CALL local_flush( 9 )
+       CALL MPI_RECV( pt(0,nys-1,nxl-1), 1, type_xy, myid, 13, comm_inter, &
+                      status, ierr )
+       WRITE ( 9, * )  '    ready'
+       CALL local_flush( 9 )
+
+!
+!--    Send the momentum flux (u) at bottom surface to the ocean model
+       WRITE ( 9, * )  '*** send usws to ocean'
+       CALL local_flush( 9 )
+       CALL MPI_SEND( usws(nys-1,nxl-1), ngp_xy, MPI_REAL, myid, 14, &
+                      comm_inter, ierr )
+       WRITE ( 9, * )  '    ready'
+       CALL local_flush( 9 )
+
+!
+!--    Send the momentum flux (v) at bottom surface to the ocean model
+       WRITE ( 9, * )  '*** send vsws to ocean'
+       CALL local_flush( 9 )
+       CALL MPI_SEND( vsws(nys-1,nxl-1), ngp_xy, MPI_REAL, myid, 15, &
+                      comm_inter, ierr )
+       WRITE ( 9, * )  '    ready'
+       CALL local_flush( 9 )
+
+    ELSEIF ( coupling_mode == 'ocean_to_atmosphere' )  THEN
+
+!
+!--    Receive heat flux at the sea surface (top) from the atmosphere model
+       WRITE ( 9, * )  '*** receive tswst from atmosphere'
+       CALL local_flush( 9 )
+       CALL MPI_RECV( tswst(nys-1,nxl-1), ngp_xy, MPI_REAL, myid, 12, &
+                      comm_inter, status, ierr )
+       WRITE ( 9, * )  '    ready'
+       CALL local_flush( 9 )
+
+!
+!--    Adjust the kinematic heat flux with respect to ocean density
+!--    (constants are the specific heat capacities for air and water)
+       tswst = tswst / rho_surface * 1005.0 / 4218.0
+
+!
+!--    Send sea surface temperature to the atmosphere model
+       WRITE ( 9, * )  '*** send pt to atmosphere'
+       CALL local_flush( 9 )
+       CALL MPI_SEND( pt(nzt,nys-1,nxl-1), 1, type_xy, myid, 13, comm_inter, &
+                      ierr )
+       WRITE ( 9, * )  '    ready'
+       CALL local_flush( 9 )
+
+!
+!--    Receive momentum flux (u) at the sea surface (top) from the atmosphere
+!--    model
+       WRITE ( 9, * )  '*** receive uswst from atmosphere'
+       CALL local_flush( 9 )
+       CALL MPI_RECV( uswst(nys-1,nxl-1), ngp_xy, MPI_REAL, myid, 14, &
+                      comm_inter, status, ierr )
+       WRITE ( 9, * )  '    ready'
+       CALL local_flush( 9 )
+
+!
+!--    Receive momentum flux (v) at the sea surface (top) from the atmosphere
+!--    model
+       WRITE ( 9, * )  '*** receive vswst from atmosphere'
+       CALL local_flush( 9 )
+       CALL MPI_RECV( vswst(nys-1,nxl-1), ngp_xy, MPI_REAL, myid, 15, &
+                      comm_inter, status, ierr )
+       WRITE ( 9, * )  '    ready'
+       CALL local_flush( 9 )
+
+!
+!--    Adjust the momentum fluxes with respect to ocean density
+       uswst = uswst / rho_surface
+       vswst = vswst / rho_surface
+
+    ENDIF
+
+    CALL cpu_log( log_point(39), 'surface_coupler', 'stop' )
+
+#endif
+
+ END SUBROUTINE surface_coupler
Index: /palm/trunk/SOURCE/swap_timelevel.f90
===================================================================
--- /palm/trunk/SOURCE/swap_timelevel.f90	(revision 101)
+++ /palm/trunk/SOURCE/swap_timelevel.f90	(revision 102)
@@ -4,5 +4,5 @@
 ! Actual revisions:
 ! -----------------
-! 
+! swaping of uswst, vswst included
 !
 ! Former revisions:
@@ -136,4 +136,6 @@
                 ENDIF
                 IF ( use_top_fluxes )  THEN
+                   uswst_m => uswst_1;  uswst => uswst_2
+                   vswst_m => vswst_1;  vswst => vswst_2
                    tswst_m => tswst_1;  tswst => tswst_2
                    IF ( humidity  .OR.  passive_scalar )  THEN
@@ -184,5 +186,7 @@
                 ENDIF
                 IF ( use_top_fluxes )  THEN
-                   tswst_m  => tswst_2;  tswst => tswst_1
+                   uswst_m => uswst_2;  uswst => uswst_1
+                   vswst_m => vswst_2;  vswst => vswst_1
+                   tswst_m => tswst_2;  tswst => tswst_1
                    IF ( humidity  .OR.  passive_scalar )  THEN
                       qswst_m => qswst_2;  qswst => qswst_1
Index: /palm/trunk/SOURCE/time_integration.f90
===================================================================
--- /palm/trunk/SOURCE/time_integration.f90	(revision 101)
+++ /palm/trunk/SOURCE/time_integration.f90	(revision 102)
@@ -4,5 +4,5 @@
 ! Actual revisions:
 ! -----------------
-! 
+! Call of new routine surface_coupler
 !
 ! Former revisions:
@@ -303,4 +303,14 @@
 
 !
+!--    Data exchange between coupled models
+       IF ( coupling_mode /= 'uncoupled' )  THEN
+          time_coupling = time_coupling + dt_3d
+          DO WHILE ( time_coupling >= dt_coupling )
+             CALL surface_coupler
+             time_coupling = time_coupling - dt_coupling
+          ENDDO
+       ENDIF
+
+!
 !--    Execute user-defined actions
        CALL user_actions( 'after_integration' )
Index: /palm/trunk/SOURCE/write_3d_binary.f90
===================================================================
--- /palm/trunk/SOURCE/write_3d_binary.f90	(revision 101)
+++ /palm/trunk/SOURCE/write_3d_binary.f90	(revision 102)
@@ -4,5 +4,5 @@
 ! Actual revisions:
 ! -----------------
-! 
+! +uswst, uswst_m, vswst, vswst_m
 !
 ! Former revisions:
@@ -201,6 +201,8 @@
     WRITE ( 14 )  'us                  ';  WRITE ( 14 )  us
     WRITE ( 14 )  'usws                ';  WRITE ( 14 )  usws
+    WRITE ( 14 )  'uswst               ';  WRITE ( 14 )  uswst
     IF ( timestep_scheme(1:5) /= 'runge' )  THEN
        WRITE ( 14 )  'usws_m              ';  WRITE ( 14 )  usws_m
+       WRITE ( 14 )  'uswst_m             ';  WRITE ( 14 )  uswst_m
     ENDIF
     IF ( ALLOCATED( us_av ) )  THEN
@@ -236,6 +238,8 @@
     ENDIF
     WRITE ( 14 )  'vsws                ';  WRITE ( 14 )  vsws
+    WRITE ( 14 )  'vswst               ';  WRITE ( 14 )  vswst
     IF ( timestep_scheme(1:5) /= 'runge' )  THEN
        WRITE ( 14 )  'vsws_m              ';  WRITE ( 14 )  vsws_m
+       WRITE ( 14 )  'vswst_m             ';  WRITE ( 14 )  vswst_m
     ENDIF
     WRITE ( 14 )  'w                   ';  WRITE ( 14 )  w
Index: /palm/trunk/SOURCE/write_var_list.f90
===================================================================
--- /palm/trunk/SOURCE/write_var_list.f90	(revision 101)
+++ /palm/trunk/SOURCE/write_var_list.f90	(revision 102)
@@ -4,5 +4,5 @@
 ! Actual revisions:
 ! -----------------
-! 
+! +top_momentumflux_u|v, time_coupling
 !
 ! Former revisions:
@@ -304,4 +304,6 @@
     WRITE ( 14 )  's_vertical_gradient_level     '
     WRITE ( 14 )  s_vertical_gradient_level
+    WRITE ( 14 )  'time_coupling                 '
+    WRITE ( 14 )  time_coupling
     WRITE ( 14 )  'time_disturb                  '
     WRITE ( 14 )  time_disturb
@@ -342,4 +344,8 @@
     WRITE ( 14 )  'top_heatflux                  '
     WRITE ( 14 )  top_heatflux
+    WRITE ( 14 )  'top_momentumflux_u            '
+    WRITE ( 14 )  top_momentumflux_u
+    WRITE ( 14 )  'top_momentumflux_v            '
+    WRITE ( 14 )  top_momentumflux_v
     WRITE ( 14 )  'top_salinityflux              '
     WRITE ( 14 )  top_salinityflux
