Changeset 151 for palm/trunk/SOURCE/init_3d_model.f90

Timestamp:

Mar 7, 2008 1:42:18 PM (16 years ago)

Author:

raasch

Message:

preliminary update for the turbulence recycling method

File:

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

palm/trunk/SOURCE/init_3d_model.f90

@@ -602 +602 @@
 
 !
-!--    apply channel flow boundary condition
+!--    Apply channel flow boundary condition
        IF ( TRIM( bc_uv_t ) == 'dirichlet_0' )  THEN
 
@@ -608 +608 @@
           v(nzt+1,:,:) = 0.0
 
-!--       for the Dirichlet condition to be correctly applied at the top, set
+!--       For the Dirichlet condition to be correctly applied at the top, set
 !--       ug and vg to zero there
           ug(nzt+1)    = 0.0
@@ -755 +755 @@
 
 !
+!--    For the moment, perturbation pressure and vertical velocity are zero
+       p = 0.0; w = 0.0
+
+!
+!--    Initialize array sums (must be defined in first call of pres)
+       sums = 0.0
+
+!
+!--    Treating cloud physics, liquid water content and precipitation amount
+!--    are zero at beginning of the simulation
+       IF ( cloud_physics )  THEN
+          ql = 0.0
+          IF ( precipitation )  precipitation_amount = 0.0
+       ENDIF
+
+!
+!--    Impose vortex with vertical axis on the initial velocity profile
+       IF ( INDEX( initializing_actions, 'initialize_vortex' ) /= 0 )  THEN
+          CALL init_rankine
+       ENDIF
+
+!
+!--    Impose temperature anomaly (advection test only)
+       IF ( INDEX( initializing_actions, 'initialize_ptanom' ) /= 0 )  THEN
+          CALL init_pt_anomaly
+       ENDIF
+
+!
+!--    If required, change the surface temperature at the start of the 3D run
+       IF ( pt_surface_initial_change /= 0.0 )  THEN
+          pt(nzb,:,:) = pt(nzb,:,:) + pt_surface_initial_change
+       ENDIF
+
+!
+!--    If required, change the surface humidity/scalar at the start of the 3D
+!--    run
+       IF ( ( humidity .OR. passive_scalar ) .AND. &
+            q_surface_initial_change /= 0.0 )  THEN
+          q(nzb,:,:) = q(nzb,:,:) + q_surface_initial_change
+       ENDIF
+
+!
+!--    Initialize the random number generator (from numerical recipes)
+       CALL random_function_ini
+
+!
+!--    Impose random perturbation on the horizontal velocity field and then
+!--    remove the divergences from the velocity field
+       IF ( create_disturbances )  THEN
+          CALL disturb_field( nzb_u_inner, tend, u )
+          CALL disturb_field( nzb_v_inner, tend, v )
+          n_sor = nsor_ini
+          CALL pres
+          n_sor = nsor
+       ENDIF
+
+!
+!--    Once again set the perturbation pressure explicitly to zero in order to
+!--    assure that it does not generate any divergences in the first time step.
+!--    At t=0 the velocity field is free of divergence (as constructed above).
+!--    Divergences being created during a time step are not yet known and thus
+!--    cannot be corrected during the time step yet.
+       p = 0.0
+
+!
+!--    Initialize old and new time levels.
+       IF ( timestep_scheme(1:5) /= 'runge' )  THEN
+          e_m = e; pt_m = pt; u_m = u; v_m = v; w_m = w; kh_m = kh; km_m = km
+       ELSE
+          te_m = 0.0; tpt_m = 0.0; tu_m = 0.0; tv_m = 0.0; tw_m = 0.0
+       ENDIF
+       e_p = e; pt_p = pt; u_p = u; v_p = v; w_p = w
+
+       IF ( humidity  .OR.  passive_scalar )  THEN
+          IF ( ASSOCIATED( q_m ) )  q_m = q
+          IF ( timestep_scheme(1:5) == 'runge' )  tq_m = 0.0
+          q_p = q
+          IF ( humidity  .AND.  ASSOCIATED( vpt_m ) )  vpt_m = vpt
+       ENDIF
+
+       IF ( ocean )  THEN
+          tsa_m = 0.0
+          sa_p  = sa
+       ENDIF
+
+
+    ELSEIF ( TRIM( initializing_actions ) == 'read_restart_data'  .OR.    &
+             TRIM( initializing_actions ) == 'read_data_for_recycling' )  &
+    THEN
+!
+!--    When reading data for initializing the recycling method, first read
+!--    some of the global variables from restart file
+       IF ( TRIM( initializing_actions ) == 'read_data_for_recycling' )  THEN
+
+          WRITE (9,*) 'before read_parts_of_var_list'
+          CALL local_flush( 9 )
+          CALL read_parts_of_var_list
+          WRITE (9,*) 'after read_parts_of_var_list'
+          CALL local_flush( 9 )
+          CALL close_file( 13 )
+!
+!--       Store temporally and horizontally averaged vertical profiles to be
+!--       used as mean inflow profiles
+          ALLOCATE( mean_inflow_profiles(nzb:nzt+1,5) )
+
+          mean_inflow_profiles(:,1) = hom_sum(:,1,0)    ! u
+          mean_inflow_profiles(:,2) = hom_sum(:,2,0)    ! v
+          mean_inflow_profiles(:,4) = hom_sum(:,4,0)    ! pt
+          mean_inflow_profiles(:,5) = hom_sum(:,8,0)    ! e
+
+!
+!--       Use these mean profiles for the inflow (provided that Dirichlet
+!--       conditions are used)
+          IF ( inflow_l )  THEN
+             DO  j = nys-1, nyn+1
+                DO  k = nzb, nzt+1
+                   u(k,j,-1)  = mean_inflow_profiles(k,1)
+                   v(k,j,-1)  = mean_inflow_profiles(k,2)
+                   w(k,j,-1)  = 0.0
+                   pt(k,j,-1) = mean_inflow_profiles(k,4)
+                   e(k,j,-1)  = mean_inflow_profiles(k,5)
+                ENDDO
+             ENDDO
+          ENDIF
+
+!
+!--       Calculate the damping factors to be used at the inflow. For a
+!--       turbulent inflow the turbulent fluctuations have to be limited
+!--       vertically because otherwise the turbulent inflow layer will grow
+!--       in time.
+          IF ( inflow_damping_height == 9999999.9 )  THEN
+!
+!--          Default: use the inversion height calculated by the prerun
+             inflow_damping_height = hom_sum(nzb+6,pr_palm,0)
+
+          ENDIF
+
+          IF ( inflow_damping_width == 9999999.9 )  THEN
+!
+!--          Default for the transition range: one tenth of the undamped layer
+             inflow_damping_width = 0.1 * inflow_damping_height
+
+          ENDIF
+
+          ALLOCATE( inflow_damping_factor(nzb:nzt+1) )
+
+          DO  k = nzb, nzt+1
+
+             IF ( zu(k) <= inflow_damping_height )  THEN
+                inflow_damping_factor(k) = 1.0
+             ELSEIF ( zu(k) <= inflow_damping_height + inflow_damping_width ) &
+             THEN
+                inflow_damping_factor(k) = 1.0 -                               &
+                                           ( zu(k) - inflow_damping_height ) / &
+                                           inflow_damping_width
+             ELSE
+                inflow_damping_factor(k) = 0.0
+             ENDIF
+
+          ENDDO
+
+       ENDIF
+
+
+!
+!--    Read binary data from restart file
+          WRITE (9,*) 'before read_3d_binary'
+          CALL local_flush( 9 )
+       CALL read_3d_binary
+          WRITE (9,*) 'after read_3d_binary'
+          CALL local_flush( 9 )
+
+!
+!--    Calculate initial temperature field and other constants used in case
+!--    of a sloping surface
+       IF ( sloping_surface )  CALL init_slope
+
+!
+!--    Initialize new time levels (only done in order to set boundary values
+!--    including ghost points)
+       e_p = e; pt_p = pt; u_p = u; v_p = v; w_p = w
+       IF ( humidity  .OR.  passive_scalar )  q_p = q
+       IF ( ocean )  sa_p = sa
+
+    ELSE
+!
+!--    Actually this part of the programm should not be reached
+       IF ( myid == 0 )  PRINT*,'+++ init_3d_model: unknown initializing ', &
+                                                    'problem'
+       CALL local_stop
+    ENDIF
+
+
+    IF (  TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
+!
+!--    Initialize old timelevels needed for radiation boundary conditions
+       IF ( outflow_l )  THEN
+          u_m_l(:,:,:) = u(:,:,1:2)
+          v_m_l(:,:,:) = v(:,:,0:1)
+          w_m_l(:,:,:) = w(:,:,0:1)
+       ENDIF
+       IF ( outflow_r )  THEN
+          u_m_r(:,:,:) = u(:,:,nx-1:nx)
+          v_m_r(:,:,:) = v(:,:,nx-1:nx)
+          w_m_r(:,:,:) = w(:,:,nx-1:nx)
+       ENDIF
+       IF ( outflow_s )  THEN
+          u_m_s(:,:,:) = u(:,0:1,:)
+          v_m_s(:,:,:) = v(:,1:2,:)
+          w_m_s(:,:,:) = w(:,0:1,:)
+       ENDIF
+       IF ( outflow_n )  THEN
+          u_m_n(:,:,:) = u(:,ny-1:ny,:)
+          v_m_n(:,:,:) = v(:,ny-1:ny,:)
+          w_m_n(:,:,:) = w(:,ny-1:ny,:)
+       ENDIF
+
+!
 !--    Calculate the initial volume flow at the right and north boundary
        IF ( conserve_volume_flow )  THEN
@@ -800 +1018 @@
        ENDIF
 
-!
-!--    For the moment, perturbation pressure and vertical velocity are zero
-       p = 0.0; w = 0.0
-
-!
-!--    Initialize array sums (must be defined in first call of pres)
-       sums = 0.0
-
-!
-!--    Treating cloud physics, liquid water content and precipitation amount
-!--    are zero at beginning of the simulation
-       IF ( cloud_physics )  THEN
-          ql = 0.0
-          IF ( precipitation )  precipitation_amount = 0.0
-       ENDIF
-
-!
-!--    Impose vortex with vertical axis on the initial velocity profile
-       IF ( INDEX( initializing_actions, 'initialize_vortex' ) /= 0 )  THEN
-          CALL init_rankine
-       ENDIF
-
-!
-!--    Impose temperature anomaly (advection test only)
-       IF ( INDEX( initializing_actions, 'initialize_ptanom' ) /= 0 )  THEN
-          CALL init_pt_anomaly
-       ENDIF
-
-!
-!--    If required, change the surface temperature at the start of the 3D run
-       IF ( pt_surface_initial_change /= 0.0 )  THEN
-          pt(nzb,:,:) = pt(nzb,:,:) + pt_surface_initial_change
-       ENDIF
-
-!
-!--    If required, change the surface humidity/scalar at the start of the 3D
-!--    run
-       IF ( ( humidity .OR. passive_scalar ) .AND. &
-            q_surface_initial_change /= 0.0 )  THEN
-          q(nzb,:,:) = q(nzb,:,:) + q_surface_initial_change
-       ENDIF
-
-!
-!--    Initialize the random number generator (from numerical recipes)
-       CALL random_function_ini
-
-!
-!--    Impose random perturbation on the horizontal velocity field and then
-!--    remove the divergences from the velocity field
-       IF ( create_disturbances )  THEN
-          CALL disturb_field( nzb_u_inner, tend, u )
-          CALL disturb_field( nzb_v_inner, tend, v )
-          n_sor = nsor_ini
-          CALL pres
-          n_sor = nsor
-       ENDIF
-
-!
-!--    Once again set the perturbation pressure explicitly to zero in order to
-!--    assure that it does not generate any divergences in the first time step.
-!--    At t=0 the velocity field is free of divergence (as constructed above).
-!--    Divergences being created during a time step are not yet known and thus
-!--    cannot be corrected during the time step yet.
-       p = 0.0
-
-!
-!--    Initialize old and new time levels.
-       IF ( timestep_scheme(1:5) /= 'runge' )  THEN
-          e_m = e; pt_m = pt; u_m = u; v_m = v; w_m = w; kh_m = kh; km_m = km
-       ELSE
-          te_m = 0.0; tpt_m = 0.0; tu_m = 0.0; tv_m = 0.0; tw_m = 0.0
-       ENDIF
-       e_p = e; pt_p = pt; u_p = u; v_p = v; w_p = w
-
-       IF ( humidity  .OR.  passive_scalar )  THEN
-          IF ( ASSOCIATED( q_m ) )  q_m = q
-          IF ( timestep_scheme(1:5) == 'runge' )  tq_m = 0.0
-          q_p = q
-          IF ( humidity  .AND.  ASSOCIATED( vpt_m ) )  vpt_m = vpt
-       ENDIF
-
-       IF ( ocean )  THEN
-          tsa_m = 0.0
-          sa_p  = sa
-       ENDIF
-
-!
-!--    Initialize old timelevels needed for radiation boundary conditions
-       IF ( outflow_l )  THEN
-          u_m_l(:,:,:) = u(:,:,1:2)
-          v_m_l(:,:,:) = v(:,:,0:1)
-          w_m_l(:,:,:) = w(:,:,0:1)
-       ENDIF
-       IF ( outflow_r )  THEN
-          u_m_r(:,:,:) = u(:,:,nx-1:nx)
-          v_m_r(:,:,:) = v(:,:,nx-1:nx)
-          w_m_r(:,:,:) = w(:,:,nx-1:nx)
-       ENDIF
-       IF ( outflow_s )  THEN
-          u_m_s(:,:,:) = u(:,0:1,:)
-          v_m_s(:,:,:) = v(:,1:2,:)
-          w_m_s(:,:,:) = w(:,0:1,:)
-       ENDIF
-       IF ( outflow_n )  THEN
-          u_m_n(:,:,:) = u(:,ny-1:ny,:)
-          v_m_n(:,:,:) = v(:,ny-1:ny,:)
-          w_m_n(:,:,:) = w(:,ny-1:ny,:)
-       ENDIF
-
-    ELSEIF ( TRIM( initializing_actions ) == 'read_restart_data'  .OR.    &
-             TRIM( initializing_actions ) == 'read_data_for_recycling' )  &
-    THEN
-!
-!--    When reading data for initializing the recycling method, first read
-!--    some of the global variables from restart file
-       IF ( TRIM( initializing_actions ) == 'read_data_for_recycling' )  THEN
-          WRITE (9,*) 'before read_parts_of_var_list'
-          CALL local_flush( 9 )
-          CALL read_parts_of_var_list
-          WRITE (9,*) 'after read_parts_of_var_list'
-          CALL local_flush( 9 )
-          CALL close_file( 13 )
-       ENDIF
-
-!
-!--    Read binary data from restart file
-          WRITE (9,*) 'before read_3d_binary'
-          CALL local_flush( 9 )
-       CALL read_3d_binary
-          WRITE (9,*) 'after read_3d_binary'
-          CALL local_flush( 9 )
-
-!
-!--    Calculate initial temperature field and other constants used in case
-!--    of a sloping surface
-       IF ( sloping_surface )  CALL init_slope
-
-!
-!--    Initialize new time levels (only done in order to set boundary values
-!--    including ghost points)
-       e_p = e; pt_p = pt; u_p = u; v_p = v; w_p = w
-       IF ( humidity  .OR.  passive_scalar )  q_p = q
-       IF ( ocean )  sa_p = sa
-
-    ELSE
-!
-!--    Actually this part of the programm should not be reached
-       IF ( myid == 0 )  PRINT*,'+++ init_3d_model: unknown initializing ', &
-                                                    'problem'
-       CALL local_stop
     ENDIF