Changeset 667 for palm/trunk/SOURCE
- Timestamp:
- Dec 23, 2010 12:06:00 PM (12 years ago)
- Location:
- palm/trunk/SOURCE
- Files:
-
- 66 edited
- 1 copied
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE
-
Property
svn:mergeinfo
set to
(toggle deleted branches)
/palm/branches/suehring 423-666 /palm/branches/letzel/masked_output/SOURCE 296-409
-
Property
svn:mergeinfo
set to
(toggle deleted branches)
-
palm/trunk/SOURCE/Makefile
r482 r667  4 4 # Current revisions: 5 5 # ------------------ 6  #  6 # +advec_ws 7 7 # 8 8 # Former revisions: … …  57 57 58 58 RCS = advec_particles.f90 advec_s_bc.f90 advec_s_pw.f90 advec_s_up.f90 \ 59  advec_s_ups.f90 advec_u_pw.f90 advec_u_up.f90 advec_u_ups.f90 \60  advec_v_pw.f90 advec_v_up.f90 advec_v_ups.f90 advec_w_pw.f90 \61  advec_w_up.f90 advec_w_ups.f90 asselin_filter.f90 average_3d_data.f90 \62  boundary_conds.f90 buoyancy.f90 calc_liquid_water_content.f90 \63  calc_precipitation.f90 calc_radiation.f90 calc_spectra.f90 \64  check_for_restart.f90 check_open.f90 check_parameters.f90 \65  close_file.f90 compute_vpt.f90 coriolis.f90 cpu_log.f90 \66  cpu_statistics.f90 data_log.f90 data_output_dvrp.f90 \67  data_output_mask.f90 data_output_profiles.f90 data_output_ptseries.f90 \68  data_output_spectra.f90 data_output_tseries.f90 data_output_2d.f90 \69  data_output_3d.f90 diffusion_e.f90 diffusion_s.f90 diffusion_u.f90 \70  diffusion_v.f90 diffusion_w.f90 diffusivities.f90 disturb_field.f90 \71  disturb_heatflux.f90 eqn_state_seawater.f90 exchange_horiz.f90 \72   59 advec_ws.f90 advec_s_ups.f90 advec_u_pw.f90 advec_u_up.f90 \  60 advec_u_ups.f90 advec_v_pw.f90 advec_v_up.f90 advec_v_ups.f90 \  61 advec_w_pw.f90 advec_w_up.f90 advec_w_ups.f90 asselin_filter.f90 \  62 average_3d_data.f90 boundary_conds.f90 buoyancy.f90 \  63 calc_liquid_water_content.f90 calc_precipitation.f90 \  64 calc_radiation.f90 calc_spectra.f90 check_for_restart.f90 \  65 check_open.f90 check_parameters.f90 close_file.f90 compute_vpt.f90 \  66 coriolis.f90 cpu_log.f90 cpu_statistics.f90 data_log.f90 \  67 data_output_dvrp.f90 data_output_mask.f90 data_output_profiles.f90 \  68 data_output_ptseries.f90 data_output_spectra.f90 data_output_tseries.f90 \  69 data_output_2d.f90 data_output_3d.f90 diffusion_e.f90 diffusion_s.f90 \  70 diffusion_u.f90 diffusion_v.f90 diffusion_w.f90 diffusivities.f90 \  71 disturb_field.f90 disturb_heatflux.f90 eqn_state_seawater.f90 \  72 exchange_horiz.f90 exchange_horiz_2d.f90 \ 73 73 fft_xy.f90 flow_statistics.f90 global_min_max.f90 \ 74 74 header.f90 impact_of_latent_heat.f90 inflow_turbulence.f90 \ … …  77 77 init_masks.f90 init_ocean.f90 init_particles.f90 init_pegrid.f90 \ 78 78 init_pt_anomaly.f90 init_rankine.f90 init_slope.f90 \ 79  interaction_droplets_ptq.f90 local_flush.f90 local_getenv.f90\80   79 interaction_droplets_ptq.f90 local_flush.f90 \  80 local_getenv.f90 local_stop.f90 local_system.f90 local_tremain.f90 \ 81 81 local_tremain_ini.f90 message.f90 modules.f90 netcdf.f90 \ 82 82 package_parin.f90 palm.f90 parin.f90 particle_boundary_conds.f90 \ … …  106 106 OBJS = advec_particles.o advec_s_bc.o advec_s_pw.o advec_s_up.o \ 107 107 advec_s_ups.o advec_u_pw.o advec_u_up.o advec_u_ups.o \ 108   108 advec_ws.o advec_v_pw.o advec_v_up.o advec_v_ups.o advec_w_pw.o \ 109 109 advec_w_up.o advec_w_ups.o asselin_filter.o average_3d_data.o \ 110 110 boundary_conds.o buoyancy.o calc_liquid_water_content.o \ … …  184 184 advec_v_up.o: modules.o 185 185 advec_v_ups.o: modules.o  186 advec_ws.o: modules.o 186 187 advec_w_pw.o: modules.o 187 188 advec_w_up.o: modules.o … …  230 231 inflow_turbulence.o: modules.o 231 232 init_1d_model.o: modules.o 232  init_3d_model.o: modules.o random_function.o  233 init_3d_model.o: modules.o random_function.o advec_ws.o 233 234 init_advec.o: modules.o 234 235 init_cloud_physics.o: modules.o … …  264 265 production_e.o: modules.o wall_fluxes.o 265 266 prognostic_equations.o: modules.o advec_s_pw.o advec_s_up.o advec_u_pw.o \  267 advec_ws.o \ 266 268 advec_u_up.o advec_v_pw.o advec_v_up.o advec_w_pw.o advec_w_up.o \ 267 269 buoyancy.o calc_precipitation.o calc_radiation.o coriolis.o \ -
palm/trunk/SOURCE/advec_particles.f90
r623 r667  4 4 ! Current revisions: 5 5 ! -----------------  6 ! Declaration of de_dx, de_dy, de_dz adapted to additional  7 ! ghost points. Furthermore the calls of exchange_horiz were modified. 6 8 ! 7 9 ! TEST: PRINT statements on unit 9 (commented out) … …  153 155 REAL :: location(1:30,1:3) 154 156 155  REAL, DIMENSION(nzb:nzt+1,nys -1:nyn+1,nxl-1:nxr+1) :: de_dx, de_dy, de_dz 157 REAL, DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: de_dx, de_dy, de_dz 156 158 157 159 REAL, DIMENSION(:,:,:), ALLOCATABLE :: trlpt, trnpt, trrpt, trspt … …  768 770 ! 769 771 !-- Lateral boundary conditions 770  CALL exchange_horiz( de_dx )771  CALL exchange_horiz( de_dy )772  CALL exchange_horiz( de_dz )773  CALL exchange_horiz( diss ) 772 CALL exchange_horiz( de_dx, nbgp )  773 CALL exchange_horiz( de_dy, nbgp )  774 CALL exchange_horiz( de_dz, nbgp )  775 CALL exchange_horiz( diss, nbgp ) 774 776 775 777 ! -
palm/trunk/SOURCE/advec_w_pw.f90
r484 r667  89 89 REAL :: gu, gv 90 90 91  92 91 gu = 2.0 * u_gtrans 93 92 gv = 2.0 * v_gtrans … …  103 102 ) 104 103 ENDDO 105  106 104 END SUBROUTINE advec_w_pw_ij 107 105 -
palm/trunk/SOURCE/asselin_filter.f90
r484 r667  4 4 ! Current revisions: 5 5 ! ----------------- 6  !  6 ! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng 7 7 ! 8 8 ! Former revisions: … …  49 49 !$OMP PARALLEL PRIVATE (i,j,k) 50 50 !$OMP DO 51  DO i = nxl -1, nxr+152  DO j = nys -1, nyn+1 51 DO i = nxlg, nxrg  52 DO j = nysg, nyng 53 53 54 54 DO k = nzb_2d(j,i), nzt+1 -
palm/trunk/SOURCE/average_3d_data.f90
r484 r667  4 4 ! Current revisions: 5 5 ! ----------------- 6  !  6 ! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng 7 7 ! 8 8 ! Former revisions: … …  58 58 59 59 CASE ( 'e' ) 60  DO i = nxl -1, nxr+161  DO j = nys -1, nyn+1 60 DO i = nxlg, nxrg  61 DO j = nysg, nyng 62 62 DO k = nzb, nzt+1 63 63 e_av(k,j,i) = e_av(k,j,i) / REAL( average_count_3d ) … …  67 67 68 68 CASE ( 'qsws*' ) 69  DO i = nxl -1, nxr+170  DO j = nys -1, nyn+1 69 DO i = nxlg, nxrg  70 DO j = nysg, nyng 71 71 qsws_av(j,i) = qsws_av(j,i) / REAL( average_count_3d ) 72 72 ENDDO … …  74 74 75 75 CASE ( 'lwp*' ) 76  DO i = nxl -1, nxr+177  DO j = nys -1, nyn+1 76 DO i = nxlg, nxrg  77 DO j = nysg, nyng 78 78 lwp_av(j,i) = lwp_av(j,i) / REAL( average_count_3d ) 79 79 ENDDO … …  81 81 82 82 CASE ( 'p' ) 83  DO i = nxl -1, nxr+184  DO j = nys -1, nyn+1 83 DO i = nxlg, nxrg  84 DO j = nysg, nyng 85 85 DO k = nzb, nzt+1 86 86 p_av(k,j,i) = p_av(k,j,i) / REAL( average_count_3d ) … …  108 108 109 109 CASE ( 'prr*' ) 110  DO i = nxl -1, nxr+1111  DO j = nys -1, nyn+1 110 DO i = nxlg, nxrg  111 DO j = nysg, nyng 112 112 precipitation_rate_av(j,i) = precipitation_rate_av(j,i) / & 113 113 REAL( average_count_3d ) … …  116 116 117 117 CASE ( 'pt' ) 118  DO i = nxl -1, nxr+1119  DO j = nys -1, nyn+1 118 DO i = nxlg, nxrg  119 DO j = nysg, nyng 120 120 DO k = nzb, nzt+1 121 121 pt_av(k,j,i) = pt_av(k,j,i) / REAL( average_count_3d ) … …  125 125 126 126 CASE ( 'q' ) 127  DO i = nxl -1, nxr+1128  DO j = nys -1, nyn+1 127 DO i = nxlg, nxrg  128 DO j = nysg, nyng 129 129 DO k = nzb, nzt+1 130 130 q_av(k,j,i) = q_av(k,j,i) / REAL( average_count_3d ) … …  134 134 135 135 CASE ( 'ql' ) 136  DO i = nxl -1, nxr+1137  DO j = nys -1, nyn+1 136 DO i = nxlg, nxrg  137 DO j = nysg, nyng 138 138 DO k = nzb, nzt+1 139 139 ql_av(k,j,i) = ql_av(k,j,i) / REAL( average_count_3d ) … …  143 143 144 144 CASE ( 'ql_c' ) 145  DO i = nxl -1, nxr+1146  DO j = nys -1, nyn+1 145 DO i = nxlg, nxrg  146 DO j = nysg, nyng 147 147 DO k = nzb, nzt+1 148 148 ql_c_av(k,j,i) = ql_c_av(k,j,i) / REAL( average_count_3d ) … …  152 152 153 153 CASE ( 'ql_v' ) 154  DO i = nxl -1, nxr+1155  DO j = nys -1, nyn+1 154 DO i = nxlg, nxrg  155 DO j = nysg, nyng 156 156 DO k = nzb, nzt+1 157 157 ql_v_av(k,j,i) = ql_v_av(k,j,i) / REAL( average_count_3d ) … …  161 161 162 162 CASE ( 'ql_vp' ) 163  DO i = nxl -1, nxr+1164  DO j = nys -1, nyn+1 163 DO i = nxlg, nxrg  164 DO j = nysg, nyng 165 165 DO k = nzb, nzt+1 166 166 ql_vp_av(k,j,i) = ql_vp_av(k,j,i) / & … …  171 171 172 172 CASE ( 'qv' ) 173  DO i = nxl -1, nxr+1174  DO j = nys -1, nyn+1 173 DO i = nxlg, nxrg  174 DO j = nysg, nyng 175 175 DO k = nzb, nzt+1 176 176 qv_av(k,j,i) = qv_av(k,j,i) / REAL( average_count_3d ) … …  180 180 181 181 CASE ( 'rho' ) 182  DO i = nxl -1, nxr+1183  DO j = nys -1, nyn+1 182 DO i = nxlg, nxrg  183 DO j = nysg, nyng 184 184 DO k = nzb, nzt+1 185 185 rho_av(k,j,i) = rho_av(k,j,i) / REAL( average_count_3d ) … …  189 189 190 190 CASE ( 's' ) 191  DO i = nxl -1, nxr+1192  DO j = nys -1, nyn+1 191 DO i = nxlg, nxrg  192 DO j = nysg, nyng 193 193 DO k = nzb, nzt+1 194 194 s_av(k,j,i) = s_av(k,j,i) / REAL( average_count_3d ) … …  198 198 199 199 CASE ( 'sa' ) 200  DO i = nxl -1, nxr+1201  DO j = nys -1, nyn+1 200 DO i = nxlg, nxrg  201 DO j = nysg, nyng 202 202 DO k = nzb, nzt+1 203 203 sa_av(k,j,i) = sa_av(k,j,i) / REAL( average_count_3d ) … …  207 207 208 208 CASE ( 'shf*' ) 209  DO i = nxl -1, nxr+1210  DO j = nys -1, nyn+1 209 DO i = nxlg, nxrg  210 DO j = nysg, nyng 211 211 shf_av(j,i) = shf_av(j,i) / REAL( average_count_3d ) 212 212 ENDDO … …  214 214 215 215 CASE ( 't*' ) 216  DO i = nxl -1, nxr+1217  DO j = nys -1, nyn+1 216 DO i = nxlg, nxrg  217 DO j = nysg, nyng 218 218 ts_av(j,i) = ts_av(j,i) / REAL( average_count_3d ) 219 219 ENDDO … …  221 221 222 222 CASE ( 'u' ) 223  DO i = nxl -1, nxr+1224  DO j = nys -1, nyn+1 223 DO i = nxlg, nxrg  224 DO j = nysg, nyng 225 225 DO k = nzb, nzt+1 226 226 u_av(k,j,i) = u_av(k,j,i) / REAL( average_count_3d ) … …  230 230 231 231 CASE ( 'u*' ) 232  DO i = nxl -1, nxr+1233  DO j = nys -1, nyn+1 232 DO i = nxlg, nxrg  233 DO j = nysg, nyng 234 234 us_av(j,i) = us_av(j,i) / REAL( average_count_3d ) 235 235 ENDDO … …  237 237 238 238 CASE ( 'v' ) 239  DO i = nxl -1, nxr+1240  DO j = nys -1, nyn+1 239 DO i = nxlg, nxrg  240 DO j = nysg, nyng 241 241 DO k = nzb, nzt+1 242 242 v_av(k,j,i) = v_av(k,j,i) / REAL( average_count_3d ) … …  246 246 247 247 CASE ( 'vpt' ) 248  DO i = nxl -1, nxr+1249  DO j = nys -1, nyn+1 248 DO i = nxlg, nxrg  249 DO j = nysg, nyng 250 250 DO k = nzb, nzt+1 251 251 vpt_av(k,j,i) = vpt_av(k,j,i) / REAL( average_count_3d ) … …  255 255 256 256 CASE ( 'w' ) 257  DO i = nxl -1, nxr+1258  DO j = nys -1, nyn+1 257 DO i = nxlg, nxrg  258 DO j = nysg, nyng 259 259 DO k = nzb, nzt+1 260 260 w_av(k,j,i) = w_av(k,j,i) / REAL( average_count_3d ) … …  264 264 265 265 CASE ( 'z0*' ) 266  DO i = nxl -1, nxr+1267  DO j = nys -1, nyn+1 266 DO i = nxlg, nxrg  267 DO j = nysg, nyng 268 268 z0_av(j,i) = z0_av(j,i) / REAL( average_count_3d ) 269 269 ENDDO -
palm/trunk/SOURCE/boundary_conds.f90
r484 r667  4 4 ! Current revisions: 5 5 ! -----------------  6 !  7 ! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng  8 ! 6 9 ! 7  !  10 ! Removed mirror boundary conditions for u and v at the bottom in case of  11 ! ibc_uv_b == 0. Instead, dirichelt boundary conditions (u=v=0) are set  12 ! in init_3d_model  13 ! 8 14 ! Former revisions: 9 15 ! ----------------- … …  70 76 IF ( range == 'main') THEN 71 77 ! 72  !-- Bottom boundary 73  IF ( ibc_uv_b == 0 ) THEN 74  ! 75  !-- Satisfying the Dirichlet condition with an extra layer below the 76  !-- surface where the u and v component change their sign 77  u_p(nzb,:,:) = -u_p(nzb+1,:,:) 78  v_p(nzb,:,:) = -v_p(nzb+1,:,:) 79  ELSE  78 !-- Bottom boundary  79 IF ( ibc_uv_b == 1 ) THEN 80 80 u_p(nzb,:,:) = u_p(nzb+1,:,:) 81 81 v_p(nzb,:,:) = v_p(nzb+1,:,:) 82 82 ENDIF 83  DO i = nxl -1, nxr+184  DO j = nys -1, nyn+1 83 DO i = nxlg, nxrg  84 DO j = nysg, nyng 85 85 w_p(nzb_w_inner(j,i),j,i) = 0.0 86 86 ENDDO … …  90 90 !-- Top boundary 91 91 IF ( ibc_uv_t == 0 ) THEN 92  u_p(nzt+1,:,:) = ug(nzt+1)93  v_p(nzt+1,:,:) = vg(nzt+1) 92 u_p(nzt+1,:,:) = ug(nzt+1)  93 v_p(nzt+1,:,:) = vg(nzt+1) 94 94 ELSE 95  u_p(nzt+1,:,:) = u_p(nzt,:,:)96  v_p(nzt+1,:,:) = v_p(nzt,:,:) 95 u_p(nzt+1,:,:) = u_p(nzt,:,:)  96 v_p(nzt+1,:,:) = v_p(nzt,:,:) 97 97 ENDIF 98 98 w_p(nzt:nzt+1,:,:) = 0.0 ! nzt is not a prognostic level (but cf. pres) … …  103 103 !-- the sea surface temperature of the coupled ocean model. 104 104 IF ( ibc_pt_b == 0 ) THEN 105  DO i = nxl -1, nxr+1106  DO j = nys -1, nyn+1 105 DO i = nxlg, nxrg  106 DO j = nysg, nyng 107 107 pt_p(nzb_s_inner(j,i),j,i) = pt(nzb_s_inner(j,i),j,i) 108 108 ENDDO 109 109 ENDDO 110 110 ELSEIF ( ibc_pt_b == 1 ) THEN 111  DO i = nxl -1, nxr+1112  DO j = nys -1, nyn+1 111 DO i = nxlg, nxrg  112 DO j = nysg, nyng 113 113 pt_p(nzb_s_inner(j,i),j,i) = pt_p(nzb_s_inner(j,i)+1,j,i) 114 114 ENDDO … …  119 119 !-- Temperature at top boundary 120 120 IF ( ibc_pt_t == 0 ) THEN 121  pt_p(nzt+1,:,:) = pt(nzt+1,:,:) 121 pt_p(nzt+1,:,:) = pt(nzt+1,:,:) 122 122 ELSEIF ( ibc_pt_t == 1 ) THEN 123  pt_p(nzt+1,:,:) = pt_p(nzt,:,:) 123 pt_p(nzt+1,:,:) = pt_p(nzt,:,:) 124 124 ELSEIF ( ibc_pt_t == 2 ) THEN 125  pt_p(nzt+1,:,:) = pt_p(nzt,:,:) + bc_pt_t_val * dzu(nzt+1) 125 pt_p(nzt+1,:,:) = pt_p(nzt,:,:) + bc_pt_t_val * dzu(nzt+1) 126 126 ENDIF 127 127 … …  130 130 !-- Generally Neumann conditions with de/dz=0 are assumed 131 131 IF ( .NOT. constant_diffusion ) THEN 132  DO i = nxl -1, nxr+1133  DO j = nys -1, nyn+1 132 DO i = nxlg, nxrg  133 DO j = nysg, nyng 134 134 e_p(nzb_s_inner(j,i),j,i) = e_p(nzb_s_inner(j,i)+1,j,i) 135 135 ENDDO … …  144 144 !-- Bottom boundary: Neumann condition because salinity flux is always 145 145 !-- given 146  DO i = nxl -1, nxr+1147  DO j = nys -1, nyn+1 146 DO i = nxlg, nxrg  147 DO j = nysg, nyng 148 148 sa_p(nzb_s_inner(j,i),j,i) = sa_p(nzb_s_inner(j,i)+1,j,i) 149 149 ENDDO … …  153 153 !-- Top boundary: Dirichlet or Neumann 154 154 IF ( ibc_sa_t == 0 ) THEN 155  sa_p(nzt+1,:,:) = sa(nzt+1,:,:) 155 sa_p(nzt+1,:,:) = sa(nzt+1,:,:) 156 156 ELSEIF ( ibc_sa_t == 1 ) THEN 157  sa_p(nzt+1,:,:) = sa_p(nzt,:,:) 157 sa_p(nzt+1,:,:) = sa_p(nzt,:,:) 158 158 ENDIF 159 159 … …  167 167 !-- Surface conditions for constant_humidity_flux 168 168 IF ( ibc_q_b == 0 ) THEN 169  DO i = nxl -1, nxr+1170  DO j = nys -1, nyn+1 169 DO i = nxlg, nxrg  170 DO j = nysg, nyng 171 171 q_p(nzb_s_inner(j,i),j,i) = q(nzb_s_inner(j,i),j,i) 172 172 ENDDO 173 173 ENDDO 174 174 ELSE 175  DO i = nxl -1, nxr+1176  DO j = nys -1, nyn+1 175 DO i = nxlg, nxrg  176 DO j = nysg, nyng 177 177 q_p(nzb_s_inner(j,i),j,i) = q_p(nzb_s_inner(j,i)+1,j,i) 178 178 ENDDO … …  182 182 !-- Top boundary 183 183 q_p(nzt+1,:,:) = q_p(nzt,:,:) + bc_q_t_val * dzu(nzt+1)  184  185 184 186 ENDIF 185 187 … …  226 228 c_max = dy / dt_3d 227 229 228  DO i = nxl -1, nxr+1 230 DO i = nxlg, nxrg 229 231 DO k = nzb+1, nzt+1 230 232 … …  299 301 !-- Bottom boundary at the outflow 300 302 IF ( ibc_uv_b == 0 ) THEN 301  u_p(nzb,-1,:) = -u_p(nzb+1,-1,:)302  v_p(nzb,0,:) = -v_p(nzb+1,0,:) 303 u_p(nzb,-1,:) = 0.0  304 v_p(nzb,0,:) = 0.0 303 305 ELSE 304 306 u_p(nzb,-1,:) = u_p(nzb+1,-1,:) … …  324 326 c_max = dy / dt_3d 325 327 326  DO i = nxl -1, nxr+1 328 DO i = nxlg, nxrg 327 329 DO k = nzb+1, nzt+1 328 330 … …  397 399 !-- Bottom boundary at the outflow 398 400 IF ( ibc_uv_b == 0 ) THEN 399  u_p(nzb,ny+1,:) = -u_p(nzb+1,ny+1,:)400  v_p(nzb,ny+1,:) = -v_p(nzb+1,ny+1,:) 401 u_p(nzb,ny+1,:) = 0.0  402 v_p(nzb,ny+1,:) = 0.0 401 403 ELSE 402 404 u_p(nzb,ny+1,:) = u_p(nzb+1,ny+1,:) … …  422 424 c_max = dx / dt_3d 423 425 424  DO j = nys -1, nyn+1 426 DO j = nysg, nyng 425 427 DO k = nzb+1, nzt+1 426 428 … …  495 497 !-- Bottom boundary at the outflow 496 498 IF ( ibc_uv_b == 0 ) THEN 497  u_p(nzb,:, -1) = -u_p(nzb+1,:,-1)498  v_p(nzb,:,-1) = -v_p(nzb+1,:,-1) 499 u_p(nzb,:,0) = 0.0  500 v_p(nzb,:,-1) = 0.0 499 501 ELSE 500  u_p(nzb,:, -1) = u_p(nzb+1,:,-1) 502 u_p(nzb,:,0) = u_p(nzb+1,:,0) 501 503 v_p(nzb,:,-1) = v_p(nzb+1,:,-1) 502 504 ENDIF … …  520 522 c_max = dx / dt_3d 521 523 522  DO j = nys -1, nyn+1 524 DO j = nysg, nyng 523 525 DO k = nzb+1, nzt+1 524 526 … …  593 595 !-- Bottom boundary at the outflow 594 596 IF ( ibc_uv_b == 0 ) THEN 595  u_p(nzb,:,nx+1) = -u_p(nzb+1,:,nx+1)596  v_p(nzb,:,nx+1) = -v_p(nzb+1,:,nx+1) 597 u_p(nzb,:,nx+1) = 0.0  598 v_p(nzb,:,nx+1) = 0.0 597 599 ELSE 598 600 u_p(nzb,:,nx+1) = u_p(nzb+1,:,nx+1) -
palm/trunk/SOURCE/calc_liquid_water_content.f90
r484 r667  4 4 ! Current revisions: 5 5 ! ----------------- 6  !  6 ! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng 7 7 ! 8 8 ! Former revisions: … …  47 47 REAL :: alpha, e_s, q_s, t_l 48 48 49  DO i = nxl -1, nxr+150  DO j = nys -1, nyn+1 49 DO i = nxlg, nxrg  50 DO j = nysg, nyng 51 51 DO k = nzb_2d(j,i)+1, nzt 52 52 -
palm/trunk/SOURCE/calc_precipitation.f90
r484 r667 Â 8 8 ! Former revisions: 9 9 ! ----------------- 10 Â ! $Id $Â 10 ! $Id: calc_precipitation.f90 484 2010-02-05 07:36:54Z raasch 11 11 ! 12 12 ! 403 2009-10-22 13:57:16Z franke -
palm/trunk/SOURCE/calc_spectra.f90
r392 r667  4 4 ! Current revisions: 5 5 ! ----------------- 6  !  6 ! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng for allocation  7 ! of tend 7 8 ! 8 9 ! Former revisions: … …  152 153 IF ( nxra > nxr .OR. nyna > nyn .OR. nza > nz ) THEN 153 154 DEALLOCATE( tend ) 154  ALLOCATE( tend(nzb:nzt+1,nys -1:nyn+1,nxl-1:nxr+1) ) 155 ALLOCATE( tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 155 156 ENDIF 156 157 -
palm/trunk/SOURCE/check_for_restart.f90
r623 r667  5 5 ! ----------------- 6 6 !  7 ! Exchange of terminate_coupled between ocean and atmosphere by PE0 7 8 ! 8 9 ! Former revisions: … …  93 94 terminate_coupled = 3 94 95 #if defined( __parallel ) 95  CALL MPI_SENDRECV( terminate_coupled, 1, MPI_INTEGER, & 96  target_id, 0, & 97  terminate_coupled_remote, 1, MPI_INTEGER, & 98  target_id, 0, & 99  comm_inter, status, ierr )  96 IF ( myid == 0 ) THEN  97 CALL MPI_SENDRECV( terminate_coupled, 1, MPI_INTEGER, &  98 target_id, 0, &  99 terminate_coupled_remote, 1, MPI_INTEGER, &  100 target_id, 0, &  101 comm_inter, status, ierr )  102 ENDIF  103 CALL MPI_BCAST( terminate_coupled_remote, 1, MPI_INTEGER, 0, comm2d, ierr) 100 104 #endif 101 105 ENDIF … …  140 144 ENDIF 141 145 #if defined( __parallel ) 142  CALL MPI_SENDRECV( terminate_coupled, 1, MPI_INTEGER, & 143  target_id, 0, & 144  terminate_coupled_remote, 1, MPI_INTEGER, & 145  target_id, 0, & 146  comm_inter, status, ierr )  146 IF ( myid == 0 ) THEN  147 CALL MPI_SENDRECV( terminate_coupled, 1, MPI_INTEGER, &  148 target_id, 0, &  149 terminate_coupled_remote, 1, MPI_INTEGER, &  150 target_id, 0, &  151 comm_inter, status, ierr )  152 ENDIF  153 CALL MPI_BCAST( terminate_coupled_remote, 1, MPI_INTEGER, 0, comm2d, ierr)  154 147 155 #endif 148 156 ENDIF -
palm/trunk/SOURCE/check_open.f90
r601 r667  4 4 ! Current revisions: 5 5 ! ----------------- 6  !  6 ! Output of total array size was adapted to nbgp. 7 7 ! 8 8 ! Former revisions: … …  278 278 !-- Output for combine_plot_fields 279 279 IF ( data_output_2d_on_each_pe .AND. myid_char /= '' ) THEN 280  WRITE (21) - 1, nx+1, -1, ny+1! total array size 280 WRITE (21) -nbgp, nx+nbgp, -nbgp, ny+nbgp ! total array size 281 281 WRITE (21) 0, nx+1, 0, ny+1 ! output part 282 282 ENDIF … …  319 319 !-- Output for combine_plot_fields 320 320 IF ( data_output_2d_on_each_pe .AND. myid_char /= '' ) THEN 321  WRITE (22) - 1, nx+1, 0, nz+1 ! total array size 321 WRITE (22) -nbgp, nx+nbgp, 0, nz+1 ! total array size 322 322 WRITE (22) 0, nx+1, 0, nz+1 ! output part 323 323 ENDIF … …  357 357 !-- Output for combine_plot_fields 358 358 IF ( data_output_2d_on_each_pe .AND. myid_char /= '' ) THEN 359  WRITE (23) - 1, ny+1, 0, nz+1 ! total array size 359 WRITE (23) -nbgp, ny+nbgp, 0, nz+1 ! total array size 360 360 WRITE (23) 0, ny+1, 0, nz+1 ! output part 361 361 ENDIF … …  392 392 !-- Specifications for combine_plot_fields 393 393 IF ( .NOT. do3d_compress ) THEN 394  WRITE ( 30 ) - 1,nx+1,-1,ny+1,0,nz_do3d 394 WRITE ( 30 ) -nbgp,nx+nbgp,-nbgp,ny+nbgp, 0 ,nz_do3d 395 395 WRITE ( 30 ) 0,nx+1,0,ny+1,0,nz_do3d 396 396 ENDIF … …  503 503 openfile(33)%opened = .TRUE. 504 504 WRITE ( 33, 3300 ) TRIM( avs_coor_file ), & 505  TRIM( avs_coor_file ), (nx+2 )*4, &506  TRIM( avs_coor_file ), (nx+2 )*4+(ny+2)*4 505 TRIM( avs_coor_file ), (nx+2*nbgp)*4, &  506 TRIM( avs_coor_file ), (nx+2*nbgp)*4+(ny+2*nbgp)*4 507 507 508 508 … …  548 548 !-- corresponding partial array of a PE only once at the top of the file 549 549 IF ( avs_output .AND. do3d_compress ) THEN 550  WRITE ( 30 ) nxl -1, nxr+1, nys-1, nyn+1, nzb, nz_do3d 550 WRITE ( 30 ) nxlg, nxrg, nysg, nyng, nzb, nz_do3d 551 551 ENDIF 552 552 … …  929 929 #endif 930 930 ENDIF 931   931 932 932 CALL handle_netcdf_error( 'check_open', 26 ) 933 933 ! … …  1325 1325 #endif 1326 1326 ENDIF 1327   1327 1328 1328 CALL handle_netcdf_error( 'check_open', 43 ) 1329 1329 -
palm/trunk/SOURCE/check_parameters.f90
r601 r667  4 4 ! Current revisions: 5 5 ! -----------------  6 ! 6 7 !  8 ! Exchange of parameters between ocean and atmosphere via PE0  9 ! Check for illegal combination of ws-scheme and timestep scheme.  10 ! Check for topography and ws-scheme.  11 ! Check for not cyclic boundary conditions in combination with ws-scheme and  12 ! loop_optimization = 'vector'.  13 ! Check for call_psolver_at_all_substeps and ws-scheme for momentum_advec.  14 !  15 ! Different processor/grid topology in atmosphere and ocean is now allowed!  16 !  17 ! Bugfixes in checking for conserve_volume_flow_mode  18 !  19 ! Adapt error messages. 7 20 ! 8 21 ! Former revisions: … …  180 193 ! 181 194 !-- Check dt_coupling, restart_time, dt_restart, end_time, dx, dy, nx and ny 182  IF ( coupling_mode /= 'uncoupled'  195 IF ( coupling_mode /= 'uncoupled') THEN 183 196 184 197 IF ( dt_coupling == 9999999.9 ) THEN … …  189 202 190 203 #if defined( __parallel ) 191  CALL MPI_SEND( dt_coupling, 1, MPI_REAL, target_id, 11, comm_inter, & 192  ierr ) 193  CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 11, comm_inter, & 194  status, ierr )  204 IF ( myid == 0 ) THEN  205 CALL MPI_SEND( dt_coupling, 1, MPI_REAL, target_id, 11, comm_inter, &  206 ierr )  207 CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 11, comm_inter, &  208 status, ierr )  209 ENDIF  210 CALL MPI_BCAST( remote, 1, MPI_REAL, 0, comm2d, ierr)  211 195 212 IF ( dt_coupling /= remote ) THEN 196 213 WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ), & … …  200 217 ENDIF 201 218 IF ( dt_coupling <= 0.0 ) THEN 202  CALL MPI_SEND( dt_max, 1, MPI_REAL, target_id, 19, comm_inter, ierr ) 203  CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 19, comm_inter, & 204  status, ierr )  219 IF ( myid == 0 ) THEN  220 CALL MPI_SEND( dt_max, 1, MPI_REAL, target_id, 19, comm_inter, ierr )  221 CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 19, comm_inter, &  222 status, ierr )  223 ENDIF  224 CALL MPI_BCAST( remote, 1, MPI_REAL, 0, comm2d, ierr)  225 205 226 dt_coupling = MAX( dt_max, remote ) 206 227 WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ), & … …  209 230 CALL message( 'check_parameters', 'PA0005', 0, 1, 0, 6, 0 ) 210 231 ENDIF 211  212  CALL MPI_SEND( restart_time, 1, MPI_REAL, target_id, 12, comm_inter, & 213  ierr ) 214  CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 12, comm_inter, & 215  status, ierr )  232 IF ( myid == 0 ) THEN  233 CALL MPI_SEND( restart_time, 1, MPI_REAL, target_id, 12, comm_inter, &  234 ierr )  235 CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 12, comm_inter, &  236 status, ierr )  237 ENDIF  238 CALL MPI_BCAST( remote, 1, MPI_REAL, 0, comm2d, ierr)  239 216 240 IF ( restart_time /= remote ) THEN 217 241 WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ), & … …  220 244 CALL message( 'check_parameters', 'PA0006', 1, 2, 0, 6, 0 ) 221 245 ENDIF 222  223  CALL MPI_SEND( dt_restart, 1, MPI_REAL, target_id, 13, comm_inter, & 224  ierr ) 225  CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 13, comm_inter, & 226  status, ierr )  246 IF ( myid == 0 ) THEN  247 CALL MPI_SEND( dt_restart, 1, MPI_REAL, target_id, 13, comm_inter, &  248 ierr )  249 CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 13, comm_inter, &  250 status, ierr )  251 ENDIF  252 CALL MPI_BCAST( remote, 1, MPI_REAL, 0, comm2d, ierr)  253 227 254 IF ( dt_restart /= remote ) THEN 228 255 WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ), & … …  233 260 234 261 simulation_time_since_reference = end_time - coupling_start_time 235  CALL MPI_SEND( simulation_time_since_reference, 1, MPI_REAL, target_id, & 236  14, comm_inter, ierr ) 237  CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 14, comm_inter, & 238  status, ierr )  262 IF ( myid == 0 ) THEN  263 CALL MPI_SEND( simulation_time_since_reference, 1, MPI_REAL, target_id, &  264 14, comm_inter, ierr )  265 CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 14, comm_inter, &  266 status, ierr )  267 ENDIF  268 CALL MPI_BCAST( remote, 1, MPI_REAL, 0, comm2d, ierr)  269 239 270 IF ( simulation_time_since_reference /= remote ) THEN 240 271 WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ), & … …  245 276 ENDIF 246 277 247  CALL MPI_SEND( dx, 1, MPI_REAL, target_id, 15, comm_inter, ierr ) 248  CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 15, comm_inter, & 249  status, ierr ) 250  IF ( dx /= remote ) THEN 251  WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ), & 252  '": dx = ', dx, '& is not equal to dx_remote = ', remote 253  CALL message( 'check_parameters', 'PA0009', 1, 2, 0, 6, 0 ) 254  ENDIF 255  256  CALL MPI_SEND( dy, 1, MPI_REAL, target_id, 16, comm_inter, ierr ) 257  CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 16, comm_inter, & 258  status, ierr ) 259  IF ( dy /= remote ) THEN 260  WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ), & 261  '": dy = ', dy, '& is not equal to dy_remote = ', remote 262  CALL message( 'check_parameters', 'PA0010', 1, 2, 0, 6, 0 ) 263  ENDIF 264  265  CALL MPI_SEND( nx, 1, MPI_INTEGER, target_id, 17, comm_inter, ierr ) 266  CALL MPI_RECV( iremote, 1, MPI_INTEGER, target_id, 17, comm_inter, & 267  status, ierr ) 268  IF ( nx /= iremote ) THEN 269  WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ), & 270  '": nx = ', nx, '& is not equal to nx_remote = ', iremote 271  CALL message( 'check_parameters', 'PA0011', 1, 2, 0, 6, 0 ) 272  ENDIF 273  274  CALL MPI_SEND( ny, 1, MPI_INTEGER, target_id, 18, comm_inter, ierr ) 275  CALL MPI_RECV( iremote, 1, MPI_INTEGER, target_id, 18, comm_inter, & 276  status, ierr ) 277  IF ( ny /= iremote ) THEN 278  WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ), & 279  '": ny = ', ny, '& is not equal to ny_remote = ', iremote 280  CALL message( 'check_parameters', 'PA0012', 1, 2, 0, 6, 0 )  278  279 IF ( myid == 0 ) THEN  280 CALL MPI_SEND( dx, 1, MPI_REAL, target_id, 15, comm_inter, ierr )  281 CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 15, comm_inter, &  282 status, ierr )  283 ENDIF  284 CALL MPI_BCAST( remote, 1, MPI_REAL, 0, comm2d, ierr)  285  286  287 IF ( coupling_mode == 'atmosphere_to_ocean') THEN  288  289 IF ( dx < remote ) THEN  290 WRITE( message_string, * ) 'coupling mode "', &  291 TRIM( coupling_mode ), &  292 '": dx in Atmosphere is not equal to or not larger then dx in ocean'  293 CALL message( 'check_parameters', 'PA0009', 1, 2, 0, 6, 0 )  294 ENDIF  295  296 IF ( (nx_a+1)*dx /= (nx_o+1)*remote ) THEN  297 WRITE( message_string, * ) 'coupling mode "', &  298 TRIM( coupling_mode ), &  299 '": Domain size in x-direction is not equal in ocean and atmosphere'  300 CALL message( 'check_parameters', 'PA0010', 1, 2, 0, 6, 0 )  301 ENDIF  302  303 ENDIF  304  305 IF ( myid == 0) THEN  306 CALL MPI_SEND( dy, 1, MPI_REAL, target_id, 16, comm_inter, ierr )  307 CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 16, comm_inter, &  308 status, ierr )  309 ENDIF  310 CALL MPI_BCAST( remote, 1, MPI_REAL, 0, comm2d, ierr)  311  312 IF ( coupling_mode == 'atmosphere_to_ocean') THEN  313  314 IF ( dy < remote ) THEN  315 WRITE( message_string, * ) 'coupling mode "', &  316 TRIM( coupling_mode ), &  317 '": dy in Atmosphere is not equal to or not larger then dy in ocean'  318 CALL message( 'check_parameters', 'PA0011', 1, 2, 0, 6, 0 )  319 ENDIF  320  321 IF ( (ny_a+1)*dy /= (ny_o+1)*remote ) THEN  322 WRITE( message_string, * ) 'coupling mode "', &  323 TRIM( coupling_mode ), &  324 '": Domain size in y-direction is not equal in ocean and atmosphere'  325 CALL message( 'check_parameters', 'PA0012', 1, 2, 0, 6, 0 )  326 ENDIF  327  328 IF ( MOD(nx_o+1,nx_a+1) /= 0 ) THEN  329 WRITE( message_string, * ) 'coupling mode "', &  330 TRIM( coupling_mode ), &  331 '": nx+1 in ocean is not divisible without remainder with nx+1 in', &  332 ' atmosphere'  333 CALL message( 'check_parameters', 'PA0339', 1, 2, 0, 6, 0 )  334 ENDIF  335  336 IF ( MOD(ny_o+1,ny_a+1) /= 0 ) THEN  337 WRITE( message_string, * ) 'coupling mode "', &  338 TRIM( coupling_mode ), &  339 '": ny+1 in ocean is not divisible without remainder with ny+1 in', &  340 ' atmosphere'  341 CALL message( 'check_parameters', 'PA0340', 1, 2, 0, 6, 0 )  342 ENDIF  343 281 344 ENDIF 282 345 #else … …  290 353 ! 291 354 !-- Exchange via intercommunicator 292  IF ( coupling_mode == 'atmosphere_to_ocean' ) THEN 355 IF ( coupling_mode == 'atmosphere_to_ocean' .AND. myid == 0 ) THEN 293 356 CALL MPI_SEND( humidity, 1, MPI_LOGICAL, target_id, 19, comm_inter, & 294 357 ierr ) 295  ELSEIF ( coupling_mode == 'ocean_to_atmosphere' ) THEN 358 ELSEIF ( coupling_mode == 'ocean_to_atmosphere' .AND. myid == 0) THEN 296 359 CALL MPI_RECV( humidity_remote, 1, MPI_LOGICAL, target_id, 19, & 297 360 comm_inter, status, ierr ) 298 361 ENDIF  362 CALL MPI_BCAST( humidity_remote, 1, MPI_LOGICAL, 0, comm2d, ierr)  363 299 364 #endif 300 365 … …  372 437 CALL message( 'check_parameters', 'PA0014', 1, 2, 0, 6, 0 ) 373 438 ENDIF  439 IF ( momentum_advec == 'ws-scheme' .OR. scalar_advec == 'ws-scheme' ) &  440 THEN  441 message_string = 'topography is still not allowed with ' // &  442 'momentum_advec = "' // TRIM( momentum_advec ) // &  443 '"or scalar_advec = "' // TRIM( scalar_advec ) //'"'  444 ! message number still needs modification  445 CALL message( 'check_parameters', 'PA0341', 1, 2, 0, 6, 0 )  446 END IF  447 374 448 ! 375 449 !-- In case of non-flat topography, check whether the convention how to … …  492 566 CALL message( 'check_parameters', 'PA0021', 1, 2, 0, 6, 0 ) 493 567 ENDIF 494   568  569 IF( momentum_advec == 'ws-scheme' .AND. &  570 call_psolver_at_all_substeps == .FALSE. ) THEN  571 message_string = 'psolver must be called at each RK3 substep when "'//&  572 TRIM(momentum_advec) // ' "is used for momentum_advec'  573 CALL message( 'check_parameters', 'PA0343', 1, 2, 0, 6, 0 )  574 END IF 495 575 ! 496 576 !-- Advection schemes: 497  IF ( momentum_advec /= 'pw-scheme' .AND. momentum_advec /= ' ups-scheme' )&498  THEN 577 IF ( momentum_advec /= 'pw-scheme' .AND. momentum_advec /= 'ws-scheme' .AND. &  578 momentum_advec /= 'ups-scheme' ) THEN 499 579 message_string = 'unknown advection scheme: momentum_advec = "' // & 500 580 TRIM( momentum_advec ) // '"' 501 581 CALL message( 'check_parameters', 'PA0022', 1, 2, 0, 6, 0 ) 502 582 ENDIF 503  IF ( ( momentum_advec == 'ups-scheme' .OR. scalar_advec == 'ups-scheme' )& 504  .AND. timestep_scheme /= 'euler' ) THEN 505  message_string = 'momentum_advec = "' // TRIM( momentum_advec ) // & 506  '" is not allowed with timestep_scheme = "' // & 507  TRIM( timestep_scheme ) // '"'  583 IF ((( momentum_advec == 'ups-scheme' .OR. scalar_advec == 'ups-scheme' )&  584 .AND. timestep_scheme /= 'euler' ) .OR. (( momentum_advec == 'ws-scheme'&  585 .OR. scalar_advec == 'ws-scheme') .AND. (timestep_scheme == 'euler' .OR. &  586 timestep_scheme == 'leapfrog+euler' .OR. timestep_scheme == 'leapfrog' &  587 .OR. timestep_scheme == 'runge-kutta-2'))) THEN  588 message_string = 'momentum_advec or scalar_advec = "' &  589 // TRIM( momentum_advec ) // '" is not allowed with timestep_scheme = "' // &  590 TRIM( timestep_scheme ) // '"' 508 591 CALL message( 'check_parameters', 'PA0023', 1, 2, 0, 6, 0 ) 509 592 ENDIF 510  511  IF ( scalar_advec /= 'pw-scheme' .AND. scalar_advec /= 'bc-scheme' .AND.& 512  scalar_advec /= 'ups-scheme' ) THEN  593 IF ( scalar_advec /= 'pw-scheme' .AND. scalar_advec /= 'ws-scheme' .AND. &  594 scalar_advec /= 'bc-scheme' .AND. scalar_advec /= 'ups-scheme' ) THEN 513 595 message_string = 'unknown advection scheme: scalar_advec = "' // & 514 596 TRIM( scalar_advec ) // '"' … …  563 645 ENDIF 564 646 565  IF ( momentum_advec /= 'pw-scheme' .AND. timestep_scheme(1:5) == 'runge') &566  THEN 647 IF ( (momentum_advec /= 'pw-scheme' .AND. momentum_advec /= 'ws-scheme') &  648 .AND. timestep_scheme(1:5) == 'runge' ) THEN 567 649 message_string = 'momentum advection scheme "' // & 568 650 TRIM( momentum_advec ) // '" & does not work with ' // & … …  712 794 ug_vertical_gradient_level_ind(1) = nzt+1 713 795 ug(nzt+1) = ug_surface 714  DO k = nzt, 0, -1 796 DO k = nzt, nzb, -1 715 797 IF ( i < 11 ) THEN 716 798 IF ( ug_vertical_gradient_level(i) > zu(k) .AND. & … …  778 860 vg_vertical_gradient_level_ind(1) = nzt+1 779 861 vg(nzt+1) = vg_surface 780  DO k = nzt, 0, -1 862 DO k = nzt, nzb, -1 781 863 IF ( i < 11 ) THEN 782 864 IF ( vg_vertical_gradient_level(i) > zu(k) .AND. & … …  1020 1102 1021 1103  1104 1022 1105 ! 1023 1106 !-- Compute Coriolis parameter … …  1159 1242 ! 1160 1243 !-- Non-cyclic lateral boundaries require the multigrid method and Piascek- 1161  !-- Willimas advection scheme. Several schemes and tools do not work with1162  !-- non-cyclic boundary conditions. 1244 !-- Willimas or Wicker - Skamarock advection scheme. Several schemes  1245 !-- and tools do not work with non-cyclic boundary conditions. 1163 1246 IF ( bc_lr /= 'cyclic' .OR. bc_ns /= 'cyclic' ) THEN 1164 1247 IF ( psolver /= 'multigrid' ) THEN … …  1167 1250 CALL message( 'check_parameters', 'PA0051', 1, 2, 0, 6, 0 ) 1168 1251 ENDIF 1169  IF ( momentum_advec /= 'pw-scheme' ) THEN  1252 IF ( momentum_advec /= 'pw-scheme' .AND. &  1253 momentum_advec /= 'ws-scheme') THEN 1170 1254 message_string = 'non-cyclic lateral boundaries do not allow ' // & 1171 1255 'momentum_advec = "' // TRIM( momentum_advec ) // '"' 1172 1256 CALL message( 'check_parameters', 'PA0052', 1, 2, 0, 6, 0 ) 1173 1257 ENDIF 1174  IF ( scalar_advec /= 'pw-scheme' ) THEN  1258 IF ( scalar_advec /= 'pw-scheme' .AND. &  1259 scalar_advec /= 'ws-scheme' ) THEN 1175 1260 message_string = 'non-cyclic lateral boundaries do not allow ' // & 1176 1261 'scalar_advec = "' // TRIM( scalar_advec ) // '"' 1177 1262 CALL message( 'check_parameters', 'PA0053', 1, 2, 0, 6, 0 ) 1178 1263 ENDIF  1264 IF ( (scalar_advec == 'ws-scheme' .OR. momentum_advec == 'ws-scheme' ) &  1265 .AND. loop_optimization == 'vector' ) THEN  1266 message_string = 'non-cyclic lateral boundaries do not allow ' // &  1267 'loop_optimization = vector and ' // &  1268 'scalar_advec = "' // TRIM( scalar_advec ) // '"'  1269 ! The error message number still needs modification.  1270 CALL message( 'check_parameters', 'PA0342', 1, 2, 0, 6, 0 )  1271 END IF 1179 1272 IF ( galilei_transformation ) THEN 1180 1273 message_string = 'non-cyclic lateral boundaries do not allow ' // & … …  1407 1500 TRIM( bc_uv_b ) // '"' 1408 1501 CALL message( 'check_parameters', 'PA0076', 1, 2, 0, 6, 0 )  1502 ENDIF  1503 !  1504 !-- In case of coupled simulations u and v at the ground in atmosphere will be  1505 !-- assigned with the u and v values of the ocean surface  1506 IF ( coupling_mode == 'atmosphere_to_ocean' ) THEN  1507 ibc_uv_b = 2 1409 1508 ENDIF 1410 1509 … …  2109 2208 hom(:,2,57,:) = SPREAD( zu, 2, statistic_regions+1 ) 2110 2209  2210 2111 2211 CASE ( 'u"pt"' ) 2112 2212 dopr_index(i) = 58 … …  2244 2344 2245 2345 END SELECT  2346 2246 2347 ! 2247 2348 !-- Check to which of the predefined coordinate systems the profile belongs … …  2584 2685 !-- Upper plot limit (grid point value) for 1D profiles 2585 2686 IF ( z_max_do1d == -1.0 ) THEN  2687 2586 2688 nz_do1d = nzt+1  2689 2587 2690 ELSE 2588 2691 DO k = nzb+1, nzt+1 … …  2737 2840 2738 2841 !  2842 2739 2843 !-- Check netcdf precison 2740 2844 ldum = .FALSE. … …  3070 3174 IF ( conserve_volume_flow ) THEN 3071 3175 IF ( TRIM( conserve_volume_flow_mode ) == 'default' ) THEN 3072  IF ( bc_lr /= 'cyclic' .OR. bc_ns /= 'cyclic' ) THEN 3073  conserve_volume_flow_mode = 'inflow_profile' 3074  ELSE 3075  conserve_volume_flow_mode = 'initial_profiles' 3076  ENDIF  3176  3177 conserve_volume_flow_mode = 'initial_profiles'  3178 3077 3179 ELSEIF ( TRIM( conserve_volume_flow_mode ) /= 'initial_profiles' .AND. & 3078 3180 TRIM( conserve_volume_flow_mode ) /= 'inflow_profile' .AND. & … …  3082 3184 CALL message( 'check_parameters', 'PA0154', 1, 2, 0, 6, 0 ) 3083 3185 ENDIF 3084  IF ( ( bc_lr /= 'cyclic' .OR. bc_ns /= 'cyclic' ).AND. &3085  TRIM( conserve_volume_flow_mode ) /= 'inflow_profile' ) THEN3086  WRITE( message_string, * ) 'non cyclic boundary conditions ', &3087  'require & conserve_volume_flow_mode = ''inflow_profile''' 3186 IF ( (bc_lr /= 'cyclic' .OR. bc_ns /= 'cyclic') .AND. &  3187 TRIM( conserve_volume_flow_mode ) == 'bulk_velocity' ) THEN  3188 WRITE( message_string, * ) 'non-cyclic boundary conditions ', &  3189 'require conserve_volume_flow_mode = ''initial_profiles''' 3088 3190 CALL message( 'check_parameters', 'PA0155', 1, 2, 0, 6, 0 ) 3089 3191 ENDIF … …  3091 3193 TRIM( conserve_volume_flow_mode ) == 'inflow_profile' ) THEN 3092 3194 WRITE( message_string, * ) 'cyclic boundary conditions ', & 3093  'require &conserve_volume_flow_mode = ''initial_profiles''', & 3195 'require conserve_volume_flow_mode = ''initial_profiles''', & 3094 3196 ' or ''bulk_velocity''' 3095 3197 CALL message( 'check_parameters', 'PA0156', 1, 2, 0, 6, 0 ) … …  3100 3202 TRIM( conserve_volume_flow_mode ) /= 'bulk_velocity' ) ) THEN 3101 3203 WRITE( message_string, * ) 'nonzero bulk velocity requires ', & 3102  'conserve_volume_flow = .T. and &', & 3204 'conserve_volume_flow = .T. and ', & 3103 3205 'conserve_volume_flow_mode = ''bulk_velocity''' 3104 3206 CALL message( 'check_parameters', 'PA0157', 1, 2, 0, 6, 0 ) … …  3139 3241 3140 3242  3243 3141 3244 END SUBROUTINE check_parameters -
palm/trunk/SOURCE/data_output_2d.f90
r623 r667  4 4 ! Current revisions: 5 5 ! ----------------- 6  !  6 ! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng in loops and  7 ! allocation of arrays local_2d and total_2d.  8 ! Calls of exchange_horiz are modiefied. 7 9 ! 8 10 ! Former revisions: … …  112 114 113 115 CASE ( 'xy' ) 114  115 116 s = 1 116  ALLOCATE( level_z( 0:nzt+1), local_2d(nxl-1:nxr+1,nys-1:nyn+1) ) 117 ALLOCATE( level_z(nzb:nzt+1), local_2d(nxlg:nxrg,nysg:nyng) ) 117 118 118 119 ! … …  130 131 IF ( iso2d_output ) CALL check_open( 21 ) 131 132 #if defined( __parallel ) 132  ALLOCATE( total_2d(- 1:nx+1,-1:ny+1) ) 133 ALLOCATE( total_2d(-nbgp:nx+nbgp,-nbgp:ny+nbgp) ) 133 134 #endif 134 135 ENDIF … …  136 137 137 138 CASE ( 'xz' ) 138  139 139 s = 2 140  ALLOCATE( local_2d(nxl -1:nxr+1,nzb:nzt+1) ) 140 ALLOCATE( local_2d(nxlg:nxrg,nzb:nzt+1) ) 141 141 142 142 ! … …  154 154 IF ( iso2d_output ) CALL check_open( 22 ) 155 155 #if defined( __parallel ) 156  ALLOCATE( total_2d(- 1:nx+1,nzb:nzt+1) ) 156 ALLOCATE( total_2d(-nbgp:nx+nbgp,nzb:nzt+1) ) 157 157 #endif 158 158 ENDIF … …  162 162 163 163 s = 3 164  ALLOCATE( local_2d(nys -1:nyn+1,nzb:nzt+1) ) 164 ALLOCATE( local_2d(nysg:nyng,nzb:nzt+1) ) 165 165 166 166 ! … …  178 178 IF ( iso2d_output ) CALL check_open( 23 ) 179 179 #if defined( __parallel ) 180  ALLOCATE( total_2d(- 1:ny+1,nzb:nzt+1) ) 180 ALLOCATE( total_2d(-nbgp:ny+nbgp,nzb:nzt+1) ) 181 181 #endif 182 182 ENDIF … …  192 192 ! 193 193 !-- Allocate a temporary array for resorting (kji -> ijk). 194  ALLOCATE( local_pf(nxl -1:nxr+1,nys-1:nyn+1,nzb:nzt+1) ) 194 ALLOCATE( local_pf(nxlg:nxrg,nysg:nyng,nzb:nzt+1) ) 195 195 196 196 ! … …  219 219 CASE ( 'lwp*_xy' ) ! 2d-array 220 220 IF ( av == 0 ) THEN 221  DO i = nxl -1, nxr+1222  DO j = nys -1, nyn+1 221 DO i = nxlg, nxrg  222 DO j = nysg, nyng 223 223 local_pf(i,j,nzb+1) = SUM( ql(nzb:nzt,j,i) * & 224 224 dzw(1:nzt+1) ) … …  226 226 ENDDO 227 227 ELSE 228  DO i = nxl -1, nxr+1229  DO j = nys -1, nyn+1 228 DO i = nxlg, nxrg  229 DO j = nysg, nyng 230 230 local_pf(i,j,nzb+1) = lwp_av(j,i) 231 231 ENDDO … …  248 248 IF ( simulated_time >= particle_advection_start ) THEN 249 249 tend = prt_count 250  CALL exchange_horiz( tend ) 250 CALL exchange_horiz( tend, nbgp ) 251 251 ELSE 252 252 tend = 0.0 253 253 ENDIF 254  DO i = nxl -1, nxr+1255  DO j = nys -1, nyn+1 254 DO i = nxlg, nxrg  255 DO j = nysg, nyng 256 256 DO k = nzb, nzt+1 257 257 local_pf(i,j,k) = tend(k,j,i) … …  261 261 resorted = .TRUE. 262 262 ELSE 263  CALL exchange_horiz( pc_av ) 263 CALL exchange_horiz( pc_av, nbgp ) 264 264 to_be_resorted => pc_av 265 265 ENDIF … …  287 287 ENDDO 288 288 ENDDO 289  CALL exchange_horiz( tend ) 289 CALL exchange_horiz( tend, nbgp ) 290 290 ELSE 291 291 tend = 0.0 292  END IF293  DO i = nxl -1, nxr+1294  DO j = nys -1, nyn+1 292 END IF  293 DO i = nxlg, nxrg  294 DO j = nysg, nyng 295 295 DO k = nzb, nzt+1 296 296 local_pf(i,j,k) = tend(k,j,i) … …  300 300 resorted = .TRUE. 301 301 ELSE 302  CALL exchange_horiz( pr_av ) 302 CALL exchange_horiz( pr_av, nbgp ) 303 303 to_be_resorted => pr_av 304 304 ENDIF … …  306 306 CASE ( 'pra*_xy' ) ! 2d-array / integral quantity => no av 307 307 CALL exchange_horiz_2d( precipitation_amount ) 308  DO i = nxl-1, nxr+1309  DO j = nys-1, nyn+1 308 DO i = nxlg, nxrg  309 DO j = nysg, nyng 310 310 local_pf(i,j,nzb+1) = precipitation_amount(j,i) 311 311 ENDDO … …  319 319 IF ( av == 0 ) THEN 320 320 CALL exchange_horiz_2d( precipitation_rate ) 321  DO i = nxl -1, nxr+1322  DO j = nys -1, nyn+1 321 DO i = nxlg, nxrg  322 DO j = nysg, nyng 323 323 local_pf(i,j,nzb+1) = precipitation_rate(j,i) 324 324 ENDDO … …  326 326 ELSE 327 327 CALL exchange_horiz_2d( precipitation_rate_av ) 328  DO i = nxl -1, nxr+1329  DO j = nys -1, nyn+1 328 DO i = nxlg, nxrg  329 DO j = nysg, nyng 330 330 local_pf(i,j,nzb+1) = precipitation_rate_av(j,i) 331 331 ENDDO … …  341 341 to_be_resorted => pt 342 342 ELSE 343  DO i = nxl-1, nxr+1344  DO j = nys-1, nyn+1 343 DO i = nxlg, nxrg  344 DO j = nysg, nyng 345 345 DO k = nzb, nzt+1 346 346 local_pf(i,j,k) = pt(k,j,i) + l_d_cp * & … …  399 399 CASE ( 'qsws*_xy' ) ! 2d-array 400 400 IF ( av == 0 ) THEN 401  DO i = nxl -1, nxr+1402  DO j = nys -1, nyn+1 401 DO i = nxlg, nxrg  402 DO j = nysg, nyng 403 403 local_pf(i,j,nzb+1) = qsws(j,i) 404 404 ENDDO 405 405 ENDDO 406 406 ELSE 407  DO i = nxl -1, nxr+1408  DO j = nys -1, nyn+1 407 DO i = nxlg, nxrg  408 DO j = nysg, nyng 409 409 local_pf(i,j,nzb+1) = qsws_av(j,i) 410 410 ENDDO … …  417 417 CASE ( 'qv_xy', 'qv_xz', 'qv_yz' ) 418 418 IF ( av == 0 ) THEN 419  DO i = nxl -1, nxr+1420  DO j = nys -1, nyn+1 419 DO i = nxlg, nxrg  420 DO j = nysg, nyng 421 421 DO k = nzb, nzt+1 422 422 local_pf(i,j,k) = q(k,j,i) - ql(k,j,i) … …  453 453 CASE ( 'shf*_xy' ) ! 2d-array 454 454 IF ( av == 0 ) THEN 455  DO i = nxl -1, nxr+1456  DO j = nys -1, nyn+1 455 DO i = nxlg, nxrg  456 DO j = nysg, nyng 457 457 local_pf(i,j,nzb+1) = shf(j,i) 458 458 ENDDO 459 459 ENDDO 460 460 ELSE 461  DO i = nxl -1, nxr+1462  DO j = nys -1, nyn+1 461 DO i = nxlg, nxrg  462 DO j = nysg, nyng 463 463 local_pf(i,j,nzb+1) = shf_av(j,i) 464 464 ENDDO … …  471 471 CASE ( 't*_xy' ) ! 2d-array 472 472 IF ( av == 0 ) THEN 473  DO i = nxl -1, nxr+1474  DO j = nys -1, nyn+1 473 DO i = nxlg, nxrg  474 DO j = nysg, nyng 475 475 local_pf(i,j,nzb+1) = ts(j,i) 476 476 ENDDO 477 477 ENDDO 478 478 ELSE 479  DO i = nxl -1, nxr+1480  DO j = nys -1, nyn+1 479 DO i = nxlg, nxrg  480 DO j = nysg, nyng 481 481 local_pf(i,j,nzb+1) = ts_av(j,i) 482 482 ENDDO … …  503 503 CASE ( 'u*_xy' ) ! 2d-array 504 504 IF ( av == 0 ) THEN 505  DO i = nxl -1, nxr+1506  DO j = nys -1, nyn+1 505 DO i = nxlg, nxrg  506 DO j = nysg, nyng 507 507 local_pf(i,j,nzb+1) = us(j,i) 508 508 ENDDO 509 509 ENDDO 510 510 ELSE 511  DO i = nxl -1, nxr+1512  DO j = nys -1, nyn+1 511 DO i = nxlg, nxrg  512 DO j = nysg, nyng 513 513 local_pf(i,j,nzb+1) = us_av(j,i) 514 514 ENDDO … …  551 551 CASE ( 'z0*_xy' ) ! 2d-array 552 552 IF ( av == 0 ) THEN 553  DO i = nxl -1, nxr+1554  DO j = nys -1, nyn+1 553 DO i = nxlg, nxrg  554 DO j = nysg, nyng 555 555 local_pf(i,j,nzb+1) = z0(j,i) 556 556 ENDDO 557 557 ENDDO 558 558 ELSE 559  DO i = nxl -1, nxr+1560  DO j = nys -1, nyn+1 559 DO i = nxlg, nxrg  560 DO j = nysg, nyng 561 561 local_pf(i,j,nzb+1) = z0_av(j,i) 562 562 ENDDO … …  593 593 !-- Resort the array to be output, if not done above 594 594 IF ( .NOT. resorted ) THEN 595  DO i = nxl -1, nxr+1596  DO j = nys -1, nyn+1 595 DO i = nxlg, nxrg  596 DO j = nysg, nyng 597 597 DO k = nzb, nzt+1 598 598 local_pf(i,j,k) = to_be_resorted(k,j,i) … …  647 647 !-- Carry out the averaging (all data are on the PE) 648 648 DO k = nzb, nzt+1 649  DO j = nys -1, nyn+1650  DO i = nxl -1, nxr+1 649 DO j = nysg, nyng  650 DO i = nxlg, nxrg 651 651 local_2d(i,j) = local_2d(i,j) + local_pf(i,j,k) 652 652 ENDDO … …  654 654 ENDDO 655 655 656  local_2d = local_2d / ( nzt -nzb + 2.0  656 local_2d = local_2d / ( nzt -nzb + 2.0) 657 657 658 658 ELSE … …  723 723 ENDIF 724 724 #endif 725  WRITE ( 21 ) nxl -1, nxr+1, nys-1, nyn+1 725 WRITE ( 21 ) nxlg, nxrg, nysg, nyng 726 726 WRITE ( 21 ) local_2d 727 727 … …  734 734 CALL MPI_BARRIER( comm2d, ierr ) 735 735 736  ngp = ( nxr -nxl+3 ) * ( nyn-nys+3) 736 ngp = ( nxrg-nxlg+1 ) * ( nyng-nysg+1 ) 737 737 IF ( myid == 0 ) THEN 738 738 ! 739 739 !-- Local array can be relocated directly. 740  total_2d(nxl -1:nxr+1,nys-1:nyn+1) = local_2d 740 total_2d(nxlg:nxrg,nysg:nyng) = local_2d 741 741 ! 742 742 !-- Receive data from all other PEs. … …  760 760 !-- Output of the total cross-section. 761 761 IF ( iso2d_output ) THEN 762  WRITE (21) total_2d( 0:nx+1,0:ny+1) 762 WRITE (21) total_2d(-nbgp:nx+nbgp,-nbgp:ny+nbgp) 763 763 ENDIF 764 764 ! 765 765 !-- Relocate the local array for the next loop increment 766 766 DEALLOCATE( local_2d ) 767  ALLOCATE( local_2d(nxl -1:nxr+1,nys-1:nyn+1) ) 767 ALLOCATE( local_2d(nxlg:nxrg,nysg:nyng) ) 768 768 769 769 #if defined( __netcdf ) … …  789 789 ! 790 790 !-- First send the local index limits to PE0 791  ind(1) = nxl -1; ind(2) = nxr+1792  ind(3) = nys -1; ind(4) = nyn+1 791 ind(1) = nxlg; ind(2) = nxrg  792 ind(3) = nysg; ind(4) = nyng 793 793 CALL MPI_SEND( ind(1), 4, MPI_INTEGER, 0, 0, & 794 794 comm2d, ierr ) 795 795 ! 796 796 !-- Send data to PE0 797  CALL MPI_SEND( local_2d(nxl -1,nys-1), ngp, & 797 CALL MPI_SEND( local_2d(nxlg,nysg), ngp, & 798 798 MPI_REAL, 0, 1, comm2d, ierr ) 799 799 ENDIF … …  882 882 883 883 ENDIF  884 884 885 ! 885 886 !-- If required, carry out averaging along y 886 887 IF ( section(is,s) == -1 ) THEN 887 888 888  ALLOCATE( local_2d_l(nxl -1:nxr+1,nzb:nzt+1) ) 889 ALLOCATE( local_2d_l(nxlg:nxrg,nzb:nzt+1) ) 889 890 local_2d_l = 0.0 890  ngp = ( nxr -nxl+3) * ( nzt-nzb+2 ) 891 ngp = ( nxrg-nxlg+1 ) * ( nzt-nzb+2 ) 891 892 ! 892 893 !-- First local averaging on the PE 893 894 DO k = nzb, nzt+1 894 895 DO j = nys, nyn 895  DO i = nxl -1, nxr+1 896 DO i = nxlg, nxrg 896 897 local_2d_l(i,k) = local_2d_l(i,k) + & 897 898 local_pf(i,j,k) … …  903 904 !-- Now do the averaging over all PEs along y 904 905 IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr ) 905  CALL MPI_ALLREDUCE( local_2d_l(nxl -1,nzb), &906  local_2d(nxl -1,nzb), ngp, MPI_REAL, & 906 CALL MPI_ALLREDUCE( local_2d_l(nxlg,nzb), &  907 local_2d(nxlg,nzb), ngp, MPI_REAL, & 907 908 MPI_SUM, comm1dy, ierr ) 908 909 #else … …  936 937 !-- BEGIN WORKAROUND--------------------------------------- 937 938 IF ( npey /= 1 .AND. section(is,s) /= -1) THEN 938  ALLOCATE( local_2d_l(nxl -1:nxr+1,nzb:nzt+1) ) 939 ALLOCATE( local_2d_l(nxlg:nxrg,nzb:nzt+1) ) 939 940 local_2d_l = 0.0 940 941 IF ( section(is,s) >= nys .AND. section(is,s) <= nyn )& … …  945 946 ! 946 947 !-- Distribute data over all PEs along y 947  ngp = ( nxr -nxl+3) * ( nzt-nzb+2 ) 948 ngp = ( nxrg-nxlg+1 ) * ( nzt-nzb+2 ) 948 949 IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr ) 949  CALL MPI_ALLREDUCE( local_2d_l(nxl -1,nzb), &950  local_2d(nxl -1,nzb), ngp, & 950 CALL MPI_ALLREDUCE( local_2d_l(nxlg,nzb), &  951 local_2d(nxlg,nzb), ngp, & 951 952 MPI_REAL, MPI_SUM, comm1dy, ierr ) 952 953 #else … …  1022 1023 ( section(is,s) == -1 .AND. nys-1 == -1 ) ) & 1023 1024 THEN 1024  WRITE (22) nxl -1, nxr+1, nzb, nzt+1 1025 WRITE (22) nxlg, nxrg, nzb, nzt+1 1025 1026 WRITE (22) local_2d 1026 1027 ELSE … …  1036 1037 CALL MPI_BARRIER( comm2d, ierr ) 1037 1038 1038  ngp = ( nxr -nxl+3) * ( nzt-nzb+2 ) 1039 ngp = ( nxrg-nxlg+1 ) * ( nzt-nzb+2 ) 1039 1040 IF ( myid == 0 ) THEN 1040 1041 ! … …  1044 1045 ( section(is,s) == -1 .AND. nys-1 == -1 ) ) & 1045 1046 THEN 1046  total_2d(nxl -1:nxr+1,nzb:nzt+1) = local_2d 1047 total_2d(nxlg:nxrg,nzb:nzt+1) = local_2d 1047 1048 ENDIF 1048 1049 ! … …  1073 1074 !-- Output of the total cross-section. 1074 1075 IF ( iso2d_output ) THEN 1075  WRITE (22) total_2d( 0:nx+1,nzb:nzt+1) 1076 WRITE (22) total_2d(-nbgp:nx+nbgp,nzb:nzt+1) 1076 1077 ENDIF 1077 1078 ! 1078 1079 !-- Relocate the local array for the next loop increment 1079 1080 DEALLOCATE( local_2d ) 1080  ALLOCATE( local_2d(nxl -1:nxr+1,nzb:nzt+1) ) 1081 ALLOCATE( local_2d(nxlg:nxrg,nzb:nzt+1) ) 1081 1082 1082 1083 #if defined( __netcdf ) … …  1099 1100 ( section(is,s) == -1 .AND. nys-1 == -1 ) ) & 1100 1101 THEN 1101  ind(1) = nxl -1; ind(2) = nxr+1 1102 ind(1) = nxlg; ind(2) = nxrg 1102 1103 ind(3) = nzb; ind(4) = nzt+1 1103 1104 ELSE … …  1110 1111 !-- If applicable, send data to PE0. 1111 1112 IF ( ind(1) /= -9999 ) THEN 1112  CALL MPI_SEND( local_2d(nxl -1,nzb), ngp, & 1113 CALL MPI_SEND( local_2d(nxlg,nzb), ngp, & 1113 1114 MPI_REAL, 0, 1, comm2d, ierr ) 1114 1115 ENDIF … …  1187 1188 IF ( section(is,s) == -1 ) THEN 1188 1189 1189  ALLOCATE( local_2d_l(nys -1:nyn+1,nzb:nzt+1) ) 1190 ALLOCATE( local_2d_l(nysg:nyng,nzb:nzt+1) ) 1190 1191 local_2d_l = 0.0 1191  ngp = ( nyn -nys+3) * ( nzt-nzb+2 ) 1192 ngp = ( nyng-nysg+1 ) * ( nzt-nzb+2 ) 1192 1193 ! 1193 1194 !-- First local averaging on the PE 1194 1195 DO k = nzb, nzt+1 1195  DO j = nys -1, nyn+1 1196 DO j = nysg, nyng 1196 1197 DO i = nxl, nxr 1197 1198 local_2d_l(j,k) = local_2d_l(j,k) + & … …  1204 1205 !-- Now do the averaging over all PEs along x 1205 1206 IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr ) 1206  CALL MPI_ALLREDUCE( local_2d_l(nys -1,nzb), &1207  local_2d(nys -1,nzb), ngp, MPI_REAL, & 1207 CALL MPI_ALLREDUCE( local_2d_l(nysg,nzb), &  1208 local_2d(nysg,nzb), ngp, MPI_REAL, & 1208 1209 MPI_SUM, comm1dx, ierr ) 1209 1210 #else … …  1237 1238 !-- BEGIN WORKAROUND--------------------------------------- 1238 1239 IF ( npex /= 1 .AND. section(is,s) /= -1) THEN 1239  ALLOCATE( local_2d_l(nys -1:nyn+1,nzb:nzt+1) ) 1240 ALLOCATE( local_2d_l(nysg:nyng,nzb:nzt+1) ) 1240 1241 local_2d_l = 0.0 1241 1242 IF ( section(is,s) >= nxl .AND. section(is,s) <= nxr )& … …  1246 1247 ! 1247 1248 !-- Distribute data over all PEs along x 1248  ngp = ( nyn -nys+3 ) * ( nzt-nzb+2 ) 1249 ngp = ( nyng-nysg+1 ) * ( nzt-nzb + 2 ) 1249 1250 IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr ) 1250  CALL MPI_ALLREDUCE( local_2d_l(nys -1,nzb), &1251  local_2d(nys -1,nzb), ngp, & 1251 CALL MPI_ALLREDUCE( local_2d_l(nysg,nzb), &  1252 local_2d(nysg,nzb), ngp, & 1252 1253 MPI_REAL, MPI_SUM, comm1dx, ierr ) 1253 1254 #else … …  1323 1324 ( section(is,s) == -1 .AND. nxl-1 == -1 ) ) & 1324 1325 THEN 1325  WRITE (23) nys -1, nyn+1, nzb, nzt+1 1326 WRITE (23) nysg, nyng, nzb, nzt+1 1326 1327 WRITE (23) local_2d 1327 1328 ELSE … …  1337 1338 CALL MPI_BARRIER( comm2d, ierr ) 1338 1339 1339  ngp = ( nyn -nys+3) * ( nzt-nzb+2 ) 1340 ngp = ( nyng-nysg+1 ) * ( nzt-nzb+2 ) 1340 1341 IF ( myid == 0 ) THEN 1341 1342 ! … …  1345 1346 ( section(is,s) == -1 .AND. nxl-1 == -1 ) ) & 1346 1347 THEN 1347  total_2d(nys -1:nyn+1,nzb:nzt+1) = local_2d 1348 total_2d(nysg:nyng,nzb:nzt+1) = local_2d 1348 1349 ENDIF 1349 1350 ! … …  1379 1380 !-- Relocate the local array for the next loop increment 1380 1381 DEALLOCATE( local_2d ) 1381  ALLOCATE( local_2d(nys -1:nyn+1,nzb:nzt+1) ) 1382 ALLOCATE( local_2d(nysg:nyng,nzb:nzt+1) ) 1382 1383 1383 1384 #if defined( __netcdf ) … …  1400 1401 ( section(is,s) == -1 .AND. nxl-1 == -1 ) ) & 1401 1402 THEN 1402  ind(1) = nys -1; ind(2) = nyn+1 1403 ind(1) = nysg; ind(2) = nyng 1403 1404 ind(3) = nzb; ind(4) = nzt+1 1404 1405 ELSE … …  1411 1412 !-- If applicable, send data to PE0. 1412 1413 IF ( ind(1) /= -9999 ) THEN 1413  CALL MPI_SEND( local_2d(nys -1,nzb), ngp, & 1414 CALL MPI_SEND( local_2d(nysg,nzb), ngp, & 1414 1415 MPI_REAL, 0, 1, comm2d, ierr ) 1415 1416 ENDIF -
palm/trunk/SOURCE/data_output_3d.f90
r647 r667  4 4 ! Current revisions: 5 5 ! ----------------- 6  !  6 ! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng in loops and  7 ! allocation of arrays. Calls of exchange_horiz are modified.  8 ! Skip-value skip_do_avs changed to a dynamic adaption of ghost points. 7 9 ! 8 10 ! Former revisions: … …  102 104 ! 103 105 !-- Allocate a temporary array with the desired output dimensions. 104  ALLOCATE( local_pf(nxl -1:nxr+1,nys-1:nyn+1,nzb:nz_do3d) ) 106 ALLOCATE( local_pf(nxlg:nxrg,nysg:nyng,nzb:nz_do3d) ) 105 107 106 108 ! … …  157 159 IF ( av == 0 ) THEN 158 160 tend = prt_count 159  CALL exchange_horiz( tend )160  DO i = nxl -1, nxr+1161  DO j = nys -1, nyn+1 161 CALL exchange_horiz( tend, nbgp )  162 DO i = nxlg, nxrg  163 DO j = nysg, nyng 162 164 DO k = nzb, nz_do3d 163 165 local_pf(i,j,k) = tend(k,j,i) … …  167 169 resorted = .TRUE. 168 170 ELSE 169  CALL exchange_horiz( pc_av ) 171 CALL exchange_horiz( pc_av, nbgp ) 170 172 to_be_resorted => pc_av 171 173 ENDIF … …  192 194 ENDDO 193 195 ENDDO 194  CALL exchange_horiz( tend )195  DO i = nxl -1, nxr+1196  DO j = nys -1, nyn+1 196 CALL exchange_horiz( tend, nbgp )  197 DO i = nxlg, nxrg  198 DO j = nysg, nyng 197 199 DO k = nzb, nzt+1 198 200 local_pf(i,j,k) = tend(k,j,i) … …  202 204 resorted = .TRUE. 203 205 ELSE 204  CALL exchange_horiz( pr_av ) 206 CALL exchange_horiz( pr_av, nbgp ) 205 207 to_be_resorted => pr_av 206 208 ENDIF … …  211 213 to_be_resorted => pt 212 214 ELSE 213  DO i = nxl -1, nxr+1214  DO j = nys -1, nyn+1 215 DO i = nxlg, nxrg  216 DO j = nysg, nyng 215 217 DO k = nzb, nz_do3d 216 218 local_pf(i,j,k) = pt(k,j,i) + l_d_cp * & … …  263 265 CASE ( 'qv' ) 264 266 IF ( av == 0 ) THEN 265  DO i = nxl -1, nxr+1266  DO j = nys -1, nyn+1 267 DO i = nxlg, nxrg  268 DO j = nysg, nyng 267 269 DO k = nzb, nz_do3d 268 270 local_pf(i,j,k) = q(k,j,i) - ql(k,j,i) … …  342 344 !-- Resort the array to be output, if not done above 343 345 IF ( .NOT. resorted ) THEN 344  DO i = nxl -1, nxr+1345  DO j = nys -1, nyn+1 346 DO i = nxlg, nxrg  347 DO j = nysg, nyng 346 348 DO k = nzb, nz_do3d 347 349 local_pf(i,j,k) = to_be_resorted(k,j,i) … …  376 378 !-- Determine the Skip-value for the next array. Record end and start 377 379 !-- require 4 byte each. 378  skip_do_avs = skip_do_avs + ( ((nx+2 )*(ny+2)*(nz_do3d+1)) * 4 + 8 ) 380 skip_do_avs = skip_do_avs + ( ((nx+2*nbgp)*(ny+2*nbgp)*(nz_do3d+1)) * 4 + 8 ) 379 381 ENDIF 380 382 … …  386 388 !-- of compressed data. 387 389 CALL write_compressed( local_pf, 30, 33, myid, nxl, nxr, nyn, nys, & 388  nzb, nz_do3d, prec ) 390 nzb, nz_do3d, prec, nbgp ) 389 391 ELSE 390 392 ! … …  400 402 WRITE ( 30 ) simulated_time, do3d_time_count(av), av 401 403 ENDIF 402  WRITE ( 30 ) nxl -1, nxr+1, nys-1, nyn+1, nzb, nz_do3d 404 WRITE ( 30 ) nxlg, nxrg, nysg, nyng, nzb, nz_do3d 403 405 WRITE ( 30 ) local_pf 404 406 … …  411 413 IF ( nxr == nx .AND. nyn /= ny ) THEN 412 414 nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_do3d(av,if), & 413  local_pf(nxl:nxr +1,nys:nyn,nzb:nz_do3d), & 415 local_pf(nxl:nxrg,nys:nyn,nzb:nz_do3d), & 414 416 start = (/ nxl+1, nys+1, nzb+1, do3d_time_count(av) /), & 415  count = (/ nxr-nxl+ 2, nyn-nys+1, nz_do3d-nzb+1, 1 /) ) 417 count = (/ nxr-nxl+1+nbgp, nyn-nys+1, nz_do3d-nzb+1, 1 /) ) 416 418 ELSEIF ( nxr /= nx .AND. nyn == ny ) THEN 417 419 nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_do3d(av,if), & 418  local_pf(nxl:nxr,nys:nyn +1,nzb:nz_do3d), & 420 local_pf(nxl:nxr,nys:nyng,nzb:nz_do3d), & 419 421 start = (/ nxl+1, nys+1, nzb+1, do3d_time_count(av) /), & 420  count = (/ nxr-nxl+1, nyn-nys+ 2, nz_do3d-nzb+1, 1 /) ) 422 count = (/ nxr-nxl+1, nyn-nys+1+nbgp, nz_do3d-nzb+1, 1 /) ) 421 423 ELSEIF ( nxr == nx .AND. nyn == ny ) THEN 422 424 nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_do3d(av,if), & 423  local_pf(nxl:nxr +1,nys:nyn+1,nzb:nz_do3d), & 425 local_pf(nxl:nxrg,nys:nyng,nzb:nz_do3d), & 424 426 start = (/ nxl+1, nys+1, nzb+1, do3d_time_count(av) /), & 425  count = (/ nxr-nxl+ 2, nyn-nys+2, nz_do3d-nzb+1, 1 /) ) 427 count = (/ nxr-nxl+1+nbgp, nyn-nys+1+nbgp, nz_do3d-nzb+1, 1 /) ) 426 428 ELSE 427 429 nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_do3d(av,if), & -
palm/trunk/SOURCE/data_output_mask.f90
r565 r667  4 4 ! Current revisions: 5 5 ! ----------------- 6  !  6 ! Calls of exchange_horiz are modified. 7 7 ! 8 8 ! Former revisions: … …  123 123 IF ( av == 0 ) THEN 124 124 tend = prt_count 125  CALL exchange_horiz( tend ) 125 CALL exchange_horiz( tend, nbgp ) 126 126 DO i = 1, mask_size_l(mid,1) 127 127 DO j = 1, mask_size_l(mid,2) … …  134 134 resorted = .TRUE. 135 135 ELSE 136  CALL exchange_horiz( pc_av ) 136 CALL exchange_horiz( pc_av, nbgp ) 137 137 to_be_resorted => pc_av 138 138 ENDIF … …  159 159 ENDDO 160 160 ENDDO 161  CALL exchange_horiz( tend ) 161 CALL exchange_horiz( tend, nbgp ) 162 162 DO i = 1, mask_size_l(mid,1) 163 163 DO j = 1, mask_size_l(mid,2) … …  170 170 resorted = .TRUE. 171 171 ELSE 172  CALL exchange_horiz( pr_av ) 172 CALL exchange_horiz( pr_av, nbgp ) 173 173 to_be_resorted => pr_av 174 174 ENDIF … …  439 439 440 440 if = if + 1  441 441 442 ENDDO 442 443 -
palm/trunk/SOURCE/diffusion_e.f90
r484 r667  4 4 ! Current revisions: 5 5 ! ----------------- 6  !  6 ! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng 7 7 ! 8 8 ! Former revisions: … …  68 68 REAL :: dvar_dz, l_stable, phi_m, var_reference 69 69 REAL :: ddzu(1:nzt+1), dd2zu(1:nzt), ddzw(1:nzt+1), & 70  l_grid(1:nzt), zu( 0:nzt+1), zw(0:nzt+1)71  REAL, DIMENSION(nzb:nzt+1,nys -1:nyn+1,nxl-1:nxr+1) :: diss, tend 70 l_grid(1:nzt), zu(nzb:nzt+1), zw(nzb:nzt+1)  71 REAL, DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: diss, tend 72 72 REAL, DIMENSION(:,:), POINTER :: rif 73 73 REAL, DIMENSION(:,:,:), POINTER :: e, km, var … …  289 289 REAL :: dvar_dz, l_stable, phi_m, var_reference 290 290 REAL :: ddzu(1:nzt+1), dd2zu(1:nzt), ddzw(1:nzt+1), & 291  l_grid(1:nzt), zu( 0:nzt+1), zw(0:nzt+1)292  REAL, DIMENSION(nzb:nzt+1,nys -1:nyn+1,nxl-1:nxr+1) ::diss, tend 291 l_grid(1:nzt), zu(nzb:nzt+1), zw(nzb:nzt+1)  292 REAL, DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: diss, tend 293 293 REAL, DIMENSION(:,:), POINTER :: rif 294 294 REAL, DIMENSION(:,:,:), POINTER :: e, km, var -
palm/trunk/SOURCE/diffusion_s.f90
r484 r667  4 4 ! Current revisions: 5 5 ! ----------------- 6  !  6 ! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng 7 7 ! 8 8 ! Former revisions: … …  64 64 REAL :: vertical_gridspace 65 65 REAL :: ddzu(1:nzt+1), ddzw(1:nzt+1) 66  REAL :: tend(nzb:nzt+1,nys -1:nyn+1,nxl-1:nxr+1) 66 REAL :: tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg) 67 67 REAL :: wall_s_flux(0:4) 68 68 REAL, DIMENSION(:,:), POINTER :: s_flux_b, s_flux_t … …  176 176 REAL :: vertical_gridspace 177 177 REAL :: ddzu(1:nzt+1), ddzw(1:nzt+1) 178  REAL :: tend(nzb:nzt+1,nys -1:nyn+1,nxl-1:nxr+1) 178 REAL :: tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg) 179 179 REAL :: wall_s_flux(0:4) 180 180 REAL, DIMENSION(:,:), POINTER :: s_flux_b, s_flux_t -
palm/trunk/SOURCE/diffusion_u.f90
r484 r667  4 4 ! Current revisions: 5 5 ! ----------------- 6  !  6 ! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng 7 7 ! 8 8 ! Former revisions: … …  72 72 INTEGER :: i, j, k 73 73 REAL :: kmym_x, kmym_y, kmyp_x, kmyp_y, kmzm, kmzp 74  REAL :: ddzu(1:nzt+1), ddzw(1:nzt+1), km_damp_y(nys -1:nyn+1)75  REAL :: tend(nzb:nzt+1,nys -1:nyn+1,nxl-1:nxr+1) 74 REAL :: ddzu(1:nzt+1), ddzw(1:nzt+1), km_damp_y(nysg:nyng)  75 REAL :: tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg) 76 76 REAL, DIMENSION(:,:), POINTER :: usws, uswst 77 77 REAL, DIMENSION(:,:,:), POINTER :: km, u, v, w … …  177 177 & - kmzm * ( ( u(k,j,i) - u(k-1,j,i) ) * ddzu(k) & 178 178 & + ( w(k-1,j,i) - w(k-1,j,i-1) ) * ddx & 179  & ) & 179 & ) & 180 180 & ) * ddzw(k) 181 181 ENDDO … …  248 248 INTEGER :: i, j, k 249 249 REAL :: kmym_x, kmym_y, kmyp_x, kmyp_y, kmzm, kmzp 250  REAL :: ddzu(1:nzt+1), ddzw(1:nzt+1), km_damp_y(nys -1:nyn+1)251  REAL :: tend(nzb:nzt+1,nys -1:nyn+1,nxl-1:nxr+1) 250 REAL :: ddzu(1:nzt+1), ddzw(1:nzt+1), km_damp_y(nysg:nyng)  251 REAL :: tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg) 252 252 REAL, DIMENSION(nzb:nzt+1) :: usvs 253 253 REAL, DIMENSION(:,:), POINTER :: usws, uswst -
palm/trunk/SOURCE/diffusion_v.f90
r484 r667  4 4 ! Current revisions: 5 5 ! ----------------- 6  !  6 ! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng 7 7 ! 8 8 ! Former revisions: … …  70 70 INTEGER :: i, j, k 71 71 REAL :: kmxm_x, kmxm_y, kmxp_x, kmxp_y, kmzm, kmzp 72  REAL :: ddzu(1:nzt+1), ddzw(1:nzt+1), km_damp_x(nxl -1:nxr+1)73  REAL :: tend(nzb:nzt+1,nys -1:nyn+1,nxl-1:nxr+1) 72 REAL :: ddzu(1:nzt+1), ddzw(1:nzt+1), km_damp_x(nxlg:nxrg)  73 REAL :: tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg) 74 74 REAL, DIMENSION(:,:), POINTER :: vsws, vswst 75 75 REAL, DIMENSION(:,:,:), POINTER :: km, u, v, w … …  246 246 INTEGER :: i, j, k 247 247 REAL :: kmxm_x, kmxm_y, kmxp_x, kmxp_y, kmzm, kmzp 248  REAL :: ddzu(1:nzt+1), ddzw(1:nzt+1), km_damp_x(nxl -1:nxr+1)249  REAL :: tend(nzb:nzt+1,nys -1:nyn+1,nxl-1:nxr+1) 248 REAL :: ddzu(1:nzt+1), ddzw(1:nzt+1), km_damp_x(nxlg:nxrg)  249 REAL :: tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg) 250 250 REAL, DIMENSION(nzb:nzt+1) :: vsus 251 251 REAL, DIMENSION(:,:), POINTER :: vsws, vswst -
palm/trunk/SOURCE/diffusion_w.f90
r484 r667  4 4 ! Current revisions: 5 5 ! ----------------- 6  !  6 ! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng 7 7 ! 8 8 ! Former revisions: … …  64 64 REAL :: kmxm_x, kmxm_z, kmxp_x, kmxp_z, kmym_y, kmym_z, kmyp_y, & 65 65 kmyp_z 66  REAL :: ddzu(1:nzt+1), ddzw(1:nzt+1), km_damp_x(nxl -1:nxr+1), &67  km_damp_y(nys -1:nyn+1)68  REAL :: tend(nzb:nzt+1,nys -1:nyn+1,nxl-1:nxr+1) 66 REAL :: ddzu(1:nzt+1), ddzw(1:nzt+1), km_damp_x(nxlg:nxrg), &  67 km_damp_y(nysg:nyng)  68 REAL :: tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg) 69 69 REAL, DIMENSION(:,:,:), POINTER :: km, u, v, w 70 70 REAL, DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) :: wsus, wsvs … …  211 211 REAL :: kmxm_x, kmxm_z, kmxp_x, kmxp_z, kmym_y, kmym_z, kmyp_y, & 212 212 kmyp_z 213  REAL :: ddzu(1:nzt+1), ddzw(1:nzt+1), km_damp_x(nxl -1:nxr+1), &214  km_damp_y(nys -1:nyn+1)215  REAL :: tend(nzb:nzt+1,nys -1:nyn+1,nxl-1:nxr+1) 213 REAL :: ddzu(1:nzt+1), ddzw(1:nzt+1), km_damp_x(nxlg:nxrg), &  214 km_damp_y(nysg:nyng)  215 REAL :: tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg) 216 216 REAL, DIMENSION(nzb:nzt+1) :: wsus, wsvs 217 217 REAL, DIMENSION(:,:,:), POINTER :: km, u, v, w -
palm/trunk/SOURCE/diffusivities.f90
r510 r667  4 4 ! Current revisions: 5 5 ! ----------------- 6  !  6 ! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng 7 7 ! 8 8 ! Former revisions: … …  54 54 REAL, SAVE :: phi_m = 1.0 55 55 56  REAL :: var(nzb:nzt+1,nys -1:nyn+1,nxl-1:nxr+1) 56 REAL :: var(nzb:nzt+1,nysg:nyng,nxlg:nxrg) 57 57 58 58 REAL, DIMENSION(1:nzt) :: l, ll, sqrt_e … …  73 73 74 74 !$OMP DO 75  DO i = nxl -1, nxr+176  DO j = nys -1, nyn+1 75 DO i = nxlg, nxrg  76 DO j = nysg, nyng 77 77 78 78 ! … …  157 157 158 158 sums_l_l(nzt+1,:,tn) = sums_l_l(nzt,:,tn) ! quasi boundary-condition for 159  ! data output 159 ! data output 160 160 161 161 !$OMP END PARALLEL … …  167 167 !-- values of the diffusivities are not needed 168 168 !$OMP PARALLEL DO 169  DO i = nxl -1, nxr+1170  DO j = nys -1, nyn+1 169 DO i = nxlg, nxrg  170 DO j = nysg, nyng 171 171 km(nzb_s_inner(j,i),j,i) = km(nzb_s_inner(j,i)+1,j,i) 172 172 km(nzt+1,j,i) = km(nzt,j,i) -
palm/trunk/SOURCE/disturb_field.f90
r484 r667  4 4 ! Current revisions: 5 5 ! ----------------- 6  !  6 ! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng  7 ! Calls of exchange_horiz are modified. 7 8 ! 8 9 ! Former revisions: … …  44 45 45 46 INTEGER :: i, j, k 46  INTEGER :: nzb_uv_inner(nys -1:nyn+1,nxl-1:nxr+1) 47 INTEGER :: nzb_uv_inner(nysg:nyng,nxlg:nxrg) 47 48 48 49 REAL :: randomnumber, & 49  dist1(nzb:nzt+1,nys -1:nyn+1,nxl-1:nxr+1), &50  field(nzb:nzt+1,nys -1:nyn+1,nxl-1:nxr+1) 50 dist1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &  51 field(nzb:nzt+1,nysg:nyng,nxlg:nxrg) 51 52 REAL, DIMENSION(:,:,:), ALLOCATABLE :: dist2 52 53 … …  57 58 !-- Create an additional temporary array and initialize the arrays needed 58 59 !-- to store the disturbance 59  ALLOCATE( dist2(nzb:nzt+1,nys -1:nyn+1,nxl-1:nxr+1) ) 60 ALLOCATE( dist2(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 60 61 dist1 = 0.0 61 62 dist2 = 0.0 … …  102 103 !-- Exchange of ghost points for the random perturbation 103 104 104  CALL exchange_horiz( dist1 ) 105   105 CALL exchange_horiz( dist1, nbgp ) 106 106 ! 107 107 !-- Applying the Shuman filter in order to smooth the perturbations. … …  128 128 !-- Exchange of ghost points for the filtered perturbation. 129 129 !-- Afterwards, filter operation and exchange of ghost points are repeated. 130  CALL exchange_horiz( dist2 ) 130 CALL exchange_horiz( dist2, nbgp ) 131 131 132 132 DO i = nxl, nxr … …  141 141 ENDDO 142 142 143  CALL exchange_horiz( dist1 ) 143 CALL exchange_horiz( dist1, nbgp ) 144 144 145 145 ! … …  148 148 !-- (diffusion criterion)) 149 149 IF ( TRIM( topography ) /= 'flat' ) THEN 150  DO i = nxl -1, nxr+1151  DO j = nys -1, nyn+1 150 DO i = nxlg, nxrg  151 DO j = nysg, nyng 152 152 dist1(nzb:nzb_uv_inner(j,i)+1,j,i) = 0.0 153 153 ENDDO … …  157 157 ! 158 158 !-- Random perturbation is added to the array to be disturbed. 159  DO i = nxl -1, nxr+1160  DO j = nys -1, nyn+1 159 DO i = nxlg, nxrg  160 DO j = nysg, nyng 161 161 DO k = disturbance_level_ind_b-2, disturbance_level_ind_t+2 162 162 field(k,j,i) = field(k,j,i) + dist1(k,j,i) -
palm/trunk/SOURCE/exchange_horiz.f90
r484 r667  1  SUBROUTINE exchange_horiz( ar  1 SUBROUTINE exchange_horiz( ar, nbgp_local) 2 2 3 3 !------------------------------------------------------------------------------! 4 4 ! Current revisions: 5 5 ! ----------------- 6  !  6 ! Dynamic exchange of ghost points with nbgp_local to ensure that no useless  7 ! ghost points exchanged in case of multigrid. type_yz(0) and type_xz(0)  8 ! used for normal grid, the remaining types used for the several grid levels.  9 ! Exchange is done via MPI-Vectors with a dynamic value of ghost points which  10 ! depend on the advection scheme. Exchange of left and right PEs is 10% faster  11 ! with MPI-Vectors than without. 7 12 ! 8 13 ! Former revisions: … …  41 46 INTEGER, DIMENSION(MPI_STATUS_SIZE,4) :: wait_stat 42 47 #endif 43  44  REAL :: ar(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1)45   48 INTEGER :: i,nbgp_local  49 REAL, DIMENSION(nzb:nzt+1,nys-nbgp_local:nyn+nbgp_local, &  50 nxl-nbgp_local:nxr+nbgp_local) :: ar 46 51 47 52 CALL cpu_log( log_point_s(2), 'exchange_horiz', 'start' ) 48 53  54 IF ( exchange_mg == .TRUE. ) THEN  55 i = grid_level  56 ELSE  57 i = 0  58 END IF 49 59 #if defined( __parallel ) 50 60 … …  56 66 !-- within the PE memory 57 67 IF ( bc_lr == 'cyclic' ) THEN 58  ar(:, nys:nyn,nxl-1) = ar(:,nys:nyn,nxr)59  ar(:, nys:nyn,nxr+1) = ar(:,nys:nyn,nxl) 68 ar(:,:,nxl-nbgp_local:nxl-1) = ar(:,:,nxr-nbgp_local+1:nxr)  69 ar(:,:,nxr+1:nxr+nbgp_local) = ar(:,:,nxl:nxl+nbgp_local-1) 60 70 ENDIF 61 71 … …  65 75 ! 66 76 !-- Send left boundary, receive right one 67  CALL MPI_ISEND( & 68  ar(nzb,nys-1,nxl), ngp_yz(grid_level), MPI_REAL, pleft, 0, & 69  comm2d, req(1), ierr ) 70  CALL MPI_IRECV( & 71  ar(nzb,nys-1,nxr+1), ngp_yz(grid_level), MPI_REAL, pright, 0, & 72  comm2d, req(2), ierr )  77 CALL MPI_ISEND(ar(nzb,nys-nbgp_local,nxl),1,type_yz(i),pleft,0,comm2d,&  78 req(1),ierr)  79 CALL MPI_IRECV(ar(nzb,nys-nbgp_local,nxr+1),1,type_yz(i),pright,0,&  80 comm2d,req(2),ierr) 73 81 ! 74 82 !-- Send right boundary, receive left one 75  CALL MPI_ISEND( & 76  ar(nzb,nys-1,nxr), ngp_yz(grid_level), MPI_REAL, pright, 1, & 77  comm2d, req(3), ierr ) 78  CALL MPI_IRECV( & 79  ar(nzb,nys-1,nxl-1), ngp_yz(grid_level), MPI_REAL, pleft, 1, & 80  comm2d, req(4), ierr )  83  84  85 CALL MPI_ISEND(ar(nzb,nys-nbgp_local,nxr+1-nbgp_local),1,type_yz(i),pright, 1, &  86 comm2d, req(3), ierr )  87 CALL MPI_IRECV(ar(nzb,nys-nbgp_local,nxl-nbgp_local),1,type_yz(i),pleft,1,&  88 comm2d,req(4), ierr)  89 81 90 CALL MPI_WAITALL( 4, req, wait_stat, ierr ) 82 91 … …  89 98 !-- within the PE memory 90 99 IF ( bc_ns == 'cyclic' ) THEN 91  ar(:,nys- 1,:) = ar(:,nyn,:)92  ar(:,nyn+1 ,:) = ar(:,nys,:) 100 ar(:,nys-nbgp_local:nys-1,:) = ar(:,nyn-nbgp_local+1:nyn,:)  101 ar(:,nyn+1:nyn+nbgp_local,:) = ar(:,nys:nys+nbgp_local-1,:) 93 102 ENDIF 94 103 … …  98 107 ! 99 108 !-- Send front boundary, receive rear one 100  CALL MPI_ISEND( ar(nzb,nys,nxl-1), 1, type_xz(grid_level), psouth, 0, &  109 !-- MPI_ISEND initial send adress changed, type_xz() is sendet nbgp times  110  111 CALL MPI_ISEND( ar(nzb,nys,nxl-nbgp_local),1, type_xz(i), psouth, 0, & 101 112 comm2d, req(1), ierr ) 102  CALL MPI_IRECV( ar(nzb,nyn+1,nxl- 1), 1, type_xz(grid_level), pnorth, 0, & 113 CALL MPI_IRECV( ar(nzb,nyn+1,nxl-nbgp_local),1, type_xz(i), pnorth, 0, & 103 114 comm2d, req(2), ierr ) 104 115 ! 105 116 !-- Send rear boundary, receive front one 106  CALL MPI_ISEND( ar(nzb,nyn ,nxl-1), 1, type_xz(grid_level), pnorth, 1, & 117 CALL MPI_ISEND( ar(nzb,nyn-nbgp_local+1,nxl-nbgp_local),1, type_xz(i), pnorth, 1, & 107 118 comm2d, req(3), ierr ) 108  CALL MPI_IRECV( ar(nzb,nys- 1,nxl-1), 1, type_xz(grid_level), psouth, 1, & 119 CALL MPI_IRECV( ar(nzb,nys-nbgp_local,nxl-nbgp_local),1, type_xz(i), psouth, 1, & 109 120 comm2d, req(4), ierr ) 110 121 call MPI_WAITALL( 4, req, wait_stat, ierr ) 111 122 112 123 ENDIF 113  114 124 115 125 #else … …  118 128 !-- Lateral boundary conditions in the non-parallel case 119 129 IF ( bc_lr == 'cyclic' ) THEN 120  ar(:,nys:nyn,nxl-1) = ar(:,nys:nyn,nxr)121  ar(:,nys:nyn,nxr+1) = ar(:,nys:nyn,nxl) 130 ar(:,:,nxl-nbgp_local:nxl-1) = ar(:,:,nxr-nbgp_local+1:nxr)  131 ar(:,:,nxr+1:nxr+nbgp_local) = ar(:,:,nxl:nxl+nbgp_local-1) 122 132 ENDIF 123 133 124 134 IF ( bc_ns == 'cyclic' ) THEN 125  ar(:,nys-1,:) = ar(:,nyn,:)126  ar(:,nyn+1,:) = ar(:,nys,:) 135 ar(:,nys-nbgp_local:nys-1,:) = ar(:,nyn-nbgp_local+1:nyn,:)  136 ar(:,nyn+1:nyn+nbgp_local,:) = ar(:,nys:nys+nbgp_local-1,:) 127 137 ENDIF 128 138 129 139 #endif 130  131 140 CALL cpu_log( log_point_s(2), 'exchange_horiz', 'stop' ) 132 141  142 133 143 END SUBROUTINE exchange_horiz -
palm/trunk/SOURCE/exchange_horiz_2d.f90
r484 r667  4 4 ! Current revisions: 5 5 ! ----------------- 6  !  6 ! Dynamic exchange of ghost points with nbgp, which depends on the advection  7 ! scheme. Exchange between left and right PEs is now done with MPI-vectors. 7 8 ! 8 9 ! Former revisions: … …  37 38 IMPLICIT NONE 38 39 39  REAL :: ar(nys-1:nyn+1,nxl-1:nxr+1)  40 REAL :: ar(nysg:nyng,nxlg:nxrg)  41 INTEGER :: i 40 42 41 43 … …  51 53 !-- One-dimensional decomposition along y, boundary values can be exchanged 52 54 !-- within the PE memory 53  ar( nys:nyn,nxl-1) = ar(nys:nyn,nxr)54  ar( nys:nyn,nxr+1) = ar(nys:nyn,nxl) 55 ar(:,nxl-nbgp:nxl-1) = ar(:,nxr-nbgp+1:nxr)  56 ar(:,nxr+1:nxr+nbgp) = ar(:,nxl:nxl+nbgp-1) 55 57 56 58 ELSE 57 59 ! 58 60 !-- Send left boundary, receive right one 59  CALL MPI_SENDRECV( ar(nys,nxl), ngp_y, MPI_REAL, pleft, 0, & 60  ar(nys,nxr+1), ngp_y, MPI_REAL, pright, 0, &  61  62 CALL MPI_SENDRECV( ar(nysg,nxl), 1, type_y, pleft, 0, &  63 ar(nysg,nxr+1), 1, type_y, pright, 0, & 61 64 comm2d, status, ierr ) 62 65 ! 63 66 !-- Send right boundary, receive left one 64  CALL MPI_SENDRECV( ar(nys ,nxr), ngp_y, MPI_REAL, pright, 1, &65  ar(nys ,nxl-1), ngp_y, MPI_REAL, pleft, 1, & 67 CALL MPI_SENDRECV( ar(nysg,nxr+1-nbgp), 1, type_y, pright, 1, &  68 ar(nysg,nxlg), 1, type_y, pleft, 1, & 66 69 comm2d, status, ierr ) 67 70 ENDIF … …  71 74 !-- One-dimensional decomposition along x, boundary values can be exchanged 72 75 !-- within the PE memory 73  ar(nys- 1,:) = ar(nyn,:)74  ar(nyn+1 ,:) = ar(nys,:) 76 ar(nys-nbgp:nys-1,:) = ar(nyn-nbgp+1:nyn,:)  77 ar(nyn+1:nyn+nbgp,:) = ar(nys:nys+nbgp-1,:) 75 78 76 79 ELSE 77 80 ! 78 81 !-- Send front boundary, receive rear one 79  CALL MPI_SENDRECV( ar(nys,nxl-1), 1, type_x, psouth, 0, & 80  ar(nyn+1,nxl-1), 1, type_x, pnorth, 0, &  82  83 CALL MPI_SENDRECV( ar(nys,nxlg), 1, type_x, psouth, 0, & !replace number of sended elements from  84 ar(nyn+1,nxlg), 1, type_x, pnorth, 0, & ! nbgp to 1 81 85 comm2d, status, ierr ) 82 86 ! 83 87 !-- Send rear boundary, receive front one 84  CALL MPI_SENDRECV( ar(nyn,nxl-1), 1, type_x, pnorth, 1, & 85  ar(nys-1,nxl-1), 1, type_x, psouth, 1, & 86  comm2d, status, ierr )  88 CALL MPI_SENDRECV( ar(nyn+1-nbgp,nxlg), 1, type_x, pnorth, 1, &  89 ar(nysg,nxlg), 1, type_x, psouth, 1, &  90 comm2d, status, ierr )  91 87 92 ENDIF 88 93 … …  92 97 !-- Lateral boundary conditions in the non-parallel case 93 98 IF ( bc_lr == 'cyclic' ) THEN 94  ar( nys:nyn,nxl-1) = ar(nys:nyn,nxr)95  ar( nys:nyn,nxr+1) = ar(nys:nyn,nxl) 99 ar(:,nxl-nbgp:nxl-1) = ar(:,nxr-nbgp+1:nxr)  100 ar(:,nxr+1:nxr+nbgp) = ar(:,nxl:nxl+nbgp-1) 96 101 ENDIF 97 102 98 103 IF ( bc_ns == 'cyclic' ) THEN 99  ar(nys-1,:) = ar(nyn,:) 100  ar(nyn+1,:) = ar(nys,:) 101  ENDIF  104 ar(nys-nbgp:nys-1,:) = ar(nyn-nbgp+1:nyn,:)  105 ar(nyn+1:nyn+nbgp,:) = ar(nys:nys+nbgp-1,:)  106 ENDIF  107 102 108 103 109 #endif … …  106 112 !-- Neumann-conditions at inflow/outflow in case of non-cyclic boundary 107 113 !-- conditions 108  IF ( inflow_l .OR. outflow_l ) ar(:,nxl-1) = ar(:,nxl) 109  IF ( inflow_r .OR. outflow_r ) ar(:,nxr+1) = ar(:,nxr) 110  IF ( inflow_s .OR. outflow_s ) ar(nys-1,:) = ar(nys,:) 111  IF ( inflow_n .OR. outflow_n ) ar(nyn+1,:) = ar(nyn,:) 112   114 IF ( inflow_l .OR. outflow_l ) THEN  115 DO i=nbgp, 1, -1  116 ar(:,nxl-i) = ar(:,nxl)  117 END DO  118 END IF  119 IF ( inflow_r .OR. outflow_r ) THEN  120 DO i=1, nbgp  121 ar(:,nxr+i) = ar(:,nxr)  122 END DO  123 END IF  124 IF ( inflow_s .OR. outflow_s ) THEN  125 DO i=nbgp, 1, -1  126 ar(nys-i,:) = ar(nys,:)  127 END DO  128 END IF  129 IF ( inflow_n .OR. outflow_n ) THEN  130 DO i=1, nbgp  131 ar(nyn+i,:) = ar(nyn,:)  132 END DO  133 END IF 113 134 CALL cpu_log( log_point_s(13), 'exchange_horiz_2d', 'stop' ) 114 135 … …  134 155 IMPLICIT NONE 135 156 136  INTEGER :: ar(nys-1:nyn+1,nxl-1:nxr+1)  157 REAL :: ar(nysg:nyng,nxlg:nxrg)  158 INTEGER :: i  159 137 160 138 161 … …  154 177 ! 155 178 !-- Send left boundary, receive right one 156  CALL MPI_SENDRECV( ar(nys ,nxl), ngp_y, MPI_INTEGER, pleft, 0, &157  ar(nys ,nxr+1), ngp_y, MPI_INTEGER, pright, 0, & 179 CALL MPI_SENDRECV( ar(nysg,nxl), 1, type_y_int, pleft, 0, &  180 ar(nysg,nxr+1), 1, type_y_int, pright, 0, & 158 181 comm2d, status, ierr ) 159 182 ! 160 183 !-- Send right boundary, receive left one 161  CALL MPI_SENDRECV( ar(nys,nxr), ngp_y, MPI_INTEGER, pright, 1, & 162  ar(nys,nxl-1), ngp_y, MPI_INTEGER, pleft, 1, & 163  comm2d, status, ierr )  184 CALL MPI_SENDRECV( ar(nysg,nxr+1-nbgp), 1, type_y_int, pright, 1, &  185 ar(nysg,nxlg), 1, type_y_int, pleft, 1, &  186 comm2d, status, ierr )  187 164 188 ENDIF 165 189 … …  168 192 !-- One-dimensional decomposition along x, boundary values can be exchanged 169 193 !-- within the PE memory 170  ar(nys-1,:) = ar(nyn,:) 171  ar(nyn+1,:) = ar(nys,:)  194 ar(nysg:nys-1,:) = ar(nyn+1-nbgp:nyn,:)  195 ar(nyn+1:nyng,:) = ar(nys:nys-1+nbgp,:)  196 172 197 173 198 ELSE 174 199 ! 175 200 !-- Send front boundary, receive rear one 176  CALL MPI_SENDRECV( ar(nys,nxl-1), 1, type_x_int, psouth, 0, & 177  ar(nyn+1,nxl-1), 1, type_x_int, pnorth, 0, & 178  comm2d, status, ierr )  201 CALL MPI_SENDRECV( ar(nys,nxlg),1, type_x_int, psouth, 0, &  202 ar(nyn+1,nxlg),1, type_x_int, pnorth, 0, &  203 comm2d, status, ierr )  204 179 205 ! 180 206 !-- Send rear boundary, receive front one 181  CALL MPI_SENDRECV( ar(nyn,nxl-1), 1, type_x_int, pnorth, 1, & 182  ar(nys-1,nxl-1), 1, type_x_int, psouth, 1, & 183  comm2d, status, ierr )  207 CALL MPI_SENDRECV( ar(nyn+1-nbgp,nxlg), nbgp, type_x_int, pnorth, 1, &  208 ar(nysg,nxlg), nbgp, type_x_int, psouth, 1, &  209 comm2d, status, ierr )  210 184 211 ENDIF 185 212 … …  194 221 195 222 IF ( bc_ns == 'cyclic' ) THEN 196  ar(nys -1,:) = ar(nyn,:)197  ar(nyn+1 ,:) = ar(nys,:) 223 ar(nysg:nys-1,:) = ar(nyn+1-nbgp:nyn,:)  224 ar(nyn+1:nyng,:) = ar(nys:nys-1+nbgp,:) 198 225 ENDIF 199 226 -
palm/trunk/SOURCE/flow_statistics.f90
r625 r667  4 4 ! Current revisions: 5 5 ! ----------------- 6  !  6 ! When advection is computed with ws-scheme, turbulent fluxes are already  7 ! computed in the respective advection routines and buffered in arrays  8 ! sums_xx_ws_l(). This is due to a consistent treatment of statistics with the  9 ! numerics and to avoid unphysical kinks near the surface.  10 ! So some if requests has to be done to dicern between fluxes from ws-scheme  11 ! other advection schemes.  12 ! Furthermore the computation of z_i is only done if the heat flux exceeds a  13 ! minimum value. This affects only simulations of a neutral boundary layer and  14 ! is due to reasons of computations in the advection scheme.  15 ! 7 16 ! 8 17 ! Former revisions: … …  102 111 REAL :: sums_ll(nzb:nzt+1,2) 103 112 104  105 113 CALL cpu_log( log_point(10), 'flow_statistics', 'start' ) 106 114 … …  133 141 sums_l(nzb+10,pr_palm,0) = sums_divnew_l(sr) ! new divergence from pres 134 142  143 !  144 !-- Copy the turbulent quantities, evaluated in the advection routines to  145 !-- the local array sums_l() for further computations  146 IF ( ws_scheme_mom ) THEN  147 !  148 !-- Boundary condition for u'u' and v'v', because below the surface no  149 !-- computation for these quantities is done.  150 DO i = nxl, nxr  151 DO j = nys, nyn  152 sums_us2_ws_l(nzb_u_inner(j,i),sr) = &  153 sums_us2_ws_l(nzb_u_inner(j,i)+1,sr)  154 sums_vs2_ws_l(nzb_v_inner(j,i),sr) = &  155 sums_vs2_ws_l(nzb_v_inner(j,i)+1,sr)  156 ENDDO  157 ENDDO  158 !  159 !-- Swap the turbulent quantities evaluated in advec_ws.  160 sums_l(:,13,0) = sums_wsus_ws_l(:,sr) ! w*u*  161 sums_l(:,15,0) = sums_wsvs_ws_l(:,sr) ! w*v*  162 sums_l(:,30,0) = sums_us2_ws_l(:,sr) ! u*2  163 sums_l(:,31,0) = sums_vs2_ws_l(:,sr) ! v*2  164 sums_l(:,32,0) = sums_ws2_ws_l(:,sr) ! w*2  165 sums_l(:,34,0) = sums_l(:,34,0) + 0.5 * &  166 (sums_us2_ws_l(:,sr) + sums_vs2_ws_l(:,sr) &  167 + sums_ws2_ws_l(:,sr)) ! e*  168 DO k = nzb, nzt  169 sums_l(nzb+5,pr_palm,0) = sums_l(nzb+5,pr_palm,0) + 0.5 * ( &  170 sums_us2_ws_l(k,sr) + sums_vs2_ws_l(k,sr) + &  171 sums_ws2_ws_l(k,sr))  172 ENDDO  173 ENDIF  174 IF ( ws_scheme_sca ) THEN  175 sums_l(:,17,0) = sums_wspts_ws_l(:,sr) ! w*pts* from advec_s_ws  176 IF ( ocean ) sums_l(:,66,0) = sums_wssas_ws_l(:,sr) ! w*sa*  177 IF ( humidity .OR. passive_scalar ) sums_l(:,49,0) = &  178 sums_wsqs_ws_l(:,sr) !w*q*  179 ENDIF 135 180 ! 136 181 !-- Horizontally averaged profiles of horizontal velocities and temperature. … …  138 183 !-- for other horizontal averages. 139 184 tn = 0  185 140 186 !$OMP PARALLEL PRIVATE( i, j, k, tn ) 141 187 #if defined( __intel_openmp_bug ) … …  215 261 ENDIF 216 262 !$OMP END PARALLEL 217  218 263 ! 219 264 !-- Summation of thread sums … …  304 349 hom(:,1,2,sr) = sums(:,2) ! v 305 350 hom(:,1,4,sr) = sums(:,4) ! pt  351 306 352 307 353 ! … …  354 400 DO j = nys, nyn 355 401 sums_l_etot = 0.0 356  sums_l_eper = 0.0357 402 DO k = nzb_s_inner(j,i), nzt+1 358  u2 = u(k,j,i)**2359  v2 = v(k,j,i)**2360  w2 = w(k,j,i)**2361  ust2 = ( u(k,j,i) - hom(k,1,1,sr) )**2362  vst2 = ( v(k,j,i) - hom(k,1,2,sr) )**2363 403 ! 364 404 !-- Prognostic and diagnostic variables … …  369 409 sums_l(k,40,tn) = sums_l(k,40,tn) + p(k,j,i) 370 410 371  !372  !-- Variances373  sums_l(k,30,tn) = sums_l(k,30,tn) + ust2 * rmask(j,i,sr)374  sums_l(k,31,tn) = sums_l(k,31,tn) + vst2 * rmask(j,i,sr)375  sums_l(k,32,tn) = sums_l(k,32,tn) + w2 * rmask(j,i,sr)376 411 sums_l(k,33,tn) = sums_l(k,33,tn) + & 377 412 ( pt(k,j,i)-hom(k,1,4,sr) )**2 * rmask(j,i,sr) … …  381 416 ( q(k,j,i)-hom(k,1,41,sr) )**2 * rmask(j,i,sr) 382 417 ENDIF 383  ! 384  !-- Higher moments 385  !-- (Computation of the skewness of w further below) 386  sums_l(k,38,tn) = sums_l(k,38,tn) + w(k,j,i) * w2 * & 387  rmask(j,i,sr) 388  ! 389  !-- Perturbation energy 390  sums_l(k,34,tn) = sums_l(k,34,tn) + 0.5 * ( ust2 + vst2 + w2 ) & 391  * rmask(j,i,sr)  418 392 419 sums_l_etot = sums_l_etot + & 393  0.5 * ( u2 + v2 + w2 ) * rmask(j,i,sr) 394  sums_l_eper = sums_l_eper + & 395  0.5 * ( ust2+vst2+w2 ) * rmask(j,i,sr)  420 0.5 * ( u(k,j,i)**2 + v(k,j,i)**2 + &  421 w(k,j,i)**2 ) * rmask(j,i,sr) 396 422 ENDDO 397 423 ! … …  401 427 !-- allow vectorization of that loop. 402 428 sums_l(nzb+4,pr_palm,tn) = sums_l(nzb+4,pr_palm,tn) + sums_l_etot 403  sums_l(nzb+5,pr_palm,tn) = sums_l(nzb+5,pr_palm,tn) + sums_l_eper404 429 ! 405 430 !-- 2D-arrays (being collected in the last column of sums_l) … …  420 445 421 446 !  447 !-- Computation of statistics when ws-scheme is not used. Else these  448 !-- quantities are evaluated in the advection routines.  449 IF ( .NOT. ws_scheme_mom ) THEN  450 !$OMP DO  451 DO i = nxl, nxr  452 DO j = nys, nyn  453 sums_l_eper = 0.0  454 DO k = nzb_s_inner(j,i), nzt+1  455 u2 = u(k,j,i)**2  456 v2 = v(k,j,i)**2  457 w2 = w(k,j,i)**2  458 ust2 = ( u(k,j,i) - hom(k,1,1,sr) )**2  459 vst2 = ( v(k,j,i) - hom(k,1,2,sr) )**2  460  461 sums_l(k,30,tn) = sums_l(k,30,tn) + ust2 * rmask(j,i,sr)  462 sums_l(k,31,tn) = sums_l(k,31,tn) + vst2 * rmask(j,i,sr)  463 sums_l(k,32,tn) = sums_l(k,32,tn) + w2 * rmask(j,i,sr)  464 !  465 !-- Higher moments  466 !-- (Computation of the skewness of w further below)  467 sums_l(k,38,tn) = sums_l(k,38,tn) + w(k,j,i) * w2 * &  468 rmask(j,i,sr)  469 !  470 !-- Perturbation energy  471  472 sums_l(k,34,tn) = sums_l(k,34,tn) + 0.5 * &  473 ( ust2 + vst2 + w2 ) * rmask(j,i,sr)  474 sums_l_eper = sums_l_eper + &  475 0.5 * ( ust2+vst2+w2 ) * rmask(j,i,sr)  476  477 ENDDO  478 sums_l(nzb+5,pr_palm,tn) = sums_l(nzb+5,pr_palm,tn) &  479 + sums_l_eper  480 ENDDO  481 ENDDO  482 ELSE  483 !$OMP DO  484 DO i = nxl, nxr  485 DO j = nys, nyn  486 DO k = nzb_s_inner(j,i), nzt + 1  487 w2 = w(k,j,i)**2  488 !  489 !-- Higher moments  490 !-- (Computation of the skewness of w further below)  491 sums_l(k,38,tn) = sums_l(k,38,tn) + w(k,j,i) * w2 * &  492 rmask(j,i,sr)  493 ENDDO  494 ENDDO  495 ENDDO  496 ENDIF  497  498 ! 422 499 !-- Horizontally averaged profiles of the vertical fluxes  500 423 501 !$OMP DO 424 502 DO i = nxl, nxr … …  587 665 pts = 0.5 * ( pt(k,j,i) - hom(k,1,4,sr) + & 588 666 pt(k+1,j,i) - hom(k+1,1,4,sr) ) 589  ! 590  !-- Momentum flux w*u* 591  sums_l(k,13,tn) = sums_l(k,13,tn) + 0.5 * & 592  ( w(k,j,i-1) + w(k,j,i) ) & 593  * ust * rmask(j,i,sr) 594  ! 595  !-- Momentum flux w*v* 596  sums_l(k,15,tn) = sums_l(k,15,tn) + 0.5 * & 597  ( w(k,j-1,i) + w(k,j,i) ) & 598  * vst * rmask(j,i,sr) 599  ! 600  !-- Heat flux w*pt* 601  !-- The following formula (comment line, not executed) does not 602  !-- work if applied to subregions 603  ! sums_l(k,17,tn) = sums_l(k,17,tn) + 0.5 * & 604  ! ( pt(k,j,i)+pt(k+1,j,i) ) & 605  ! * w(k,j,i) * rmask(j,i,sr) 606  sums_l(k,17,tn) = sums_l(k,17,tn) + pts * w(k,j,i) * & 607  rmask(j,i,sr) 608  !  667 609 668 !-- Higher moments 610 669 sums_l(k,35,tn) = sums_l(k,35,tn) + pts * w(k,j,i)**2 * & … …  617 676 !-- but so far there is no other suitable place to calculate) 618 677 IF ( ocean ) THEN 619  pts = 0.5 * ( sa(k,j,i) - hom(k,1,23,sr) + &  678 IF( .NOT. ws_scheme_sca ) THEN  679 pts = 0.5 * ( sa(k,j,i) - hom(k,1,23,sr) + & 620 680 sa(k+1,j,i) - hom(k+1,1,23,sr) ) 621  sums_l(k,66,tn) = sums_l(k,66,tn) + pts * w(k,j,i) * & 681 sums_l(k,66,tn) = sums_l(k,66,tn) + pts * w(k,j,i) * & 622 682 rmask(j,i,sr)  683 ENDIF 623 684 sums_l(k,64,tn) = sums_l(k,64,tn) + rho(k,j,i) * & 624 685 rmask(j,i,sr) … …  635 696 sums_l(k,46,tn) = sums_l(k,46,tn) + pts * w(k,j,i) * & 636 697 rmask(j,i,sr) 637  pts = 0.5 * ( q(k,j,i) - hom(k,1,41,sr) + & 638  q(k+1,j,i) - hom(k+1,1,41,sr) ) 639  sums_l(k,49,tn) = sums_l(k,49,tn) + pts * w(k,j,i) * & 640  rmask(j,i,sr)  698 641 699 IF ( cloud_physics .OR. cloud_droplets ) THEN 642 700 pts = 0.5 * & … …  652 710 ! 653 711 !-- Passive scalar flux 654  IF ( passive_scalar ) THEN 712 IF ( passive_scalar .AND. ( .NOT. ws_scheme_sca )) THEN 655 713 pts = 0.5 * ( q(k,j,i) - hom(k,1,41,sr) + & 656 714 q(k+1,j,i) - hom(k+1,1,41,sr) ) … …  661 719 ! 662 720 !-- Energy flux w*e* 663  sums_l(k,37,tn) = sums_l(k,37,tn) + w(k,j,i) * 0.5 * & 664  ( ust**2 + vst**2 + w(k,j,i)**2 )&665  * rmask(j,i,sr)666   721 !-- has to be adjusted  722 sums_l(k,37,tn) = sums_l(k,37,tn) + w(k,j,i) * 0.5 * &  723 ( ust**2 + vst**2 + w(k,j,i)**2 )&  724 * rmask(j,i,sr) 667 725 ENDDO 668 726 ENDDO 669 727 ENDDO  728 !- for reasons of speed optimization the loop is splitted, to avoid if-else  729 !- statements inside the loop  730 !- Fluxes which have been computed in part directly inside the advection routines  731 !- treated seperatly.  732 !- First treat the momentum fluxes  733 IF ( .NOT. ws_scheme_mom ) THEN  734 !$OMP DO  735 DO i = nxl, nxr  736 DO j = nys, nyn  737 DO k = nzb_diff_s_inner(j,i)-1, nzt_diff  738 ust = 0.5 * ( u(k,j,i) - hom(k,1,1,sr) + &  739 u(k+1,j,i) - hom(k+1,1,1,sr) )  740 vst = 0.5 * ( v(k,j,i) - hom(k,1,2,sr) + &  741 v(k+1,j,i) - hom(k+1,1,2,sr) )  742 !  743 !-- Momentum flux w*u*  744 sums_l(k,13,tn) = sums_l(k,13,tn) + 0.5 * &  745 ( w(k,j,i-1) + w(k,j,i) ) &  746 * ust * rmask(j,i,sr)  747 !  748 !-- Momentum flux w*v*  749 sums_l(k,15,tn) = sums_l(k,15,tn) + 0.5 * &  750 ( w(k,j-1,i) + w(k,j,i) ) &  751 * vst * rmask(j,i,sr)  752 ENDDO  753 ENDDO  754 ENDDO  755  756 ENDIF  757 IF ( .NOT. ws_scheme_sca ) THEN  758 !$OMP DO  759 DO i = nxl, nxr  760 DO j = nys, nyn  761 DO k = nzb_diff_s_inner(j,i) - 1, nzt_diff  762 !- vertical heat flux  763 sums_l(k,17,tn) = sums_l(k,17,tn) + 0.5 * &  764 ( pt(k,j,i) - hom(k,1,4,sr) + &  765 pt(k+1,j,i) - hom(k+1,1,4,sr) ) &  766 * w(k,j,i) * rmask(j,i,sr)  767 IF ( humidity ) THEN  768 pts = 0.5 * ( q(k,j,i) - hom(k,1,41,sr) + &  769 q(k+1,j,i) - hom(k+1,1,41,sr) )  770 sums_l(k,49,tn) = sums_l(k,49,tn) + pts * w(k,j,i) * &  771 rmask(j,i,sr)  772 ENDIF  773 ENDDO  774 ENDDO  775 ENDDO  776  777 ENDIF  778 670 779 671 780 ! … …  814 923 815 924 #if defined( __parallel )  925 816 926 ! 817 927 !-- Compute total sum from local sums … …  843 953 sums(k,70:pr_palm-2) = sums(k,70:pr_palm-2)/ ngp_2dh_s_inner(k,sr) 844 954 ENDDO  955 845 956 !-- Upstream-parts 846 957 sums(nzb:nzb+11,pr_palm-1) = sums(nzb:nzb+11,pr_palm-1) / ngp_3d(sr) … …  868 979 ENDDO 869 980 ENDIF 870  871 981 ! 872 982 !-- Collect horizontal average in hom. … …  934 1044 ! upstream-parts u_x, u_y, u_z, v_x, 935 1045 ! v_y, usw. (in last but one profile) 936  hom(:,1,pr_palm,sr) = sums(:,pr_palm)  1046 hom(:,1,pr_palm,sr) = sums(:,pr_palm) 937 1047 ! u*, w'u', w'v', t* (in last profile) 938 1048 … …  951 1061 z_i(1) = 0.0 952 1062 first = .TRUE.  1063 953 1064 IF ( ocean ) THEN 954 1065 DO k = nzt, nzb+1, -1 955  IF ( first .AND. hom(k,1,18,sr) < 0.0 ) THEN  1066 IF ( first .AND. hom(k,1,18,sr) < 0.0 &  1067 .AND. abs(hom(k,1,18,sr)) > 1.0E-8) THEN 956 1068 first = .FALSE. 957 1069 height = zw(k) 958 1070 ENDIF 959 1071 IF ( hom(k,1,18,sr) < 0.0 .AND. &  1072 abs(hom(k,1,18,sr)) > 1.0E-8 .AND. & 960 1073 hom(k-1,1,18,sr) > hom(k,1,18,sr) ) THEN 961 1074 IF ( zw(k) < 1.5 * height ) THEN … …  969 1082 ELSE 970 1083 DO k = nzb, nzt-1 971  IF ( first .AND. hom(k,1,18,sr) < 0.0 ) THEN  1084 IF ( first .AND. hom(k,1,18,sr) < 0.0 &  1085 .AND. abs(hom(k,1,18,sr)) > 1.0E-8 ) THEN 972 1086 first = .FALSE. 973 1087 height = zw(k) 974 1088 ENDIF 975  IF ( hom(k,1,18,sr) < 0.0 .AND. &  1089 IF ( hom(k,1,18,sr) < 0.0 .AND. &  1090 abs(hom(k,1,18,sr)) > 1.0E-8 .AND. & 976 1091 hom(k+1,1,18,sr) > hom(k,1,18,sr) ) THEN 977 1092 IF ( zw(k) < 1.5 * height ) THEN … …  1026 1141 !-- the characteristic convective boundary layer temperature. 1027 1142 !-- The horizontal average at nzb+1 is input for the average temperature. 1028  IF ( hom(nzb,1,18,sr) > 0.0 .AND. z_i(1) /= 0.0 ) THEN  1143 IF ( hom(nzb,1,18,sr) > 0.0 .AND. abs(hom(nzb,1,18,sr)) > 1.0E-8 &  1144 .AND. z_i(1) /= 0.0 ) THEN 1029 1145 hom(nzb+8,1,pr_palm,sr) = ( g / hom(nzb+1,1,4,sr) * & 1030 1146 hom(nzb,1,18,sr) * & -
palm/trunk/SOURCE/global_min_max.f90
r623 r667  5 5 ! Current revisions: 6 6 ! ----------------- 7  !  7 ! Adapting of the index arrays, because MINLOC assumes lowerbound at 1 and not  8 ! at nbgp. 8 9 ! 9 10 ! Former revisions: … …  59 60 !-- Determine the local minimum 60 61 fmin_ijk_l = MINLOC( ar ) 61  fmin_ijk_l(1) = i1 + fmin_ijk_l(1) - 1! MINLOC assumes lowerbound = 162  fmin_ijk_l(2) = j1 + fmin_ijk_l(2) - 1 62 fmin_ijk_l(1) = i1 + fmin_ijk_l(1) - nbgp ! MINLOC assumes lowerbound = 1  63 fmin_ijk_l(2) = j1 + fmin_ijk_l(2) - nbgp 63 64 fmin_ijk_l(3) = k1 + fmin_ijk_l(3) - 1 64 65 fmin_l(1) = ar(fmin_ijk_l(1),fmin_ijk_l(2),fmin_ijk_l(3)) … …  100 101 !-- Determine the local maximum 101 102 fmax_ijk_l = MAXLOC( ar ) 102  fmax_ijk_l(1) = i1 + fmax_ijk_l(1) - 1! MAXLOC assumes lowerbound = 1103  fmax_ijk_l(2) = j1 + fmax_ijk_l(2) - 1 103 fmax_ijk_l(1) = i1 + fmax_ijk_l(1) - nbgp ! MAXLOC assumes lowerbound = 1  104 fmax_ijk_l(2) = j1 + fmax_ijk_l(2) - nbgp 104 105 fmax_ijk_l(3) = k1 + fmax_ijk_l(3) - 1 105 106 fmax_l(1) = ar(fmax_ijk_l(1),fmax_ijk_l(2),fmax_ijk_l(3)) … …  221 222 IF ( fmax_ijk(1) < 0 ) THEN 222 223 value = -value 223  value_ijk(1) = -value_ijk(1) - 10  224 value_ijk(1) = -value_ijk(1) - 10 !??? 224 225 ENDIF 225 226 … …  228 229 ! 229 230 !-- Limit index values to the range 0..nx, 0..ny 230  IF ( value_ijk(3) == -1 ) value_ijk(3) = nx231  IF ( value_ijk(3) == nx+1 ) value_ijk(3) = 0232  IF ( value_ijk(2) == -1 ) value_ijk(2) = ny233  IF ( value_ijk(2) == ny+1 ) value_ijk(2) = 0 231 IF ( value_ijk(3) < 0 ) value_ijk(3) = nx +1 + value_ijk(3)  232 IF ( value_ijk(3) > nx ) value_ijk(3) = value_ijk(3) - (nx+1)  233 IF ( value_ijk(2) < 0 ) value_ijk(2) = ny +1 + value_ijk(2)  234 IF ( value_ijk(2) > ny ) value_ijk(2) = value_ijk(2) - (ny+1) 234 235 235 236 -
palm/trunk/SOURCE/header.f90
r581 r667  4 4 ! Current revisions: 5 5 ! ----------------- 6  !  6 ! Output of advection scheme.  7 ! Modified output of Prandtl-layer height. 7 8 ! 8 9 ! Former revisions: … …  250 251 IF ( momentum_advec == 'pw-scheme' ) THEN 251 252 WRITE ( io, 113 ) 252  ELSE  253 ELSEIF (momentum_advec == 'ws-scheme' ) THEN  254 WRITE ( io, 503 )  255 ELSEIF (momentum_advec == 'ups-scheme' ) THEN 253 256 WRITE ( io, 114 ) 254 257 IF ( cut_spline_overshoot ) WRITE ( io, 124 ) … …  267 270 IF ( scalar_advec == 'pw-scheme' ) THEN 268 271 WRITE ( io, 116 )  272 ELSEIF ( scalar_advec == 'ws-scheme' ) THEN  273 WRITE ( io, 504 ) 269 274 ELSEIF ( scalar_advec == 'ups-scheme' ) THEN 270 275 WRITE ( io, 117 ) … …  575 580 ELSEIF( ibc_pt_t == 2 ) THEN 576 581 roben = TRIM( roben ) // ' pt(nzt+1) = pt(nzt) + dpt/dz_ini'  582 577 583 ENDIF 578 584 … …  662 668 663 669 IF ( prandtl_layer ) THEN 664  WRITE ( io, 305 ) 0.5 *(zu(1)-zu(0)), roughness_length, kappa, & 670 WRITE ( io, 305 ) (zu(1)-zu(0)), roughness_length, kappa, & 665 671 rif_min, rif_max 666 672 IF ( .NOT. constant_heatflux ) WRITE ( io, 308 ) … …  1981 1987 ' Dissipation calculation: ',A/) 1982 1988 502 FORMAT (' Damping layer starts from ',F7.1,' m (GP ',I4,')'/)  1989 503 FORMAT (' --> Momentum advection via Wicker-Skamarock-Scheme 5th order')  1990 504 FORMAT (' --> Scalar advection via Wicker-Skamarock-Scheme 5th order') 1983 1991 1984 1992 -
palm/trunk/SOURCE/inflow_turbulence.f90
r623 r667  4 4 ! Current revisions: 5 5 ! ----------------- 6  !  6 ! Using nbgp recycling planes for a better resolution of the turbulent flow  7 ! near the inflow. 7 8 ! 8 9 ! Former revisions: … …  35 36 IMPLICIT NONE 36 37 37  INTEGER :: i, imax, j, k, ngp_ifd, ngp_pr 38 INTEGER :: i, imax, j, k, l, ngp_ifd, ngp_pr 38 39 39 40 REAL, DIMENSION(1:2) :: volume_flow_l, volume_flow_offset 40  REAL, DIMENSION(nzb:nzt+1,5 ) :: avpr, avpr_l41  REAL, DIMENSION(nzb:nzt+1,nys -1:nyn+1,5) :: inflow_dist 41 REAL, DIMENSION(nzb:nzt+1,5,nbgp) :: avpr, avpr_l  42 REAL, DIMENSION(nzb:nzt+1,nysg:nyng,5,nbgp) :: inflow_dist 42 43 43 44 CALL cpu_log( log_point(40), 'inflow_turbulence', 'start' ) 44 45 45 46 ! 46  !-- Carry out horizontalaveraging in the recycling plane 47 !-- Carry out spanwise averaging in the recycling plane 47 48 avpr_l = 0.0 48  ngp_pr = ( nzt - nzb + 2 ) * 5 49  ngp_ifd = ngp_pr * ( nyn - nys + 3) 49 ngp_pr = ( nzt - nzb + 2 ) * 5 * nbgp  50 ngp_ifd = ngp_pr * ( nyn - nys + 1 + 2 * nbgp ) 50 51 51 52 ! 52 53 !-- First, local averaging within the recycling domain 53  IF ( recycling_plane >= nxl ) THEN 54  55  imax = MIN( nxr, recycling_plane ) 56  57  DO i = nxl, imax  54  55 i = recycling_plane  56  57 #if defined( __parallel )  58 IF ( myidx == id_recycling ) THEN  59  60 DO l = 1, nbgp 58 61 DO j = nys, nyn 59  DO k = nzb, nzt +160  61  avpr_l(k,1 ) = avpr_l(k,1) + u(k,j,i)62  avpr_l(k,2 ) = avpr_l(k,2) + v(k,j,i)63  avpr_l(k,3 ) = avpr_l(k,3) + w(k,j,i)64  avpr_l(k,4 ) = avpr_l(k,4) + pt(k,j,i)65  avpr_l(k,5 ) = avpr_l(k,5) + e(k,j,i) 62 DO k = nzb, nzt + 1  63  64 avpr_l(k,1,l) = avpr_l(k,1,l) + u(k,j,i)  65 avpr_l(k,2,l) = avpr_l(k,2,l) + v(k,j,i)  66 avpr_l(k,3,l) = avpr_l(k,3,l) + w(k,j,i)  67 avpr_l(k,4,l) = avpr_l(k,4,l) + pt(k,j,i)  68 avpr_l(k,5,l) = avpr_l(k,5,l) + e(k,j,i) 66 69 67 70 ENDDO 68 71 ENDDO 69  ENDDO 70  71  ENDIF 72  73  ! WRITE (9,*) '*** averaged profiles avpr_l' 74  ! DO k = nzb, nzt+1 75  ! WRITE (9,'(F5.1,1X,F5.1,1X,F5.1,1X,F6.1,1X,F7.2)') avpr_l(k,1),avpr_l(k,2),avpr_l(k,3),avpr_l(k,4),avpr_l(k,5) 76  ! ENDDO 77  ! WRITE (9,*) ' ' 78  79  #if defined( __parallel )  72 i = i + 1  73 ENDDO  74  75 ENDIF 80 76 ! 81 77 !-- Now, averaging over all PEs 82 78 IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr ) 83  CALL MPI_ALLREDUCE( avpr_l(nzb,1), avpr(nzb,1), ngp_pr, MPI_REAL, MPI_SUM, & 84  comm2d, ierr )  79 CALL MPI_ALLREDUCE( avpr_l(nzb,1,1), avpr(nzb,1,1), ngp_pr, &  80 MPI_REAL, MPI_SUM, comm2d, ierr )  81 85 82 #else  83 DO l = 1, nbgp  84 DO j = nys, nyn  85 DO k = nzb, nzt + 1  86  87 avpr_l(k,1,l) = avpr_l(k,1,l) + u(k,j,i)  88 avpr_l(k,2,l) = avpr_l(k,2,l) + v(k,j,i)  89 avpr_l(k,3,l) = avpr_l(k,3,l) + w(k,j,i)  90 avpr_l(k,4,l) = avpr_l(k,4,l) + pt(k,j,i)  91 avpr_l(k,5,l) = avpr_l(k,5,l) + e(k,j,i)  92  93 ENDDO  94 ENDDO  95 i = i + 1  96 ENDDO  97 86 98 avpr = avpr_l 87 99 #endif 88 100 89  avpr = avpr / ( ( ny + 1 ) * ( recycling_plane + 1 ) ) 90  91  ! WRITE (9,*) '*** averaged profiles' 92  ! DO k = nzb, nzt+1 93  ! WRITE (9,'(F5.1,1X,F5.1,1X,F5.1,1X,F6.1,1X,F7.2)') avpr(k,1),avpr(k,2),avpr(k,3),avpr(k,4),avpr(k,5) 94  ! ENDDO 95  ! WRITE (9,*) ' ' 96   101 avpr = avpr / ( ny + 1 ) 97 102 ! 98 103 !-- Calculate the disturbances at the recycling plane … …  101 106 #if defined( __parallel ) 102 107 IF ( myidx == id_recycling ) THEN 103  104  DO j = nys-1, nyn+1  108 DO l = 1, nbgp  109 DO j = nysg, nyng  110 DO k = nzb, nzt + 1  111  112 inflow_dist(k,j,1,l) = u(k,j,i+1) - avpr(k,1,l)  113 inflow_dist(k,j,2,l) = v(k,j,i) - avpr(k,2,l)  114 inflow_dist(k,j,3,l) = w(k,j,i) - avpr(k,3,l)  115 inflow_dist(k,j,4,l) = pt(k,j,i) - avpr(k,4,l)  116 inflow_dist(k,j,5,l) = e(k,j,i) - avpr(k,5,l)  117  118 ENDDO  119 ENDDO  120 i = i + 1  121 ENDDO  122  123 ENDIF  124 #else  125 DO l = 1, nbgp  126 DO j = nysg, nyng 105 127 DO k = nzb, nzt+1 106 128 107  inflow_dist(k,j,1) = u(k,j,i+1) - avpr(k,1) 108  inflow_dist(k,j,2) = v(k,j,i) - avpr(k,2) 109  inflow_dist(k,j,3) = w(k,j,i) - avpr(k,3) 110  inflow_dist(k,j,4) = pt(k,j,i) - avpr(k,4) 111  inflow_dist(k,j,5) = e(k,j,i) - avpr(k,5) 112  113  ENDDO 114  ENDDO 115  116  ENDIF 117  #else 118  DO j = nys-1, nyn+1 119  DO k = nzb, nzt+1 120  121  inflow_dist(k,j,1) = u(k,j,i+1) - avpr(k,1) 122  inflow_dist(k,j,2) = v(k,j,i) - avpr(k,2) 123  inflow_dist(k,j,3) = w(k,j,i) - avpr(k,3) 124  inflow_dist(k,j,4) = pt(k,j,i) - avpr(k,4) 125  inflow_dist(k,j,5) = e(k,j,i) - avpr(k,5) 126  127  ENDDO  129 inflow_dist(k,j,1,l) = u(k,j,i+1) - avpr(k,1,l)  130 inflow_dist(k,j,2,l) = v(k,j,i) - avpr(k,2,l)  131 inflow_dist(k,j,3,l) = w(k,j,i) - avpr(k,3,l)  132 inflow_dist(k,j,4,l) = pt(k,j,i) - avpr(k,4,l)  133 inflow_dist(k,j,5,l) = e(k,j,i) - avpr(k,5,l)  134  135 ENDDO  136 ENDDO  137 i = i + 1 128 138 ENDDO 129 139 #endif … …  134 144 IF ( myidx == id_recycling .AND. myidx /= id_inflow ) THEN 135 145 136  CALL MPI_SEND( inflow_dist(nzb,nys -1,1), ngp_ifd, MPI_REAL, & 146 CALL MPI_SEND( inflow_dist(nzb,nysg,1,1), ngp_ifd, MPI_REAL, & 137 147 id_inflow, 1, comm1dx, ierr ) 138 148 … …  140 150 141 151 inflow_dist = 0.0 142  CALL MPI_RECV( inflow_dist(nzb,nys -1,1), ngp_ifd, MPI_REAL, & 152 CALL MPI_RECV( inflow_dist(nzb,nysg,1,1), ngp_ifd, MPI_REAL, & 143 153 id_recycling, 1, comm1dx, status, ierr ) 144 154 … …  150 160 IF ( nxl == 0 ) THEN 151 161 152  DO j = nys-1, nyn+1 153  DO k = nzb, nzt+1 154  155  ! WRITE (9,*) 'j=',j,' k=',k 156  ! WRITE (9,*) 'mean_u = ', mean_inflow_profiles(k,1), ' dist_u = ',& 157  ! inflow_dist(k,j,1) 158  ! WRITE (9,*) 'mean_v = ', mean_inflow_profiles(k,2), ' dist_v = ',& 159  ! inflow_dist(k,j,2) 160  ! WRITE (9,*) 'mean_w = 0.0', ' dist_w = ',& 161  ! inflow_dist(k,j,3) 162  ! WRITE (9,*) 'mean_pt = ', mean_inflow_profiles(k,4), ' dist_pt = ',& 163  ! inflow_dist(k,j,4) 164  ! WRITE (9,*) 'mean_e = ', mean_inflow_profiles(k,5), ' dist_e = ',& 165  ! inflow_dist(k,j,5) 166  u(k,j,0) = mean_inflow_profiles(k,1) + & 167  inflow_dist(k,j,1) * inflow_damping_factor(k) 168  v(k,j,-1) = mean_inflow_profiles(k,2) + & 169  inflow_dist(k,j,2) * inflow_damping_factor(k) 170  w(k,j,-1) = inflow_dist(k,j,3) * inflow_damping_factor(k) 171  pt(k,j,-1) = mean_inflow_profiles(k,4) + & 172  inflow_dist(k,j,4) * inflow_damping_factor(k) 173  e(k,j,-1) = mean_inflow_profiles(k,5) + & 174  inflow_dist(k,j,5) * inflow_damping_factor(k) 175  e(k,j,-1) = MAX( e(k,j,-1), 0.0 )  162 DO j = nysg, nyng  163 DO k = nzb, nzt + 1  164  165 u(k,j,-nbgp+1:0) = mean_inflow_profiles(k,1) + &  166 inflow_dist(k,j,1,1:nbgp) * inflow_damping_factor(k)  167 v(k,j,-nbgp:-1) = mean_inflow_profiles(k,2) + &  168 inflow_dist(k,j,2,1:nbgp) * inflow_damping_factor(k)  169 w(k,j,-nbgp:-1) = inflow_dist(k,j,3,1:nbgp) * inflow_damping_factor(k)  170 pt(k,j,-nbgp:-1) = mean_inflow_profiles(k,4) + &  171 inflow_dist(k,j,4,1:nbgp) * inflow_damping_factor(k)  172 e(k,j,-nbgp:-1) = mean_inflow_profiles(k,5) + &  173 inflow_dist(k,j,5,1:nbgp) * inflow_damping_factor(k)  174 e(k,j,-nbgp:-1) = MAX( e(k,j,-nbgp:-1), 0.0 ) 176 175 177 176 ENDDO -
palm/trunk/SOURCE/init_1d_model.f90
r392 r667  63 63 l1d_m(nzb:nzt+1), rif1d(nzb:nzt+1), te_e(nzb:nzt+1), & 64 64 te_em(nzb:nzt+1), te_u(nzb:nzt+1), te_um(nzb:nzt+1), & 65  te_v(nzb:nzt+1), te_vm(nzb:nzt+1), u1d(nzb:nzt+1), & 65 te_v(nzb:nzt+1), te_vm(nzb:nzt+1), u1d(nzb:nzt+1), & 66 66 u1d_m(nzb:nzt+1), u1d_p(nzb:nzt+1), v1d(nzb:nzt+1), & 67 67 v1d_m(nzb:nzt+1), v1d_p(nzb:nzt+1) ) … …  385 385 !-- boundary condition applies to u and v. 386 386 !-- The boundary condition for e is set further below ( (u*/cm)**2 ). 387  u1d_p(nzb) = -u1d_p(nzb+1) 388  v1d_p(nzb) = -v1d_p(nzb+1)  387 ! u1d_p(nzb) = -u1d_p(nzb+1)  388 ! v1d_p(nzb) = -v1d_p(nzb+1)  389  390 u1d_p(nzb) = 0.0  391 v1d_p(nzb) = 0.0 389 392 390 393 ! -
palm/trunk/SOURCE/init_3d_model.f90
r631 r667  7 7 ! Current revisions: 8 8 ! ----------------- 9  ! Bugfix: type conversion of '1' to 64bit for the MAX function (ngp_3d_inner)  9 ! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng in loops and  10 ! allocation of arrays. Calls of exchange_horiz are modified.  11 ! Call ws_init to initialize arrays needed for statistical evaluation and  12 ! optimization when ws-scheme is used.  13 !  14 !  15 ! Initial volume flow is now calculated by using the variable hom_sum.  16 ! Therefore the correction of initial volume flow for non-flat topography  17 ! removed (removed u_nzb_p1_for_vfc and v_nzb_p1_for_vfc)  18 !  19 ! Changed surface boundary conditions for u and v in case of ibc_uv_b == 0 from  20 ! mirror bc to dirichlet boundary conditions (u=v=0), so that k=nzb is  21 ! representative for the height z0  22 !  23 ! Bugfix: type conversion of '1' to 64bit for the MAX function (ngp_3d_inner) 10 24 ! 11 25 ! Former revisions: … …  101 115 !------------------------------------------------------------------------------! 102 116  117 USE advec_ws 103 118 USE arrays_3d 104 119 USE averaging … …  147 162 ngp_2dh_s_inner(nzb:nzt+1,0:statistic_regions), & 148 163 ngp_2dh_s_inner_l(nzb:nzt+1,0:statistic_regions), & 149  rmask(nys -1:nyn+1,nxl-1:nxr+1,0:statistic_regions), & 164 rmask(nysg:nyng,nxlg:nxrg,0:statistic_regions), & 150 165 sums(nzb:nzt+1,pr_palm+max_pr_user), & 151 166 sums_l(nzb:nzt+1,pr_palm+max_pr_user,0:threads_per_task-1), & … …  154 169 sums_wsts_bc_l(nzb:nzt+1,0:statistic_regions), & 155 170 ts_value(dots_max,0:statistic_regions) ) 156  ALLOCATE( km_damp_x(nxl -1:nxr+1), km_damp_y(nys-1:nyn+1) )157  158  ALLOCATE( rif_1(nys -1:nyn+1,nxl-1:nxr+1), shf_1(nys-1:nyn+1,nxl-1:nxr+1), &159  ts(nys -1:nyn+1,nxl-1:nxr+1), tswst_1(nys-1:nyn+1,nxl-1:nxr+1), &160  us(nys -1:nyn+1,nxl-1:nxr+1), usws_1(nys-1:nyn+1,nxl-1:nxr+1), &161  uswst_1(nys -1:nyn+1,nxl-1:nxr+1), &162  vsws_1(nys -1:nyn+1,nxl-1:nxr+1), &163  vswst_1(nys -1:nyn+1,nxl-1:nxr+1), z0(nys-1:nyn+1,nxl-1:nxr+1) ) 171 ALLOCATE( km_damp_x(nxlg:nxrg), km_damp_y(nysg:nyng) )  172  173 ALLOCATE( rif_1(nysg:nyng,nxlg:nxrg), shf_1(nysg:nyng,nxlg:nxrg), &  174 ts(nysg:nyng,nxlg:nxrg), tswst_1(nysg:nyng,nxlg:nxrg), &  175 us(nysg:nyng,nxlg:nxrg), usws_1(nysg:nyng,nxlg:nxrg), &  176 uswst_1(nysg:nyng,nxlg:nxrg), &  177 vsws_1(nysg:nyng,nxlg:nxrg), &  178 vswst_1(nysg:nyng,nxlg:nxrg), z0(nysg:nyng,nxlg:nxrg) ) 164 179 165 180 IF ( timestep_scheme(1:5) /= 'runge' ) THEN 166 181 ! 167 182 !-- Leapfrog scheme needs two timelevels of diffusion quantities 168  ALLOCATE( rif_2(nys -1:nyn+1,nxl-1:nxr+1), &169  shf_2(nys -1:nyn+1,nxl-1:nxr+1), &170  tswst_2(nys -1:nyn+1,nxl-1:nxr+1), &171  usws_2(nys -1:nyn+1,nxl-1:nxr+1), &172  uswst_2(nys -1:nyn+1,nxl-1:nxr+1), &173  vswst_2(nys -1:nyn+1,nxl-1:nxr+1), &174  vsws_2(nys -1:nyn+1,nxl-1:nxr+1) ) 183 ALLOCATE( rif_2(nysg:nyng,nxlg:nxrg), &  184 shf_2(nysg:nyng,nxlg:nxrg), &  185 tswst_2(nysg:nyng,nxlg:nxrg), &  186 usws_2(nysg:nyng,nxlg:nxrg), &  187 uswst_2(nysg:nyng,nxlg:nxrg), &  188 vswst_2(nysg:nyng,nxlg:nxrg), &  189 vsws_2(nysg:nyng,nxlg:nxrg) ) 175 190 ENDIF 176 191 177 192 ALLOCATE( d(nzb+1:nzta,nys:nyna,nxl:nxra), & 178  e_1(nzb:nzt+1,nys -1:nyn+1,nxl-1:nxr+1), &179  e_2(nzb:nzt+1,nys -1:nyn+1,nxl-1:nxr+1), &180  e_3(nzb:nzt+1,nys -1:nyn+1,nxl-1:nxr+1), &181  kh_1(nzb:nzt+1,nys -1:nyn+1,nxl-1:nxr+1), &182  km_1(nzb:nzt+1,nys -1:nyn+1,nxl-1:nxr+1), &183  p(nzb:nzt+1,nys -1:nyn+1,nxl-1:nxr+1), &184  pt_1(nzb:nzt+1,nys -1:nyn+1,nxl-1:nxr+1), &185  pt_2(nzb:nzt+1,nys -1:nyn+1,nxl-1:nxr+1), &186  pt_3(nzb:nzt+1,nys -1:nyn+1,nxl-1:nxr+1), &187  tend(nzb:nzt+1,nys -1:nyn+1,nxl-1:nxr+1), &188  u_1(nzb:nzt+1,nys -1:nyn+1,nxl-1:nxr+1), &189  u_2(nzb:nzt+1,nys -1:nyn+1,nxl-1:nxr+1), &190  u_3(nzb:nzt+1,nys -1:nyn+1,nxl-1:nxr+1), &191  v_1(nzb:nzt+1,nys -1:nyn+1,nxl-1:nxr+1), &192  v_2(nzb:nzt+1,nys -1:nyn+1,nxl-1:nxr+1), &193  v_3(nzb:nzt+1,nys -1:nyn+1,nxl-1:nxr+1), &194  w_1(nzb:nzt+1,nys -1:nyn+1,nxl-1:nxr+1), &195  w_2(nzb:nzt+1,nys -1:nyn+1,nxl-1:nxr+1), &196  w_3(nzb:nzt+1,nys -1:nyn+1,nxl-1:nxr+1) ) 193 e_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &  194 e_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &  195 e_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &  196 kh_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &  197 km_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &  198 p(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &  199 pt_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &  200 pt_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &  201 pt_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &  202 tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &  203 u_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &  204 u_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &  205 u_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &  206 v_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &  207 v_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &  208 v_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &  209 w_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &  210 w_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &  211 w_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 197 212 198 213 IF ( timestep_scheme(1:5) /= 'runge' ) THEN 199  ALLOCATE( kh_2(nzb:nzt+1,nys -1:nyn+1,nxl-1:nxr+1), &200  km_2(nzb:nzt+1,nys -1:nyn+1,nxl-1:nxr+1) ) 214 ALLOCATE( kh_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &  215 km_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 201 216 ENDIF 202 217 … …  204 219 ! 205 220 !-- 2D-humidity/scalar arrays 206  ALLOCATE ( qs(nys -1:nyn+1,nxl-1:nxr+1), &207  qsws_1(nys -1:nyn+1,nxl-1:nxr+1), &208  qswst_1(nys -1:nyn+1,nxl-1:nxr+1) ) 221 ALLOCATE ( qs(nysg:nyng,nxlg:nxrg), &  222 qsws_1(nysg:nyng,nxlg:nxrg), &  223 qswst_1(nysg:nyng,nxlg:nxrg) ) 209 224 210 225 IF ( timestep_scheme(1:5) /= 'runge' ) THEN 211  ALLOCATE( qsws_2(nys -1:nyn+1,nxl-1:nxr+1), &212  qswst_2(nys -1:nyn+1,nxl-1:nxr+1) ) 226 ALLOCATE( qsws_2(nysg:nyng,nxlg:nxrg), &  227 qswst_2(nysg:nyng,nxlg:nxrg) ) 213 228 ENDIF 214 229 ! 215 230 !-- 3D-humidity/scalar arrays 216  ALLOCATE( q_1(nzb:nzt+1,nys -1:nyn+1,nxl-1:nxr+1), &217  q_2(nzb:nzt+1,nys -1:nyn+1,nxl-1:nxr+1), &218  q_3(nzb:nzt+1,nys -1:nyn+1,nxl-1:nxr+1) ) 231 ALLOCATE( q_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &  232 q_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &  233 q_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 219 234 220 235 ! 221 236 !-- 3D-arrays needed for humidity only 222 237 IF ( humidity ) THEN 223  ALLOCATE( vpt_1(nzb:nzt+1,nys -1:nyn+1,nxl-1:nxr+1) ) 238 ALLOCATE( vpt_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 224 239 225 240 IF ( timestep_scheme(1:5) /= 'runge' ) THEN 226  ALLOCATE( vpt_2(nzb:nzt+1,nys -1:nyn+1,nxl-1:nxr+1) ) 241 ALLOCATE( vpt_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 227 242 ENDIF 228 243 … …  230 245 ! 231 246 !-- Liquid water content 232  ALLOCATE ( ql_1(nzb:nzt+1,nys -1:nyn+1,nxl-1:nxr+1) ) 247 ALLOCATE ( ql_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 233 248 ! 234 249 !-- Precipitation amount and rate (only needed if output is switched) 235  ALLOCATE( precipitation_amount(nys -1:nyn+1,nxl-1:nxr+1), &236  precipitation_rate(nys -1:nyn+1,nxl-1:nxr+1) ) 250 ALLOCATE( precipitation_amount(nysg:nyng,nxlg:nxrg), &  251 precipitation_rate(nysg:nyng,nxlg:nxrg) ) 237 252 ENDIF 238 253 … …  241 256 !-- Liquid water content, change in liquid water content, 242 257 !-- real volume of particles (with weighting), volume of particles 243  ALLOCATE ( ql_1(nzb:nzt+1,nys -1:nyn+1,nxl-1:nxr+1), &244  ql_2(nzb:nzt+1,nys -1:nyn+1,nxl-1:nxr+1), &245  ql_v(nzb:nzt+1,nys -1:nyn+1,nxl-1:nxr+1), &246  ql_vp(nzb:nzt+1,nys -1:nyn+1,nxl-1:nxr+1) ) 258 ALLOCATE ( ql_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &  259 ql_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &  260 ql_v(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &  261 ql_vp(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 247 262 ENDIF 248 263 … …  252 267 253 268 IF ( ocean ) THEN 254  ALLOCATE( saswsb_1(nys -1:nyn+1,nxl-1:nxr+1), &255  saswst_1(nys -1:nyn+1,nxl-1:nxr+1) )256  ALLOCATE( prho_1(nzb:nzt+1,nys -1:nyn+1,nxl-1:nxr+1), &257  rho_1(nzb:nzt+1,nys -1:nyn+1,nxl-1:nxr+1), &258  sa_1(nzb:nzt+1,nys -1:nyn+1,nxl-1:nxr+1), &259  sa_2(nzb:nzt+1,nys -1:nyn+1,nxl-1:nxr+1), &260  sa_3(nzb:nzt+1,nys -1:nyn+1,nxl-1:nxr+1) ) 269 ALLOCATE( saswsb_1(nysg:nyng,nxlg:nxrg), &  270 saswst_1(nysg:nyng,nxlg:nxrg) )  271 ALLOCATE( prho_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &  272 rho_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &  273 sa_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &  274 sa_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &  275 sa_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 261 276 prho => prho_1 262 277 rho => rho_1 ! routines calc_mean_profile and diffusion_e require 263 278 ! density to be apointer 264 279 IF ( humidity_remote ) THEN 265  ALLOCATE( qswst_remote(nys -1:nyn+1,nxl-1:nxr+1)) 280 ALLOCATE( qswst_remote(nysg:nyng,nxlg:nxrg)) 266 281 qswst_remote = 0.0 267 282 ENDIF … …  272 287 !-- particle velocities 273 288 IF ( use_sgs_for_particles ) THEN 274  ALLOCATE ( diss(nzb:nzt+1,nys -1:nyn+1,nxl-1:nxr+1) ) 289 ALLOCATE ( diss(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 275 290 ELSE 276 291 ALLOCATE ( diss(2,2,2) ) ! required because diss is used as a … …  288 303 !-- 3D-arrays for the leaf area density and the canopy drag coefficient 289 304 IF ( plant_canopy ) THEN 290  ALLOCATE ( lad_s(nzb:nzt+1,nys -1:nyn+1,nxl-1:nxr+1), &291  lad_u(nzb:nzt+1,nys -1:nyn+1,nxl-1:nxr+1), &292  lad_v(nzb:nzt+1,nys -1:nyn+1,nxl-1:nxr+1), &293  lad_w(nzb:nzt+1,nys -1:nyn+1,nxl-1:nxr+1), &294  cdc(nzb:nzt+1,nys -1:nyn+1,nxl-1:nxr+1) ) 305 ALLOCATE ( lad_s(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &  306 lad_u(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &  307 lad_v(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &  308 lad_w(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &  309 cdc(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 295 310 296 311 IF ( passive_scalar ) THEN 297  ALLOCATE ( sls(nzb:nzt+1,nys -1:nyn+1,nxl-1:nxr+1), &298  sec(nzb:nzt+1,nys -1:nyn+1,nxl-1:nxr+1) ) 312 ALLOCATE ( sls(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &  313 sec(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 299 314 ENDIF 300 315 301 316 IF ( cthf /= 0.0 ) THEN 302  ALLOCATE ( lai(nzb:nzt+1,nys -1:nyn+1,nxl-1:nxr+1), &303  canopy_heat_flux(nzb:nzt+1,nys -1:nyn+1,nxl-1:nxr+1) ) 317 ALLOCATE ( lai(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &  318 canopy_heat_flux(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 304 319 ENDIF 305 320 … …  309 324 !-- 4D-array for storing the Rif-values at vertical walls 310 325 IF ( topography /= 'flat' ) THEN 311  ALLOCATE( rif_wall(nzb:nzt+1,nys -1:nyn+1,nxl-1:nxr+1,1:4) ) 326 ALLOCATE( rif_wall(nzb:nzt+1,nysg:nyng,nxlg:nxrg,1:4) ) 312 327 rif_wall = 0.0 313  ENDIF314  315  !316  !-- Velocities at nzb+1 needed for volume flow control317  IF ( conserve_volume_flow ) THEN318  ALLOCATE( u_nzb_p1_for_vfc(nys:nyn), v_nzb_p1_for_vfc(nxl:nxr) )319  u_nzb_p1_for_vfc = 0.0320  v_nzb_p1_for_vfc = 0.0321 328 ENDIF 322 329 … …  325 332 !-- are needed for radiation boundary conditions 326 333 IF ( outflow_l ) THEN 327  ALLOCATE( u_m_l(nzb:nzt+1,nys -1:nyn+1,1:2), &328  v_m_l(nzb:nzt+1,nys -1:nyn+1,0:1), &329  w_m_l(nzb:nzt+1,nys -1:nyn+1,0:1) ) 334 ALLOCATE( u_m_l(nzb:nzt+1,nysg:nyng,1:2), &  335 v_m_l(nzb:nzt+1,nysg:nyng,0:1), &  336 w_m_l(nzb:nzt+1,nysg:nyng,0:1) ) 330 337 ENDIF 331 338 IF ( outflow_r ) THEN 332  ALLOCATE( u_m_r(nzb:nzt+1,nys -1:nyn+1,nx-1:nx), &333  v_m_r(nzb:nzt+1,nys -1:nyn+1,nx-1:nx), &334  w_m_r(nzb:nzt+1,nys -1:nyn+1,nx-1:nx) ) 339 ALLOCATE( u_m_r(nzb:nzt+1,nysg:nyng,nx-1:nx), &  340 v_m_r(nzb:nzt+1,nysg:nyng,nx-1:nx), &  341 w_m_r(nzb:nzt+1,nysg:nyng,nx-1:nx) ) 335 342 ENDIF 336 343 IF ( outflow_l .OR. outflow_r ) THEN 337  ALLOCATE( c_u(nzb:nzt+1,nys -1:nyn+1), c_v(nzb:nzt+1,nys-1:nyn+1), &338  c_w(nzb:nzt+1,nys -1:nyn+1) ) 344 ALLOCATE( c_u(nzb:nzt+1,nysg:nyng), c_v(nzb:nzt+1,nysg:nyng), &  345 c_w(nzb:nzt+1,nysg:nyng) ) 339 346 ENDIF 340 347 IF ( outflow_s ) THEN 341  ALLOCATE( u_m_s(nzb:nzt+1,0:1,nxl -1:nxr+1), &342  v_m_s(nzb:nzt+1,1:2,nxl -1:nxr+1), &343  w_m_s(nzb:nzt+1,0:1,nxl -1:nxr+1) ) 348 ALLOCATE( u_m_s(nzb:nzt+1,0:1,nxlg:nxrg), &  349 v_m_s(nzb:nzt+1,1:2,nxlg:nxrg), &  350 w_m_s(nzb:nzt+1,0:1,nxlg:nxrg) ) 344 351 ENDIF 345 352 IF ( outflow_n ) THEN 346  ALLOCATE( u_m_n(nzb:nzt+1,ny-1:ny,nxl -1:nxr+1), &347  v_m_n(nzb:nzt+1,ny-1:ny,nxl -1:nxr+1), &348  w_m_n(nzb:nzt+1,ny-1:ny,nxl -1:nxr+1) ) 353 ALLOCATE( u_m_n(nzb:nzt+1,ny-1:ny,nxlg:nxrg), &  354 v_m_n(nzb:nzt+1,ny-1:ny,nxlg:nxrg), &  355 w_m_n(nzb:nzt+1,ny-1:ny,nxlg:nxrg) ) 349 356 ENDIF 350 357 IF ( outflow_s .OR. outflow_n ) THEN 351  ALLOCATE( c_u(nzb:nzt+1,nxl -1:nxr+1), c_v(nzb:nzt+1,nxl-1:nxr+1), &352  c_w(nzb:nzt+1,nxl -1:nxr+1) ) 358 ALLOCATE( c_u(nzb:nzt+1,nxlg:nxrg), c_v(nzb:nzt+1,nxlg:nxrg), &  359 c_w(nzb:nzt+1,nxlg:nxrg) ) 353 360 ENDIF 354 361 … …  435 442 ! 436 443 !-- Transfer initial profiles to the arrays of the 3D model 437  DO i = nxl -1, nxr+1438  DO j = nys -1, nyn+1 444 DO i = nxlg, nxrg  445 DO j = nysg, nyng 439 446 e(:,j,i) = e1d 440 447 kh(:,j,i) = kh1d … …  447 454 448 455 IF ( humidity .OR. passive_scalar ) THEN 449  DO i = nxl -1, nxr+1450  DO j = nys -1, nyn+1 456 DO i = nxlg, nxrg  457 DO j = nysg, nyng 451 458 q(:,j,i) = q_init 452 459 ENDDO … …  455 462 456 463 IF ( .NOT. constant_diffusion ) THEN 457  DO i = nxl -1, nxr+1458  DO j = nys -1, nyn+1 464 DO i = nxlg, nxrg  465 DO j = nysg, nyng 459 466 e(:,j,i) = e1d 460 467 ENDDO … …  505 512 ENDDO 506 513 ENDDO 507  IF ( conserve_volume_flow ) THEN 508  IF ( nxr == nx ) THEN 509  DO j = nys, nyn 510  DO k = nzb + 1, nzb_u_inner(j,nx) 511  u_nzb_p1_for_vfc(j) = u_nzb_p1_for_vfc(j) + & 512  u1d(k) * dzu(k) 513  ENDDO 514  ENDDO 515  ENDIF 516  IF ( nyn == ny ) THEN 517  DO i = nxl, nxr 518  DO k = nzb + 1, nzb_v_inner(ny,i) 519  v_nzb_p1_for_vfc(i) = v_nzb_p1_for_vfc(i) + & 520  v1d(k) * dzu(k) 521  ENDDO 522  ENDDO 523  ENDIF 524  ENDIF  514 525 515 ! 526 516 !-- WARNING: The extra boundary conditions set after running the … …  530 520 !-- --------- advection scheme: keep u and v zero one layer below 531 521 !-- the topography. 532  IF ( ibc_uv_b == 0 ) THEN 533  ! 534  !-- Satisfying the Dirichlet condition with an extra layer below 535  !-- the surface where the u and v component change their sign. 536  DO i = nxl-1, nxr+1 537  DO j = nys-1, nyn+1 538  IF ( nzb_u_inner(j,i) == 0 ) u(0,j,i) = -u(1,j,i) 539  IF ( nzb_v_inner(j,i) == 0 ) v(0,j,i) = -v(1,j,i) 540  ENDDO 541  ENDDO 542  543  ELSE  522 !  523 !-- Following was removed, because mirror boundary condition are  524 !-- replaced by dirichlet boundary conditions  525 !  526 ! IF ( ibc_uv_b == 0 ) THEN  527 !!  528 !!-- Satisfying the Dirichlet condition with an extra layer below  529 !!-- the surface where the u and v component change their sign.  530 ! DO i = nxl-1, nxr+1  531 ! DO j = nys-1, nyn+1  532 ! IF ( nzb_u_inner(j,i) == 0 ) u(0,j,i) = -u(1,j,i)  533 ! IF ( nzb_v_inner(j,i) == 0 ) v(0,j,i) = -v(1,j,i)  534 ! ENDDO  535 ! ENDDO  536 !  537 ! ELSE  538 IF ( ibc_uv_b == 1 ) THEN 544 539 ! 545 540 !-- Neumann condition … …  560 555 !-- Use constructed initial profiles (velocity constant with height, 561 556 !-- temperature profile with constant gradient) 562  DO i = nxl -1, nxr+1563  DO j = nys -1, nyn+1 557 DO i = nxlg, nxrg  558 DO j = nysg, nyng 564 559 pt(:,j,i) = pt_init 565 560 u(:,j,i) = u_init … …  574 569 !-- in the limiting formula!). The original values are stored to be later 575 570 !-- used for volume flow control. 576  DO i = nxl -1, nxr+1577  DO j = nys -1, nyn+1 571 DO i = nxlg, nxrg  572 DO j = nysg, nyng 578 573 u(nzb:nzb_u_inner(j,i)+1,j,i) = 0.0 579 574 v(nzb:nzb_v_inner(j,i)+1,j,i) = 0.0 580 575 ENDDO 581 576 ENDDO 582  IF ( conserve_volume_flow ) THEN583  IF ( nxr == nx ) THEN584  DO j = nys, nyn585  DO k = nzb + 1, nzb_u_inner(j,nx) + 1586  u_nzb_p1_for_vfc(j) = u_nzb_p1_for_vfc(j) + &587  u_init(k) * dzu(k)588  ENDDO589  ENDDO590  ENDIF591  IF ( nyn == ny ) THEN592  DO i = nxl, nxr593  DO k = nzb + 1, nzb_v_inner(ny,i) + 1594  v_nzb_p1_for_vfc(i) = v_nzb_p1_for_vfc(i) + &595  v_init(k) * dzu(k)596  ENDDO597  ENDDO598  ENDIF599  ENDIF600 577 601 578 IF ( humidity .OR. passive_scalar ) THEN 602  DO i = nxl -1, nxr+1603  DO j = nys -1, nyn+1 579 DO i = nxlg, nxrg  580 DO j = nysg, nyng 604 581 q(:,j,i) = q_init 605 582 ENDDO … …  608 585 609 586 IF ( ocean ) THEN 610  DO i = nxl -1, nxr+1611  DO j = nys -1, nyn+1 587 DO i = nxlg, nxrg  588 DO j = nysg, nyng 612 589 sa(:,j,i) = sa_init 613 590 ENDDO … …  660 637 661 638 ENDIF  639 !  640 !-- Bottom boundary  641 IF ( ibc_uv_b == 0 .OR. ibc_uv_b == 2 ) THEN  642 u(nzb,:,:) = 0.0  643 v(nzb,:,:) = 0.0  644 ENDIF 662 645 663 646 ! … …  683 666 hom(:,1,5,:) = SPREAD( u(:,nys,nxl), 2, statistic_regions+1 ) 684 667 hom(:,1,6,:) = SPREAD( v(:,nys,nxl), 2, statistic_regions+1 ) 685  IF ( ibc_uv_b == 0 ) THEN 686  hom(nzb,1,5,:) = -hom(nzb+1,1,5,:) ! due to satisfying the Dirichlet 687  hom(nzb,1,6,:) = -hom(nzb+1,1,6,:) ! condition with an extra layer 688  ! below the surface where the u and v component change their sign  668 IF ( ibc_uv_b == 0 .OR. ibc_uv_b == 2) THEN  669 hom(nzb,1,5,:) = 0.0  670 hom(nzb,1,6,:) = 0.0 689 671 ENDIF 690 672 hom(:,1,7,:) = SPREAD( pt(:,nys,nxl), 2, statistic_regions+1 ) … …  733 715 !-- Over topography surface_heatflux is replaced by wall_heatflux(0) 734 716 IF ( TRIM( topography ) /= 'flat' ) THEN 735  DO i = nxl -1, nxr+1736  DO j = nys -1, nyn+1 717 DO i = nxlg, nxrg  718 DO j = nysg, nyng 737 719 IF ( nzb_s_inner(j,i) /= 0 ) THEN 738 720 shf(j,i) = wall_heatflux(0) … …  755 737 IF ( TRIM( topography ) /= 'flat' ) THEN 756 738 wall_qflux = wall_humidityflux 757  DO i = nxl -1, nxr+1758  DO j = nys -1, nyn+1 739 DO i = nxlg, nxrg  740 DO j = nysg, nyng 759 741 IF ( nzb_s_inner(j,i) /= 0 ) THEN 760 742 qsws(j,i) = wall_qflux(0) … …  827 809 ENDIF 828 810 829  !830  !-- Calculate the initial volume flow at the right and north boundary831  IF ( conserve_volume_flow ) THEN832  833  volume_flow_initial_l = 0.0834  volume_flow_area_l = 0.0835  836  IF ( nxr == nx ) THEN837  DO j = nys, nyn838  DO k = nzb_2d(j,nx) + 1, nzt839  volume_flow_initial_l(1) = volume_flow_initial_l(1) + &840  u(k,j,nx) * dzu(k)841  volume_flow_area_l(1) = volume_flow_area_l(1) + dzu(k)842  ENDDO843  !844  !-- Correction if velocity at nzb+1 has been set zero further above845  volume_flow_initial_l(1) = volume_flow_initial_l(1) + &846  u_nzb_p1_for_vfc(j)847  ENDDO848  ENDIF849  850  IF ( nyn == ny ) THEN851  DO i = nxl, nxr852  DO k = nzb_2d(ny,i) + 1, nzt853  volume_flow_initial_l(2) = volume_flow_initial_l(2) + &854  v(k,ny,i) * dzu(k)855  volume_flow_area_l(2) = volume_flow_area_l(2) + dzu(k)856  ENDDO857  !858  !-- Correction if velocity at nzb+1 has been set zero further above859  volume_flow_initial_l(2) = volume_flow_initial_l(2) + &860  v_nzb_p1_for_vfc(i)861  ENDDO862  ENDIF863  864  #if defined( __parallel )865  IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )866  CALL MPI_ALLREDUCE( volume_flow_initial_l(1), volume_flow_initial(1),&867  2, MPI_REAL, MPI_SUM, comm2d, ierr )868  IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )869  CALL MPI_ALLREDUCE( volume_flow_area_l(1), volume_flow_area(1), &870  2, MPI_REAL, MPI_SUM, comm2d, ierr )871  #else872  volume_flow_initial = volume_flow_initial_l873  volume_flow_area = volume_flow_area_l874  #endif875  !876  !-- In case of 'bulk_velocity' mode, volume_flow_initial is overridden877  !-- and calculated from u|v_bulk instead.878  IF ( TRIM( conserve_volume_flow_mode ) == 'bulk_velocity' ) THEN879  volume_flow_initial(1) = u_bulk * volume_flow_area(1)880  volume_flow_initial(2) = v_bulk * volume_flow_area(2)881  ENDIF882  883  ENDIF884 811 885 812 ! … …  968 895 sa_p = sa 969 896 ENDIF 970   897 971 898 972 899 ELSEIF ( TRIM( initializing_actions ) == 'read_restart_data' .OR. & 973   900 TRIM( initializing_actions ) == 'cyclic_fill' ) & 974 901 THEN 975 902 ! … …  978 905 IF ( TRIM( initializing_actions ) == 'cyclic_fill' ) THEN 979 906 980  WRITE (9,*) 'before read_parts_of_var_list'981  CALL local_flush( 9 )982 907 CALL read_parts_of_var_list 983  WRITE (9,*) 'after read_parts_of_var_list'984  CALL local_flush( 9 )985 908 CALL close_file( 13 ) 986 909 … …  1002 925 !-- conditions are used) 1003 926 IF ( inflow_l ) THEN 1004  DO j = nys -1, nyn+1 927 DO j = nysg, nyng 1005 928 DO k = nzb, nzt+1 1006  u(k,j, -1) = mean_inflow_profiles(k,1)1007  v(k,j, -1) = mean_inflow_profiles(k,2)1008  w(k,j, -1) = 0.01009  pt(k,j, -1) = mean_inflow_profiles(k,4)1010  e(k,j, -1) = mean_inflow_profiles(k,5) 929 u(k,j,nxlg:-1) = mean_inflow_profiles(k,1)  930 v(k,j,nxlg:-1) = mean_inflow_profiles(k,2)  931 w(k,j,nxlg:-1) = 0.0  932 pt(k,j,nxlg:-1) = mean_inflow_profiles(k,4)  933 e(k,j,nxlg:-1) = mean_inflow_profiles(k,5) 1011 934 ENDDO 1012 935 ENDDO … …  1064 987 ! 1065 988 !-- Read binary data from restart file 1066  WRITE (9,*) 'before read_3d_binary' 1067  CALL local_flush( 9 )  989 1068 990 CALL read_3d_binary 1069  WRITE (9,*) 'after read_3d_binary'1070  CALL local_flush( 9 )1071 991 1072 992 ! … …  1074 994 IF ( TRIM( initializing_actions ) == 'cyclic_fill' .AND. & 1075 995 topography /= 'flat' ) THEN 1076  !1077  !-- Correction of initial volume flow1078  IF ( conserve_volume_flow ) THEN1079  IF ( nxr == nx ) THEN1080  DO j = nys, nyn1081  DO k = nzb + 1, nzb_u_inner(j,nx)1082  u_nzb_p1_for_vfc(j) = u_nzb_p1_for_vfc(j) + &1083  u(k,j,nx) * dzu(k)1084  ENDDO1085  ENDDO1086  ENDIF1087  IF ( nyn == ny ) THEN1088  DO i = nxl, nxr1089  DO k = nzb + 1, nzb_v_inner(ny,i)1090  v_nzb_p1_for_vfc(i) = v_nzb_p1_for_vfc(i) + &1091  v(k,ny,i) * dzu(k)1092  ENDDO1093  ENDDO1094  ENDIF1095  ENDIF1096  1097 996 ! 1098 997 !-- Inside buildings set velocities and TKE back to zero. … …  1100 999 !-- maybe revise later. 1101 1000 IF ( timestep_scheme(1:5) == 'runge' ) THEN 1102  DO i = nxl -1, nxr+11103  DO j = nys -1, nyn+1 1001 DO i = nxlg, nxrg  1002 DO j = nysg, nyng 1104 1003 u (nzb:nzb_u_inner(j,i),j,i) = 0.0 1105 1004 v (nzb:nzb_v_inner(j,i),j,i) = 0.0 … …  1118 1017 ENDDO 1119 1018 ELSE 1120  DO i = nxl -1, nxr+11121  DO j = nys -1, nyn+1 1019 DO i = nxlg, nxrg  1020 DO j = nysg, nyng 1122 1021 u (nzb:nzb_u_inner(j,i),j,i) = 0.0 1123 1022 v (nzb:nzb_v_inner(j,i),j,i) = 0.0 … …  1139 1038 1140 1039 ! 1141  !-- Calculate the initial volume flow at the right and north boundary1142  IF ( conserve_volume_flow .AND. &1143  TRIM( initializing_actions ) == 'cyclic_fill' ) THEN1144  1145  volume_flow_initial_l = 0.01146  volume_flow_area_l = 0.01147  1148  IF ( nxr == nx ) THEN1149  DO j = nys, nyn1150  DO k = nzb_2d(j,nx) + 1, nzt1151  volume_flow_initial_l(1) = volume_flow_initial_l(1) + &1152  u(k,j,nx) * dzu(k)1153  volume_flow_area_l(1) = volume_flow_area_l(1) + dzu(k)1154  ENDDO1155  !1156  !-- Correction if velocity inside buildings has been set to zero1157  !-- further above1158  volume_flow_initial_l(1) = volume_flow_initial_l(1) + &1159  u_nzb_p1_for_vfc(j)1160  ENDDO1161  ENDIF1162  1163  IF ( nyn == ny ) THEN1164  DO i = nxl, nxr1165  DO k = nzb_2d(ny,i) + 1, nzt1166  volume_flow_initial_l(2) = volume_flow_initial_l(2) + &1167  v(k,ny,i) * dzu(k)1168  volume_flow_area_l(2) = volume_flow_area_l(2) + dzu(k)1169  ENDDO1170  !1171  !-- Correction if velocity inside buildings has been set to zero1172  !-- further above1173  volume_flow_initial_l(2) = volume_flow_initial_l(2) + &1174  v_nzb_p1_for_vfc(i)1175  ENDDO1176  ENDIF1177  1178  #if defined( __parallel )1179  IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )1180  CALL MPI_ALLREDUCE( volume_flow_initial_l(1), volume_flow_initial(1),&1181  2, MPI_REAL, MPI_SUM, comm2d, ierr )1182  IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )1183  CALL MPI_ALLREDUCE( volume_flow_area_l(1), volume_flow_area(1), &1184  2, MPI_REAL, MPI_SUM, comm2d, ierr )1185  #else1186  volume_flow_initial = volume_flow_initial_l1187  volume_flow_area = volume_flow_area_l1188  #endif1189  ENDIF1190  1191  1192  !1193 1040 !-- Calculate initial temperature field and other constants used in case 1194 1041 !-- of a sloping surface … …  1243 1090 w_m_n(:,:,:) = w(:,ny-1:ny,:) 1244 1091 ENDIF 1245  1246  ENDIF 1247   1092  1093 ENDIF  1094 !  1095 !-- Calculate the initial volume flow at the right and north boundary  1096 IF ( conserve_volume_flow ) THEN  1097  1098 volume_flow_initial_l = 0.0  1099 volume_flow_area_l = 0.0  1100  1101 IF ( TRIM( initializing_actions ) == 'cyclic_fill' ) THEN  1102  1103 IF ( nxr == nx ) THEN  1104 DO j = nys, nyn  1105 DO k = nzb_2d(j,nx) + 1, nzt  1106 volume_flow_initial_l(1) = volume_flow_initial_l(1) + &  1107 hom_sum(k,1,0) * dzw(k)  1108 volume_flow_area_l(1) = volume_flow_area_l(1) + dzw(k)  1109 ENDDO  1110 ENDDO  1111 ENDIF  1112  1113 IF ( nyn == ny ) THEN  1114 DO i = nxl, nxr  1115 DO k = nzb_2d(ny,i) + 1, nzt  1116 volume_flow_initial_l(2) = volume_flow_initial_l(2) + &  1117 hom_sum(k,2,0) * dzw(k)  1118 volume_flow_area_l(2) = volume_flow_area_l(2) + dzw(k)  1119 ENDDO  1120 ENDDO  1121 ENDIF  1122  1123 ELSEIF ( TRIM( initializing_actions ) /= 'read_restart_data' ) THEN  1124  1125 IF ( nxr == nx ) THEN  1126 DO j = nys, nyn  1127 DO k = nzb_2d(j,nx) + 1, nzt  1128 volume_flow_initial_l(1) = volume_flow_initial_l(1) + &  1129 u(k,j,nx) * dzw(k)  1130 volume_flow_area_l(1) = volume_flow_area_l(1) + dzw(k)  1131 ENDDO  1132 ENDDO  1133 ENDIF  1134  1135 IF ( nyn == ny ) THEN  1136 DO i = nxl, nxr  1137 DO k = nzb_2d(ny,i) + 1, nzt  1138 volume_flow_initial_l(2) = volume_flow_initial_l(2) + &  1139 v(k,ny,i) * dzw(k)  1140 volume_flow_area_l(2) = volume_flow_area_l(2) + dzw(k)  1141 ENDDO  1142 ENDDO  1143 ENDIF  1144  1145 ENDIF  1146  1147 #if defined( __parallel )  1148 CALL MPI_ALLREDUCE( volume_flow_initial_l(1), volume_flow_initial(1),&  1149 2, MPI_REAL, MPI_SUM, comm2d, ierr )  1150 CALL MPI_ALLREDUCE( volume_flow_area_l(1), volume_flow_area(1), &  1151 2, MPI_REAL, MPI_SUM, comm2d, ierr )  1152  1153 CALL MPI_ALLREDUCE( volume_flow_initial_l(2), volume_flow_initial(2),&  1154 2, MPI_REAL, MPI_SUM, comm2d, ierr )  1155 CALL MPI_ALLREDUCE( volume_flow_area_l(2), volume_flow_area(2), &  1156 2, MPI_REAL, MPI_SUM, comm2d, ierr )  1157  1158 #else  1159 volume_flow_initial = volume_flow_initial_l  1160 volume_flow_area = volume_flow_area_l  1161 #endif  1162  1163 !  1164 !-- In case of 'bulk_velocity' mode, volume_flow_initial is overridden  1165 !-- and calculated from u|v_bulk instead.  1166 IF ( TRIM( conserve_volume_flow_mode ) == 'bulk_velocity' ) THEN  1167 volume_flow_initial(1) = u_bulk * volume_flow_area(1)  1168 volume_flow_initial(2) = v_bulk * volume_flow_area(2)  1169 ENDIF  1170  1171 ENDIF 1248 1172 ! 1249 1173 !-- Initialization of the leaf area density … …  1254 1178 CASE( 'block' ) 1255 1179 1256  DO i = nxl -1, nxr+11257  DO j = nys -1, nyn+1 1180 DO i = nxlg, nxrg  1181 DO j = nysg, nyng 1258 1182 lad_s(:,j,i) = lad(:) 1259 1183 cdc(:,j,i) = drag_coefficient … …  1277 1201 END SELECT 1278 1202 1279  CALL exchange_horiz( lad_s )1280  CALL exchange_horiz( cdc ) 1203 CALL exchange_horiz( lad_s, nbgp )  1204 CALL exchange_horiz( cdc, nbgp ) 1281 1205 1282 1206 IF ( passive_scalar ) THEN 1283  CALL exchange_horiz( sls )1284  CALL exchange_horiz( sec ) 1207 CALL exchange_horiz( sls, nbgp )  1208 CALL exchange_horiz( sec, nbgp ) 1285 1209 ENDIF 1286 1210 … …  1311 1235 lad_w(nzt+1,:,:) = lad_w(nzt,:,:) 1312 1236 1313  CALL exchange_horiz( lad_u )1314  CALL exchange_horiz( lad_v )1315  CALL exchange_horiz( lad_w ) 1237 CALL exchange_horiz( lad_u, nbgp )  1238 CALL exchange_horiz( lad_v, nbgp )  1239 CALL exchange_horiz( lad_w, nbgp ) 1316 1240 1317 1241 ! … …  1322 1246 !-- integration of the leaf area density 1323 1247 lai(:,:,:) = 0.0 1324  DO i = nxl -1, nxr+11325  DO j = nys -1, nyn+1 1248 DO i = nxlg, nxrg  1249 DO j = nysg, nyng 1326 1250 DO k = pch_index-1, 0, -1 1327 1251 lai(k,j,i) = lai(k+1,j,i) + & … …  1339 1263 !-- Evaluation of the upward kinematic vertical heat flux within the 1340 1264 !-- canopy 1341  DO i = nxl -1, nxr+11342  DO j = nys -1, nyn+1 1265 DO i = nxlg, nxrg  1266 DO j = nysg, nyng 1343 1267 DO k = 0, pch_index 1344 1268 canopy_heat_flux(k,j,i) = cthf * & … …  1384 1308 !-- Initialize quantities for special advections schemes 1385 1309 CALL init_advec  1310 IF ( momentum_advec == 'ws-scheme' .OR. &  1311 scalar_advec == 'ws-scheme' ) CALL ws_init 1386 1312 1387 1313 ! … …  1439 1365 IF ( bc_lr == 'dirichlet/radiation' ) THEN 1440 1366 1441  DO i = nxl -1, nxr+1 1367 DO i = nxlg, nxrg 1442 1368 IF ( i >= nx - outflow_damping_width ) THEN 1443 1369 km_damp_x(i) = km_damp_max * MIN( 1.0, & … …  1452 1378 ELSEIF ( bc_lr == 'radiation/dirichlet' ) THEN 1453 1379 1454  DO i = nxl -1, nxr+1 1380 DO i = nxlg, nxrg 1455 1381 IF ( i <= outflow_damping_width ) THEN 1456 1382 km_damp_x(i) = km_damp_max * MIN( 1.0, & … …  1467 1393 IF ( bc_ns == 'dirichlet/radiation' ) THEN 1468 1394 1469  DO j = nys -1, nyn+1 1395 DO j = nysg, nyng 1470 1396 IF ( j >= ny - outflow_damping_width ) THEN 1471 1397 km_damp_y(j) = km_damp_max * MIN( 1.0, & … …  1480 1406 ELSEIF ( bc_ns == 'radiation/dirichlet' ) THEN 1481 1407 1482  DO j = nys -1, nyn+1 1408 DO j = nysg, nyng 1483 1409 IF ( j <= outflow_damping_width ) THEN 1484 1410 km_damp_y(j) = km_damp_max * MIN( 1.0, & … …  1594 1520 !-- buoyancy, etc. A zero value will occur for cases where all grid points of 1595 1521 !-- the respective subdomain lie below the surface topography 1596  ngp_2dh_outer = MAX( 1, ngp_2dh_outer(:,:) )  1522 ngp_2dh_outer = MAX( 1, ngp_2dh_outer(:,:) ) 1597 1523 ngp_3d_inner = MAX( INT(1, KIND = SELECTED_INT_KIND( 18 )), & 1598 1524 ngp_3d_inner(:) ) 1599  ngp_2dh_s_inner = MAX( 1, ngp_2dh_s_inner(:,:) )  1525 ngp_2dh_s_inner = MAX( 1, ngp_2dh_s_inner(:,:) ) 1600 1526 1601 1527 DEALLOCATE( ngp_2dh_l, ngp_2dh_outer_l, ngp_3d_inner_l, ngp_3d_inner_tmp ) -
palm/trunk/SOURCE/init_advec.f90
r484 r667 Â 8 8 ! Former revisions: 9 9 ! ----------------- 10 Â ! $Id $Â 10 ! $Id: init_advec.f90 484 2010-02-05 07:36:54Z raasch 11 11 ! RCS Log replace by Id keyword, revision history cleaned up 12 12 ! -
palm/trunk/SOURCE/init_coupling.f90
r484 r667  4 4 ! Current revisions: 5 5 ! ----------------- 6  !  6 !  7 ! determination of target_id's moved to init_pegrid  8 ! 7 9 ! 8 10 ! Former revisions: … …  24 26 USE pegrid 25 27 USE control_parameters  28 USE indices 26 29 27 30 IMPLICIT NONE … …  48 51 !-- ATTENTION: numprocs will be reset according to the new communicators 49 52 #if defined ( __parallel )  53  54 !myid_absolut = myid 50 55 IF ( myid == 0 ) THEN 51 56 READ (*,*,ERR=10,END=10) coupling_mode, bc_data(1), bc_data(2) … …  94 99 95 100 IF ( myid < bc_data(1) ) THEN 96  inter_color = 0 97  numprocs = bc_data(1)  101 inter_color = 0  102 numprocs = bc_data(1)  103 coupling_mode = 'atmosphere_to_ocean' 98 104 ELSE 99  inter_color = 1 100  numprocs = bc_data(2)  105 inter_color = 1  106 numprocs = bc_data(2)  107 coupling_mode = 'ocean_to_atmosphere' 101 108 ENDIF 102  ! 103  !-- Calculate the target PE for coupling and set up the new communicators. 104  !-- Currently only 1:1 topologies are supported. 105  target_id = myid - ISIGN( numprocs, inter_color - 1 ) 106  IF ( inter_color == 0 ) THEN 107  coupling_mode = 'atmosphere_to_ocean' 108  ELSE 109  coupling_mode = 'ocean_to_atmosphere' 110  ENDIF  109 111 110 CALL MPI_COMM_SPLIT( MPI_COMM_WORLD, inter_color, 0, comm_palm, ierr ) 112 111 comm2d = comm_palm -
palm/trunk/SOURCE/init_grid.f90
r559 r667  4 4 ! Current revisions: 5 5 ! ----------------- 6  !  6 ! Definition of new array bounds nxlg, nxrg, nysg, nyng on each PE.  7 ! Furthermore the allocation of arrays and steering of loops is done with these  8 ! parameters. Call of exchange_horiz are modified.  9 ! In case of dirichlet bounday condition at the bottom zu(0)=0.0  10 ! dzu_mg has to be set explicitly for a equally spaced grid near bottom.  11 ! ddzu_pres added to use a equally spaced grid near bottom. 7 12 ! 8 13 ! Former revisions: … …  76 81 77 82 REAL, DIMENSION(:,:,:), ALLOCATABLE :: distance 78   83  84 !  85 ! Computation of the array bounds.  86 nxlg = nxl - nbgp  87 nxrg = nxr + nbgp  88 nysg = nys - nbgp  89 nyng = nyn + nbgp 79 90 ! 80 91 !-- Allocate grid arrays 81 92 ALLOCATE( ddzu(1:nzt+1), ddzw(1:nzt+1), dd2zu(1:nzt), dzu(1:nzt+1), & 82  dzw(1:nzt+1), l_grid(1:nzt), zu( 0:nzt+1), zw(0:nzt+1) ) 93 dzw(1:nzt+1), l_grid(1:nzt), zu(nzb:nzt+1), zw(nzb:nzt+1) ) 83 94 84 95 ! … …  97 108 ! 98 109 !-- Grid for atmosphere with surface at z=0 (k=0, w-grid). 99  !-- Since the w-level lies on the surface, the first u-level (staggered!)100  !-- lies below the surface (used for "mirror" boundary condition).101 110 !-- The first u-level above the surface corresponds to the top of the 102 111 !-- Prandtl-layer. 103  zu(0) = - dz * 0.5  112  113 IF ( ibc_uv_b == 0 .OR. ibc_uv_b == 2 ) THEN  114 zu(0) = 0.0  115 ! zu(0) = - dz * 0.5  116 ELSE  117 zu(0) = - dz * 0.5  118 ENDIF 104 119 zu(1) = dz * 0.5 105 120 … …  173 188 dd2zu(k) = 1.0 / ( dzu(k) + dzu(k+1) ) 174 189 ENDDO  190  191 !  192 !-- In case of FFT method or SOR swap inverse grid lenght ddzu to ddzu_fft and  193 !-- modify the lowest entry. This is necessary for atmosphere runs, because  194 !-- zu(0) and so ddzu(1) changed. Accompanied with this modified arrays  195 !-- pressure solver uses wrong values and this causes kinks in the profiles  196 !-- of turbulent quantities.  197 IF ( psolver /= 'multigrid' ) THEN  198 ALLOCATE( ddzu_pres(1:nzt+1) )  199 ddzu_pres = ddzu  200 IF( .NOT. ocean ) ddzu_pres(1) = ddzu_pres(2)  201 ENDIF 175 202 176 203 ! … …  187 214 188 215 dzu_mg(:,maximum_grid_level) = dzu  216 !  217 !-- To ensure a equally spaced grid. For ocean runs this is not necessary,  218 !-- because zu(0) does not changed so far. Also this would cause errors  219 !-- if a vertical stretching for ocean runs is used.  220 IF ( .NOT. ocean ) dzu_mg(1,maximum_grid_level) = dzu(2) 189 221 dzw_mg(:,maximum_grid_level) = dzw 190 222 nzt_l = nzt … …  239 271 !-- the flag arrays needed for the multigrid method 240 272 gls = 2**( maximum_grid_level )  273 241 274 ALLOCATE( corner_nl(nys:nyn,nxl:nxr), corner_nr(nys:nyn,nxl:nxr), & 242 275 corner_sl(nys:nyn,nxl:nxr), corner_sr(nys:nyn,nxl:nxr), & 243  nzb_local(-gls:ny+gls,-gls:nx+gls), nzb_tmp(-1:ny+1,-1:nx+1), &  276 nzb_local(-gls:ny+gls,-gls:nx+gls), &  277 nzb_tmp(-nbgp:ny+nbgp,-nbgp:nx+nbgp), & 244 278 wall_l(nys:nyn,nxl:nxr), wall_n(nys:nyn,nxl:nxr), & 245 279 wall_r(nys:nyn,nxl:nxr), wall_s(nys:nyn,nxl:nxr) ) 246  ALLOCATE( fwxm(nys-1:nyn+1,nxl-1:nxr+1), fwxp(nys-1:nyn+1,nxl-1:nxr+1), & 247  fwym(nys-1:nyn+1,nxl-1:nxr+1), fwyp(nys-1:nyn+1,nxl-1:nxr+1), & 248  fxm(nys-1:nyn+1,nxl-1:nxr+1), fxp(nys-1:nyn+1,nxl-1:nxr+1), & 249  fym(nys-1:nyn+1,nxl-1:nxr+1), fyp(nys-1:nyn+1,nxl-1:nxr+1), & 250  nzb_s_inner(nys-1:nyn+1,nxl-1:nxr+1), & 251  nzb_s_outer(nys-1:nyn+1,nxl-1:nxr+1), & 252  nzb_u_inner(nys-1:nyn+1,nxl-1:nxr+1), & 253  nzb_u_outer(nys-1:nyn+1,nxl-1:nxr+1), & 254  nzb_v_inner(nys-1:nyn+1,nxl-1:nxr+1), & 255  nzb_v_outer(nys-1:nyn+1,nxl-1:nxr+1), & 256  nzb_w_inner(nys-1:nyn+1,nxl-1:nxr+1), & 257  nzb_w_outer(nys-1:nyn+1,nxl-1:nxr+1), & 258  nzb_diff_s_inner(nys-1:nyn+1,nxl-1:nxr+1), & 259  nzb_diff_s_outer(nys-1:nyn+1,nxl-1:nxr+1), & 260  nzb_diff_u(nys-1:nyn+1,nxl-1:nxr+1), & 261  nzb_diff_v(nys-1:nyn+1,nxl-1:nxr+1), & 262  nzb_2d(nys-1:nyn+1,nxl-1:nxr+1), & 263  wall_e_x(nys-1:nyn+1,nxl-1:nxr+1), & 264  wall_e_y(nys-1:nyn+1,nxl-1:nxr+1), & 265  wall_u(nys-1:nyn+1,nxl-1:nxr+1), & 266  wall_v(nys-1:nyn+1,nxl-1:nxr+1), & 267  wall_w_x(nys-1:nyn+1,nxl-1:nxr+1), & 268  wall_w_y(nys-1:nyn+1,nxl-1:nxr+1) ) 269  270  ALLOCATE( l_wall(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )  280 ALLOCATE( fwxm(nysg:nyng,nxlg:nxrg), fwxp(nysg:nyng,nxlg:nxrg), &  281 fwym(nysg:nyng,nxlg:nxrg), fwyp(nysg:nyng,nxlg:nxrg), &  282 fxm(nysg:nyng,nxlg:nxrg), fxp(nysg:nyng,nxlg:nxrg), &  283 fym(nysg:nyng,nxlg:nxrg), fyp(nysg:nyng,nxlg:nxrg), &  284 nzb_s_inner(nysg:nyng,nxlg:nxrg), &  285 nzb_s_outer(nysg:nyng,nxlg:nxrg), &  286 nzb_u_inner(nysg:nyng,nxlg:nxrg), &  287 nzb_u_outer(nysg:nyng,nxlg:nxrg), &  288 nzb_v_inner(nysg:nyng,nxlg:nxrg), &  289 nzb_v_outer(nysg:nyng,nxlg:nxrg), &  290 nzb_w_inner(nysg:nyng,nxlg:nxrg), &  291 nzb_w_outer(nysg:nyng,nxlg:nxrg), &  292 nzb_diff_s_inner(nysg:nyng,nxlg:nxrg), &  293 nzb_diff_s_outer(nysg:nyng,nxlg:nxrg), &  294 nzb_diff_u(nysg:nyng,nxlg:nxrg), &  295 nzb_diff_v(nysg:nyng,nxlg:nxrg), &  296 nzb_2d(nysg:nyng,nxlg:nxrg), &  297 wall_e_x(nysg:nyng,nxlg:nxrg), &  298 wall_e_y(nysg:nyng,nxlg:nxrg), &  299 wall_u(nysg:nyng,nxlg:nxrg), &  300 wall_v(nysg:nyng,nxlg:nxrg), &  301 wall_w_x(nysg:nyng,nxlg:nxrg), &  302 wall_w_y(nysg:nyng,nxlg:nxrg) )  303  304  305  306 ALLOCATE( l_wall(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 271 307 272 308 nzb_s_inner = nzb; nzb_s_outer = nzb … …  327 363 vertical_influence(0) = vertical_influence(1) 328 364 329  DO i = nxl -1, nxr+1330  DO j = nys -1, nyn+1 365 DO i = nxlg, nxrg  366 DO j = nysg, nyng 331 367 DO k = nzb_s_inner(j,i) + 1, & 332 368 nzb_s_inner(j,i) + vertical_influence(nzb_s_inner(j,i)) … …  477 513 nzb_local(:,-gls:-1) = nzb_local(:,nx-gls+1:nx) 478 514 nzb_local(:,nx+1:nx+gls) = nzb_local(:,0:gls-1)  515  516 479 517 480 518 GOTO 12 … …  588 626 !-- nzb_s_outer: 589 627 !-- extend nzb_local east-/westwards first, then north-/southwards 590  nzb_tmp = nzb_local(- 1:ny+1,-1:nx+1) 628 nzb_tmp = nzb_local(-nbgp:ny+nbgp,-nbgp:nx+nbgp) 591 629 DO j = -1, ny + 1 592 630 DO i = 0, nx … …  620 658 !-- nzb_u_inner: 621 659 !-- extend nzb_local rightwards only 622  nzb_tmp = nzb_local(- 1:ny+1,-1:nx+1) 660 nzb_tmp = nzb_local(-nbgp:ny+nbgp,-nbgp:nx+nbgp) 623 661 DO j = -1, ny + 1 624 662 DO i = 0, nx + 1 … …  626 664 ENDDO 627 665 ENDDO 628  nzb_u_inner = nzb_tmp(nys -1:nyn+1,nxl-1:nxr+1) 666 nzb_u_inner = nzb_tmp(nysg:nyng,nxlg:nxrg) 629 667 630 668 ! … …  652 690 !-- nzb_v_inner: 653 691 !-- extend nzb_local northwards only 654  nzb_tmp = nzb_local(- 1:ny+1,-1:nx+1) 692 nzb_tmp = nzb_local(-nbgp:ny+nbgp,-nbgp:nx+nbgp) 655 693 DO i = -1, nx + 1 656 694 DO j = 0, ny + 1 … …  658 696 ENDDO 659 697 ENDDO 660  nzb_v_inner = nzb_tmp(nys- 1:nyn+1,nxl-1:nxr+1) 698 nzb_v_inner = nzb_tmp(nys-nbgp:nyn+nbgp,nxl-nbgp:nxr+nbgp) 661 699 662 700 ! … …  1096 1134 ! 1097 1135 !-- Need to set lateral boundary conditions for l_wall 1098  CALL exchange_horiz( l_wall )  1136  1137 CALL exchange_horiz( l_wall, nbgp ) 1099 1138 1100 1139 DEALLOCATE( corner_nl, corner_nr, corner_sl, corner_sr, nzb_local, & -
palm/trunk/SOURCE/init_particles.f90
r623 r667  4 4 ! Current revisions: 5 5 ! ----------------- 6  !  6 ! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng for allocation  7 ! of arrays. 7 8 ! 8 9 ! Former revisions: … …  185 186 ENDIF 186 187 187  ALLOCATE( prt_count(nzb:nzt+1,nys -1:nyn+1,nxl-1:nxr+1), &188  prt_start_index(nzb:nzt+1,nys -1:nyn+1,nxl-1:nxr+1), & 188 ALLOCATE( prt_count(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &  189 prt_start_index(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 189 190 particle_mask(maximum_number_of_particles), & 190 191 particles(maximum_number_of_particles) ) … …  209 210 !-- particles, which can be also periodically released at later times. 210 211 !-- Also allocate array for particle tail coordinates, if needed. 211  ALLOCATE( prt_count(nzb:nzt+1,nys -1:nyn+1,nxl-1:nxr+1), &212  prt_start_index(nzb:nzt+1,nys -1:nyn+1,nxl-1:nxr+1), & 212 ALLOCATE( prt_count(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &  213 prt_start_index(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 213 214 particle_mask(maximum_number_of_particles), & 214 215 particles(maximum_number_of_particles) ) -
palm/trunk/SOURCE/init_pegrid.f90
r647 r667  4 4 ! Current revisions: 5 5 ! -----------------  6 !  7 ! Moved determination of target_id's from init_coupling  8 !  9 ! Determination of parameters needed for coupling (coupling_topology, ngp_a, ngp_o)  10 ! with different grid/processor-topology in ocean and atmosphere  11 !  12 !  13 ! Adaption of ngp_xy, ngp_y to a dynamic number of ghost points.  14 ! The maximum_grid_level changed from 1 to 0. 0 is the normal grid, 1 to  15 ! maximum_grid_level the grids for multigrid, in which 0 and 1 are normal grids.  16 ! This distinction is due to reasons of data exchange and performance for the  17 ! normal grid and grids in poismg.  18 ! The definition of MPI-Vectors adapted to a dynamic numer of ghost points.  19 ! New MPI-Vectors for data exchange between left and right boundaries added.  20 ! This is due to reasons of performance (10% faster). 6 21 ! 7 22 ! ATTENTION: nnz_x undefined problem still has to be solved!!!!!!!! … …  79 94 80 95  96 81 97 IMPLICIT NONE 82 98 … …  88 104 89 105 INTEGER, DIMENSION(:), ALLOCATABLE :: ind_all, nxlf, nxrf, nynf, nysf  106  107 INTEGER, DIMENSION(2) :: pdims_remote 90 108 91 109 LOGICAL :: found … …  103 121 104 122 #if defined( __parallel )  123 105 124 ! 106 125 !-- Determine the processor topology or check it, if prescribed by the user … …  624 643 #endif 625 644  645 !  646 !-- Determine the number of ghost points  647 IF (scalar_advec == 'ws-scheme' .OR. momentum_advec == 'ws-scheme') THEN  648 nbgp = 3  649 ELSE  650 nbgp = 1  651 END IF  652 626 653 ! 627 654 !-- In case of coupled runs, create a new MPI derived datatype for the 628 655 !-- exchange of surface (xy) data . 629 656 !-- Gridpoint number for the exchange of ghost points (xy-plane) 630  ngp_xy = ( nxr - nxl + 3 ) * ( nyn - nys + 3 )  657  658 ngp_xy = ( nxr - nxl + 1 + 2 * nbgp ) * ( nyn - nys + 1 + 2 * nbgp ) 631 659 632 660 ! … …  635 663 CALL MPI_TYPE_VECTOR( ngp_xy, 1, nzt-nzb+2, MPI_REAL, type_xy, ierr ) 636 664 CALL MPI_TYPE_COMMIT( type_xy, ierr )  665  666  667 IF ( TRIM( coupling_mode ) .NE. 'uncoupled' ) THEN  668  669 !  670 !-- Pass the number of grid points of the atmosphere model to  671 !-- the ocean model and vice versa  672 IF ( coupling_mode == 'atmosphere_to_ocean' ) THEN  673  674 nx_a = nx  675 ny_a = ny  676  677 IF ( myid == 0 ) THEN  678 CALL MPI_SEND( nx_a, 1, MPI_INTEGER, numprocs, 1, &  679 comm_inter, ierr )  680 CALL MPI_SEND( ny_a, 1, MPI_INTEGER, numprocs, 2, &  681 comm_inter, ierr )  682 CALL MPI_SEND( pdims, 2, MPI_INTEGER, numprocs, 3, &  683 comm_inter, ierr )  684 CALL MPI_RECV( nx_o, 1, MPI_INTEGER, numprocs, 4, &  685 comm_inter, status, ierr )  686 CALL MPI_RECV( ny_o, 1, MPI_INTEGER, numprocs, 5, &  687 comm_inter, status, ierr )  688 CALL MPI_RECV( pdims_remote, 2, MPI_INTEGER, numprocs, 6, &  689 comm_inter, status, ierr )  690 ENDIF  691  692 CALL MPI_BCAST( nx_o, 1, MPI_INTEGER, 0, comm2d, ierr)  693 CALL MPI_BCAST( ny_o, 1, MPI_INTEGER, 0, comm2d, ierr)  694 CALL MPI_BCAST( pdims_remote, 2, MPI_INTEGER, 0, comm2d, ierr)  695  696 ELSEIF ( coupling_mode == 'ocean_to_atmosphere' ) THEN  697  698 nx_o = nx  699 ny_o = ny  700  701 IF ( myid == 0 ) THEN  702 CALL MPI_RECV( nx_a, 1, MPI_INTEGER, 0, 1, &  703 comm_inter, status, ierr )  704 CALL MPI_RECV( ny_a, 1, MPI_INTEGER, 0, 2, &  705 comm_inter, status, ierr )  706 CALL MPI_RECV( pdims_remote, 2, MPI_INTEGER, 0, 3, &  707 comm_inter, status, ierr )  708 CALL MPI_SEND( nx_o, 1, MPI_INTEGER, 0, 4, &  709 comm_inter, ierr )  710 CALL MPI_SEND( ny_o, 1, MPI_INTEGER, 0, 5, &  711 comm_inter, ierr )  712 CALL MPI_SEND( pdims, 2, MPI_INTEGER, 0, 6, &  713 comm_inter, ierr )  714 ENDIF  715  716 CALL MPI_BCAST( nx_a, 1, MPI_INTEGER, 0, comm2d, ierr)  717 CALL MPI_BCAST( ny_a, 1, MPI_INTEGER, 0, comm2d, ierr)  718 CALL MPI_BCAST( pdims_remote, 2, MPI_INTEGER, 0, comm2d, ierr)  719  720 ENDIF  721  722 ngp_a = (nx_a+1+2*nbgp)*(ny_a+1+2*nbgp)  723 ngp_o = (nx_o+1+2*nbgp)*(ny_o+1+2*nbgp)  724  725 !  726 !-- determine if the horizontal grid and the number of PEs  727 !-- in ocean and atmosphere is same or not  728 !-- (different number of PEs still not implemented)  729 IF ( nx_o == nx_a .AND. ny_o == ny_a .AND. &  730 pdims(1) == pdims_remote(1) .AND. pdims(2) == pdims_remote(2) ) &  731 THEN  732 coupling_topology = 0  733 ELSE  734 coupling_topology = 1  735 ENDIF  736  737 !  738 !-- Determine the target PEs for the exchange between ocean and  739 !-- atmosphere (comm2d)  740 IF ( coupling_topology == 0) THEN  741 IF ( TRIM( coupling_mode ) .EQ. 'atmosphere_to_ocean' ) THEN  742 target_id = myid + numprocs  743 ELSE  744 target_id = myid  745 ENDIF  746  747 ELSE  748  749 !  750 !-- In case of nonequivalent topology in ocean and atmosphere only for  751 !-- PE0 in ocean and PE0 in atmosphere a target_id is needed, since  752 !-- data echxchange between ocean and atmosphere will be done only by  753 !-- those PEs.  754 IF ( myid == 0 ) THEN  755 IF ( TRIM( coupling_mode ) .EQ. 'atmosphere_to_ocean' ) THEN  756 target_id = numprocs  757 ELSE  758 target_id = 0  759 ENDIF  760 print*, coupling_mode, myid, " -> ", target_id, "numprocs: ", numprocs  761 ENDIF  762 ENDIF  763  764 ENDIF  765  766 637 767 #endif 638 768 … …  854 984 ELSE 855 985 856  maximum_grid_level = 1 986 maximum_grid_level = 0 857 987 858 988 ENDIF … …  863 993 ! 864 994 !-- Gridpoint number for the exchange of ghost points (y-line for 2D-arrays) 865  ngp_y = nyn - nys + 1  995 ngp_y = nyn - nys + 1 + 2 * nbgp 866 996 867 997 ! 868 998 !-- Define a new MPI derived datatype for the exchange of ghost points in 869 999 !-- y-direction for 2D-arrays (line) 870 Â