Changeset 336 for palm/trunk/SOURCE/init_ocean.f90

Timestamp:

Jun 10, 2009 11:19:35 AM (15 years ago)

Author:

raasch

Message:

several small bugfixes; some more dvrp changes

File:

• palm/trunk/SOURCE/init_ocean.f90 (modified) (4 diffs)

palm/trunk/SOURCE/init_ocean.f90

@@ -24 +24 @@
     USE control_parameters
     USE eqn_state_seawater_mod
+    USE pegrid
     USE grid_variables
     USE indices
@@ -29 +30 @@
     IMPLICIT NONE
 
-    INTEGER ::  k
+    INTEGER ::  k, n
 
     REAL    ::  sa_l, pt_l, rho_l
+
+    REAL, DIMENSION(nzb:nzt+1) ::  rho_init
 
     ALLOCATE( hyp(nzb:nzt+1) )
@@ -45 +48 @@
 
     hyp(nzt)      = hyp(nzt+1) + rho_surface * g * 0.5 * dzu(nzt+1)
-    rho_reference = rho_surface * 0.5 * dzu(nzt+1)
+    rho_init(nzt) = rho_surface
 
     DO  k = nzt-1, 0, -1
+       hyp(k) = hyp(k+1) + rho_surface * g * dzu(k)
+    ENDDO
 
-       sa_l = 0.5 * ( sa_init(k) + sa_init(k+1) )
-       pt_l = 0.5 * ( pt_init(k) + pt_init(k+1) )
+    IF ( myid == 0 )  THEN
+       print*,'hydro pres using rho_surface'
+       DO  k = nzt+1, 0, -1
+          print*, 'k = ', k, ' hyp = ', hyp(k)
+       ENDDO
+       print*, ' '
+    ENDIF
 
-       rho_l = eqn_state_seawater_func( hyp(k+1), pt_l, sa_l )
+    DO  n = 1, 5
 
-       hyp(k)        = hyp(k+1) + rho_l * g * dzu(k+1)
-       rho_reference = rho_reference + rho_l * dzu(k+1)
+       rho_reference = rho_surface * 0.5 * dzu(nzt+1)
+
+       DO  k = nzt-1, 0, -1
+
+          sa_l = 0.5 * ( sa_init(k) + sa_init(k+1) )
+          pt_l = 0.5 * ( pt_init(k) + pt_init(k+1) )
+
+          rho_init(k) = eqn_state_seawater_func( hyp(k), pt_l, sa_l )
+
+          rho_reference = rho_reference + rho_init(k) * dzu(k+1)
+
+       ENDDO
+
+       rho_reference = rho_reference / ( zw(nzt) - zu(nzb) )
+
+       DO  k = nzt-1, 0, -1
+          hyp(k) = hyp(k+1) + g * 0.5 * ( rho_init(k) + rho_init(k+1 ) ) * &
+                              dzu(k+1)
+       ENDDO
+
+       IF ( myid == 0 )  THEN
+          print*,'hydro pres / rho  n = ', n
+          DO  k = nzt+1, 0, -1
+             print*, 'k = ', k, ' hyp = ', hyp(k), ' rho = ', rho_init(k)
+          ENDDO
+          print*, ' '
+       ENDIF
 
     ENDDO
-
-    rho_reference = rho_reference / ( zw(nzt) - zu(nzb) )
 
 !
@@ -70 +103 @@
 
        prho_reference = prho_reference + dzu(k+1) * &
-                        eqn_state_seawater_func( hyp(0), pt_l, sa_l )
+                        eqn_state_seawater_func( 0.0, pt_l, sa_l )
 
     ENDDO