Ignore:
Timestamp:
Jun 10, 2009 11:19:35 AM (12 years ago)
Author:
raasch
Message:

several small bugfixes; some more dvrp changes

File:
1 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SOURCE/init_ocean.f90

    r139 r336  
    2424    USE control_parameters
    2525    USE eqn_state_seawater_mod
     26    USE pegrid
    2627    USE grid_variables
    2728    USE indices
     
    2930    IMPLICIT NONE
    3031
    31     INTEGER ::  k
     32    INTEGER ::  k, n
    3233
    3334    REAL    ::  sa_l, pt_l, rho_l
     35
     36    REAL, DIMENSION(nzb:nzt+1) ::  rho_init
    3437
    3538    ALLOCATE( hyp(nzb:nzt+1) )
     
    4548
    4649    hyp(nzt)      = hyp(nzt+1) + rho_surface * g * 0.5 * dzu(nzt+1)
    47     rho_reference = rho_surface * 0.5 * dzu(nzt+1)
     50    rho_init(nzt) = rho_surface
    4851
    4952    DO  k = nzt-1, 0, -1
     53       hyp(k) = hyp(k+1) + rho_surface * g * dzu(k)
     54    ENDDO
    5055
    51        sa_l = 0.5 * ( sa_init(k) + sa_init(k+1) )
    52        pt_l = 0.5 * ( pt_init(k) + pt_init(k+1) )
     56    IF ( myid == 0 )  THEN
     57       print*,'hydro pres using rho_surface'
     58       DO  k = nzt+1, 0, -1
     59          print*, 'k = ', k, ' hyp = ', hyp(k)
     60       ENDDO
     61       print*, ' '
     62    ENDIF
    5363
    54        rho_l = eqn_state_seawater_func( hyp(k+1), pt_l, sa_l )
     64    DO  n = 1, 5
    5565
    56        hyp(k)        = hyp(k+1) + rho_l * g * dzu(k+1)
    57        rho_reference = rho_reference + rho_l * dzu(k+1)
     66       rho_reference = rho_surface * 0.5 * dzu(nzt+1)
     67
     68       DO  k = nzt-1, 0, -1
     69
     70          sa_l = 0.5 * ( sa_init(k) + sa_init(k+1) )
     71          pt_l = 0.5 * ( pt_init(k) + pt_init(k+1) )
     72
     73          rho_init(k) = eqn_state_seawater_func( hyp(k), pt_l, sa_l )
     74
     75          rho_reference = rho_reference + rho_init(k) * dzu(k+1)
     76
     77       ENDDO
     78
     79       rho_reference = rho_reference / ( zw(nzt) - zu(nzb) )
     80
     81       DO  k = nzt-1, 0, -1
     82          hyp(k) = hyp(k+1) + g * 0.5 * ( rho_init(k) + rho_init(k+1 ) ) * &
     83                              dzu(k+1)
     84       ENDDO
     85
     86       IF ( myid == 0 )  THEN
     87          print*,'hydro pres / rho  n = ', n
     88          DO  k = nzt+1, 0, -1
     89             print*, 'k = ', k, ' hyp = ', hyp(k), ' rho = ', rho_init(k)
     90          ENDDO
     91          print*, ' '
     92       ENDIF
    5893
    5994    ENDDO
    60 
    61     rho_reference = rho_reference / ( zw(nzt) - zu(nzb) )
    6295
    6396!
     
    70103
    71104       prho_reference = prho_reference + dzu(k+1) * &
    72                         eqn_state_seawater_func( hyp(0), pt_l, sa_l )
     105                        eqn_state_seawater_func( 0.0, pt_l, sa_l )
    73106
    74107    ENDDO
Note: See TracChangeset for help on using the changeset viewer.