Index: palm/trunk/DOC/app/chapter_4.5.7.html
===================================================================
--- palm/trunk/DOC/app/chapter_4.5.7.html	(revision 377)
+++ palm/trunk/DOC/app/chapter_4.5.7.html	(revision 388)
@@ -1,10 +1,10 @@
 <!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.0 Transitional//EN">
 <html><head>
+
 <meta http-equiv="CONTENT-TYPE" content="text/html; charset=windows-1252"><title>PALM chapter 4.5.7</title> <meta name="GENERATOR" content="StarOffice 7 (Win32)"> <meta name="AUTHOR" content="Siegfried Raasch"> <meta name="CREATED" content="20041029;14344622"> <meta name="CHANGED" content="20050119;9531085"> <meta name="KEYWORDS" content="parallel LES model"> <style>
 <!--
 @page { size: 21cm 29.7cm }
 -->
-</style></head>
-<body style="direction: ltr;" lang="en-US"><h4 style="line-height: 100%;"><font size="4">4.5.7
+</style></head><body style="direction: ltr;" lang="en-US"><h4 style="line-height: 100%;"><font size="4">4.5.7
 Plots of
 isosurfaces, 2d cross sections and particles with dvrp</font></h4>
@@ -14,5 +14,5 @@
 system. If you are interested in using <span style="font-weight: bold;">dvrp</span> on your system, please contact the PALM developers.</p><br><span style="font-weight: bold;">General remarks:</span><p style="line-height: 100%;">The <span style="font-weight: bold;">dvrp</span>
 software was originally developed by the RRZN
-(Stephan Olbrich, Carsten Chmielewski) and is meanwhile continuosly developed and improved under the name <span style="font-weight: bold;">dsvr</span> by the University of D&uuml;sseldorf (Prof. Stephan Olbrich, see webpage of the <a href="http://www.dsvr-software.de/">dsvr-software</a>). It allows to create 3d-animations with PALM,
+(Stephan Olbrich, Carsten Chmielewski) and is meanwhile continuosly developed and improved under the name <span style="font-weight: bold;">dsvr</span> by the University of Düsseldorf (Prof. Stephan Olbrich, see webpage of the <a href="http://www.dsvr-software.de/">dsvr-software</a>). It allows to create 3d-animations with PALM,
 which can be displayed via a special plugin for internet browsers. With
 suitable graphic hardware (e.g. NVIDIA quattro FX cards) even stereoscopic views are
@@ -50,5 +50,5 @@
 Methoden zur visuellen und haptischen tele-immersiven
 Exploration
-komplexer Volumen- und Str&ouml;mungsdaten aus
+komplexer Volumen- und Strömungsdaten aus
 parallelisierten,
 dynamischen 3D-Simulationen" (2005-2007). </p>
@@ -57,5 +57,5 @@
 PALM software package (see chapter <a href="chapter_3.7.html">3.7</a>).
 To use this package, the additional option<tt style="font-family: monospace;"> </tt><font style="font-family: Courier New,Courier,monospace;" size="2">-p
-<span style="font-family: Helvetica,Arial,sans-serif;">&ldquo;</span>dvrp_graphics<span style="font-family: Helvetica,Arial,sans-serif;">&rdquo;</span></font><span style="font-family: monospace;">
+<span style="font-family: Helvetica,Arial,sans-serif;"></span>dvrp_graphics<span style="font-family: Helvetica,Arial,sans-serif;"></span></font><span style="font-family: monospace;">
 </span>has to be given in the <b>mrun</b> call. This
 automatically links the
@@ -75,5 +75,7 @@
 parallel</span><br><br><br></li><li>Add the dvrp-steering parameters to your NAMELIST-parameter file, e.g.<br><br><span style="font-family: Courier New,Courier,monospace;">&nbsp;&amp;d3par &nbsp;end_time = 3600.0,<br>&nbsp; ... &nbsp; &nbsp; /<br><br>&nbsp;&amp;dvrp_graphics_par &nbsp; dvrp_username = '<span style="font-weight: bold;">&lt;UNAME&gt;</span>',<br>&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; dvrp_host = '130.75.105.2',<br>&nbsp;
 &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;
-dvrp_directory = '<span style="font-weight: bold;">&lt;UNAME&gt;</span>/<span style="font-weight: bold;">&lt;UDIR&gt;</span>', ... /<br></span><br>For other dvrp-parameters see <a href="chapter_4.2.html#Paketparameter">chapter 4.2</a>. An example parameter file can be found in directory <span style="font-family: Courier New,Courier,monospace;">..../trunk/EXAMPLES/dvr &nbsp;</span>.<span style="font-family: Courier New,Courier,monospace;"><br><br></span></li><li>On the IMUK-system, create a subdirectory in which the dvrp-data are stored.<br><br><span style="font-family: Courier New,Courier,monospace;">&nbsp; &nbsp;mkdir /data/raasch/Dvrp_daten/<span style="font-weight: bold;">&lt;UNAME&gt;</span></span><br><br>This directory has to be created only once, before dvrp is used for the first time.<br></li><li>Submit the job&nbsp; with <span style="font-weight: bold;">mrun</span>-command<br><br><span style="font-family: Courier New,Courier,monospace;">&nbsp; &nbsp;mrun .... -p dvrp_graphics ....</span><br><br></li><li><span style="font-family: Courier New,Courier,monospace;"><span style="font-family: Times New Roman,Times,serif;">After the job has finished, the dvrp-output can be found on the IMUK-cluster in a subdirectory under <span style="font-family: Courier New,Courier,monospace;">/data/raasch/Dvrp_daten</span>. The name of the subdirectory is determined by the dvrp-parameter <a href="chapter_4.2.html#dvrp_directory">dvrp_directory</a> (see above), i.e. if the user has set <span style="font-family: Courier New,Courier,monospace;">dvrp_directory</span> = <span style="font-style: italic;">'<span style="font-weight: bold;">&lt;UNAME&gt;</span>/movie_1'</span>, the dvrp-data are stored under&nbsp;<span style="font-family: Courier New,Courier,monospace;">/data/raasch/Dvrp_daten/<span style="font-weight: bold;">&lt;UNAME&gt;</span>/movie_1</span>.<br><br></span></span></li><li><span style="font-family: Courier New,Courier,monospace;"><span style="font-family: Times New Roman,Times,serif;">Change to this subdirectory (e.g. </span></span><span style="font-family: Courier New,Courier,monospace;"><span style="font-family: Times New Roman,Times,serif;"><span style="font-family: Courier New,Courier,monospace;">/data/raasch/Dvrp_daten/&lt;UNAME&gt;/movie_1</span></span></span><span style="font-family: Courier New,Courier,monospace;"><span style="font-family: Times New Roman,Times,serif;">) and enter the command<br><br><span style="font-family: Courier New,Courier,monospace;">&nbsp; &nbsp;process_dvr_output</span><br><br>It will create a file with name <span style="font-family: Courier New,Courier,monospace;">all_streams_streaming.html</span>.<br><br></span></span></li><li><span style="font-family: Courier New,Courier,monospace;"><span style="font-family: Times New Roman,Times,serif;">Call the opera-browser (i.e. enter the command <span style="font-family: Courier New,Courier,monospace;">opera</span>) and open the file </span></span><span style="font-family: Courier New,Courier,monospace;"><span style="font-family: Times New Roman,Times,serif;"><span style="font-family: Courier New,Courier,monospace;">all_streams_streaming.html</span>.</span></span><span style="font-family: Courier New,Courier,monospace;"><span style="font-family: Times New Roman,Times,serif;"></span><br></span><br></li></ol><span style="font-weight: bold;">Further features:<br></span>The script<span style="font-family: Courier New,Courier,monospace;"> process_dvr_output </span>has additional options:<br><br><span style="font-weight: bold; font-family: Courier New,Courier,monospace;">-s</span> &nbsp; &nbsp; : <span style="font-weight: bold;">create sequence output</span>.
+dvrp_directory = '<span style="font-weight: bold;">&lt;UNAME&gt;</span>/<span style="font-weight: bold;">&lt;UDIR&gt;</span>', ... /<br></span><br>For other dvrp-parameters see <a href="chapter_4.2.html#Paketparameter">chapter 4.2</a>. An example parameter file can be found in directory <span style="font-family: Courier New,Courier,monospace;">..../trunk/EXAMPLES/dvr &nbsp;</span>.<span style="font-family: Courier New,Courier,monospace;"><br><br></span></li><li>On the IMUK-system, create a subdirectory in which the dvrp-data are stored and give it access permit for group <span style="font-style: italic;">palm</span>.<br><br><span style="font-family: Courier New,Courier,monospace;">&nbsp; &nbsp;mkdir /data/raasch/Dvrp_daten/<span style="font-weight: bold;">&lt;UNAME&gt;<br>
+</span>&nbsp;&nbsp; chmod g+rwx</span><span style="font-family: Courier New,Courier,monospace;"> /data/raasch/Dvrp_daten/<span style="font-weight: bold;">&lt;UNAME&gt;</span></span><br><br>This directory has to be created only once, before dvrp is used for the first time.<br>
+<br></li><li>Submit the job&nbsp; with <span style="font-weight: bold;">mrun</span>-command<br><br><span style="font-family: Courier New,Courier,monospace;">&nbsp; &nbsp;mrun .... -p dvrp_graphics ....</span><br><br></li><li><span style="font-family: Courier New,Courier,monospace;"><span style="font-family: Times New Roman,Times,serif;">After the job has finished, the dvrp-output can be found on the IMUK-cluster in a subdirectory under <span style="font-family: Courier New,Courier,monospace;">/data/raasch/Dvrp_daten</span>. The name of the subdirectory is determined by the dvrp-parameter <a href="chapter_4.2.html#dvrp_directory">dvrp_directory</a> (see above), i.e. if the user has set <span style="font-family: Courier New,Courier,monospace;">dvrp_directory</span> = <span style="font-style: italic;">'<span style="font-weight: bold;">&lt;UNAME&gt;</span>/movie_1'</span>, the dvrp-data are stored under&nbsp;<span style="font-family: Courier New,Courier,monospace;">/data/raasch/Dvrp_daten/<span style="font-weight: bold;">&lt;UNAME&gt;</span>/movie_1</span>.<br><br></span></span></li><li><span style="font-family: Courier New,Courier,monospace;"><span style="font-family: Times New Roman,Times,serif;">Change to this subdirectory (e.g. </span></span><span style="font-family: Courier New,Courier,monospace;"><span style="font-family: Times New Roman,Times,serif;"><span style="font-family: Courier New,Courier,monospace;">/data/raasch/Dvrp_daten/&lt;UNAME&gt;/movie_1</span></span></span><span style="font-family: Courier New,Courier,monospace;"><span style="font-family: Times New Roman,Times,serif;">) and enter the command<br><br><span style="font-family: Courier New,Courier,monospace;">&nbsp; &nbsp;process_dvr_output</span><br><br>It will create a file with name <span style="font-family: Courier New,Courier,monospace;">all_streams_streaming.html</span>.<br><br></span></span></li><li><span style="font-family: Courier New,Courier,monospace;"><span style="font-family: Times New Roman,Times,serif;">Call the opera-browser (i.e. enter the command <span style="font-family: Courier New,Courier,monospace;">opera</span>) and open the file </span></span><span style="font-family: Courier New,Courier,monospace;"><span style="font-family: Times New Roman,Times,serif;"><span style="font-family: Courier New,Courier,monospace;">all_streams_streaming.html</span>.</span></span><span style="font-family: Courier New,Courier,monospace;"><span style="font-family: Times New Roman,Times,serif;"></span><br></span><br></li></ol><span style="font-weight: bold;">Further features:<br></span>The script<span style="font-family: Courier New,Courier,monospace;"> process_dvr_output </span>has additional options:<br><br><span style="font-weight: bold; font-family: Courier New,Courier,monospace;">-s</span> &nbsp; &nbsp; : <span style="font-weight: bold;">create sequence output</span>.
 All streams (and static scenes) are collected to one sequence (one
 static scene), which can be displayed with a browser without using the
@@ -84,5 +86,5 @@
 all_streams_sequence.html to a directory on your local computer and
 open the file all_streams_sequence.html with your browser. For this,
-you will need to install the the dvr-plugin (see webpage of the <a href="http://www.dsvr-software.de/">dsvr-software</a>) on your local computer, which is also available for Windows.<br><br><span style="font-weight: bold;">-a</span> &nbsp; &nbsp; &nbsp;: <span style="font-weight: bold;">acceleration factor</span>. If, in case of sequence mode, the performance of the animation is to slow, you can accelerate it. E.g., by using <span style="font-family: Courier New,Courier,monospace;">"-a 2</span>", only every second frame of the original streams will be used for the sequence.<br><br><br><br><br><span style="font-weight: bold;">Current limitations (May 09):</span><br>Only a special opera-version on host "bora" can be used. Log-in on "bora" and call &nbsp;<span style="font-family: Courier New,Courier,monospace;">/usr/local/bin/opera</span>. <span style="font-weight: bold;">Before that(!!!)</span>, enter the command "<span style="font-family: Courier New,Courier,monospace;">export LD_PRELOAD=libXm.so</span>" !
+you will need to install the the dvr-plugin (see webpage of the <a href="http://www.dsvr-software.de/">dsvr-software</a>) on your local computer, which is also available for Windows.<br><br><span style="font-weight: bold;">-a</span> &nbsp; &nbsp; &nbsp;: <span style="font-weight: bold;">acceleration factor</span>. If, in case of sequence mode, the performance of the animation is to slow, you can accelerate it. E.g., by using <span style="font-family: Courier New,Courier,monospace;">"-a 2</span>", only every second frame of the original streams will be used for the sequence.<br><br><br><br><br><span style="font-weight: bold;">Current limitations (May 09):</span><br>Only a special opera-version on host "bora" can be used. Log-in on "bora" and call &nbsp;<span style="font-family: Courier New,Courier,monospace;">/usr/bin/opera</span>. <span style="font-weight: bold;">Before that(!!!)</span>, enter the command "<span style="font-family: Courier New,Courier,monospace;">export LD_PRELOAD=libXm.so</span>" !
 <hr><p style="line-height: 100%;"><br>
 <font color="#000080"><font color="#000080"><a href="chapter_4.5.6.html"><font color="#000080"><img src="left.gif" name="Grafik1" align="bottom" border="2" height="32" width="32"></font></a><a href="index.html"><font color="#000080"><img src="up.gif" name="Grafik2" align="bottom" border="2" height="32" width="32"></font></a><a href="chapter_4.6.html"><font color="#000080"><img src="right.gif" name="Grafik3" align="bottom" border="2" height="32" width="32"></font></a></font></font></p><p style="line-height: 100%;">&nbsp;<i>Last
Index: palm/trunk/SCRIPTS/mrun
===================================================================
--- palm/trunk/SCRIPTS/mrun	(revision 377)
+++ palm/trunk/SCRIPTS/mrun	(revision 388)
@@ -3187,4 +3187,7 @@
                       export MPI_CONNECTIONS_THRESHOLD=8192
                       mpiexec_mpt -np $ii   ./a.out  $ROPTS  < runfile_atmos
+
+                          # next is test for openmp usage
+                      # mpiexec -n $ii -pernode  ./a.out  $ROPTS  < runfile_atmos
                    elif [[ $( echo $mpilib | cut -c1-3 ) = mva ]]
                    then
Index: palm/trunk/SOURCE/CURRENT_MODIFICATIONS
===================================================================
--- palm/trunk/SOURCE/CURRENT_MODIFICATIONS	(revision 377)
+++ palm/trunk/SOURCE/CURRENT_MODIFICATIONS	(revision 388)
@@ -125,4 +125,10 @@
 Errors:
 ------
+Bugfix: Initial hydrostatic pressure profile in case of ocean runs is now
+calculated in 5 iteration steps. (init_ocean)
+
+Bugfix: wrong sign in buoyancy production of ocean part in case of not using
+the reference density (only in 3D routine production_e) (production_e)
+
 Bugfix: output of averaged 2d/3d quantities requires that an avaraging
 interval has been set, respective error message is included (check_parameters)
@@ -183,3 +189,3 @@
 Bugfix: initial setting of time_coupling in coupled restart runs (time_integration)
 
-advec_particles, check_parameters, cpu_log, data_output_2d, data_output_3d, header, init_3d_model, init_particles, modules, netcdf, prandtl_fluxes, production_e, read_var_list, time_integration, user_last_actions, write_var_list
+advec_particles, check_parameters, cpu_log, data_output_2d, data_output_3d, header, init_3d_model, init_particles, init_ocean, modules, netcdf, prandtl_fluxes, production_e, read_var_list, time_integration, user_last_actions, write_var_list
Index: palm/trunk/SOURCE/check_parameters.f90
===================================================================
--- palm/trunk/SOURCE/check_parameters.f90	(revision 377)
+++ palm/trunk/SOURCE/check_parameters.f90	(revision 388)
@@ -4,4 +4,5 @@
 ! Actual revisions:
 ! -----------------
+! Check profiles fpr prho and hyp.
 ! Bugfix: output of averaged 2d/3d quantities requires that an avaraging
 ! interval has been set, respective error message is included
@@ -933,5 +934,6 @@
 
 !
-!--    If required compute the profile of leaf area density used in the plant canopy model 
+!--    If required compute the profile of leaf area density used in the plant
+!--    canopy model 
        IF ( plant_canopy ) THEN
        
@@ -967,5 +969,6 @@
 
 !
-!--       In case of no given leaf area density gradients, choose a vanishing gradient
+!--       In case of no given leaf area density gradients, choose a vanishing
+!--       gradient
           IF ( lad_vertical_gradient_level(1) == -9999999.9 ) THEN
              lad_vertical_gradient_level(1) = 0.0
@@ -2089,7 +2092,14 @@
 
           CASE ( 'rho' )
-             dopr_index(i) = 64
-             dopr_unit(i)  = 'kg/m3'
-             hom(:,2,64,:) = SPREAD( zu, 2, statistic_regions+1 )
+             IF ( .NOT. ocean ) THEN
+                message_string = 'data_output_pr = ' // &
+                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
+                                 'lemented for ocean = .FALSE.'
+                CALL message( 'check_parameters', 'PA0091', 1, 2, 0, 6, 0 )
+             ELSE
+                dopr_index(i) = 64
+                dopr_unit(i)  = 'kg/m3'
+                hom(:,2,64,:) = SPREAD( zu, 2, statistic_regions+1 )
+             ENDIF
 
           CASE ( 'w"sa"' )
@@ -2150,4 +2160,21 @@
                 hom(:,2,70,:) = SPREAD( zu, 2, statistic_regions+1 )
              ENDIF
+
+          CASE ( 'prho' )
+             IF ( .NOT. ocean ) THEN
+                message_string = 'data_output_pr = ' // &
+                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
+                                 'lemented for ocean = .FALSE.'
+                CALL message( 'check_parameters', 'PA0091', 1, 2, 0, 6, 0 )
+             ELSE
+                dopr_index(i) = 71
+                dopr_unit(i)  = 'kg/m3'
+                hom(:,2,71,:) = SPREAD( zu, 2, statistic_regions+1 )
+             ENDIF
+
+          CASE ( 'hyp' )
+             dopr_index(i) = 72
+             dopr_unit(i)  = 'kPa'
+             hom(:,2,72,:) = SPREAD( zu, 2, statistic_regions+1 )
 
           CASE DEFAULT
@@ -2929,6 +2956,6 @@
 !-- Check pressure gradient conditions
     IF ( dp_external .AND. conserve_volume_flow )  THEN
-       WRITE( message_string, * )  'Both dp_external and conserve_volume_flow', &
-            ' are .TRUE. but one of them must be .FALSE.'
+       WRITE( message_string, * )  'Both dp_external and conserve_volume_flo', &
+            'w are .TRUE. but one of them must be .FALSE.'
        CALL message( 'check_parameters', 'PA0150', 1, 2, 0, 6, 0 )
     ENDIF
@@ -2940,6 +2967,6 @@
        ENDIF
        IF ( .NOT. ANY( dpdxy /= 0.0 ) )  THEN
-          WRITE( message_string, * )  'dp_external is .TRUE. but dpdxy is zero',&
-               ', i.e. the external pressure gradient & will not be applied'
+          WRITE( message_string, * )  'dp_external is .TRUE. but dpdxy is ze', &
+               'ro, i.e. the external pressure gradient & will not be applied'
           CALL message( 'check_parameters', 'PA0152', 0, 1, 0, 6, 0 )
        ENDIF
Index: palm/trunk/SOURCE/eqn_state_seawater.f90
===================================================================
--- palm/trunk/SOURCE/eqn_state_seawater.f90	(revision 377)
+++ palm/trunk/SOURCE/eqn_state_seawater.f90	(revision 388)
@@ -4,5 +4,6 @@
 ! Actual revisions:
 ! -----------------
-! First constant in array den also defined as type double
+! Potential density is additionally calculated in eqn_state_seawater,
+! first constant in array den also defined as type double
 !
 ! Former revisions:
@@ -69,5 +70,5 @@
        INTEGER ::  i, j, k
 
-       REAL ::  p1, p2, p3, pt1, pt2, pt3, pt4, sa1, sa15, sa2
+       REAL ::  pden, pnom, p1, p2, p3, pt1, pt2, pt3, pt4, sa1, sa15, sa2
 
        DO  i = nxl, nxr
@@ -76,6 +77,4 @@
 !
 !--             Pressure is needed in dbar
-!                p1 = hyp(0) * 1E-4
-!                p1 = 0.0
                 p1 = hyp(k) * 1E-4
                 p2 = p1 * p1
@@ -93,22 +92,35 @@
                 sa2  = sa1 * sa1
 
-                rho(k,j,i) =                                                   &
-                        ( nom(1)           + nom(2)*pt1     + nom(3)*pt2     + &
-                          nom(4)*pt3       + nom(5)*sa1     + nom(6)*sa1*pt1 + &
-                          nom(7)*sa2       + nom(8)*p1      + nom(9)*p1*pt2  + &
-                          nom(10)*p1*sa1   + nom(11)*p2     + nom(12)*p2*pt2   &
-                        ) /                                                    &
-                        ( den(1)           + den(2)*pt1     + den(3)*pt2     + &
-                          den(4)*pt3       + den(5)*pt4     + den(6)*sa1     + &
-                          den(7)*sa1*pt1   + den(8)*sa1*pt3 + den(9)*sa15    + &
-                          den(10)*sa15*pt2 + den(11)*p1     + den(12)*p2*pt3 + &
-                          den(13)*p3*pt1                                       &
-                        )
+                pnom = nom(1)           + nom(2)*pt1     + nom(3)*pt2     + &
+                       nom(4)*pt3       + nom(5)*sa1     + nom(6)*sa1*pt1 + &
+                       nom(7)*sa2
+
+                pden = den(1)           + den(2)*pt1     + den(3)*pt2     + &
+                       den(4)*pt3       + den(5)*pt4     + den(6)*sa1     + &
+                       den(7)*sa1*pt1   + den(8)*sa1*pt3 + den(9)*sa15    + &
+                       den(10)*sa15*pt2
+
+!
+!--             Potential density (without pressure terms)
+                prho(k,j,i) = pnom / pden
+
+                pnom = pnom +             nom(8)*p1      + nom(9)*p1*pt2  + &
+                       nom(10)*p1*sa1   + nom(11)*p2     + nom(12)*p2*pt2
+
+                pden = pden +             den(11)*p1     + den(12)*p2*pt3 + &
+                       den(13)*p3*pt1
+
+!
+!--             In-situ density
+                rho(k,j,i) = pnom / pden
 
              ENDDO
 !
 !--          Neumann conditions are assumed at bottom and top boundary
-             rho(nzt+1,j,i)            = rho(nzt,j,i)
-             rho(nzb_s_inner(j,i),j,i) = rho(nzb_s_inner(j,i)+1,j,i)
+             prho(nzt+1,j,i)            = prho(nzt,j,i)
+             prho(nzb_s_inner(j,i),j,i) = prho(nzb_s_inner(j,i)+1,j,i)
+             rho(nzt+1,j,i)             = rho(nzt,j,i)
+             rho(nzb_s_inner(j,i),j,i)  = rho(nzb_s_inner(j,i)+1,j,i)
+
           ENDDO
        ENDDO
@@ -129,11 +141,9 @@
        INTEGER ::  i, j, k
 
-       REAL ::  p1, p2, p3, pt1, pt2, pt3, pt4, sa1, sa15, sa2
+       REAL ::  pden, pnom, p1, p2, p3, pt1, pt2, pt3, pt4, sa1, sa15, sa2
 
        DO  k = nzb_s_inner(j,i)+1, nzt
 !
 !--       Pressure is needed in dbar
-!          p1 = hyp(0) * 1E-4
-!          p1 = 0.0
           p1 = hyp(k) * 1E-4
           p2 = p1 * p1
@@ -151,20 +161,35 @@
           sa2  = sa1 * sa1
 
-          rho(k,j,i) = ( nom(1)           + nom(2)*pt1     + nom(3)*pt2     + &
-                         nom(4)*pt3       + nom(5)*sa1     + nom(6)*sa1*pt1 + &
-                         nom(7)*sa2       + nom(8)*p1      + nom(9)*p1*pt2  + &
-                         nom(10)*p1*sa1   + nom(11)*p2     + nom(12)*p2*pt2   &
-                       ) /                                                    &
-                       ( den(1)           + den(2)*pt1     + den(3)*pt2     + &
-                         den(4)*pt3       + den(5)*pt4     + den(6)*sa1     + &
-                         den(7)*sa1*pt1   + den(8)*sa1*pt3 + den(9)*sa15    + &
-                         den(10)*sa15*pt2 + den(11)*p1     + den(12)*p2*pt3 + &
-                         den(13)*p3*pt1                                       &
-                       )
+          pnom = nom(1)           + nom(2)*pt1     + nom(3)*pt2     + &
+                 nom(4)*pt3       + nom(5)*sa1     + nom(6)*sa1*pt1 + &
+                 nom(7)*sa2
+
+          pden = den(1)           + den(2)*pt1     + den(3)*pt2     + &
+                 den(4)*pt3       + den(5)*pt4     + den(6)*sa1     + &
+                 den(7)*sa1*pt1   + den(8)*sa1*pt3 + den(9)*sa15    + &
+                 den(10)*sa15*pt2
+
+!
+!--       Potential density (without pressure terms)
+          prho(k,j,i) = pnom / pden
+
+          pnom = pnom +             nom(8)*p1      + nom(9)*p1*pt2  + &
+                 nom(10)*p1*sa1   + nom(11)*p2     + nom(12)*p2*pt2
+          pden = pden +             den(11)*p1     + den(12)*p2*pt3 + &
+                 den(13)*p3*pt1
+
+!
+!--       In-situ density
+          rho(k,j,i) = pnom / pden
+
+
        ENDDO
+
 !
 !--    Neumann conditions are assumed at bottom and top boundary
-       rho(nzt+1,j,i)            = rho(nzt,j,i)
-       rho(nzb_s_inner(j,i),j,i) = rho(nzb_s_inner(j,i)+1,j,i)
+       prho(nzt+1,j,i)            = prho(nzt,j,i)
+       prho(nzb_s_inner(j,i),j,i) = prho(nzb_s_inner(j,i)+1,j,i)
+       rho(nzt+1,j,i)             = rho(nzt,j,i)
+       rho(nzb_s_inner(j,i),j,i)  = rho(nzb_s_inner(j,i)+1,j,i)
 
     END SUBROUTINE eqn_state_seawater_ij
Index: palm/trunk/SOURCE/flow_statistics.f90
===================================================================
--- palm/trunk/SOURCE/flow_statistics.f90	(revision 377)
+++ palm/trunk/SOURCE/flow_statistics.f90	(revision 388)
@@ -4,4 +4,6 @@
 ! Current revisions:
 ! -----------------
+! Vertical profiles of potential density and hydrostatic pressure are
+! calculated.
 ! Added missing timeseries calculation of w"q"(0), moved timeseries q* to the
 ! end.
@@ -531,5 +533,5 @@
                 IF ( humidity )  THEN
                    sums_l(nzt,48,tn) = sums_l(nzt,48,tn) + &
-                                       qswst(j,i) * rmask(j,i,sr)  ! w"q" (w"qv")
+                                       qswst(j,i) * rmask(j,i,sr) ! w"q" (w"qv")
                    IF ( cloud_physics )  THEN
                       sums_l(nzt,45,tn) = sums_l(nzt,45,tn) + (           &
@@ -546,5 +548,5 @@
                 IF ( passive_scalar )  THEN
                    sums_l(nzt,48,tn) = sums_l(nzt,48,tn) + &
-                                       qswst(j,i) * rmask(j,i,sr)  ! w"q" (w"qv")
+                                       qswst(j,i) * rmask(j,i,sr) ! w"q" (w"qv")
                 ENDIF
              ENDIF
@@ -597,4 +599,6 @@
                                                        rmask(j,i,sr)
                    sums_l(k,64,tn) = sums_l(k,64,tn) + rho(k,j,i) * &
+                                                       rmask(j,i,sr)
+                   sums_l(k,71,tn) = sums_l(k,71,tn) + prho(k,j,i) * &
                                                        rmask(j,i,sr)
                 ENDIF
@@ -644,5 +648,8 @@
 !
 !--    Density at top follows Neumann condition
-       IF ( ocean )  sums_l(nzt+1,64,tn) = sums_l(nzt,64,tn)
+       IF ( ocean )  THEN
+          sums_l(nzt+1,64,tn) = sums_l(nzt,64,tn)
+          sums_l(nzt+1,71,tn) = sums_l(nzt,71,tn)
+       ENDIF
 
 !
@@ -897,4 +904,6 @@
        hom(:,1,69,sr) = sums(:,69)     ! w"e + w"p"/rho
        hom(:,1,70,sr) = sums(:,70)     ! q*2
+       hom(:,1,71,sr) = sums(:,71)     ! prho
+       hom(:,1,72,sr) = hyp * 1E-4     ! hyp in kPa
 
        hom(:,1,pr_palm-1,sr) = sums(:,pr_palm-1)
@@ -916,5 +925,4 @@
 !--    is less than 1.5 times the height where the heat flux becomes negative
 !--    (positive) for the first time.
-!--    NOTE: This criterion is still capable of improving!
        z_i(1) = 0.0
        first = .TRUE.
Index: palm/trunk/SOURCE/init_3d_model.f90
===================================================================
--- palm/trunk/SOURCE/init_3d_model.f90	(revision 377)
+++ palm/trunk/SOURCE/init_3d_model.f90	(revision 388)
@@ -7,4 +7,5 @@
 ! Current revisions:
 ! -----------------
+! Initialization of prho added.
 ! bugfix: correction of initial volume flow for non-flat topography
 ! bugfix: zero initialization of arrays within buildings for 'cyclic_fill'
@@ -234,10 +235,12 @@
        ALLOCATE( saswsb_1(nys-1:nyn+1,nxl-1:nxr+1), &
                  saswst_1(nys-1:nyn+1,nxl-1:nxr+1) )
-       ALLOCATE( rho_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
-                 sa_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1),  &
-                 sa_2(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1),  &
+       ALLOCATE( prho_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
+                 rho_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1),  &
+                 sa_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1),   &
+                 sa_2(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1),   &
                  sa_3(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
-       rho => rho_1  ! routine calc_mean_profile requires density to be a
-                     ! pointer
+       prho => prho_1
+       rho  => rho_1  ! routines calc_mean_profile and diffusion_e require
+                      ! density to be apointer
        IF ( humidity_remote )  THEN
           ALLOCATE( qswst_remote(nys-1:nyn+1,nxl-1:nxr+1) )
@@ -1328,4 +1331,5 @@
 !--    Initialize quantities needed for the ocean model
        CALL init_ocean
+
     ELSE
 !
Index: palm/trunk/SOURCE/init_ocean.f90
===================================================================
--- palm/trunk/SOURCE/init_ocean.f90	(revision 377)
+++ palm/trunk/SOURCE/init_ocean.f90	(revision 388)
@@ -4,5 +4,6 @@
 ! Actual revisions:
 ! -----------------
-! Bugfix: First calculation of hyp(0) changed
+! Bugfix: Initial profiles of hydrostatic pressure and density are calculated
+! iteratively. First calculation of hyp(0) changed.
 !
 ! Former revisions:
@@ -32,5 +33,5 @@
     INTEGER ::  k, n
 
-    REAL    ::  sa_l, pt_l, rho_l
+    REAL    ::  sa_l, pt_l
 
     REAL, DIMENSION(nzb:nzt+1) ::  rho_init
@@ -45,4 +46,5 @@
 !-- Calculate initial vertical profile of hydrostatic pressure (in Pa)
 !-- and the reference density (used later in buoyancy term)
+!-- First step: Calculate pressure using reference density
     hyp(nzt+1) = surface_pressure * 100.0
 
@@ -55,12 +57,7 @@
     hyp(0) = hyp(1) + rho_surface * g * dzu(1)
 
-    IF ( myid == 0 )  THEN
-       print*,'hydro pres using rho_surface'
-       DO  k = nzt+1, 0, -1
-          print*, 'k = ', k, ' hyp = ', hyp(k)
-       ENDDO
-       print*, ' '
-    ENDIF
-
+!
+!-- Second step: Iteratively calculate in situ density (based on presssure)
+!-- and pressure (based on in situ density)
     DO  n = 1, 5
 
@@ -85,12 +82,4 @@
        ENDDO
 
-       IF ( myid == 0 )  THEN
-          print*,'hydro pres / rho  n = ', n
-          DO  k = nzt+1, 0, -1
-             print*, 'k = ', k, ' hyp = ', hyp(k), ' rho = ', rho_init(k)
-          ENDDO
-          print*, ' '
-       ENDIF
-
     ENDDO
 
@@ -111,6 +100,6 @@
 
 !
-!-- Calculate the initial potential density, based on the initial
-!-- temperature and salinity profile
+!-- Calculate the 3d array of initial in situ and potential density,
+!-- based on the initial temperature and salinity profile
     CALL eqn_state_seawater
 
Index: palm/trunk/SOURCE/modules.f90
===================================================================
--- palm/trunk/SOURCE/modules.f90	(revision 377)
+++ palm/trunk/SOURCE/modules.f90	(revision 388)
@@ -5,4 +5,5 @@
 ! Current revisions:
 ! -----------------
+! +prho, prho_1
 ! +bc_lr_cyc, bc_ns_cyc
 ! +output_for_t0 
@@ -169,12 +170,12 @@
 
     REAL, DIMENSION(:,:,:), ALLOCATABLE, TARGET ::                             &
-          e_1, e_2, e_3, kh_1, kh_2, km_1, km_2, p, pt_1, pt_2, pt_3, q_1,     &
-          q_2, q_3, ql_1, ql_2, rho_1, sa_1, sa_2, sa_3, u_1, u_2, u_3, v_1,   &
-          v_2, v_3, vpt_1, vpt_2, w_1, w_2, w_3
+          e_1, e_2, e_3, kh_1, kh_2, km_1, km_2, p, prho_1, pt_1, pt_2, pt_3,  &
+          q_1, q_2, q_3, ql_1, ql_2, rho_1, sa_1, sa_2, sa_3, u_1, u_2, u_3,   &
+          v_1, v_2, v_3, vpt_1, vpt_2, w_1, w_2, w_3
 
     REAL, DIMENSION(:,:,:), POINTER ::                                         &
-          e, e_m, e_p, kh, kh_m, km, km_m, pt, pt_m, pt_p, q, q_m, q_p, ql,    &
-          ql_c, rho, sa, sa_p, te_m, tpt_m, tq_m, tsa_m, tu_m, tv_m, tw_m, u,  &
-          u_m, u_p, v, v_m, v_p, vpt, vpt_m, w, w_m, w_p
+          e, e_m, e_p, kh, kh_m, km, km_m, prho, pt, pt_m, pt_p, q, q_m, q_p,  &
+          ql, ql_c, rho, sa, sa_p, te_m, tpt_m, tq_m, tsa_m, tu_m, tv_m, tw_m, &
+          u, u_m, u_p, v, v_m, v_p, vpt, vpt_m, w, w_m, w_p
 
     REAL, DIMENSION(:,:,:,:), ALLOCATABLE ::  rif_wall
Index: palm/trunk/SOURCE/production_e.f90
===================================================================
--- palm/trunk/SOURCE/production_e.f90	(revision 377)
+++ palm/trunk/SOURCE/production_e.f90	(revision 388)
@@ -4,4 +4,6 @@
 ! Actual revisions:
 ! -----------------
+! Bugfix: wrong sign in buoyancy production of ocean part in case of not using
+!         the reference density (only in 3D routine production_e)
 ! Bugfix to avoid zero division by km_neutral
 !
@@ -450,6 +452,6 @@
                    DO  j = nys, nyn
                       DO  k = nzb_s_inner(j,i)+1, nzt
-                         tend(k,j,i) = tend(k,j,i) +                    &
-                                       kh(k,j,i) * g / prho_reference * &
+                         tend(k,j,i) = tend(k,j,i) +                   &
+                                       kh(k,j,i) * g / rho_reference * &
                                        ( rho(k+1,j,i)-rho(k-1,j,i) ) * dd2zu(k)
                       ENDDO
@@ -487,5 +489,5 @@
                    DO  j = nys, nyn
                       DO  k = nzb_s_inner(j,i)+1, nzt
-                         tend(k,j,i) = tend(k,j,i) -                &
+                         tend(k,j,i) = tend(k,j,i) +                &
                                        kh(k,j,i) * g / rho(k,j,i) * &
                                        ( rho(k+1,j,i)-rho(k-1,j,i) ) * dd2zu(k)
@@ -946,5 +948,5 @@
 !--             bottom and top surface layer
                 DO  k = nzb_s_inner(j,i)+1, nzt
-                   tend(k,j,i) = tend(k,j,i) + kh(k,j,i) * g / prho_reference * &
+                   tend(k,j,i) = tend(k,j,i) + kh(k,j,i) * g / rho_reference * &
                                       ( rho(k+1,j,i) - rho(k-1,j,i) ) * dd2zu(k)
                 ENDDO
Index: palm/trunk/SOURCE/prognostic_equations.f90
===================================================================
--- palm/trunk/SOURCE/prognostic_equations.f90	(revision 377)
+++ palm/trunk/SOURCE/prognostic_equations.f90	(revision 388)
@@ -4,4 +4,5 @@
 ! Actual revisions:
 ! -----------------
+! prho is used instead of rho in diffusion_e,
 ! external pressure gradient
 !
@@ -330,5 +331,5 @@
           CALL coriolis( i, j, 3 )
           IF ( ocean )  THEN
-             CALL buoyancy( i, j, rho, prho_reference, 3, 64 )
+             CALL buoyancy( i, j, rho, rho_reference, 3, 64 )
           ELSE
              IF ( .NOT. humidity )  THEN
@@ -725,10 +726,10 @@
                 IF ( .NOT. humidity )  THEN
                    IF ( ocean )  THEN
-                      CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, km,   &
-                                        l_grid, rho, prho_reference, rif, tend, &
-                                        zu, zw )
+                      CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, km,  &
+                                        l_grid, prho, prho_reference, rif,     &
+                                        tend, zu, zw )
                    ELSE
-                      CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, km, &
-                                        l_grid, pt, pt_reference, rif, tend,  &
+                      CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, km,  &
+                                        l_grid, pt, pt_reference, rif, tend,   &
                                         zu, zw )
                    ENDIF
@@ -769,5 +770,5 @@
                       IF ( ocean )  THEN
                          CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e,  &
-                                           km, l_grid, rho, prho_reference,   &
+                                           km, l_grid, prho, prho_reference,  &
                                            rif, tend, zu, zw )
                       ELSE
@@ -1013,5 +1014,5 @@
           CALL coriolis( i, j, 3 )
           IF ( ocean )  THEN
-             CALL buoyancy( i, j, rho, prho_reference, 3, 64 )
+             CALL buoyancy( i, j, rho, rho_reference, 3, 64 )
           ELSE
              IF ( .NOT. humidity )  THEN
@@ -1263,6 +1264,6 @@
                 IF ( .NOT. humidity )  THEN
                    IF ( ocean )  THEN
-                      CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, &
-                                        km, l_grid, rho, prho_reference,  &
+                      CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e,  &
+                                        km, l_grid, prho, prho_reference,  &
                                         rif, tend, zu, zw )
                    ELSE
@@ -1552,5 +1553,5 @@
     CALL coriolis( 3 )
     IF ( ocean )  THEN
-       CALL buoyancy( rho, prho_reference, 3, 64 )
+       CALL buoyancy( rho, rho_reference, 3, 64 )
     ELSE
        IF ( .NOT. humidity )  THEN
@@ -1967,6 +1968,6 @@
           IF ( .NOT. humidity )  THEN
              IF ( ocean )  THEN
-                CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, rho, &
-                                  prho_reference, rif, tend, zu, zw )
+                CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, &
+                                  prho, prho_reference, rif, tend, zu, zw )
              ELSE
                 CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, pt, &
@@ -2004,5 +2005,5 @@
                 IF ( ocean )  THEN
                    CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, &
-                                     rho, prho_reference, rif, tend, zu, zw )
+                                     prho, prho_reference, rif, tend, zu, zw )
                 ELSE
                    CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, &
Index: palm/trunk/SOURCE/time_integration.f90
===================================================================
--- palm/trunk/SOURCE/time_integration.f90	(revision 377)
+++ palm/trunk/SOURCE/time_integration.f90	(revision 388)
@@ -4,4 +4,5 @@
 ! Actual revisions:
 ! -----------------
+! Using prho instead of rho in diffusvities.
 ! Coupling with independent precursor runs.
 ! Bugfix: output of particle time series only if particle advection is switched
@@ -284,5 +285,5 @@
              IF ( .NOT. humidity ) THEN
                 IF ( ocean )  THEN
-                   CALL diffusivities( rho, prho_reference )
+                   CALL diffusivities( prho, prho_reference )
                 ELSE
                    CALL diffusivities( pt, pt_reference )
