Changeset 97 for palm/trunk/SOURCE


Ignore:
Timestamp:
Jun 21, 2007 8:23:15 AM (17 years ago)
Author:
raasch
Message:

New:
---
ocean version including prognostic equation for salinity and equation of state for seawater. Routine buoyancy can be used with both temperature and density.
+ inipar-parameters bc_sa_t, bottom_salinityflux, ocean, sa_surface, sa_vertical_gradient, sa_vertical_gradient_level, top_salinityflux

advec_s_bc, average_3d_data, boundary_conds, buoyancy, check_parameters, data_output_2d, data_output_3d, diffusion_e, flow_statistics, header, init_grid, init_3d_model, modules, netcdf, parin, production_e, prognostic_equations, read_var_list, sum_up_3d_data, swap_timelevel, time_integration, user_interface, write_var_list, write_3d_binary

New:
eqn_state_seawater, init_ocean

Changed:


inipar-parameter use_pt_reference renamed use_reference

hydro_press renamed hyp, routine calc_mean_pt_profile renamed calc_mean_profile

format adjustments for the ocean version (run_control)

advec_particles, buoyancy, calc_liquid_water_content, check_parameters, diffusion_e, diffusivities, header, init_cloud_physics, modules, production_e, prognostic_equations, run_control

Errors:


Bugfix: height above topography instead of height above level k=0 is used for calculating the mixing length (diffusion_e and diffusivities).

Bugfix: error in boundary condition for TKE removed (advec_s_bc)

advec_s_bc, diffusion_e, prognostic_equations

Location:
palm/trunk/SOURCE
Files:
18 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SOURCE/CURRENT_MODIFICATIONS

    r96 r97  
    11New:
    22---
    3 ocean version including prognostic equation for salinity and equation of state for seawater
    4 + inipar-parameters bc_sa_t, bottom_salinityflux, ocean, sa_surface, sa_vertical_gradient, sa_vertical_gradient_level, top_salinityflux
    53
    6 average_3d_data, boundary_conds, buoyancy, check_parameters, data_output_2d, data_output_3d, diffusion_e, flow_statistics, header, init_grid, init_3d_model, modules, parin, prognostic_equations, read_var_list, sum_up_3d_data, swap_timelevel, time_integration, user_interface, write_var_list, write_3d_binary
    7 
    8 New:
    9 eqn_state_seawater, init_ocean
    104
    115Changed:
    126-------
    13 hydro_press renamed hyp, routine calc_mean_pt_profile renamed calc_mean_profile
    14 
    15 advec_particles, buoyancy, calc_liquid_water_content, init_cloud_physics, modules, prognostic_equations
    167
    178
     
    1910------
    2011
    21 Bugfix: height above topography instead of height above level k=0 is used for calculating the mixing length (diffusion_e and diffusivities).
    2212
    23 diffusion_e, prognostic_equations
    24 
    25 
    26 
    27 
    28 To be completed (adjusted to ocean version):
    29 ----------------
    30 
    31 boundary conditions for salinity (check_parameters)
    32 output of salinity profile (header)
    33 
    34 initialization of ug,vg profiles
    35 
    36 calculation of w*
    37 output of theta* and z_i (e.g. run_control file)
    38 
    39 default disturbance level
    40 
    41 Bott-Chlond scheme to be extended for salinity
    42 
    43 Upstream-spline scheme to be extended for salinity?????
  • palm/trunk/SOURCE/advec_s_bc.f90

    r77 r97  
    44! Actual revisions:
    55! -----------------
    6 !
     6! Advection of salinity included
     7! Bugfix: Error in boundary condition for TKE removed
    78!
    89! Former revisions:
     
    755756       ENDIF
    756757
     758    ELSEIF ( sk_char == 'sa' )  THEN
     759
     760!
     761!--    Salinity boundary condition at the bottom boundary.
     762!--    So far, always Neumann (i.e. here zero gradient) is used
     763       DO  i = nxl, nxr
     764          DO  j = nys, nyn
     765             sk_p(nzb-1,j,i) = sk_p(nzb,j,i)
     766             sk_p(nzb-2,j,i) = sk_p(nzb,j,i)
     767          ENDDO
     768       ENDDO
     769
     770!
     771!--    Salinity boundary condition at the top boundary.
     772!--    Dirichlet or Neumann (zero gradient)
     773       DO  i = nxl, nxr
     774          DO  j = nys, nyn
     775             sk_p(nzt+2,j,i)   = sk_p(nzt+1,j,i)
     776             sk_p(nzt+3,j,i)   = sk_p(nzt+1,j,i)
     777          ENDDO
     778       ENDDO
     779
    757780    ELSEIF ( sk_char == 'q' )  THEN
    758781
    759782!
    760 !--    Specific humidity boundary condition at the bottom boundary
    761        IF ( ibc_q_b == 0 )  THEN
    762 !
    763 !--       Dirichlet (fixed surface humidity)
    764           DO  i = nxl, nxr
    765              DO  j = nys, nyn
    766                 sk_p(nzb-1,j,i) = sk_p(nzb,j,i)
    767                 sk_p(nzb-2,j,i) = sk_p(nzb,j,i)
    768              ENDDO
    769           ENDDO
    770 
    771        ELSE
    772 !
    773 !--       Neumann (i.e. here zero gradient)
    774           DO  i = nxl, nxr
    775              DO  j = nys, nyn
    776                 sk_p(nzb-1,j,i) = sk_p(nzb,j,i)
    777                 sk_p(nzb-2,j,i) = sk_p(nzb,j,i)
    778              ENDDO
    779           ENDDO
    780 
    781        ENDIF
     783!--    Specific humidity boundary condition at the bottom boundary.
     784!--    Dirichlet (fixed surface humidity) or Neumann (i.e. zero gradient)
     785       DO  i = nxl, nxr
     786          DO  j = nys, nyn
     787             sk_p(nzb-1,j,i) = sk_p(nzb,j,i)
     788             sk_p(nzb-2,j,i) = sk_p(nzb,j,i)
     789          ENDDO
     790       ENDDO
    782791
    783792!
     
    809818!
    810819!--    TKE boundary condition at bottom and top boundary (generally Neumann)
    811        sk_p(nzb-1,j,i) = sk_p(nzb,j,i)
    812        sk_p(nzb-2,j,i) = sk_p(nzb,j,i)
    813        sk_p(nzt+2,j,i) = sk_p(nzt+1,j,i)
    814        sk_p(nzt+3,j,i) = sk_p(nzt+1,j,i)
     820       DO  i = nxl, nxr
     821          DO  j = nys, nyn
     822             sk_p(nzb-1,j,i) = sk_p(nzb,j,i)
     823             sk_p(nzb-2,j,i) = sk_p(nzb,j,i)
     824             sk_p(nzt+2,j,i) = sk_p(nzt+1,j,i)
     825             sk_p(nzt+3,j,i) = sk_p(nzt+1,j,i)
     826          ENDDO
     827       ENDDO
    815828
    816829    ELSE
  • palm/trunk/SOURCE/buoyancy.f90

    r96 r97  
    44! Actual revisions:
    55! -----------------
     6! Routine reneralized to be used with temperature AND density:
     7! argument theta renamed var, new argument var_reference,
     8! use_pt_reference renamed use_reference,
    69! calc_mean_pt_profile renamed calc_mean_profile
    710!
     
    4649! Call for all grid points
    4750!------------------------------------------------------------------------------!
    48     SUBROUTINE buoyancy( theta, wind_component, pr )
     51    SUBROUTINE buoyancy( var, var_reference, wind_component, pr )
    4952
    5053       USE arrays_3d
     
    5760
    5861       INTEGER ::  i, j, k, pr, wind_component
    59        REAL, DIMENSION(:,:,:), POINTER  ::  theta
     62       REAL    ::  var_reference
     63       REAL, DIMENSION(:,:,:), POINTER  ::  var
    6064
    6165
     
    6367!
    6468!--       Normal case: horizontal surface
    65           IF ( use_pt_reference )  THEN
    66              DO  i = nxl, nxr
    67                 DO  j = nys, nyn
    68                    DO  k = nzb_s_inner(j,i)+1, nzt-1
    69                       tend(k,j,i) = tend(k,j,i) + g * 0.5 * (                 &
    70                         ( theta(k,j,i)   - hom(k,1,pr,0)   ) / pt_reference + &
    71                         ( theta(k+1,j,i) - hom(k+1,1,pr,0) ) / pt_reference   &
     69          IF ( use_reference )  THEN
     70             DO  i = nxl, nxr
     71                DO  j = nys, nyn
     72                   DO  k = nzb_s_inner(j,i)+1, nzt-1
     73                      tend(k,j,i) = tend(k,j,i) + atmos_ocean_sign * g * 0.5 * &
     74                                                            (                  &
     75                          ( var(k,j,i)   - hom(k,1,pr,0)   ) / var_reference + &
     76                          ( var(k+1,j,i) - hom(k+1,1,pr,0) ) / var_reference   &
    7277                                                            )
    7378                   ENDDO
     
    7883                DO  j = nys, nyn
    7984                   DO  k = nzb_s_inner(j,i)+1, nzt-1
    80                       tend(k,j,i) = tend(k,j,i) + g * 0.5 * (                  &
    81                         ( theta(k,j,i)   - hom(k,1,pr,0)   ) / hom(k,1,pr,0) + &
    82                         ( theta(k+1,j,i) - hom(k+1,1,pr,0) ) / hom(k+1,1,pr,0) &
     85                      tend(k,j,i) = tend(k,j,i) + atmos_ocean_sign * g * 0.5 * &
     86                                                            (                  &
     87                          ( var(k,j,i)   - hom(k,1,pr,0)   ) / hom(k,1,pr,0) + &
     88                          ( var(k+1,j,i) - hom(k+1,1,pr,0) ) / hom(k+1,1,pr,0) &
    8389                                                            )
    8490                   ENDDO
     
    136142! Call for grid point i,j
    137143!------------------------------------------------------------------------------!
    138     SUBROUTINE buoyancy_ij( i, j, theta, wind_component, pr )
     144    SUBROUTINE buoyancy_ij( i, j, var, var_reference, wind_component, pr )
    139145
    140146       USE arrays_3d
     
    147153
    148154       INTEGER ::  i, j, k, pr, wind_component
    149        REAL, DIMENSION(:,:,:), POINTER  ::  theta
     155       REAL    ::  var_reference
     156       REAL, DIMENSION(:,:,:), POINTER  ::  var
    150157
    151158
     
    153160!
    154161!--       Normal case: horizontal surface
    155           IF ( use_pt_reference )  THEN
    156              DO  k = nzb_s_inner(j,i)+1, nzt-1
    157                  tend(k,j,i) = tend(k,j,i) + g * 0.5 * (                      &
    158                         ( theta(k,j,i)   - hom(k,1,pr,0)   ) / pt_reference + &
    159                         ( theta(k+1,j,i) - hom(k+1,1,pr,0) ) / pt_reference   &
    160                                                        )
    161              ENDDO
    162           ELSE
    163              DO  k = nzb_s_inner(j,i)+1, nzt-1
    164                  tend(k,j,i) = tend(k,j,i) + g * 0.5 * (                       &
    165                         ( theta(k,j,i)   - hom(k,1,pr,0)   ) / hom(k,1,pr,0) + &
    166                         ( theta(k+1,j,i) - hom(k+1,1,pr,0) ) / hom(k+1,1,pr,0) &
    167                                                        )
     162          IF ( use_reference )  THEN
     163             DO  k = nzb_s_inner(j,i)+1, nzt-1
     164                 tend(k,j,i) = tend(k,j,i) + atmos_ocean_sign * g * 0.5 * (   &
     165                         ( var(k,j,i)   - hom(k,1,pr,0)   ) / var_reference + &
     166                         ( var(k+1,j,i) - hom(k+1,1,pr,0) ) / var_reference   &
     167                                                                          )
     168             ENDDO
     169          ELSE
     170             DO  k = nzb_s_inner(j,i)+1, nzt-1
     171                 tend(k,j,i) = tend(k,j,i) + atmos_ocean_sign * g * 0.5 * (    &
     172                          ( var(k,j,i)   - hom(k,1,pr,0)   ) / hom(k,1,pr,0) + &
     173                          ( var(k+1,j,i) - hom(k+1,1,pr,0) ) / hom(k+1,1,pr,0) &
     174                                                                          )
    168175             ENDDO
    169176          ENDIF
  • palm/trunk/SOURCE/check_parameters.f90

    r96 r97  
    44! Actual revisions:
    55! -----------------
    6 ! Initial salinity profile is calculated.
     6! Initial salinity profile is calculated, salinity boundary conditions are
     7! checked,
    78! z_max_do1d is checked only in case of ocean = .f.,
    8 ! +initial temperature profile for the ocean version,
     9! +initial temperature and geostrophic velocity profiles for the ocean version,
     10! use_pt_reference renamed use_reference
    911!
    1012! Former revisions:
     
    163165          WRITE( action, '(A,A)' )  'timestep_scheme = ', timestep_scheme
    164166       ENDIF
     167       IF ( momentum_advec == 'ups-scheme' )  THEN
     168          WRITE( action, '(A,A)' )  'momentum_advec = ', momentum_advec
     169       ENDIF
    165170       IF ( action /= ' ' )  THEN
    166171          IF ( myid == 0 )  THEN
     
    447452       i = 1
    448453       gradient = 0.0
    449        ug_vertical_gradient_level_ind(1) = 0
    450        ug(0) = ug_surface
    451        DO  k = 1, nzt+1
    452           IF ( ug_vertical_gradient_level(i) < zu(k)  .AND. &
    453                ug_vertical_gradient_level(i) >= 0.0 )  THEN
    454              gradient = ug_vertical_gradient(i) / 100.0
    455              ug_vertical_gradient_level_ind(i) = k - 1
    456              i = i + 1
    457              IF ( i > 10 )  THEN
    458                 IF ( myid == 0 )  THEN
    459                    PRINT*, '+++ check_parameters: upper bound 10 of array', &
    460                            ' "ug_vertical_gradient_level_ind" exceeded'
    461                 ENDIF
    462                 CALL local_stop
    463              ENDIF
    464           ENDIF
    465           IF ( gradient /= 0.0 )  THEN
    466              IF ( k /= 1 )  THEN
    467                 ug(k) = ug(k-1) + dzu(k) * gradient
     454
     455       IF ( .NOT. ocean )  THEN
     456
     457          ug_vertical_gradient_level_ind(1) = 0
     458          ug(0) = ug_surface
     459          DO  k = 1, nzt+1
     460             IF ( ug_vertical_gradient_level(i) < zu(k)  .AND. &
     461                  ug_vertical_gradient_level(i) >= 0.0 )  THEN
     462                gradient = ug_vertical_gradient(i) / 100.0
     463                ug_vertical_gradient_level_ind(i) = k - 1
     464                i = i + 1
     465                IF ( i > 10 )  THEN
     466                   IF ( myid == 0 )  THEN
     467                      PRINT*, '+++ check_parameters: upper bound 10 of array', &
     468                              ' "ug_vertical_gradient_level_ind" exceeded'
     469                   ENDIF
     470                   CALL local_stop
     471                ENDIF
     472             ENDIF
     473             IF ( gradient /= 0.0 )  THEN
     474                IF ( k /= 1 )  THEN
     475                   ug(k) = ug(k-1) + dzu(k) * gradient
     476                ELSE
     477                   ug(k) = ug_surface + 0.5 * dzu(k) * gradient
     478                ENDIF
    468479             ELSE
    469                 ug(k) = ug_surface + 0.5 * dzu(k) * gradient
    470              ENDIF
    471           ELSE
    472              ug(k) = ug(k-1)
    473           ENDIF
    474        ENDDO
     480                ug(k) = ug(k-1)
     481             ENDIF
     482          ENDDO
     483
     484       ELSE
     485
     486          ug_vertical_gradient_level_ind(1) = nzt+1
     487          DO  k = nzt, 0, -1
     488             IF ( ug_vertical_gradient_level(i) > zu(k)  .AND. &
     489                  ug_vertical_gradient_level(i) <= 0.0 )  THEN
     490                gradient = ug_vertical_gradient(i) / 100.0
     491                ug_vertical_gradient_level_ind(i) = k + 1
     492                i = i + 1
     493                IF ( i > 10 )  THEN
     494                   IF ( myid == 0 )  THEN
     495                      PRINT*, '+++ check_parameters: upper bound 10 of array', &
     496                              ' "ug_vertical_gradient_level_ind" exceeded'
     497                   ENDIF
     498                   CALL local_stop
     499                ENDIF
     500             ENDIF
     501             IF ( gradient /= 0.0 )  THEN
     502                IF ( k /= nzt )  THEN
     503                   ug(k) = ug(k+1) - dzu(k+1) * gradient
     504                ELSE
     505                   ug(k)   = ug_surface - 0.5 * dzu(k+1) * gradient
     506                   ug(k+1) = ug_surface + 0.5 * dzu(k+1) * gradient
     507                ENDIF
     508             ELSE
     509                ug(k) = ug(k+1)
     510             ENDIF
     511          ENDDO
     512
     513       ENDIF
    475514
    476515       u_init = ug
     
    478517!
    479518!--    In case of no given gradients for ug, choose a vanishing gradient
    480        IF ( ug_vertical_gradient_level(1) == -1.0 )  THEN
     519       IF ( ug_vertical_gradient_level(1) == -9999999.9 )  THEN
    481520          ug_vertical_gradient_level(1) = 0.0
    482521       ENDIF 
     
    488527       i = 1
    489528       gradient = 0.0
    490        vg_vertical_gradient_level_ind(1) = 0
    491        vg(0) = vg_surface
    492        DO  k = 1, nzt+1
    493           IF ( vg_vertical_gradient_level(i) < zu(k)  .AND. &
    494                vg_vertical_gradient_level(i) >= 0.0 )  THEN
    495              gradient = vg_vertical_gradient(i) / 100.0
    496              vg_vertical_gradient_level_ind(i) = k - 1
    497              i = i + 1
    498              IF ( i > 10 )  THEN
    499                 IF ( myid == 0 )  THEN
    500                    PRINT*, '+++ check_parameters: upper bound 10 of array', &
    501                            ' "vg_vertical_gradient_level_ind" exceeded'
    502                 ENDIF
    503                 CALL local_stop
    504              ENDIF
    505           ENDIF
    506           IF ( gradient /= 0.0 )  THEN
    507              IF ( k /= 1 )  THEN
    508                 vg(k) = vg(k-1) + dzu(k) * gradient
     529
     530       IF ( .NOT. ocean )  THEN
     531
     532          vg_vertical_gradient_level_ind(1) = 0
     533          vg(0) = vg_surface
     534          DO  k = 1, nzt+1
     535             IF ( vg_vertical_gradient_level(i) < zu(k)  .AND. &
     536                  vg_vertical_gradient_level(i) >= 0.0 )  THEN
     537                gradient = vg_vertical_gradient(i) / 100.0
     538                vg_vertical_gradient_level_ind(i) = k - 1
     539                i = i + 1
     540                IF ( i > 10 )  THEN
     541                   IF ( myid == 0 )  THEN
     542                      PRINT*, '+++ check_parameters: upper bound 10 of array', &
     543                              ' "vg_vertical_gradient_level_ind" exceeded'
     544                   ENDIF
     545                   CALL local_stop
     546                ENDIF
     547             ENDIF
     548             IF ( gradient /= 0.0 )  THEN
     549                IF ( k /= 1 )  THEN
     550                   vg(k) = vg(k-1) + dzu(k) * gradient
     551                ELSE
     552                   vg(k) = vg_surface + 0.5 * dzu(k) * gradient
     553                ENDIF
    509554             ELSE
    510                 vg(k) = vg_surface + 0.5 * dzu(k) * gradient
    511              ENDIF
    512           ELSE
    513              vg(k) = vg(k-1)
    514           ENDIF
    515        ENDDO
     555                vg(k) = vg(k-1)
     556             ENDIF
     557          ENDDO
     558
     559       ELSE
     560
     561          vg_vertical_gradient_level_ind(1) = 0
     562          DO  k = nzt, 0, -1
     563             IF ( vg_vertical_gradient_level(i) > zu(k)  .AND. &
     564                  vg_vertical_gradient_level(i) <= 0.0 )  THEN
     565                gradient = vg_vertical_gradient(i) / 100.0
     566                vg_vertical_gradient_level_ind(i) = k + 1
     567                i = i + 1
     568                IF ( i > 10 )  THEN
     569                   IF ( myid == 0 )  THEN
     570                      PRINT*, '+++ check_parameters: upper bound 10 of array', &
     571                              ' "vg_vertical_gradient_level_ind" exceeded'
     572                   ENDIF
     573                   CALL local_stop
     574                ENDIF
     575             ENDIF
     576             IF ( gradient /= 0.0 )  THEN
     577                IF ( k /= nzt )  THEN
     578                   vg(k) = vg(k+1) - dzu(k+1) * gradient
     579                ELSE
     580                   vg(k)   = vg_surface - 0.5 * dzu(k+1) * gradient
     581                   vg(k+1) = vg_surface + 0.5 * dzu(k+1) * gradient
     582                ENDIF
     583             ELSE
     584                vg(k) = vg(k+1)
     585             ENDIF
     586          ENDDO
     587
     588       ENDIF
    516589
    517590       v_init = vg
     
    519592!
    520593!--    In case of no given gradients for vg, choose a vanishing gradient
    521        IF ( vg_vertical_gradient_level(1) == -1.0 )  THEN
     594       IF ( vg_vertical_gradient_level(1) == -9999999.9 )  THEN
    522595          vg_vertical_gradient_level(1) = 0.0
    523596       ENDIF
     
    707780
    708781!
    709 !-- Reference temperature to be used in buoyancy terms
    710     IF ( pt_reference /= 9999999.9 )  use_pt_reference = .TRUE.
     782!-- Ocean runs always use reference values in the buoyancy term. Therefore
     783!-- set the reference temperature equal to the surface temperature.
     784    IF ( ocean  .AND.  pt_reference == 9999999.9 )  pt_reference = pt_surface
     785
     786!
     787!-- Reference value has to be used in buoyancy terms
     788    IF ( pt_reference /= 9999999.9 )  use_reference = .TRUE.
     789
     790!
     791!-- Sign of buoyancy/stability terms
     792    IF ( ocean )  atmos_ocean_sign = -1.0
     793
     794!
     795!-- Ocean version is using flux boundary conditions at the top
     796    IF ( ocean )  use_top_fluxes = .TRUE.
    711797
    712798!
     
    10201106
    10211107       IF ( top_salinityflux == 9999999.9 )  constant_top_salinityflux = .FALSE.
     1108       IF ( ibc_sa_t == 1  .AND.   top_salinityflux == 9999999.9 )  THEN
     1109          IF ( myid == 0 )  THEN
     1110             PRINT*, '+++ check_parameters:'
     1111             PRINT*, '    boundary_condition: bc_sa_t = ', bc_sa_t
     1112             PRINT*, '    requires to set top_salinityflux '
     1113          ENDIF
     1114          CALL local_stop
     1115       ENDIF
    10221116
    10231117!
     
    24742568!
    24752569!-- Determine upper and lower hight level indices for random perturbations
    2476     IF ( disturbance_level_b == -1.0 )  THEN
    2477        disturbance_level_b     = zu(nzb+3)
    2478        disturbance_level_ind_b = nzb + 3
     2570    IF ( disturbance_level_b == -9999999.9 )  THEN
     2571       IF ( ocean ) THEN
     2572          disturbance_level_b     = zu((nzt*2)/3)
     2573          disturbance_level_ind_b = ( nzt * 2 ) / 3
     2574       ELSE
     2575          disturbance_level_b     = zu(nzb+3)
     2576          disturbance_level_ind_b = nzb + 3
     2577       ENDIF
    24792578    ELSEIF ( disturbance_level_b < zu(3) )  THEN
    24802579       IF ( myid == 0 )  PRINT*, '+++ check_parameters:  disturbance_level_b=',&
     
    24962595    ENDIF
    24972596
    2498     IF ( disturbance_level_t == -1.0 )  THEN
    2499        disturbance_level_t     = zu(nzt/3)
    2500        disturbance_level_ind_t = nzt / 3
     2597    IF ( disturbance_level_t == -9999999.9 )  THEN
     2598       IF ( ocean )  THEN
     2599          disturbance_level_t     = zu(nzt-3)
     2600          disturbance_level_ind_t = nzt - 3
     2601       ELSE
     2602          disturbance_level_t     = zu(nzt/3)
     2603          disturbance_level_ind_t = nzt / 3
     2604       ENDIF
    25012605    ELSEIF ( disturbance_level_t > zu(nzt-2) )  THEN
    25022606       IF ( myid == 0 )  PRINT*, '+++ check_parameters:  disturbance_level_t=',&
  • palm/trunk/SOURCE/diffusion_e.f90

    r94 r97  
    88! This is also a bugfix, because the height above the topography is now
    99! used instead of the height above level k=0.
     10! theta renamed var, dpt_dz renamed dvar_dz, +new argument var_reference
     11! use_pt_reference renamed use_reference
    1012!
    1113! Former revisions:
     
    5052! Call for all grid points
    5153!------------------------------------------------------------------------------!
    52     SUBROUTINE diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, theta, &
    53                             rif, tend, zu, zw )
     54    SUBROUTINE diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, var, &
     55                            var_reference, rif, tend, zu, zw )
    5456
    5557       USE control_parameters
     
    6163
    6264       INTEGER ::  i, j, k
    63        REAL            ::  dpt_dz, l_stable, phi_m
     65       REAL            ::  dvar_dz, l_stable, phi_m, var_reference
    6466       REAL            ::  ddzu(1:nzt+1), dd2zu(1:nzt), ddzw(1:nzt+1), &
    6567                           l_grid(1:nzt), zu(0:nzt+1), zw(0:nzt+1)
    6668       REAL, DIMENSION(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) :: diss, tend
    6769       REAL, DIMENSION(:,:), POINTER   ::  rif
    68        REAL, DIMENSION(:,:,:), POINTER ::  e, km, theta
     70       REAL, DIMENSION(:,:,:), POINTER ::  e, km, var
    6971       REAL, DIMENSION(nzb+1:nzt,nys:nyn) ::  dissipation, l, ll
    7072 
     
    7375!--    This if clause must be outside the k-loop because otherwise
    7476!--    runtime errors occur with -C hopt on NEC
    75        IF ( use_pt_reference )  THEN
     77       IF ( use_reference )  THEN
    7678
    7779          DO  i = nxl, nxr
     
    9193!
    9294!--                Calculate the mixing length (for dissipation)
    93                    dpt_dz = ( theta(k+1,j,i) - theta(k-1,j,i) ) * dd2zu(k)
    94                    IF ( dpt_dz > 0.0 ) THEN
     95                   dvar_dz = atmos_ocean_sign * &
     96                             ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
     97                   IF ( dvar_dz > 0.0 ) THEN
    9598                      l_stable = 0.76 * SQRT( e(k,j,i) ) / &
    96                                         SQRT( g / pt_reference * dpt_dz ) + 1E-5
     99                                 SQRT( g / var_reference * dvar_dz ) + 1E-5
    97100                   ELSE
    98101                      l_stable = l_grid(k)
     
    180183!
    181184!--                Calculate the mixing length (for dissipation)
    182                    dpt_dz = ( theta(k+1,j,i) - theta(k-1,j,i) ) * dd2zu(k)
    183                    IF ( dpt_dz > 0.0 ) THEN
     185                   dvar_dz = atmos_ocean_sign * &
     186                             ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
     187                   IF ( dvar_dz > 0.0 ) THEN
    184188                      l_stable = 0.76 * SQRT( e(k,j,i) ) / &
    185                                         SQRT( g / theta(k,j,i) * dpt_dz ) + 1E-5
     189                                        SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5
    186190                   ELSE
    187191                      l_stable = l_grid(k)
     
    270274!------------------------------------------------------------------------------!
    271275    SUBROUTINE diffusion_e_ij( i, j, ddzu, dd2zu, ddzw, diss, e, km, l_grid, &
    272                                theta, rif, tend, zu, zw )
     276                               var, var_reference, rif, tend, zu, zw )
    273277
    274278       USE control_parameters
     
    280284
    281285       INTEGER         ::  i, j, k
    282        REAL            ::  dpt_dz, l_stable, phi_m
     286       REAL            ::  dvar_dz, l_stable, phi_m, var_reference
    283287       REAL            ::  ddzu(1:nzt+1), dd2zu(1:nzt), ddzw(1:nzt+1), &
    284288                           l_grid(1:nzt), zu(0:nzt+1), zw(0:nzt+1)
    285289       REAL, DIMENSION(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) ::  diss, tend
    286290       REAL, DIMENSION(:,:), POINTER   ::  rif
    287        REAL, DIMENSION(:,:,:), POINTER ::  e, km, theta
     291       REAL, DIMENSION(:,:,:), POINTER ::  e, km, var
    288292       REAL, DIMENSION(nzb+1:nzt)    ::  dissipation, l, ll
    289293
     
    303307!--    Calculate the mixing length (for dissipation)
    304308       DO  k = nzb_s_inner(j,i)+1, nzt
    305           dpt_dz = ( theta(k+1,j,i) - theta(k-1,j,i) ) * dd2zu(k)
    306           IF ( dpt_dz > 0.0 ) THEN
    307              IF ( use_pt_reference )  THEN
     309          dvar_dz = atmos_ocean_sign * &
     310                    ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
     311          IF ( dvar_dz > 0.0 ) THEN
     312             IF ( use_reference )  THEN
    308313                l_stable = 0.76 * SQRT( e(k,j,i) ) / &
    309                                   SQRT( g / pt_reference * dpt_dz ) + 1E-5
     314                                  SQRT( g / var_reference * dvar_dz ) + 1E-5
    310315             ELSE
    311316                l_stable = 0.76 * SQRT( e(k,j,i) ) / &
    312                                   SQRT( g / theta(k,j,i) * dpt_dz ) + 1E-5
     317                                  SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5
    313318             ENDIF
    314319          ELSE
  • palm/trunk/SOURCE/diffusivities.f90

    r94 r97  
    1  SUBROUTINE diffusivities( theta )
     1 SUBROUTINE diffusivities( var, var_reference )
    22
    33!------------------------------------------------------------------------------!
     
    77! This is also a bugfix, because the height above the topography is now
    88! used instead of the height above level k=0.
     9! theta renamed var, dpt_dz renamed dvar_dz, +new argument var_reference
     10! use_pt_reference renamed use_reference
    911!
    1012! Former revisions:
     
    4143    INTEGER ::  i, j, k, omp_get_thread_num, sr, tn
    4244
    43     REAL    ::  dpt_dz, l_stable, phi_m = 1.0
     45    REAL    ::  dvar_dz, l_stable, var_reference
    4446
    45     REAL    ::  theta(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1)
     47    REAL, SAVE ::  phi_m = 1.0
     48
     49    REAL    ::  var(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1)
    4650
    4751    REAL, DIMENSION(1:nzt) ::  l, ll, sqrt_e
     
    5862!
    5963!-- Compute the turbulent diffusion coefficient for momentum
    60     !$OMP PARALLEL PRIVATE (dpt_dz,i,j,k,l,ll,l_stable,phi_m,sqrt_e,sr,tn)
     64    !$OMP PARALLEL PRIVATE (dvar_dz,i,j,k,l,ll,l_stable,phi_m,sqrt_e,sr,tn)
    6165!$  tn = omp_get_thread_num()
    6266
     
    9498!--       Determine the mixing length
    9599          DO  k = nzb_s_inner(j,i)+1, nzt
    96              dpt_dz = ( theta(k+1,j,i) - theta(k-1,j,i) ) * dd2zu(k)
    97              IF ( dpt_dz > 0.0 ) THEN
    98                 IF ( use_pt_reference )  THEN
     100             dvar_dz = atmos_ocean_sign * &  ! inverse effect of pt/rho gradient
     101                       ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
     102             IF ( dvar_dz > 0.0 ) THEN
     103                IF ( use_reference )  THEN
    99104                   l_stable = 0.76 * sqrt_e(k) / &
    100                                      SQRT( g / pt_reference * dpt_dz ) + 1E-5
     105                                     SQRT( g / var_reference * dvar_dz ) + 1E-5
    101106                ELSE
    102107                   l_stable = 0.76 * sqrt_e(k) / &
    103                                      SQRT( g / theta(k,j,i) * dpt_dz ) + 1E-5
     108                                     SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5
    104109                ENDIF
    105110             ELSE
  • palm/trunk/SOURCE/eqn_state_seawater.f90

    r96 r97  
    88! Former revisions:
    99! -----------------
    10 ! $Id:$
     10! $Id$
    1111!
    1212! Initial revision
     
    1818! salinity, and pressure.
    1919! For coefficients see Jackett et al., 2006: J. Atm. Ocean Tech.
     20! eqn_state_seawater calculates the potential density referred at hyp(0).
     21! eqn_state_seawater_func calculates density.
    2022!------------------------------------------------------------------------------!
    2123
     
    7072       DO  i = nxl, nxr
    7173          DO  j = nys, nyn
    72              DO  k = nzb_u_inner(j,i)+1, nzt
     74             DO  k = nzb_s_inner(j,i)+1, nzt
    7375!
    7476!--             Pressure is needed in dbar
    75                 p1 = hyp(k) * 1E-4
     77                p1 = hyp(0) * 1E-4
    7678                p2 = p1 * p1
    7779                p3 = p2 * p1
     
    102104
    103105             ENDDO
     106!
     107!--          Neumann conditions are assumed at bottom and top boundary
     108             rho(nzt+1,j,i)            = rho(nzt,j,i)
     109             rho(nzb_s_inner(j,i),j,i) = rho(nzb_s_inner(j,i)+1,j,i)
    104110          ENDDO
    105111       ENDDO
     
    122128       REAL ::  p1, p2, p3, pt1, pt2, pt3, pt4, sa1, sa15, sa2
    123129
    124        DO  k = nzb_u_inner(j,i)+1, nzt
     130       DO  k = nzb_s_inner(j,i)+1, nzt
    125131!
    126132!--       Pressure is needed in dbar
    127           p1 = hyp(k) * 1E-4
     133          p1 = hyp(0) * 1E-4
    128134          p2 = p1 * p1
    129135          p3 = p2 * p1
     
    152158                       )
    153159       ENDDO
     160!
     161!--    Neumann conditions are assumed at bottom and top boundary
     162       rho(nzt+1,j,i)            = rho(nzt,j,i)
     163       rho(nzb_s_inner(j,i),j,i) = rho(nzb_s_inner(j,i)+1,j,i)
    154164
    155165    END SUBROUTINE eqn_state_seawater_ij
  • palm/trunk/SOURCE/flow_statistics.f90

    r96 r97  
    55! -----------------
    66! Statistics for ocean version (salinity, density) added,
    7 ! calculation of Deardorff velocity scale adjusted to be used with the ocean
    8 ! version (HAS STILL TO BE COMPLETED!!!)
     7! calculation of z_i and Deardorff velocity scale adjusted to be used with
     8! the ocean version
    99!
    1010! Former revisions:
     
    558558!
    559559!--             Salinity flux and density (density does not belong to here,
    560 !--             but so far there is no suitable place to calculate)
     560!--             but so far there is no other suitable place to calculate)
    561561                IF ( ocean )  THEN
    562562                   pts = 0.5 * ( sa(k,j,i)   - hom(k,1,23,sr) + &
     
    611611
    612612!
     613!--    Density at top follows Neumann condition
     614       IF ( ocean )  sums_l(nzt+1,64,tn) = sums_l(nzt,64,tn)
     615
     616!
    613617!--    Divergence of vertical flux of resolved scale energy and pressure
    614618!--    fluctuations. First calculate the products, then the divergence.
     
    709713!
    710714!--       Fluxes at the surface must be zero (e.g. due to the Prandtl-layer)
    711           sums(nzb,58) = 0.0
    712           sums(nzb,59) = 0.0
    713           sums(nzb,60) = 0.0
    714           sums(nzb,61) = 0.0
    715           sums(nzb,62) = 0.0
    716           sums(nzb,63) = 0.0
     715          sums_l(nzb,58,tn) = 0.0
     716          sums_l(nzb,59,tn) = 0.0
     717          sums_l(nzb,60,tn) = 0.0
     718          sums_l(nzb,61,tn) = 0.0
     719          sums_l(nzb,62,tn) = 0.0
     720          sums_l(nzb,63,tn) = 0.0
    717721
    718722       ENDIF
     
    857861       z_i(1) = 0.0
    858862       first = .TRUE.
    859 !       IF ( .NOT. ocean )  THEN
     863       IF ( ocean )  THEN
     864          DO  k = nzt, nzb+1, -1
     865             IF ( first .AND. hom(k,1,18,sr) < 0.0 )  THEN
     866                first = .FALSE.
     867                height = zw(k)
     868             ENDIF
     869             IF ( hom(k,1,18,sr) < 0.0  .AND. &
     870                  hom(k-1,1,18,sr) > hom(k,1,18,sr) )  THEN
     871                IF ( zw(k) < 1.5 * height )  THEN
     872                   z_i(1) = zw(k)
     873                ELSE
     874                   z_i(1) = height
     875                ENDIF
     876                EXIT
     877             ENDIF
     878          ENDDO
     879       ELSE
    860880          DO  k = nzb, nzt-1
    861881             IF ( first .AND. hom(k,1,18,sr) < 0.0 )  THEN
     
    873893             ENDIF
    874894          ENDDO
    875 !       ELSE
    876 !       ENDIF
    877 
    878 !
    879 !--    Second scheme: Starting from the top model boundary, look for the first
    880 !--    characteristic kink in the temperature profile, where the originally
    881 !--    stable stratification notably weakens.
     895       ENDIF
     896
     897!
     898!--    Second scheme: Starting from the top/bottom model boundary, look for
     899!--    the first characteristic kink in the temperature profile, where the
     900!--    originally stable stratification notably weakens.
    882901       z_i(2) = 0.0
    883        DO  k = nzt-1, nzb+1, -1
    884           IF ( ( hom(k+1,1,4,sr) - hom(k,1,4,sr) ) > &
    885                2.0 * ( hom(k,1,4,sr) - hom(k-1,1,4,sr) ) )  THEN
    886              z_i(2) = zu(k)
    887              EXIT
    888           ENDIF
    889        ENDDO
     902       IF ( ocean )  THEN
     903          DO  k = nzb+1, nzt-1
     904             IF ( ( hom(k,1,4,sr) - hom(k-1,1,4,sr) ) > &
     905                  2.0 * ( hom(k+1,1,4,sr) - hom(k,1,4,sr) ) )  THEN
     906                z_i(2) = zu(k)
     907                EXIT
     908             ENDIF
     909          ENDDO
     910       ELSE
     911          DO  k = nzt-1, nzb+1, -1
     912             IF ( ( hom(k+1,1,4,sr) - hom(k,1,4,sr) ) > &
     913                  2.0 * ( hom(k,1,4,sr) - hom(k-1,1,4,sr) ) )  THEN
     914                z_i(2) = zu(k)
     915                EXIT
     916             ENDIF
     917          ENDDO
     918       ENDIF
    890919
    891920       hom(nzb+6,1,pr_palm,sr) = z_i(1)
  • palm/trunk/SOURCE/header.f90

    r94 r97  
    44! Actual revisions:
    55! -----------------
    6 ! Output of model height, stretch level, Prandtl-layer height and initial
    7 ! temperature profile adjusted to be used also with the ocean version.
     6! Adjustments for the ocean version.
     7! use_pt_reference renamed use_reference
    88!
    99! Former revisions:
     
    103103          PRINT*,'+++ header:  unknown action(s): ',initializing_actions
    104104       ENDIF
     105    ENDIF
     106    IF ( ocean )  THEN
     107       run_classification = 'ocean - ' // run_classification
     108    ELSE
     109       run_classification = 'atmosphere - ' // run_classification
    105110    ENDIF
    106111
     
    391396       roben = 'e(nzt+1) = e(nzt) = e(nzt-1)'
    392397
    393        WRITE ( io, 301 )  runten, roben       
    394 
    395     ENDIF
    396 
    397     IF ( humidity  .OR.  passive_scalar )  THEN
    398        IF ( humidity )  THEN
    399           IF ( ibc_q_b == 0 )  THEN
    400              runten = 'q(0)     = q_surface'
    401           ELSE
    402              runten = 'q(0)     = q(1)'
    403           ENDIF
    404           IF ( ibc_q_t == 0 )  THEN
    405              roben =  'q(nzt)   = q_top'
    406           ELSE
    407              roben =  'q(nzt)   = q(nzt-1) + dq/dz'
    408           ENDIF
    409        ELSE
    410           IF ( ibc_q_b == 0 )  THEN
    411              runten = 's(0)     = s_surface'
    412           ELSE
    413              runten = 's(0)     = s(1)'
    414           ENDIF
    415           IF ( ibc_q_t == 0 )  THEN
    416              roben =  's(nzt)   = s_top'
    417           ELSE
    418              roben =  's(nzt)   = s(nzt-1) + ds/dz'
    419           ENDIF
    420        ENDIF
    421 
    422        WRITE ( io, 302 ) runten, roben
    423 
     398       WRITE ( io, 301 )  'e', runten, roben       
     399
     400    ENDIF
     401
     402    IF ( ocean )  THEN
     403       runten = 'sa(0)    = sa(1)'
     404       IF ( ibc_sa_t == 0 )  THEN
     405          roben =  'sa(nzt+1) = sa_surface'
     406       ELSE
     407          roben =  'sa(nzt+1) = sa(nzt)'
     408       ENDIF
     409       WRITE ( io, 301 ) 'sa', runten, roben
     410    ENDIF
     411
     412    IF ( humidity )  THEN
     413       IF ( ibc_q_b == 0 )  THEN
     414          runten = 'q(0)     = q_surface'
     415       ELSE
     416          runten = 'q(0)     = q(1)'
     417       ENDIF
     418       IF ( ibc_q_t == 0 )  THEN
     419          roben =  'q(nzt)   = q_top'
     420       ELSE
     421          roben =  'q(nzt)   = q(nzt-1) + dq/dz'
     422       ENDIF
     423       WRITE ( io, 301 ) 'q', runten, roben
     424    ENDIF
     425
     426    IF ( passive_scalar )  THEN
     427       IF ( ibc_q_b == 0 )  THEN
     428          runten = 's(0)     = s_surface'
     429       ELSE
     430          runten = 's(0)     = s(1)'
     431       ENDIF
     432       IF ( ibc_q_t == 0 )  THEN
     433          roben =  's(nzt)   = s_top'
     434       ELSE
     435          roben =  's(nzt)   = s(nzt-1) + ds/dz'
     436       ENDIF
     437       WRITE ( io, 301 ) 's', runten, roben
    424438    ENDIF
    425439
     
    442456       IF ( constant_top_heatflux )  THEN
    443457          WRITE ( io, 306 )  top_heatflux
     458       ENDIF
     459       IF ( ocean  .AND.  constant_top_salinityflux )  THEN
     460          WRITE ( io, 309 )  top_salinityflux
    444461       ENDIF
    445462       IF ( humidity  .OR.  passive_scalar )  THEN
     
    887904!-- Other quantities
    888905    WRITE ( io, 411 )  g
    889     IF ( use_pt_reference )  WRITE ( io, 412 )  pt_reference
     906    IF ( use_reference )  THEN
     907       IF ( ocean )  THEN
     908          WRITE ( io, 412 )  prho_reference
     909       ELSE
     910          WRITE ( io, 413 )  pt_reference
     911       ENDIF
     912    ENDIF
    890913
    891914!
     
    10111034                             TRIM( gradients ), TRIM( slices )
    10121035       ENDIF
     1036    ENDIF
     1037
     1038!
     1039!-- Initial salinity profile
     1040!-- Building output strings, starting with surface salinity
     1041    IF ( ocean )  THEN
     1042       WRITE ( temperatures, '(F6.2)' )  sa_surface
     1043       gradients = '------'
     1044       slices = '     0'
     1045       coordinates = '   0.0'
     1046       i = 1
     1047       DO  WHILE ( sa_vertical_gradient_level_ind(i) /= -9999 )
     1048
     1049          WRITE (coor_chr,'(F7.2)')  sa_init(sa_vertical_gradient_level_ind(i))
     1050          temperatures = TRIM( temperatures ) // ' ' // TRIM( coor_chr )
     1051
     1052          WRITE (coor_chr,'(F7.2)')  sa_vertical_gradient(i)
     1053          gradients = TRIM( gradients ) // ' ' // TRIM( coor_chr )
     1054
     1055          WRITE (coor_chr,'(I7)')  sa_vertical_gradient_level_ind(i)
     1056          slices = TRIM( slices ) // ' ' // TRIM( coor_chr )
     1057
     1058          WRITE (coor_chr,'(F7.1)')  sa_vertical_gradient_level(i)
     1059          coordinates = TRIM( coordinates ) // ' '  // TRIM( coor_chr )
     1060
     1061          i = i + 1
     1062       ENDDO
     1063
     1064       WRITE ( io, 425 )  TRIM( coordinates ), TRIM( temperatures ), &
     1065                          TRIM( gradients ), TRIM( slices )
    10131066    ENDIF
    10141067
     
    11351188
    11361189 99 FORMAT (1X,78('-'))
    1137 100 FORMAT (/1X,'*************************',11X,28('-')/        &
     1190100 FORMAT (/1X,'*************************',11X,42('-')/        &
    11381191            1X,'* ',A,' *',11X,A/                               &
    1139             1X,'*************************',11X,28('-')//        &
     1192            1X,'*************************',11X,42('-')//        &
    11401193            ' Date:            ',A8,11X,'Run:       ',A20/      &
    11411194            ' Time:            ',A8,11X,'Run-No.:   ',I2.2/     &
     
    12561309             ' B. bound.: ',A/ &
    12571310             ' T. bound.: ',A)
    1258 301 FORMAT (/'                     e'// &
    1259              ' B. bound.: ',A/ &
    1260              ' T. bound.: ',A)
    1261 302 FORMAT (/'                     q'// &
     1311301 FORMAT (/'                     ',A// &
    12621312             ' B. bound.: ',A/ &
    12631313             ' T. bound.: ',A)
     
    12681318             '       zp = ',F6.2,' m   z0 = ',F6.4,' m   kappa = ',F4.2/ &
    12691319             '       Rif value range:   ',F6.2,' <= rif <=',F6.2)
    1270 306 FORMAT ('       Predefined constant heatflux:   ',F6.3,' K m/s')
     1320306 FORMAT ('       Predefined constant heatflux:   ',F9.6,' K m/s')
    12711321307 FORMAT ('       Heatflux has a random normal distribution')
    12721322308 FORMAT ('       Predefined surface temperature')
     1323309 FORMAT ('       Predefined constant salinityflux:   ',F9.6,' psu m/s')
    12731324310 FORMAT (//'    1D-Model:'// &
    12741325             '       Rif value range:   ',F6.2,' <= rif <=',F6.2)
     
    13521403            '                            f*    = ',F9.6,' 1/s')
    13531404411 FORMAT (/'    Gravity             :   g     = ',F4.1,' m/s**2')
    1354 412 FORMAT (/'    Reference temperature in buoyancy terms: ',F8.4,' K')
     1405412 FORMAT (/'    Reference density in buoyancy terms: ',F8.3,' kg/m**3')
     1406413 FORMAT (/'    Reference temperature in buoyancy terms: ',F8.4,' K')
    13551407415 FORMAT (/'    Cloud physics parameters:'/ &
    13561408             '    ------------------------'/)
     
    13821434424 FORMAT (/'    Characteristic levels of the geo. wind component vg:'// &
    13831435            '       Height:      ',A,'  m'/ &
    1384             '       vg:          ',A,'  m/S'/ &
     1436            '       vg:          ',A,'  m/s'/ &
    13851437            '       Gradient:    ',A,'  1/100s'/ &
    13861438            '       Gridpoint:   ',A)
     1439425 FORMAT (/'    Characteristic levels of the initial salinity profile:'// &
     1440            '       Height:     ',A,'  m'/ &
     1441            '       Salinity:   ',A,'  psu'/ &
     1442            '       Gradient:   ',A,'  psu/100m'/ &
     1443            '       Gridpoint:  ',A)
    13871444450 FORMAT (//' LES / Turbulence quantities:'/ &
    13881445              ' ---------------------------'/)
  • palm/trunk/SOURCE/init_3d_model.f90

    r96 r97  
    529529       hom(:,1,24,:) = SPREAD( kh(:,nys,nxl), 2, statistic_regions+1 )
    530530
     531       IF ( ocean )  THEN
     532!
     533!--       Store initial salinity profile
     534          hom(:,1,26,:)  = SPREAD( sa(:,nys,nxl), 2, statistic_regions+1 )
     535       ENDIF
    531536
    532537       IF ( humidity )  THEN
  • palm/trunk/SOURCE/init_ocean.f90

    r96 r97  
    4040    hyp(nzt+1) = surface_pressure * 100.0
    4141
    42     hyp(nzt)   = hyp(nzt+1) + rho_surface * g * 0.5 * dzu(nzt+1)
    43     rho_ref    = rho_surface * 0.5 * dzu(nzt+1)
     42    hyp(nzt)      = hyp(nzt+1) + rho_surface * g * 0.5 * dzu(nzt+1)
     43    rho_reference = rho_surface * 0.5 * dzu(nzt+1)
    4444
    4545    DO  k = nzt-1, 0, -1
     
    4848       pt_l = 0.5 * ( pt_init(k) + pt_init(k+1) )
    4949
    50        rho_l   = eqn_state_seawater_func( hyp(k+1), pt_l, sa_l )
     50       rho_l = eqn_state_seawater_func( hyp(k+1), pt_l, sa_l )
    5151
    52        hyp(k)  = hyp(k+1) + rho_l * g * dzu(k+1)
    53        rho_ref = rho_ref + rho_l * dzu(k+1)
     52       hyp(k)        = hyp(k+1) + rho_l * g * dzu(k+1)
     53       rho_reference = rho_reference + rho_l * dzu(k+1)
    5454
    5555    ENDDO
    5656
    57     rho_ref = rho_ref / ( zw(nzt) - zu(nzb) )
    58     print*, '*** rho_ref = ', rho_ref
     57    rho_reference = rho_reference / ( zw(nzt) - zu(nzb) )
     58
     59!
     60!-- Calculate the reference potential density
     61    prho_reference = 0.0
     62    DO  k = 0, nzt
     63
     64       sa_l = 0.5 * ( sa_init(k) + sa_init(k+1) )
     65       pt_l = 0.5 * ( pt_init(k) + pt_init(k+1) )
     66
     67       prho_reference = prho_reference + dzu(k+1) * &
     68                        eqn_state_seawater_func( hyp(0), pt_l, sa_l )
     69
     70    ENDDO
     71
     72    prho_reference = prho_reference / ( zu(nzt) - zu(nzb) )
    5973
    6074
  • palm/trunk/SOURCE/modules.f90

    r96 r97  
    55! Actual revisions:
    66! -----------------
    7 ! +ocean, r, + salinity variables
     7! +atmos_ocean_sign, ocean, r, + salinity variables
    88! defaults of .._vertical_gradient_levels changed from -1.0 to -9999999.9
    9 ! hydro_press renamed hyp
     9! hydro_press renamed hyp, use_pt_reference renamed use_reference
    1010!
    1111! Former revisions:
     
    324324                sloping_surface = .FALSE., stop_dt = .FALSE., &
    325325                terminate_run = .FALSE., use_prior_plot1d_parameters = .FALSE.,&
    326                 use_pt_reference = .FALSE., use_surface_fluxes = .FALSE., &
     326                use_reference = .FALSE., use_surface_fluxes = .FALSE., &
    327327                use_top_fluxes = .FALSE., use_ug_for_galilei_tr = .TRUE., &
    328328                use_upstream_for_tke = .FALSE., wall_adjustment = .TRUE.
     
    333333    REAL ::  advected_distance_x = 0.0, advected_distance_y = 0.0, &
    334334             alpha_surface = 0.0, asselin_filter_factor = 0.1, &
     335             atmos_ocean_sign = 1.0, &
    335336             averaging_interval = 0.0, averaging_interval_pr = 9999999.9, &
    336337             averaging_interval_sp = 9999999.9, bc_pt_t_val, bc_q_t_val, &
     
    340341             building_wall_south = 9999999.9, cfl_factor = -1.0, &
    341342             cos_alpha_surface, disturbance_amplitude = 0.25, &
    342              disturbance_energy_limit = 0.01, disturbance_level_b = -1.0, &
    343              disturbance_level_t = -1.0, dt = -1.0, dt_averaging_input = 0.0, &
     343             disturbance_energy_limit = 0.01, &
     344             disturbance_level_b = -9999999.9, &
     345             disturbance_level_t = -9999999.9, &
     346             dt = -1.0, dt_averaging_input = 0.0, &
    344347             dt_averaging_input_pr = 9999999.9, dt_data_output = 9999999.9, &
    345348             dt_data_output_av = 9999999.9, dt_disturb = 9999999.9, &
     
    360363             overshoot_limit_w = 0.0, particle_maximum_age = 9999999.9, &
    361364             phi = 55.0, prandtl_number = 1.0, &
    362              precipitation_amount_interval = 9999999.9, &
    363              pt_reference = 9999999.9, &
    364              pt_slope_offset = 0.0, pt_surface = 300.0, &
    365              pt_surface_initial_change = 0.0, q_surface = 0.0, &
    366              q_surface_initial_change = 0.0, rayleigh_damping_factor = -1.0, &
    367              rayleigh_damping_height = -1.0, residual_limit = 1.0E-4, &
    368              restart_time = 9999999.9, rho_ref, rho_surface, rif_max = 1.0, &
     365             precipitation_amount_interval = 9999999.9, prho_reference, &
     366             pt_reference = 9999999.9, pt_slope_offset = 0.0, &
     367             pt_surface = 300.0, pt_surface_initial_change = 0.0, &
     368             q_surface = 0.0, q_surface_initial_change = 0.0, &
     369             rayleigh_damping_factor = -1.0, rayleigh_damping_height = -1.0, &
     370             residual_limit = 1.0E-4, restart_time = 9999999.9, rho_reference, &
     371             rho_surface, rif_max = 1.0, &
    369372             rif_min = -5.0, roughness_length = 0.1, sa_surface = 35.0, &
    370373             simulated_time = 0.0, simulated_time_at_begin, sin_alpha_surface, &
     
    401404             tsc(10) = (/ 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 /), &
    402405             ug_vertical_gradient(10) = 0.0, &
    403              ug_vertical_gradient_level(10) = -1.0, &
     406             ug_vertical_gradient_level(10) = -9999999.9, &
    404407             vg_vertical_gradient(10) = 0.0, &
    405              vg_vertical_gradient_level(10) = -1.0, &
     408             vg_vertical_gradient_level(10) = -9999999.9, &
    406409             volume_flow(1:2) = 0.0, volume_flow_area(1:2) = 0.0, &
    407410             volume_flow_initial(1:2) = 0.0, wall_heatflux(0:4) = 0.0
  • palm/trunk/SOURCE/netcdf.f90

    r77 r97  
    77! Current revisions:
    88! ------------------
    9 !
     9! Grids defined for rho and sa
    1010!
    1111! Former revisions:
     
    345345!--             Most variables are defined on the scalar grid
    346346                CASE ( 'e', 'p', 'pc', 'pr', 'pt', 'q', 'ql', 'ql_c', 'ql_v', &
    347                        'ql_vp', 'qv', 's', 'vpt' )
     347                       'ql_vp', 'qv', 'rho', 's', 'sa', 'vpt' )
    348348
    349349                   grid_x = 'x'
     
    870870                      CASE ( 'e_xy', 'p_xy', 'pc_xy', 'pr_xy', 'pt_xy', 'q_xy',&
    871871                             'ql_xy', 'ql_c_xy', 'ql_v_xy', 'ql_vp_xy',        &
    872                              'qv_xy', 's_xy', 'vpt_xy' )
     872                             'qv_xy', 'rho_xy', 's_xy', 'sa_xy', 'vpt_xy' )
    873873
    874874                         grid_x = 'x'
     
    14141414                   CASE ( 'e_xz', 'p_xz', 'pc_xz', 'pr_xz', 'pt_xz', 'q_xz',  &
    14151415                          'ql_xz', 'ql_c_xz', 'ql_v_xz', 'ql_vp_xz', 'qv_xz', &
    1416                           's_xz', 'vpt_xz' )
     1416                          'rho_xz', 's_xz', 'sa_xz', 'vpt_xz' )
    14171417
    14181418                      grid_x = 'x'
     
    19251925                   CASE ( 'e_yz', 'p_yz', 'pc_yz', 'pr_yz', 'pt_yz', 'q_yz',  &
    19261926                          'ql_yz', 'ql_c_yz', 'ql_v_yz', 'ql_vp_yz', 'qv_yz', &
    1927                           's_yz', 'vpt_yz' )
     1927                          'rho_yz', 's_yz', 'sa_yz', 'vpt_yz' )
    19281928
    19291929                      grid_x = 'x'
  • palm/trunk/SOURCE/palm.f90

    r90 r97  
    6363    INTEGER           ::  i, run_description_header_i(80)
    6464
    65     version = 'PALM 3.2b'
     65    version = 'PALM 3.3'
    6666
    6767#if defined( __parallel )
  • palm/trunk/SOURCE/production_e.f90

    r77 r97  
    44! Actual revisions:
    55! -----------------
    6 !
     6! energy production by density flux (in ocean) added
     7! use_pt_reference renamed use_reference
    78!
    89! Former revisions:
     
    364365          IF ( .NOT. humidity )  THEN
    365366
    366              IF ( use_pt_reference )  THEN
    367 
    368                 DO  j = nys, nyn
    369                    DO  k = nzb_diff_s_inner(j,i), nzt_diff
    370                       tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g/pt_reference * &
     367             IF ( use_reference )  THEN
     368
     369                IF ( ocean )  THEN
     370!
     371!--                So far in the ocean no special treatment of density flux in
     372!--                the bottom and top surface layer
     373                   DO  j = nys, nyn
     374                      DO  k = nzb_s_inner(j,i), nzt
     375                         tend(k,j,i) = tend(k,j,i) +                    &
     376                                       kh(k,j,i) * g / prho_reference * &
     377                                       ( rho(k+1,j,i)-rho(k-1,j,i) ) * dd2zu(k)
     378                      ENDDO
     379                   ENDDO
     380
     381                ELSE
     382
     383                   DO  j = nys, nyn
     384                      DO  k = nzb_diff_s_inner(j,i), nzt_diff
     385                         tend(k,j,i) = tend(k,j,i) +                  &
     386                                       kh(k,j,i) * g / pt_reference * &
    371387                                       ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k)
     388                      ENDDO
     389
     390                      IF ( use_surface_fluxes )  THEN
     391                         k = nzb_diff_s_inner(j,i)-1
     392                         tend(k,j,i) = tend(k,j,i) + g / pt_reference * shf(j,i)
     393                      ENDIF
     394
     395                      IF ( use_top_fluxes )  THEN
     396                         k = nzt
     397                         tend(k,j,i) = tend(k,j,i) + g / pt_reference * &
     398                                                     tswst(j,i)
     399                      ENDIF
    372400                   ENDDO
    373401
    374                    IF ( use_surface_fluxes )  THEN
    375                       k = nzb_diff_s_inner(j,i)-1
    376                       tend(k,j,i) = tend(k,j,i) + g / pt_reference * shf(j,i)
    377                    ENDIF
    378 
    379                    IF ( use_top_fluxes )  THEN
    380                       k = nzt
    381                       tend(k,j,i) = tend(k,j,i) + g / pt_reference * tswst(j,i)
    382                    ENDIF
    383                 ENDDO
     402                ENDIF
    384403
    385404             ELSE
    386405
    387                 DO  j = nys, nyn
    388                    DO  k = nzb_diff_s_inner(j,i), nzt_diff
    389                       tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / pt(k,j,i) * &
     406                IF ( ocean )  THEN
     407!
     408!--                So far in the ocean no special treatment of density flux in
     409!--                the bottom and top surface layer
     410                   DO  j = nys, nyn
     411                      DO  k = nzb_s_inner(j,i), nzt
     412                         tend(k,j,i) = tend(k,j,i) -                &
     413                                       kh(k,j,i) * g / rho(k,j,i) * &
     414                                       ( rho(k+1,j,i)-rho(k-1,j,i) ) * dd2zu(k)
     415                      ENDDO
     416                   ENDDO
     417
     418                ELSE
     419
     420                   DO  j = nys, nyn
     421                      DO  k = nzb_diff_s_inner(j,i), nzt_diff
     422                         tend(k,j,i) = tend(k,j,i) -               &
     423                                       kh(k,j,i) * g / pt(k,j,i) * &
    390424                                       ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k)
     425                      ENDDO
     426
     427                      IF ( use_surface_fluxes )  THEN
     428                         k = nzb_diff_s_inner(j,i)-1
     429                         tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * shf(j,i)
     430                      ENDIF
     431
     432                      IF ( use_top_fluxes )  THEN
     433                         k = nzt
     434                         tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * tswst(j,i)
     435                      ENDIF
    391436                   ENDDO
    392437
    393                    IF ( use_surface_fluxes )  THEN
    394                       k = nzb_diff_s_inner(j,i)-1
    395                       tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * shf(j,i)
    396                    ENDIF
    397 
    398                    IF ( use_top_fluxes )  THEN
    399                       k = nzt
    400                       tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * tswst(j,i)
    401                    ENDIF
    402                 ENDDO
     438                ENDIF
    403439
    404440             ENDIF
     
    762798       IF ( .NOT. humidity )  THEN
    763799
    764           IF ( use_pt_reference )  THEN
    765 
    766              DO  k = nzb_diff_s_inner(j,i), nzt_diff
    767                 tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / pt_reference * &
     800          IF ( use_reference )  THEN
     801
     802             IF ( ocean )  THEN
     803!
     804!--             So far in the ocean no special treatment of density flux in the
     805!--             bottom and top surface layer
     806                DO  k = nzb_s_inner(j,i), nzt
     807                   tend(k,j,i) = tend(k,j,i) + kh(k,j,i) * g / prho_reference * &
     808                                      ( rho(k+1,j,i) - rho(k-1,j,i) ) * dd2zu(k)
     809                ENDDO
     810
     811             ELSE
     812
     813                DO  k = nzb_diff_s_inner(j,i), nzt_diff
     814                   tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / pt_reference * &
    768815                                       ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k)
    769              ENDDO
    770 
    771              IF ( use_surface_fluxes )  THEN
    772                 k = nzb_diff_s_inner(j,i)-1
    773                 tend(k,j,i) = tend(k,j,i) + g / pt_reference * shf(j,i)
     816                ENDDO
     817
     818                IF ( use_surface_fluxes )  THEN
     819                   k = nzb_diff_s_inner(j,i)-1
     820                   tend(k,j,i) = tend(k,j,i) + g / pt_reference * shf(j,i)
     821                ENDIF
     822
     823                IF ( use_top_fluxes )  THEN
     824                   k = nzt
     825                   tend(k,j,i) = tend(k,j,i) + g / pt_reference * tswst(j,i)
     826                ENDIF
     827
    774828             ENDIF
    775829
    776              IF ( use_top_fluxes )  THEN
    777                 k = nzt
    778                 tend(k,j,i) = tend(k,j,i) + g / pt_reference * tswst(j,i)
     830          ELSE
     831
     832             IF ( ocean )  THEN
     833!
     834!--             So far in the ocean no special treatment of density flux in the
     835!--             bottom and top surface layer
     836                DO  k = nzb_s_inner(j,i), nzt
     837                   tend(k,j,i) = tend(k,j,i) + kh(k,j,i) * g / rho(k,j,i) * &
     838                                      ( rho(k+1,j,i) - rho(k-1,j,i) ) * dd2zu(k)
     839                ENDDO
     840
     841             ELSE
     842
     843                DO  k = nzb_diff_s_inner(j,i), nzt_diff
     844                   tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / pt(k,j,i) * &
     845                                       ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k)
     846                ENDDO
     847
     848                IF ( use_surface_fluxes )  THEN
     849                   k = nzb_diff_s_inner(j,i)-1
     850                   tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * shf(j,i)
     851                ENDIF
     852
     853                IF ( use_top_fluxes )  THEN
     854                   k = nzt
     855                   tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * tswst(j,i)
     856                ENDIF
     857
    779858             ENDIF
    780859
    781           ELSE
    782 
    783              DO  k = nzb_diff_s_inner(j,i), nzt_diff
    784                 tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / pt(k,j,i) * &
    785                                        ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k)
    786              ENDDO
    787 
    788              IF ( use_surface_fluxes )  THEN
    789                 k = nzb_diff_s_inner(j,i)-1
    790                 tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * shf(j,i)
    791              ENDIF
    792 
    793              IF ( use_top_fluxes )  THEN
    794                 k = nzt
    795                 tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * tswst(j,i)
    796              ENDIF
    797 
    798          ENDIF
     860          ENDIF
    799861
    800862       ELSE
  • palm/trunk/SOURCE/prognostic_equations.f90

    r96 r97  
    55! -----------------
    66! prognostic equation for salinity, density is calculated from equation of
    7 ! state for seawater, +eqn_state_seawater_mod
    8 ! new argument zw in calls of diffusion_e, calc_mean_pt_profile renamed
     7! state for seawater and is used for calculation of buoyancy,
     8! +eqn_state_seawater_mod
     9! diffusion_e is called with argument rho in case of ocean runs,
     10! new argument zw in calls of diffusion_e, new argument pt_/prho_reference
     11! in calls of buoyancy and diffusion_e, calc_mean_pt_profile renamed
    912! calc_mean_profile
    1013!
     
    148151          ENDIF
    149152          CALL coriolis( i, j, 1 )
    150           IF ( sloping_surface )  CALL buoyancy( i, j, pt, 1, 4 )
     153          IF ( sloping_surface )  CALL buoyancy( i, j, pt, pt_reference, 1, 4 )
    151154          CALL user_actions( i, j, 'u-tendency' )
    152155
     
    283286          ENDIF
    284287          CALL coriolis( i, j, 3 )
    285           IF ( .NOT. humidity )  THEN
    286              CALL buoyancy( i, j, pt, 3, 4 )
     288          IF ( ocean )  THEN
     289             CALL buoyancy( i, j, rho, prho_reference, 3, 64 )
    287290          ELSE
    288              CALL buoyancy( i, j, vpt, 3, 44 )
     291             IF ( .NOT. humidity )  THEN
     292                CALL buoyancy( i, j, pt, pt_reference, 3, 4 )
     293             ELSE
     294                CALL buoyancy( i, j, vpt, pt_reference, 3, 44 )
     295             ENDIF
    289296          ENDIF
    290297          CALL user_actions( i, j, 'w-tendency' )
     
    654661                  .NOT. use_upstream_for_tke )  THEN
    655662                IF ( .NOT. humidity )  THEN
    656                    CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, km, &
    657                                      l_grid, pt, rif, tend, zu, zw )
     663                   IF ( ocean )  THEN
     664                      CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, km,   &
     665                                        l_grid, rho, prho_reference, rif, tend, &
     666                                        zu, zw )
     667                   ELSE
     668                      CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, km, &
     669                                        l_grid, pt, pt_reference, rif, tend,  &
     670                                        zu, zw )
     671                   ENDIF
    658672                ELSE
    659                    CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, km, &
    660                                      l_grid, vpt, rif, tend, zu, zw )
     673                   CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, km,     &
     674                                     l_grid, vpt, pt_reference, rif, tend, zu, &
     675                                     zw )
    661676                ENDIF
    662677             ELSE
     
    680695                   IF ( .NOT. humidity )  THEN
    681696                      CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e_m, &
    682                                         km_m, l_grid, pt_m, rif_m, tend, zu, &
    683                                         zw )
     697                                        km_m, l_grid, pt_m, pt_reference,  &
     698                                        rif_m, tend, zu, zw )
    684699                   ELSE
    685700                      CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e_m, &
    686                                         km_m, l_grid, vpt_m, rif_m, tend, zu, &
    687                                         zw )
     701                                        km_m, l_grid, vpt_m, pt_reference, &
     702                                        rif_m, tend, zu, zw )
    688703                   ENDIF
    689704                ELSE
    690705                   IF ( .NOT. humidity )  THEN
    691                       CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, km, &
    692                                         l_grid, pt, rif, tend, zu, zw )
     706                      IF ( ocean )  THEN
     707                         CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e,  &
     708                                           km, l_grid, rho, prho_reference,   &
     709                                           rif, tend, zu, zw )
     710                      ELSE
     711                         CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e,  &
     712                                           km, l_grid, pt, pt_reference, rif, &
     713                                           tend, zu, zw )
     714                      ENDIF
    693715                   ELSE
    694716                      CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, km, &
    695                                         l_grid, vpt, rif, tend, zu, zw )
     717                                        l_grid, vpt, pt_reference, rif, tend, &
     718                                        zu, zw )
    696719                   ENDIF
    697720                ENDIF
     
    797820             ENDIF
    798821             CALL coriolis( i, j, 1 )
    799              IF ( sloping_surface )  CALL buoyancy( i, j, pt, 1, 4 )
     822             IF ( sloping_surface )  CALL buoyancy( i, j, pt, pt_reference, 1, &
     823                                                    4 )
    800824             CALL user_actions( i, j, 'u-tendency' )
    801825
     
    896920             ENDIF
    897921             CALL coriolis( i, j, 3 )
    898              IF ( .NOT. humidity )  THEN
    899                 CALL buoyancy( i, j, pt, 3, 4 )
     922             IF ( ocean )  THEN
     923                CALL buoyancy( i, j, rho, prho_reference, 3, 64 )
    900924             ELSE
    901                 CALL buoyancy( i, j, vpt, 3, 44 )
     925                IF ( .NOT. humidity )  THEN
     926                   CALL buoyancy( i, j, pt, pt_reference, 3, 4 )
     927                ELSE
     928                   CALL buoyancy( i, j, vpt, pt_reference, 3, 44 )
     929                ENDIF
    902930             ENDIF
    903931             CALL user_actions( i, j, 'w-tendency' )
     
    11141142                   IF ( .NOT. humidity )  THEN
    11151143                      CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e_m, &
    1116                                         km_m, l_grid, pt_m, rif_m, tend, zu, &
    1117                                         zw )
     1144                                        km_m, l_grid, pt_m, pt_reference,  &
     1145                                        rif_m, tend, zu, zw )
    11181146                   ELSE
    11191147                      CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e_m, &
    1120                                         km_m, l_grid, vpt_m, rif_m, tend, zu, &
    1121                                         zw )
     1148                                        km_m, l_grid, vpt_m, pt_reference, &
     1149                                        rif_m, tend, zu, zw )
    11221150                   ENDIF
    11231151                ELSE
    11241152                   IF ( .NOT. humidity )  THEN
    1125                       CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, km, &
    1126                                         l_grid, pt, rif, tend, zu, zw )
     1153                      IF ( ocean )  THEN
     1154                         CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, &
     1155                                           km, l_grid, rho, prho_reference,  &
     1156                                           rif, tend, zu, zw )
     1157                      ELSE
     1158                         CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e,  &
     1159                                           km, l_grid, pt, pt_reference, rif, &
     1160                                           tend, zu, zw )
     1161                      ENDIF
    11271162                   ELSE
    11281163                      CALL diffusion_e( i, j, ddzu, dd2zu, ddzw, diss, e, km, &
    1129                                         l_grid, vpt, rif, tend, zu, zw )
     1164                                        l_grid, vpt, pt_reference, rif, tend, &
     1165                                        zu, zw )
    11301166                   ENDIF
    11311167                ENDIF
     
    12241260    ENDIF
    12251261    CALL coriolis( 1 )
    1226     IF ( sloping_surface )  CALL buoyancy( pt, 1, 4 )
     1262    IF ( sloping_surface )  CALL buoyancy( pt, pt_reference, 1, 4 )
    12271263    CALL user_actions( 'u-tendency' )
    12281264
     
    13671403    ENDIF
    13681404    CALL coriolis( 3 )
    1369     IF ( .NOT. humidity )  THEN
    1370        CALL buoyancy( pt, 3, 4 )
     1405    IF ( ocean )  THEN
     1406       CALL buoyancy( rho, prho_reference, 3, 64 )
    13711407    ELSE
    1372        CALL buoyancy( vpt, 3, 44 )
     1408       IF ( .NOT. humidity )  THEN
     1409          CALL buoyancy( pt, pt_reference, 3, 4 )
     1410       ELSE
     1411          CALL buoyancy( vpt, pt_reference, 3, 44 )
     1412       ENDIF
    13731413    ENDIF
    13741414    CALL user_actions( 'w-tendency' )
     
    17531793       THEN
    17541794          IF ( .NOT. humidity )  THEN
    1755              CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, pt, &
    1756                                rif, tend, zu, zw )
     1795             IF ( ocean )  THEN
     1796                CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, rho, &
     1797                                  prho_reference, rif, tend, zu, zw )
     1798             ELSE
     1799                CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, pt, &
     1800                                  pt_reference, rif, tend, zu, zw )
     1801             ENDIF
    17571802          ELSE
    17581803             CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, vpt, &
    1759                                rif, tend, zu, zw )
     1804                               pt_reference, rif, tend, zu, zw )
    17601805          ENDIF
    17611806       ELSE
     
    17771822             IF ( .NOT. humidity )  THEN
    17781823                CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e_m, km_m, l_grid, &
    1779                                   pt_m, rif_m, tend, zu, zw )
     1824                                  pt_m, pt_reference, rif_m, tend, zu, zw )
    17801825             ELSE
    17811826                CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e_m, km_m, l_grid, &
    1782                                   vpt_m, rif_m, tend, zu, zw )
     1827                                  vpt_m, pt_reference, rif_m, tend, zu, zw )
    17831828             ENDIF
    17841829          ELSE
    17851830             IF ( .NOT. humidity )  THEN
    1786                 CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, pt, &
    1787                                   rif, tend, zu, zw )
     1831                IF ( ocean )  THEN
     1832                   CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, &
     1833                                     rho, prho_reference, rif, tend, zu, zw )
     1834                ELSE
     1835                   CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, &
     1836                                     pt, pt_reference, rif, tend, zu, zw )
     1837                ENDIF
    17881838             ELSE
    17891839                CALL diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, vpt, &
    1790                                   rif, tend, zu, zw )
     1840                                  pt_reference, rif, tend, zu, zw )
    17911841             ENDIF
    17921842          ENDIF
  • palm/trunk/SOURCE/run_control.f90

    r90 r97  
    44! Actual revisions:
    55! -----------------
    6 !
     6! Timestep format changed
    77!
    88! Former revisions:
     
    113113100 FORMAT (///'Run-control output:'/ &
    114114              &'------------------'// &
    115            &'RUN  ITER. HH:MM:SS.SS   DT(E)     UMAX     VMAX     WMAX     U*', &
    116            &'    W*   THETA*   Z_I     ENERG.   DISTENERG    DIVOLD     DIVNE', &
    117            &'W     UMAX(KJI)    VMAX(KJI)    WMAX(KJI)   ADVECX   ADVECY   MG', &
    118            &'CYC'/ &
    119            &'----------------------------------------------------------------', &
    120            &'----------------------------------------------------------------', &
    121            &'----------------------------------------------------------------', &
    122            &'--')
    123 101 FORMAT (I3,1X,I6,1X,A8,F3.2,1X,F7.4,A1,A1,F8.4,A1,F8.4,A1,F8.4,2X,F5.3,2X, &
     115          &'RUN  ITER. HH:MM:SS.SS    DT(E)     UMAX     VMAX     WMAX     U', &
     116          &'*    W*   THETA*    Z_I     ENERG.   DISTENERG    DIVOLD     DIV', &
     117          &'NEW     UMAX(KJI)    VMAX(KJI)    WMAX(KJI)   ADVECX   ADVECY   ', &
     118          &'MGCYC'/                                                            &
     119          &'----------------------------------------------------------------', &
     120          &'----------------------------------------------------------------', &
     121          &'----------------------------------------------------------------', &
     122          &'-----')
     123101 FORMAT (I3,1X,I6,1X,A8,F3.2,1X,F8.4,A1,A1,F8.4,A1,F8.4,A1,F8.4,2X,F5.3,2X, &
    124124            F4.2, &
    125             2X,F6.3,2X,F5.0,1X,4(E10.3,1X),3(3(I4),1X),F8.3,1X,F8.3,5X,I3)
     125            2X,F6.3,2X,F6.0,1X,4(E10.3,1X),3(3(I4),1X),F8.3,1X,F8.3,5X,I3)
    126126
    127127 END SUBROUTINE run_control
  • palm/trunk/SOURCE/time_integration.f90

    r95 r97  
    44! Actual revisions:
    55! -----------------
    6 ! Ghostpoint exchange for salinity and density
     6! diffusivities is called with argument rho in case of ocean runs,
     7! new argument pt_/prho_reference in calls of diffusivities,
     8! ghostpoint exchange for salinity and density
    79!
    810! Former revisions:
     
    249251             CALL cpu_log( log_point(17), 'diffusivities', 'start' )
    250252             IF ( .NOT. humidity ) THEN
    251                 CALL diffusivities( pt )
     253                IF ( ocean )  THEN
     254                   CALL diffusivities( rho, prho_reference )
     255                ELSE
     256                   CALL diffusivities( pt, pt_reference )
     257                ENDIF
    252258             ELSE
    253                 CALL diffusivities( vpt )
     259                CALL diffusivities( vpt, pt_reference )
    254260             ENDIF
    255261             CALL cpu_log( log_point(17), 'diffusivities', 'stop' )
Note: See TracChangeset for help on using the changeset viewer.