Changeset 388


Ignore:
Timestamp:
Sep 23, 2009 9:40:33 AM (12 years ago)
Author:
raasch
Message:

in-situ AND potential density are calculated and used in the ocean version

Location:
palm/trunk
Files:
12 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/DOC/app/chapter_4.5.7.html

    r337 r388  
    11<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.0 Transitional//EN">
    22<html><head>
     3
    34<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>
    45<!--
    56@page { size: 21cm 29.7cm }
    67-->
    7 </style></head>
    8 <body style="direction: ltr;" lang="en-US"><h4 style="line-height: 100%;"><font size="4">4.5.7
     8</style></head><body style="direction: ltr;" lang="en-US"><h4 style="line-height: 100%;"><font size="4">4.5.7
    99Plots of
    1010isosurfaces, 2d cross sections and particles with dvrp</font></h4>
     
    1414system. 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>
    1515software was originally developed by the RRZN
    16 (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,
     16(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,
    1717which can be displayed via a special plugin for internet browsers. With
    1818suitable graphic hardware (e.g. NVIDIA quattro FX cards) even stereoscopic views are
     
    5050Methoden zur visuellen und haptischen tele-immersiven
    5151Exploration
    52 komplexer Volumen- und Str&ouml;mungsdaten aus
     52komplexer Volumen- und Strömungsdaten aus
    5353parallelisierten,
    5454dynamischen 3D-Simulationen" (2005-2007). </p>
     
    5757PALM software package (see chapter <a href="chapter_3.7.html">3.7</a>).
    5858To use this package, the additional option<tt style="font-family: monospace;"> </tt><font style="font-family: Courier New,Courier,monospace;" size="2">-p
    59 <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;">
     59<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;">
    6060</span>has to be given in the <b>mrun</b> call. This
    6161automatically links the
     
    7575parallel</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;
    7676&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp;
    77 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>.
     77dvrp_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>
     78</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>
     79<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>.
    7880All streams (and static scenes) are collected to one sequence (one
    7981static scene), which can be displayed with a browser without using the
     
    8486all_streams_sequence.html to a directory on your local computer and
    8587open the file all_streams_sequence.html with your browser. For this,
    86 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>" !
     88you 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>" !
    8789<hr><p style="line-height: 100%;"><br>
    8890<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
  • palm/trunk/SCRIPTS/mrun

    r377 r388  
    31873187                      export MPI_CONNECTIONS_THRESHOLD=8192
    31883188                      mpiexec_mpt -np $ii   ./a.out  $ROPTS  < runfile_atmos
     3189
     3190                          # next is test for openmp usage
     3191                      # mpiexec -n $ii -pernode  ./a.out  $ROPTS  < runfile_atmos
    31893192                   elif [[ $( echo $mpilib | cut -c1-3 ) = mva ]]
    31903193                   then
  • palm/trunk/SOURCE/CURRENT_MODIFICATIONS

    r377 r388  
    125125Errors:
    126126------
     127Bugfix: Initial hydrostatic pressure profile in case of ocean runs is now
     128calculated in 5 iteration steps. (init_ocean)
     129
     130Bugfix: wrong sign in buoyancy production of ocean part in case of not using
     131the reference density (only in 3D routine production_e) (production_e)
     132
    127133Bugfix: output of averaged 2d/3d quantities requires that an avaraging
    128134interval has been set, respective error message is included (check_parameters)
     
    183189Bugfix: initial setting of time_coupling in coupled restart runs (time_integration)
    184190
    185 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
     191advec_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
  • palm/trunk/SOURCE/check_parameters.f90

    r376 r388  
    44! Actual revisions:
    55! -----------------
     6! Check profiles fpr prho and hyp.
    67! Bugfix: output of averaged 2d/3d quantities requires that an avaraging
    78! interval has been set, respective error message is included
     
    933934
    934935!
    935 !--    If required compute the profile of leaf area density used in the plant canopy model
     936!--    If required compute the profile of leaf area density used in the plant
     937!--    canopy model
    936938       IF ( plant_canopy ) THEN
    937939       
     
    967969
    968970!
    969 !--       In case of no given leaf area density gradients, choose a vanishing gradient
     971!--       In case of no given leaf area density gradients, choose a vanishing
     972!--       gradient
    970973          IF ( lad_vertical_gradient_level(1) == -9999999.9 ) THEN
    971974             lad_vertical_gradient_level(1) = 0.0
     
    20892092
    20902093          CASE ( 'rho' )
    2091              dopr_index(i) = 64
    2092              dopr_unit(i)  = 'kg/m3'
    2093              hom(:,2,64,:) = SPREAD( zu, 2, statistic_regions+1 )
     2094             IF ( .NOT. ocean ) THEN
     2095                message_string = 'data_output_pr = ' // &
     2096                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
     2097                                 'lemented for ocean = .FALSE.'
     2098                CALL message( 'check_parameters', 'PA0091', 1, 2, 0, 6, 0 )
     2099             ELSE
     2100                dopr_index(i) = 64
     2101                dopr_unit(i)  = 'kg/m3'
     2102                hom(:,2,64,:) = SPREAD( zu, 2, statistic_regions+1 )
     2103             ENDIF
    20942104
    20952105          CASE ( 'w"sa"' )
     
    21502160                hom(:,2,70,:) = SPREAD( zu, 2, statistic_regions+1 )
    21512161             ENDIF
     2162
     2163          CASE ( 'prho' )
     2164             IF ( .NOT. ocean ) THEN
     2165                message_string = 'data_output_pr = ' // &
     2166                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
     2167                                 'lemented for ocean = .FALSE.'
     2168                CALL message( 'check_parameters', 'PA0091', 1, 2, 0, 6, 0 )
     2169             ELSE
     2170                dopr_index(i) = 71
     2171                dopr_unit(i)  = 'kg/m3'
     2172                hom(:,2,71,:) = SPREAD( zu, 2, statistic_regions+1 )
     2173             ENDIF
     2174
     2175          CASE ( 'hyp' )
     2176             dopr_index(i) = 72
     2177             dopr_unit(i)  = 'kPa'
     2178             hom(:,2,72,:) = SPREAD( zu, 2, statistic_regions+1 )
    21522179
    21532180          CASE DEFAULT
     
    29292956!-- Check pressure gradient conditions
    29302957    IF ( dp_external .AND. conserve_volume_flow )  THEN
    2931        WRITE( message_string, * )  'Both dp_external and conserve_volume_flow', &
    2932             ' are .TRUE. but one of them must be .FALSE.'
     2958       WRITE( message_string, * )  'Both dp_external and conserve_volume_flo', &
     2959            'w are .TRUE. but one of them must be .FALSE.'
    29332960       CALL message( 'check_parameters', 'PA0150', 1, 2, 0, 6, 0 )
    29342961    ENDIF
     
    29402967       ENDIF
    29412968       IF ( .NOT. ANY( dpdxy /= 0.0 ) )  THEN
    2942           WRITE( message_string, * )  'dp_external is .TRUE. but dpdxy is zero',&
    2943                ', i.e. the external pressure gradient & will not be applied'
     2969          WRITE( message_string, * )  'dp_external is .TRUE. but dpdxy is ze', &
     2970               'ro, i.e. the external pressure gradient & will not be applied'
    29442971          CALL message( 'check_parameters', 'PA0152', 0, 1, 0, 6, 0 )
    29452972       ENDIF
  • palm/trunk/SOURCE/eqn_state_seawater.f90

    r336 r388  
    44! Actual revisions:
    55! -----------------
    6 ! First constant in array den also defined as type double
     6! Potential density is additionally calculated in eqn_state_seawater,
     7! first constant in array den also defined as type double
    78!
    89! Former revisions:
     
    6970       INTEGER ::  i, j, k
    7071
    71        REAL ::  p1, p2, p3, pt1, pt2, pt3, pt4, sa1, sa15, sa2
     72       REAL ::  pden, pnom, p1, p2, p3, pt1, pt2, pt3, pt4, sa1, sa15, sa2
    7273
    7374       DO  i = nxl, nxr
     
    7677!
    7778!--             Pressure is needed in dbar
    78 !                p1 = hyp(0) * 1E-4
    79 !                p1 = 0.0
    8079                p1 = hyp(k) * 1E-4
    8180                p2 = p1 * p1
     
    9392                sa2  = sa1 * sa1
    9493
    95                 rho(k,j,i) =                                                   &
    96                         ( nom(1)           + nom(2)*pt1     + nom(3)*pt2     + &
    97                           nom(4)*pt3       + nom(5)*sa1     + nom(6)*sa1*pt1 + &
    98                           nom(7)*sa2       + nom(8)*p1      + nom(9)*p1*pt2  + &
    99                           nom(10)*p1*sa1   + nom(11)*p2     + nom(12)*p2*pt2   &
    100                         ) /                                                    &
    101                         ( den(1)           + den(2)*pt1     + den(3)*pt2     + &
    102                           den(4)*pt3       + den(5)*pt4     + den(6)*sa1     + &
    103                           den(7)*sa1*pt1   + den(8)*sa1*pt3 + den(9)*sa15    + &
    104                           den(10)*sa15*pt2 + den(11)*p1     + den(12)*p2*pt3 + &
    105                           den(13)*p3*pt1                                       &
    106                         )
     94                pnom = nom(1)           + nom(2)*pt1     + nom(3)*pt2     + &
     95                       nom(4)*pt3       + nom(5)*sa1     + nom(6)*sa1*pt1 + &
     96                       nom(7)*sa2
     97
     98                pden = den(1)           + den(2)*pt1     + den(3)*pt2     + &
     99                       den(4)*pt3       + den(5)*pt4     + den(6)*sa1     + &
     100                       den(7)*sa1*pt1   + den(8)*sa1*pt3 + den(9)*sa15    + &
     101                       den(10)*sa15*pt2
     102
     103!
     104!--             Potential density (without pressure terms)
     105                prho(k,j,i) = pnom / pden
     106
     107                pnom = pnom +             nom(8)*p1      + nom(9)*p1*pt2  + &
     108                       nom(10)*p1*sa1   + nom(11)*p2     + nom(12)*p2*pt2
     109
     110                pden = pden +             den(11)*p1     + den(12)*p2*pt3 + &
     111                       den(13)*p3*pt1
     112
     113!
     114!--             In-situ density
     115                rho(k,j,i) = pnom / pden
    107116
    108117             ENDDO
    109118!
    110119!--          Neumann conditions are assumed at bottom and top boundary
    111              rho(nzt+1,j,i)            = rho(nzt,j,i)
    112              rho(nzb_s_inner(j,i),j,i) = rho(nzb_s_inner(j,i)+1,j,i)
     120             prho(nzt+1,j,i)            = prho(nzt,j,i)
     121             prho(nzb_s_inner(j,i),j,i) = prho(nzb_s_inner(j,i)+1,j,i)
     122             rho(nzt+1,j,i)             = rho(nzt,j,i)
     123             rho(nzb_s_inner(j,i),j,i)  = rho(nzb_s_inner(j,i)+1,j,i)
     124
    113125          ENDDO
    114126       ENDDO
     
    129141       INTEGER ::  i, j, k
    130142
    131        REAL ::  p1, p2, p3, pt1, pt2, pt3, pt4, sa1, sa15, sa2
     143       REAL ::  pden, pnom, p1, p2, p3, pt1, pt2, pt3, pt4, sa1, sa15, sa2
    132144
    133145       DO  k = nzb_s_inner(j,i)+1, nzt
    134146!
    135147!--       Pressure is needed in dbar
    136 !          p1 = hyp(0) * 1E-4
    137 !          p1 = 0.0
    138148          p1 = hyp(k) * 1E-4
    139149          p2 = p1 * p1
     
    151161          sa2  = sa1 * sa1
    152162
    153           rho(k,j,i) = ( nom(1)           + nom(2)*pt1     + nom(3)*pt2     + &
    154                          nom(4)*pt3       + nom(5)*sa1     + nom(6)*sa1*pt1 + &
    155                          nom(7)*sa2       + nom(8)*p1      + nom(9)*p1*pt2  + &
    156                          nom(10)*p1*sa1   + nom(11)*p2     + nom(12)*p2*pt2   &
    157                        ) /                                                    &
    158                        ( den(1)           + den(2)*pt1     + den(3)*pt2     + &
    159                          den(4)*pt3       + den(5)*pt4     + den(6)*sa1     + &
    160                          den(7)*sa1*pt1   + den(8)*sa1*pt3 + den(9)*sa15    + &
    161                          den(10)*sa15*pt2 + den(11)*p1     + den(12)*p2*pt3 + &
    162                          den(13)*p3*pt1                                       &
    163                        )
     163          pnom = nom(1)           + nom(2)*pt1     + nom(3)*pt2     + &
     164                 nom(4)*pt3       + nom(5)*sa1     + nom(6)*sa1*pt1 + &
     165                 nom(7)*sa2
     166
     167          pden = den(1)           + den(2)*pt1     + den(3)*pt2     + &
     168                 den(4)*pt3       + den(5)*pt4     + den(6)*sa1     + &
     169                 den(7)*sa1*pt1   + den(8)*sa1*pt3 + den(9)*sa15    + &
     170                 den(10)*sa15*pt2
     171
     172!
     173!--       Potential density (without pressure terms)
     174          prho(k,j,i) = pnom / pden
     175
     176          pnom = pnom +             nom(8)*p1      + nom(9)*p1*pt2  + &
     177                 nom(10)*p1*sa1   + nom(11)*p2     + nom(12)*p2*pt2
     178          pden = pden +             den(11)*p1     + den(12)*p2*pt3 + &
     179                 den(13)*p3*pt1
     180
     181!
     182!--       In-situ density
     183          rho(k,j,i) = pnom / pden
     184
     185
    164186       ENDDO
     187
    165188!
    166189!--    Neumann conditions are assumed at bottom and top boundary
    167        rho(nzt+1,j,i)            = rho(nzt,j,i)
    168        rho(nzb_s_inner(j,i),j,i) = rho(nzb_s_inner(j,i)+1,j,i)
     190       prho(nzt+1,j,i)            = prho(nzt,j,i)
     191       prho(nzb_s_inner(j,i),j,i) = prho(nzb_s_inner(j,i)+1,j,i)
     192       rho(nzt+1,j,i)             = rho(nzt,j,i)
     193       rho(nzb_s_inner(j,i),j,i)  = rho(nzb_s_inner(j,i)+1,j,i)
    169194
    170195    END SUBROUTINE eqn_state_seawater_ij
  • palm/trunk/SOURCE/flow_statistics.f90

    r343 r388  
    44! Current revisions:
    55! -----------------
     6! Vertical profiles of potential density and hydrostatic pressure are
     7! calculated.
    68! Added missing timeseries calculation of w"q"(0), moved timeseries q* to the
    79! end.
     
    531533                IF ( humidity )  THEN
    532534                   sums_l(nzt,48,tn) = sums_l(nzt,48,tn) + &
    533                                        qswst(j,i) * rmask(j,i,sr)  ! w"q" (w"qv")
     535                                       qswst(j,i) * rmask(j,i,sr) ! w"q" (w"qv")
    534536                   IF ( cloud_physics )  THEN
    535537                      sums_l(nzt,45,tn) = sums_l(nzt,45,tn) + (           &
     
    546548                IF ( passive_scalar )  THEN
    547549                   sums_l(nzt,48,tn) = sums_l(nzt,48,tn) + &
    548                                        qswst(j,i) * rmask(j,i,sr)  ! w"q" (w"qv")
     550                                       qswst(j,i) * rmask(j,i,sr) ! w"q" (w"qv")
    549551                ENDIF
    550552             ENDIF
     
    597599                                                       rmask(j,i,sr)
    598600                   sums_l(k,64,tn) = sums_l(k,64,tn) + rho(k,j,i) * &
     601                                                       rmask(j,i,sr)
     602                   sums_l(k,71,tn) = sums_l(k,71,tn) + prho(k,j,i) * &
    599603                                                       rmask(j,i,sr)
    600604                ENDIF
     
    644648!
    645649!--    Density at top follows Neumann condition
    646        IF ( ocean )  sums_l(nzt+1,64,tn) = sums_l(nzt,64,tn)
     650       IF ( ocean )  THEN
     651          sums_l(nzt+1,64,tn) = sums_l(nzt,64,tn)
     652          sums_l(nzt+1,71,tn) = sums_l(nzt,71,tn)
     653       ENDIF
    647654
    648655!
     
    897904       hom(:,1,69,sr) = sums(:,69)     ! w"e + w"p"/rho
    898905       hom(:,1,70,sr) = sums(:,70)     ! q*2
     906       hom(:,1,71,sr) = sums(:,71)     ! prho
     907       hom(:,1,72,sr) = hyp * 1E-4     ! hyp in kPa
    899908
    900909       hom(:,1,pr_palm-1,sr) = sums(:,pr_palm-1)
     
    916925!--    is less than 1.5 times the height where the heat flux becomes negative
    917926!--    (positive) for the first time.
    918 !--    NOTE: This criterion is still capable of improving!
    919927       z_i(1) = 0.0
    920928       first = .TRUE.
  • palm/trunk/SOURCE/init_3d_model.f90

    r359 r388  
    77! Current revisions:
    88! -----------------
     9! Initialization of prho added.
    910! bugfix: correction of initial volume flow for non-flat topography
    1011! bugfix: zero initialization of arrays within buildings for 'cyclic_fill'
     
    234235       ALLOCATE( saswsb_1(nys-1:nyn+1,nxl-1:nxr+1), &
    235236                 saswst_1(nys-1:nyn+1,nxl-1:nxr+1) )
    236        ALLOCATE( rho_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
    237                  sa_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1),  &
    238                  sa_2(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1),  &
     237       ALLOCATE( prho_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1), &
     238                 rho_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1),  &
     239                 sa_1(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1),   &
     240                 sa_2(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1),   &
    239241                 sa_3(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
    240        rho => rho_1  ! routine calc_mean_profile requires density to be a
    241                      ! pointer
     242       prho => prho_1
     243       rho  => rho_1  ! routines calc_mean_profile and diffusion_e require
     244                      ! density to be apointer
    242245       IF ( humidity_remote )  THEN
    243246          ALLOCATE( qswst_remote(nys-1:nyn+1,nxl-1:nxr+1) )
     
    13281331!--    Initialize quantities needed for the ocean model
    13291332       CALL init_ocean
     1333
    13301334    ELSE
    13311335!
  • palm/trunk/SOURCE/init_ocean.f90

    r366 r388  
    44! Actual revisions:
    55! -----------------
    6 ! Bugfix: First calculation of hyp(0) changed
     6! Bugfix: Initial profiles of hydrostatic pressure and density are calculated
     7! iteratively. First calculation of hyp(0) changed.
    78!
    89! Former revisions:
     
    3233    INTEGER ::  k, n
    3334
    34     REAL    ::  sa_l, pt_l, rho_l
     35    REAL    ::  sa_l, pt_l
    3536
    3637    REAL, DIMENSION(nzb:nzt+1) ::  rho_init
     
    4546!-- Calculate initial vertical profile of hydrostatic pressure (in Pa)
    4647!-- and the reference density (used later in buoyancy term)
     48!-- First step: Calculate pressure using reference density
    4749    hyp(nzt+1) = surface_pressure * 100.0
    4850
     
    5557    hyp(0) = hyp(1) + rho_surface * g * dzu(1)
    5658
    57     IF ( myid == 0 )  THEN
    58        print*,'hydro pres using rho_surface'
    59        DO  k = nzt+1, 0, -1
    60           print*, 'k = ', k, ' hyp = ', hyp(k)
    61        ENDDO
    62        print*, ' '
    63     ENDIF
    64 
     59!
     60!-- Second step: Iteratively calculate in situ density (based on presssure)
     61!-- and pressure (based on in situ density)
    6562    DO  n = 1, 5
    6663
     
    8582       ENDDO
    8683
    87        IF ( myid == 0 )  THEN
    88           print*,'hydro pres / rho  n = ', n
    89           DO  k = nzt+1, 0, -1
    90              print*, 'k = ', k, ' hyp = ', hyp(k), ' rho = ', rho_init(k)
    91           ENDDO
    92           print*, ' '
    93        ENDIF
    94 
    9584    ENDDO
    9685
     
    111100
    112101!
    113 !-- Calculate the initial potential density, based on the initial
    114 !-- temperature and salinity profile
     102!-- Calculate the 3d array of initial in situ and potential density,
     103!-- based on the initial temperature and salinity profile
    115104    CALL eqn_state_seawater
    116105
  • palm/trunk/SOURCE/modules.f90

    r367 r388  
    55! Current revisions:
    66! -----------------
     7! +prho, prho_1
    78! +bc_lr_cyc, bc_ns_cyc
    89! +output_for_t0
     
    169170
    170171    REAL, DIMENSION(:,:,:), ALLOCATABLE, TARGET ::                             &
    171           e_1, e_2, e_3, kh_1, kh_2, km_1, km_2, p, pt_1, pt_2, pt_3, q_1,     &
    172           q_2, q_3, ql_1, ql_2, rho_1, sa_1, sa_2, sa_3, u_1, u_2, u_3, v_1,   &
    173           v_2, v_3, vpt_1, vpt_2, w_1, w_2, w_3
     172          e_1, e_2, e_3, kh_1, kh_2, km_1, km_2, p, prho_1, pt_1, pt_2, pt_3,  &
     173          q_1, q_2, q_3, ql_1, ql_2, rho_1, sa_1, sa_2, sa_3, u_1, u_2, u_3,   &
     174          v_1, v_2, v_3, vpt_1, vpt_2, w_1, w_2, w_3
    174175
    175176    REAL, DIMENSION(:,:,:), POINTER ::                                         &
    176           e, e_m, e_p, kh, kh_m, km, km_m, pt, pt_m, pt_p, q, q_m, q_p, ql,    &
    177           ql_c, rho, sa, sa_p, te_m, tpt_m, tq_m, tsa_m, tu_m, tv_m, tw_m, u, &
    178           u_m, u_p, v, v_m, v_p, vpt, vpt_m, w, w_m, w_p
     177          e, e_m, e_p, kh, kh_m, km, km_m, prho, pt, pt_m, pt_p, q, q_m, q_p,  &
     178          ql, ql_c, rho, sa, sa_p, te_m, tpt_m, tq_m, tsa_m, tu_m, tv_m, tw_m, &
     179          u, u_m, u_p, v, v_m, v_p, vpt, vpt_m, w, w_m, w_p
    179180
    180181    REAL, DIMENSION(:,:,:,:), ALLOCATABLE ::  rif_wall
  • palm/trunk/SOURCE/production_e.f90

    r364 r388  
    44! Actual revisions:
    55! -----------------
     6! Bugfix: wrong sign in buoyancy production of ocean part in case of not using
     7!         the reference density (only in 3D routine production_e)
    68! Bugfix to avoid zero division by km_neutral
    79!
     
    450452                   DO  j = nys, nyn
    451453                      DO  k = nzb_s_inner(j,i)+1, nzt
    452                          tend(k,j,i) = tend(k,j,i) +                    &
    453                                        kh(k,j,i) * g / prho_reference * &
     454                         tend(k,j,i) = tend(k,j,i) +                   &
     455                                       kh(k,j,i) * g / rho_reference * &
    454456                                       ( rho(k+1,j,i)-rho(k-1,j,i) ) * dd2zu(k)
    455457                      ENDDO
     
    487489                   DO  j = nys, nyn
    488490                      DO  k = nzb_s_inner(j,i)+1, nzt
    489                          tend(k,j,i) = tend(k,j,i) -                &
     491                         tend(k,j,i) = tend(k,j,i) +                &
    490492                                       kh(k,j,i) * g / rho(k,j,i) * &
    491493                                       ( rho(k+1,j,i)-rho(k-1,j,i) ) * dd2zu(k)
     
    946948!--             bottom and top surface layer
    947949                DO  k = nzb_s_inner(j,i)+1, nzt
    948                    tend(k,j,i) = tend(k,j,i) + kh(k,j,i) * g / prho_reference * &
     950                   tend(k,j,i) = tend(k,j,i) + kh(k,j,i) * g / rho_reference * &
    949951                                      ( rho(k+1,j,i) - rho(k-1,j,i) ) * dd2zu(k)
    950952                ENDDO
  • palm/trunk/SOURCE/prognostic_equations.f90

    r240 r388  
    44! Actual revisions:
    55! -----------------
     6! prho is used instead of rho in diffusion_e,
    67! external pressure gradient
    78!
     
    330331          CALL coriolis( i, j, 3 )
    331332          IF ( ocean )  THEN
    332              CALL buoyancy( i, j, rho, prho_reference, 3, 64 )
     333             CALL buoyancy( i, j, rho, rho_reference, 3, 64 )
    333334          ELSE
    334335             IF ( .NOT. humidity )  THEN
     
    725726                IF ( .NOT. humidity )  THEN
    726727                   IF ( ocean )  THEN
    727                       CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, km,   &
    728                                         l_grid, rho, prho_reference, rif, tend, &
    729                                         zu, zw )
     728                      CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, km,  &
     729                                        l_grid, prho, prho_reference, rif,    &
     730                                        tend, zu, zw )
    730731                   ELSE
    731                       CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, km, &
    732                                         l_grid, pt, pt_reference, rif, tend,  &
     732                      CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, km,  &
     733                                        l_grid, pt, pt_reference, rif, tend,   &
    733734                                        zu, zw )
    734735                   ENDIF
     
    769770                      IF ( ocean )  THEN
    770771                         CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e,  &
    771                                            km, l_grid, rho, prho_reference,   &
     772                                           km, l_grid, prho, prho_reference,  &
    772773                                           rif, tend, zu, zw )
    773774                      ELSE
     
    10131014          CALL coriolis( i, j, 3 )
    10141015          IF ( ocean )  THEN
    1015              CALL buoyancy( i, j, rho, prho_reference, 3, 64 )
     1016             CALL buoyancy( i, j, rho, rho_reference, 3, 64 )
    10161017          ELSE
    10171018             IF ( .NOT. humidity )  THEN
     
    12631264                IF ( .NOT. humidity )  THEN
    12641265                   IF ( ocean )  THEN
    1265                       CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, &
    1266                                         km, l_grid, rho, prho_reference,  &
     1266                      CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e,  &
     1267                                        km, l_grid, prho, prho_reference,  &
    12671268                                        rif, tend, zu, zw )
    12681269                   ELSE
     
    15521553    CALL coriolis( 3 )
    15531554    IF ( ocean )  THEN
    1554        CALL buoyancy( rho, prho_reference, 3, 64 )
     1555       CALL buoyancy( rho, rho_reference, 3, 64 )
    15551556    ELSE
    15561557       IF ( .NOT. humidity )  THEN
     
    19671968          IF ( .NOT. humidity )  THEN
    19681969             IF ( ocean )  THEN
    1969                 CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, rho, &
    1970                                   prho_reference, rif, tend, zu, zw )
     1970                CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, &
     1971                                  prho, prho_reference, rif, tend, zu, zw )
    19711972             ELSE
    19721973                CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, pt, &
     
    20042005                IF ( ocean )  THEN
    20052006                   CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, &
    2006                                      rho, prho_reference, rif, tend, zu, zw )
     2007                                     prho, prho_reference, rif, tend, zu, zw )
    20072008                ELSE
    20082009                   CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, &
  • palm/trunk/SOURCE/time_integration.f90

    r348 r388  
    44! Actual revisions:
    55! -----------------
     6! Using prho instead of rho in diffusvities.
    67! Coupling with independent precursor runs.
    78! Bugfix: output of particle time series only if particle advection is switched
     
    284285             IF ( .NOT. humidity ) THEN
    285286                IF ( ocean )  THEN
    286                    CALL diffusivities( rho, prho_reference )
     287                   CALL diffusivities( prho, prho_reference )
    287288                ELSE
    288289                   CALL diffusivities( pt, pt_reference )
Note: See TracChangeset for help on using the changeset viewer.