Context Navigation

start_em.F

Last change on this file was 1, checked in by lfita, 10 years ago

-- --- Opening of the WRF+LMDZ coupling repository --- -- -

WRF: version v3.3
LMDZ: version v1818

More details in:

File size: 56.2 KB

Line
1	!-------------------------------------------------------------------
2
3	SUBROUTINE start_domain_em ( grid, allowed_to_read &
4	! Actual arguments generated from Registry
5	# include "dummy_new_args.inc"
6	!
7	)
8
9	USE module_domain, ONLY : domain, wrfu_timeinterval, get_ijk_from_grid, &
10	domain_setgmtetc
11	USE module_state_description
12	USE module_model_constants
13	USE module_bc, ONLY : boundary_condition_check, set_physical_bc2d
14	USE module_bc_em
15	USE module_configure, ONLY : grid_config_rec_type
16	USE module_tiles, ONLY : set_tiles
17	#ifdef DM_PARALLEL
18	USE module_dm, ONLY : wrf_dm_min_real, wrf_dm_maxval, ntasks_x, ntasks_y, &
19	local_communicator_periodic, local_communicator, mytask, ntasks
20	#else
21	USE module_dm, ONLY : wrf_dm_min_real
22	#endif
23	USE module_comm_dm
24
25	USE module_physics_init
26	USE module_fr_sfire_driver_wrf, ONLY : sfire_driver_em_init
27	USE module_stoch, ONLY : SETUP_STOCH,rand_seed, update_stoch
28	#ifdef WRF_CHEM
29	USE module_aerosols_sorgam, ONLY: sum_pm_sorgam
30	USE module_gocart_aerosols, ONLY: sum_pm_gocart
31	USE module_mosaic_driver, ONLY: sum_pm_mosaic
32	USE module_input_tracer, ONLY: initialize_tracer
33	#endif
34
35	!!debug
36	!USE module_compute_geop
37
38	USE module_model_constants
39	USE module_avgflx_em, ONLY : zero_avgflx
40
41	IMPLICIT NONE
42	! Input data.
43	TYPE (domain) :: grid
44
45	LOGICAL , INTENT(IN) :: allowed_to_read
46
47	! Definitions of dummy arguments to this routine (generated from Registry).
48	# include "dummy_new_decl.inc"
49
50	! Structure that contains run-time configuration (namelist) data for domain
51	TYPE (grid_config_rec_type) :: config_flags
52
53	! Local data
54	INTEGER :: &
55	ids, ide, jds, jde, kds, kde, &
56	ims, ime, jms, jme, kms, kme, &
57	ips, ipe, jps, jpe, kps, kpe, &
58	its, ite, jts, jte, kts, kte, &
59	ij,i,j,k,ii,jj,kk,loop,error,l
60
61	INTEGER :: imsx, imex, jmsx, jmex, kmsx, kmex, &
62	ipsx, ipex, jpsx, jpex, kpsx, kpex, &
63	imsy, imey, jmsy, jmey, kmsy, kmey, &
64	ipsy, ipey, jpsy, jpey, kpsy, kpey
65
66	INTEGER :: i_m
67
68	REAL :: p00, t00, a, tiso, p_surf, pd_surf, temp
69	#ifdef WRF_CHEM
70	REAL RGASUNIV ! universal gas constant [ J/mol-K ]
71	PARAMETER ( RGASUNIV = 8.314510 )
72	REAL,DIMENSION(grid%sm31:grid%em31,grid%sm32:grid%em32,grid%sm33:grid%em33) :: &
73	z_at_w,convfac
74	REAL :: tempfac
75	#endif
76
77	REAL :: qvf1, qvf2, qvf
78	REAL :: MPDT
79	REAL :: spongeweight
80	LOGICAL :: first_trip_for_this_domain, start_of_simulation, fill_w_flag
81	LOGICAL, EXTERNAL :: wrf_dm_on_monitor
82
83	#ifndef WRF_CHEM
84	REAL,ALLOCATABLE,DIMENSION(:,:,:) :: cldfra_old
85	#endif
86
87	REAL :: lat1 , lat2 , lat3 , lat4
88	REAL :: lon1 , lon2 , lon3 , lon4
89	INTEGER :: num_points_lat_lon , iloc , jloc
90	CHARACTER (LEN=132) :: message
91	TYPE(WRFU_TimeInterval) :: stepTime
92	REAL, DIMENSION(:,:), ALLOCATABLE :: clat_glob
93	logical :: f_flux ! flag for computing averaged fluxes in cu_gd
94
95	INTEGER :: idex, jdex
96	INTEGER :: im1,ip1,jm1,jp1
97	REAL :: hx,hy,pi
98
99	CALL get_ijk_from_grid ( grid , &
100	ids, ide, jds, jde, kds, kde, &
101	ims, ime, jms, jme, kms, kme, &
102	ips, ipe, jps, jpe, kps, kpe, &
103	imsx, imex, jmsx, jmex, kmsx, kmex, &
104	ipsx, ipex, jpsx, jpex, kpsx, kpex, &
105	imsy, imey, jmsy, jmey, kmsy, kmey, &
106	ipsy, ipey, jpsy, jpey, kpsy, kpey )
107
108	kts = kps ; kte = kpe ! note that tile is entire patch
109	its = ips ; ite = ipe ! note that tile is entire patch
110	jts = jps ; jte = jpe ! note that tile is entire patch
111	#ifndef WRF_CHEM
112	ALLOCATE(CLDFRA_OLD(IMS:IME,KMS:KME,JMS:JME),STAT=I) ; CLDFRA_OLD = 0.
113	#endif
114	CALL model_to_grid_config_rec ( grid%id , model_config_rec , config_flags )
115
116	IF ( ( MOD (ide-ids,config_flags%parent_grid_ratio) .NE. 0 ) .OR. &
117	( MOD (jde-jds,config_flags%parent_grid_ratio) .NE. 0 ) ) THEN
118	WRITE(message, FMT='("Nested dimensions are illegal for domain ",I2,": Both &
119	&MOD(",I4,"-",I1,",",I2,") and MOD(",I4,"-",I1,",",I2,") must = 0" )') &
120	grid%id,ide,ids,config_flags%parent_grid_ratio,jde,jds,config_flags%parent_grid_ratio
121	CALL wrf_error_fatal ( message )
122	END IF
123
124	! --- SETUP AND INITIALIZE STOCHASTIC KINETIC ENERGY BACKSCATTER SCHEME ---
125
126	IF ( (grid%id.eq.1) .AND. ( .NOT. grid%did_stoch ) ) THEN
127	grid%did_stoch = .TRUE.
128	stoch_force_select: SELECT CASE(config_flags%stoch_force_opt)
129
130	CASE (STOCH_BACKSCATTER)
131
132	IF ( wrf_dm_on_monitor () ) THEN
133	CALL rand_seed ( config_flags, grid%SEED1, grid%SEED2, grid%NENS )
134	ENDIF
135	#ifdef DM_PARALLEL
136	CALL wrf_dm_bcast_bytes ( grid%SEED1, IWORDSIZE )
137	CALL wrf_dm_bcast_bytes ( grid%SEED2, IWORDSIZE )
138	#endif
139
140	call SETUP_STOCH(grid%VERTSTRUCC,grid%VERTSTRUCS, &
141	grid%SPT_AMP,grid%SPSTREAM_AMP, &
142	grid%stoch_vertstruc_opt, &
143	grid%SEED1,grid%SEED2,grid%time_step, &
144	grid%DX,grid%DY, &
145	grid%TOT_BACKSCAT_PSI,grid%TOT_BACKSCAT_T, &
146	ids, ide, jds, jde, kds, kde, &
147	ims, ime, jms, jme, kms, kme, &
148	its, ite, jts, jte, kts, kte )
149
150	END SELECT stoch_force_select
151	END IF
152	! --- END SETUP STOCHASTIC KINETIC ENERGY BACKSCATTER SCHEME ----------
153
154	IF ( config_flags%polar ) THEN
155	!write(0,*)__FILE__,__LINE__,' clat ',ips,ipe,jps,jpe
156	!do j = jps,jpe
157	!write(0,*)__FILE__,__LINE__,' clat ',ids,j,grid%clat(ips,j)
158	!enddo
159
160	#ifdef DM_PARALLEL
161	! WARNING: this might present scaling issues on very large numbers of processors
162	ALLOCATE( clat_glob(ids:ide,jds:jde) )
163
164	CALL wrf_patch_to_global_real ( grid%clat, clat_glob, grid%domdesc, 'xy', 'xy', &
165	ids, ide, jds, jde, 1, 1, &
166	ims, ime, jms, jme, 1, 1, &
167	its, ite, jts, jte, 1, 1 )
168
169	CALL wrf_dm_bcast_real ( clat_glob , (ide-ids+1)*(jde-jds+1) )
170
171	grid%clat_xxx(ipsx:ipex,jpsx:jpex) = clat_glob(ipsx:ipex,jpsx:jpex)
172
173	DEALLOCATE( clat_glob )
174	#endif
175	ENDIF
176
177	! here we check to see if the boundary conditions are set properly
178
179	CALL boundary_condition_check( config_flags, bdyzone, error, grid%id )
180
181	!kludge - need to stop CG from resetting precip and phys tendencies to zero
182	! when we are in here due to a nest being spawned, we want to still
183	! recompute the base state, but that is about it
184	! This is temporary and will need to be changed when grid%itimestep is removed.
185
186	IF ( grid%itimestep .EQ. 0 ) THEN
187	first_trip_for_this_domain = .TRUE.
188	ELSE
189	first_trip_for_this_domain = .FALSE.
190	END IF
191
192	IF ( .not. ( config_flags%restart .or. grid%moved ) ) THEN
193	grid%itimestep=0
194	ENDIF
195
196	IF ( config_flags%restart .or. grid%moved ) THEN
197	first_trip_for_this_domain = .TRUE.
198	ENDIF
199
200	! wig: Add a combined exponential+linear weight on the mother boundaries
201	! following code changes by Ruby Leung. For the nested grid, there
202	! appears to be some problems when a sponge is used. The points where
203	! processors meet have problematic values.
204
205	CALL lbc_fcx_gcx ( grid%fcx , grid%gcx , grid%spec_bdy_width , &
206	grid%spec_zone , grid%relax_zone , grid%dt , config_flags%spec_exp , &
207	config_flags%specified , config_flags%nested )
208
209	IF ( config_flags%nested ) THEN
210	grid%dtbc = 0.
211	ENDIF
212
213	IF ( ( grid%id .NE. 1 ) .AND. ( .NOT. config_flags%input_from_file ) ) THEN
214
215	! Every time a domain starts or every time a domain moves, this routine is called. We want
216	! the center (middle) lat/lon of the grid for the metacode. The lat/lon values are
217	! defined at mass points. Depending on the even/odd points in the SN and WE directions,
218	! we end up with the middle point as either 1 point or an average of either 2 or 4 points.
219	! Add to this, the need to make sure that we are on the correct patch to retrieve the
220	! value of the lat/lon, AND that the lat/lons (for an average) may not all be on the same
221	! patch. Once we find the correct value for lat lon, we need to keep it around on all patches,
222	! which is where the wrf_dm_min_real calls come in.
223	! If this is the most coarse domain, we do not go in here. Also, if there is an input file
224	! (which has the right values for the middle lat/lon) we do not go in this IF test.
225
226	IF ( ( MOD(ide,2) .EQ. 0 ) .AND. ( MOD(jde,2) .EQ. 0 ) ) THEN
227	num_points_lat_lon = 1
228	iloc = ide/2
229	jloc = jde/2
230	IF ( ( ips .LE. iloc ) .AND. ( ipe .GE. iloc ) .AND. &
231	( jps .LE. jloc ) .AND. ( jpe .GE. jloc ) ) THEN
232	lat1 = grid%xlat (iloc,jloc)
233	lon1 = grid%xlong(iloc,jloc)
234	ELSE
235	lat1 = 99999.
236	lon1 = 99999.
237	END IF
238	lat1 = wrf_dm_min_real ( lat1 )
239	lon1 = wrf_dm_min_real ( lon1 )
240	CALL nl_set_cen_lat ( grid%id , lat1 )
241	CALL nl_set_cen_lon ( grid%id , lon1 )
242	ELSE IF ( ( MOD(ide,2) .NE. 0 ) .AND. ( MOD(jde,2) .EQ. 0 ) ) THEN
243	num_points_lat_lon = 2
244	iloc = (ide-1)/2
245	jloc = jde /2
246	IF ( ( ips .LE. iloc ) .AND. ( ipe .GE. iloc ) .AND. &
247	( jps .LE. jloc ) .AND. ( jpe .GE. jloc ) ) THEN
248	lat1 = grid%xlat (iloc,jloc)
249	lon1 = grid%xlong(iloc,jloc)
250	ELSE
251	lat1 = 99999.
252	lon1 = 99999.
253	END IF
254	lat1 = wrf_dm_min_real ( lat1 )
255	lon1 = wrf_dm_min_real ( lon1 )
256
257	iloc = (ide+1)/2
258	jloc = jde /2
259	IF ( ( ips .LE. iloc ) .AND. ( ipe .GE. iloc ) .AND. &
260	( jps .LE. jloc ) .AND. ( jpe .GE. jloc ) ) THEN
261	lat2 = grid%xlat (iloc,jloc)
262	lon2 = grid%xlong(iloc,jloc)
263	ELSE
264	lat2 = 99999.
265	lon2 = 99999.
266	END IF
267	lat2 = wrf_dm_min_real ( lat2 )
268	lon2 = wrf_dm_min_real ( lon2 )
269
270	CALL nl_set_cen_lat ( grid%id , ( lat1 + lat2 ) * 0.50 )
271	CALL nl_set_cen_lon ( grid%id , ( lon1 + lon2 ) * 0.50 )
272	ELSE IF ( ( MOD(ide,2) .EQ. 0 ) .AND. ( MOD(jde,2) .NE. 0 ) ) THEN
273	num_points_lat_lon = 2
274	iloc = ide /2
275	jloc = (jde-1)/2
276	IF ( ( ips .LE. iloc ) .AND. ( ipe .GE. iloc ) .AND. &
277	( jps .LE. jloc ) .AND. ( jpe .GE. jloc ) ) THEN
278	lat1 = grid%xlat (iloc,jloc)
279	lon1 = grid%xlong(iloc,jloc)
280	ELSE
281	lat1 = 99999.
282	lon1 = 99999.
283	END IF
284	lat1 = wrf_dm_min_real ( lat1 )
285	lon1 = wrf_dm_min_real ( lon1 )
286
287	iloc = ide /2
288	jloc = (jde+1)/2
289	IF ( ( ips .LE. iloc ) .AND. ( ipe .GE. iloc ) .AND. &
290	( jps .LE. jloc ) .AND. ( jpe .GE. jloc ) ) THEN
291	lat2 = grid%xlat (iloc,jloc)
292	lon2 = grid%xlong(iloc,jloc)
293	ELSE
294	lat2 = 99999.
295	lon2 = 99999.
296	END IF
297	lat2 = wrf_dm_min_real ( lat2 )
298	lon2 = wrf_dm_min_real ( lon2 )
299
300	CALL nl_set_cen_lat ( grid%id , ( lat1 + lat2 ) * 0.50 )
301	CALL nl_set_cen_lon ( grid%id , ( lon1 + lon2 ) * 0.50 )
302	ELSE IF ( ( MOD(ide,2) .NE. 0 ) .AND. ( MOD(jde,2) .NE. 0 ) ) THEN
303	num_points_lat_lon = 4
304	iloc = (ide-1)/2
305	jloc = (jde-1)/2
306	IF ( ( ips .LE. iloc ) .AND. ( ipe .GE. iloc ) .AND. &
307	( jps .LE. jloc ) .AND. ( jpe .GE. jloc ) ) THEN
308	lat1 = grid%xlat (iloc,jloc)
309	lon1 = grid%xlong(iloc,jloc)
310	ELSE
311	lat1 = 99999.
312	lon1 = 99999.
313	END IF
314	lat1 = wrf_dm_min_real ( lat1 )
315	lon1 = wrf_dm_min_real ( lon1 )
316
317	iloc = (ide+1)/2
318	jloc = (jde-1)/2
319	IF ( ( ips .LE. iloc ) .AND. ( ipe .GE. iloc ) .AND. &
320	( jps .LE. jloc ) .AND. ( jpe .GE. jloc ) ) THEN
321	lat2 = grid%xlat (iloc,jloc)
322	lon2 = grid%xlong(iloc,jloc)
323	ELSE
324	lat2 = 99999.
325	lon2 = 99999.
326	END IF
327	lat2 = wrf_dm_min_real ( lat2 )
328	lon2 = wrf_dm_min_real ( lon2 )
329
330	iloc = (ide-1)/2
331	jloc = (jde+1)/2
332	IF ( ( ips .LE. iloc ) .AND. ( ipe .GE. iloc ) .AND. &
333	( jps .LE. jloc ) .AND. ( jpe .GE. jloc ) ) THEN
334	lat3 = grid%xlat (iloc,jloc)
335	lon3 = grid%xlong(iloc,jloc)
336	ELSE
337	lat3 = 99999.
338	lon3 = 99999.
339	END IF
340	lat3 = wrf_dm_min_real ( lat3 )
341	lon3 = wrf_dm_min_real ( lon3 )
342
343	iloc = (ide+1)/2
344	jloc = (jde+1)/2
345	IF ( ( ips .LE. iloc ) .AND. ( ipe .GE. iloc ) .AND. &
346	( jps .LE. jloc ) .AND. ( jpe .GE. jloc ) ) THEN
347	lat4 = grid%xlat (iloc,jloc)
348	lon4 = grid%xlong(iloc,jloc)
349	ELSE
350	lat4 = 99999.
351	lon4 = 99999.
352	END IF
353	lat4 = wrf_dm_min_real ( lat4 )
354	lon4 = wrf_dm_min_real ( lon4 )
355
356	CALL nl_set_cen_lat ( grid%id , ( lat1 + lat2 + lat3 + lat4 ) * 0.25 )
357	CALL nl_set_cen_lon ( grid%id , ( lon1 + lon2 + lon3 + lon4 ) * 0.25 )
358	END IF
359	END IF
360
361	IF ( .NOT. config_flags%restart .AND. &
362	(( config_flags%input_from_hires ) .OR. ( config_flags%input_from_file ))) THEN
363
364	IF ( config_flags%map_proj .EQ. 0 ) THEN
365	CALL wrf_error_fatal ( 'start_domain: Idealized case cannot have a separate nested input file' )
366	END IF
367
368	IF ( config_flags%use_baseparam_fr_nml ) then
369	CALL nl_get_base_pres ( 1 , p00 )
370	CALL nl_get_base_temp ( 1 , t00 )
371	CALL nl_get_base_lapse ( 1 , a )
372	CALL nl_get_iso_temp ( 1 , tiso )
373
374	ELSE
375	! get these constants from model data
376
377	t00 = grid%t00
378	p00 = grid%p00
379	a = grid%tlp
380	tiso = grid%tiso
381
382	IF (t00 .LT. 100. .or. p00 .LT. 10000.) THEN
383	WRITE(wrf_err_message,*)&
384	'start_em: did not find base state parameters in wrfinput. Add use_baseparam_fr_nml = .t. in &dynamics and rerun'
385	CALL wrf_error_fatal(TRIM(wrf_err_message))
386	ENDIF
387
388	ENDIF
389
390	! Base state potential temperature and inverse density (alpha = 1/rho) from
391	! the half eta levels and the base-profile surface pressure. Compute 1/rho
392	! from equation of state. The potential temperature is a perturbation from t0.
393
394	DO j = jts, MIN(jte,jde-1)
395	DO i = its, MIN(ite,ide-1)
396
397	! Base state pressure is a function of eta level and terrain, only, plus
398	! the hand full of constants: p00 (sea level pressure, Pa), t00 (sea level
399	! temperature, K), and A (temperature difference, from 1000 mb to 300 mb, K).
400
401	p_surf = p00 * EXP ( -t00/a + ( (t00/a)*2 - 2.ggrid%ht(i,j)/a/r_d ) *0.5 )
402
403	DO k = 1, kte-1
404	grid%pb(i,k,j) = grid%znu(k)*(p_surf - grid%p_top) + grid%p_top
405	temp = MAX ( tiso, t00 + A*LOG(grid%pb(i,k,j)/p00) )
406	grid%t_init(i,k,j) = temp(p00/grid%pb(i,k,j))*(r_d/cp) - t0
407	! grid%t_init(i,k,j) = (t00 + ALOG(grid%pb(i,k,j)/p00))(p00/grid%pb(i,k,j))**(r_d/cp) - t0
408	grid%alb(i,k,j) = (r_d/p1000mb)(grid%t_init(i,k,j)+t0)(grid%pb(i,k,j)/p1000mb)**cvpm
409	END DO
410
411	! Base state mu is defined as base state surface pressure minus grid%p_top
412
413	grid%mub(i,j) = p_surf - grid%p_top
414
415	! Integrate base geopotential, starting at terrain elevation. This assures that
416	! the base state is in exact hydrostatic balance with respect to the model equations.
417	! This field is on full levels.
418
419	grid%phb(i,1,j) = grid%ht(i,j) * g
420	DO k = 2,kte
421	grid%phb(i,k,j) = grid%phb(i,k-1,j) - grid%dnw(k-1)grid%mub(i,j)grid%alb(i,k-1,j)
422	END DO
423	END DO
424	END DO
425
426	ENDIF
427
428	IF(.not.config_flags%restart)THEN
429
430	! if this is for a nested domain, the defined/interpolated fields are the _2
431
432	IF ( first_trip_for_this_domain ) THEN
433
434	! data that is expected to be zero must be explicitly initialized as such
435	! grid%h_diabatic = 0.
436
437	DO j = jts,min(jte,jde-1)
438	DO k = kts,kte-1
439	DO i = its, min(ite,ide-1)
440	IF ( grid%imask_nostag(i,j) .EQ. 1 ) THEN
441	grid%t_1(i,k,j)=grid%t_2(i,k,j)
442	ENDIF
443	ENDDO
444	ENDDO
445	ENDDO
446
447	DO j = jts,min(jte,jde-1)
448	DO k = kts,kte
449	DO i = its, min(ite,ide-1)
450	grid%ph_1(i,k,j)=grid%ph_2(i,k,j)
451	ENDDO
452	ENDDO
453	ENDDO
454
455	DO j = jts,min(jte,jde-1)
456	DO i = its, min(ite,ide-1)
457	IF ( grid%imask_nostag(i,j) .EQ. 1 ) THEN
458	grid%mu_1(i,j)=grid%mu_2(i,j)
459	ENDIF
460	ENDDO
461	ENDDO
462	END IF
463
464	! reconstitute base-state fields
465
466	IF(config_flags%max_dom .EQ. 1)THEN
467	! with single domain, grid%t_init from wrfinput is OK to use
468	DO j = jts,min(jte,jde-1)
469	DO k = kts,kte-1
470	DO i = its, min(ite,ide-1)
471	IF ( grid%imask_nostag(i,j) .EQ. 1 ) THEN
472	grid%pb(i,k,j) = grid%znu(k)*grid%mub(i,j)+grid%p_top
473	grid%alb(i,k,j) = (r_d/p1000mb)(grid%t_init(i,k,j)+t0)(grid%pb(i,k,j)/p1000mb)**cvpm
474	ENDIF
475	ENDDO
476	ENDDO
477	ENDDO
478	ELSE
479	! with nests, grid%t_init generally needs recomputations (since it is not interpolated)
480	DO j = jts,min(jte,jde-1)
481	DO k = kts,kte-1
482	DO i = its, min(ite,ide-1)
483	IF ( grid%imask_nostag(i,j) .EQ. 1 ) THEN
484	grid%pb(i,k,j) = grid%znu(k)*grid%mub(i,j)+grid%p_top
485	grid%alb(i,k,j) = -grid%rdnw(k)*(grid%phb(i,k+1,j)-grid%phb(i,k,j))/grid%mub(i,j)
486	grid%t_init(i,k,j) = grid%alb(i,k,j)(p1000mb/r_d)/((grid%pb(i,k,j)/p1000mb)*cvpm) - t0
487	ENDIF
488	ENDDO
489	ENDDO
490	ENDDO
491	ENDIF
492
493	! Use equations from calc_p_rho_phi to derive p and al from ph
494
495	DO j=jts,min(jte,jde-1)
496	DO k=kts,kte-1
497	DO i=its,min(ite,ide-1)
498	qvf = 1.+rvovrd*moist(i,k,j,P_QV)
499	grid%al(i,k,j)=-1./(grid%mub(i,j)+grid%mu_1(i,j))(grid%alb(i,k,j)grid%mu_1(i,j) &
500	+grid%rdnw(k)*(grid%ph_1(i,k+1,j)-grid%ph_1(i,k,j)))
501	grid%p(i,k,j)=p1000mb( (r_d(t0+grid%t_1(i,k,j))*qvf)/ &
502	(p1000mb(grid%al(i,k,j)+grid%alb(i,k,j))) )*cpovcv &
503	-grid%pb(i,k,j)
504	grid%p_hyd(i,k,j) = grid%p(i,k,j) + grid%pb(i,k,j)
505	grid%alt(i,k,j) = grid%al(i,k,j) + grid%alb(i,k,j)
506	ENDDO
507	ENDDO
508	ENDDO
509
510	#if 0
511	DO j = jts,min(jte,jde-1)
512
513	k = kte-1
514	DO i = its, min(ite,ide-1)
515	IF ( grid%imask_nostag(i,j) .EQ. 1 ) THEN
516	qvf1 = 0.5*(moist(i,k,j,P_QV)+moist(i,k,j,P_QV))
517	qvf2 = 1./(1.+qvf1)
518	qvf1 = qvf1*qvf2
519	grid%p(i,k,j) = - 0.5(grid%mu_1(i,j)+qvf1grid%mub(i,j))/grid%rdnw(k)/qvf2
520	grid%p_hyd(i,k,j) = grid%p(i,k,j) + grid%pb(i,k,j)
521	qvf = 1. + rvovrd*moist(i,k,j,P_QV)
522	grid%alt(i,k,j) = (r_d/p1000mb)(grid%t_1(i,k,j)+t0)qvf(((grid%p(i,k,j)+grid%pb(i,k,j))/p1000mb)*cvpm)
523	grid%al(i,k,j) = grid%alt(i,k,j) - grid%alb(i,k,j)
524	ENDIF
525	ENDDO
526
527	DO k = kte-2, 1, -1
528	DO i = its, min(ite,ide-1)
529	IF ( grid%imask_nostag(i,j) .EQ. 1 ) THEN
530	qvf1 = 0.5*(moist(i,k,j,P_QV)+moist(i,k+1,j,P_QV))
531	qvf2 = 1./(1.+qvf1)
532	qvf1 = qvf1*qvf2
533	grid%p(i,k,j) = grid%p(i,k+1,j) - (grid%mu_1(i,j) + qvf1*grid%mub(i,j))/qvf2/grid%rdn(k+1)
534	grid%p_hyd(i,k,j) = grid%p(i,k,j) + grid%pb(i,k,j)
535	qvf = 1. + rvovrd*moist(i,k,j,P_QV)
536	grid%alt(i,k,j) = (r_d/p1000mb)(grid%t_1(i,k,j)+t0)qvf* &
537	(((grid%p(i,k,j)+grid%pb(i,k,j))/p1000mb)**cvpm)
538	grid%al(i,k,j) = grid%alt(i,k,j) - grid%alb(i,k,j)
539	ENDIF
540	ENDDO
541	ENDDO
542
543	ENDDO
544	#endif
545
546	ENDIF
547
548	IF ( grid%press_adj .and. ( grid%id .NE. 1 ) .AND. .NOT. ( config_flags%restart ) .AND. &
549	( ( config_flags%input_from_hires ) .OR. ( config_flags%input_from_file ) ) ) THEN
550	DO j = jts, MIN(jte,jde-1)
551	DO i = its, MIN(ite,ide-1)
552	grid%mu_2(i,j) = grid%mu_2(i,j) + grid%al(i,1,j) / ( grid%alt(i,1,j) * grid%alb(i,1,j) ) * &
553	g * ( grid%ht(i,j) - grid%ht_fine(i,j) )
554	END DO
555	END DO
556	DO j = jts,min(jte,jde-1)
557	DO i = its, min(ite,ide-1)
558	grid%mu_1(i,j)=grid%mu_2(i,j)
559	ENDDO
560	ENDDO
561
562	END IF
563
564	IF ( first_trip_for_this_domain ) THEN
565
566	CALL wrf_debug ( 100 , 'start_domain_em: Before call to phy_init' )
567
568	! namelist MPDT does not exist yet, so set it here
569	! MPDT is the call frequency for microphysics in minutes (0 means every step)
570	MPDT = 0.
571
572	! set GMT outside of phy_init because phy_init may not be called on this
573	! process if, for example, it is a moving nest and if this part of the domain is not
574	! being initialized (not the leading edge).
575	CALL domain_setgmtetc( grid, start_of_simulation )
576
577	!-----------------------------------------------------------------------------
578	! Adaptive time step: Added by T. Hutchinson, WSI 11/6/07
579	!
580	!
581
582	IF ( ( grid%use_adaptive_time_step ) .AND. &
583	( ( grid%dfi_opt .EQ. DFI_NODFI ) .OR. ( grid%dfi_stage .EQ. DFI_FST ) ) ) THEN
584
585	! Calculate any variables that were not set
586
587	if (grid%starting_time_step == -1) then
588	grid%starting_time_step = NINT(6 * MIN(grid%dx,grid%dy) / 1000)
589	endif
590
591	if (grid%max_time_step == -1) then
592	grid%max_time_step = 3*grid%starting_time_step
593	endif
594
595	if (grid%min_time_step == -1) then
596	grid%min_time_step = 0.5*grid%starting_time_step
597	endif
598
599	! Set a starting timestep.
600
601	grid%dt = grid%starting_time_step
602
603	! Initialize max cfl values
604
605	grid%last_max_vert_cfl = 0
606	grid%last_max_horiz_cfl = 0
607
608	! Check to assure that time_step_sound is to be dynamically set.
609
610	CALL nl_set_time_step_sound ( 1 , 0 )
611	grid%time_step_sound = 0
612
613	grid%max_msftx=MAXVAL(grid%msftx)
614	grid%max_msfty=MAXVAL(grid%msfty)
615	#ifdef DM_PARALLEL
616	CALL wrf_dm_maxval(grid%max_msftx, idex, jdex)
617	CALL wrf_dm_maxval(grid%max_msfty, idex, jdex)
618	#endif
619
620	! This first call just initializes variables.
621	! If a restart, get initialized variables from restart file
622
623	IF ( .NOT. ( config_flags%restart ) ) then
624	CALL adapt_timestep(grid, config_flags)
625	END IF
626
627
628	END IF
629
630	! End of adaptive time step modifications
631	!-----------------------------------------------------------------------------
632
633
634	CALL set_tiles ( grid , grid%imask_nostag, ims, ime, jms, jme, ips, ipe, jps, jpe )
635	!
636	! Phy init can do reads and broadcasts when initializing physics -- landuse for example. However, if
637	! we're running on a reduced mesh (that is, some tasks don't have any work) we have to at least let them
638	! pass through this code so the broadcasts don't hang on the other, active tasks. Set the number of
639	! tiles to a minimum of 1 and assume that the backwards patch ranges (ips=0, ipe=-1) will prevent
640	! anything else from happening on the blank tasks. JM 20080605
641	!
642	if ( allowed_to_read ) grid%num_tiles = max(1,grid%num_tiles)
643	!
644	! Phy_init is not necessarily thread-safe; do not multi-thread this loop.
645	! The tiling is to handle the fact that we may be masking off part of the computation.
646	!
647	DO ij = 1, grid%num_tiles
648
649	!tgs do not need physics initialization for backward DFI integration
650	IF ( grid%dfi_opt .NE. DFI_NODFI .and. grid%dfi_stage .EQ. DFI_FST) THEN !tgs
651	grid%stepra=nint(grid%RADT*60./grid%DT)
652	grid%stepra=max(grid%stepra,1)
653	grid%stepbl=nint(grid%BLDT*60./grid%DT)
654	grid%stepbl=max(grid%stepbl,1)
655	grid%stepcu=nint(grid%CUDT*60./grid%DT)
656	grid%stepcu=max(grid%stepcu,1)
657	grid%stepfg=nint(grid%FGDT*60./grid%DT)
658	grid%stepfg=max(grid%stepfg,1)
659	ENDIF
660
661	IF ( ( grid%dfi_opt .EQ. DFI_NODFI ) .or. &
662	( ( grid%dfi_stage .NE. DFI_BCK ) .and. &
663	( grid%dfi_stage .NE. DFI_STARTBCK ) ) ) THEN !tgs, mods by tah
664
665	CALL phy_init ( grid%id , config_flags, grid%DT, grid%RESTART, grid%znw, grid%znu, &
666	grid%p_top, grid%tsk, grid%RADT,grid%BLDT,grid%CUDT, MPDT, &
667	grid%rucuten, grid%rvcuten, grid%rthcuten, &
668	grid%rqvcuten, grid%rqrcuten, grid%rqccuten, &
669	grid%rqscuten, grid%rqicuten, &
670	grid%rushten, grid%rvshten, grid%rthshten, &
671	grid%rqvshten, grid%rqrshten, grid%rqcshten, &
672	grid%rqsshten, grid%rqishten, grid%rqgshten, &
673	grid%rublten,grid%rvblten,grid%rthblten, &
674	grid%rqvblten,grid%rqcblten,grid%rqiblten, &
675	grid%rthraten,grid%rthratenlw,grid%rthratensw, &
676	grid%stepbl,grid%stepra,grid%stepcu, &
677	grid%w0avg, grid%rainnc, grid%rainc, grid%raincv, grid%rainncv, &
678	grid%snownc, grid%snowncv, grid%graupelnc, grid%graupelncv, &
679	grid%nca,grid%swrad_scat, &
680	grid%cldefi,grid%lowlyr, &
681	grid%mass_flux, &
682	grid%rthften, grid%rqvften, &
683	grid%cldfra, &
684	#ifdef WRF_CHEM
685	grid%cldfra_old, &
686	#endif
687	#ifndef WRF_CHEM
688	cldfra_old, &
689	#endif
690	grid%glw,grid%gsw,grid%emiss,grid%embck, &
691	grid%lu_index, &
692	grid%landuse_ISICE, grid%landuse_LUCATS, &
693	grid%landuse_LUSEAS, grid%landuse_ISN, &
694	grid%lu_state, &
695	grid%xlat,grid%xlong,grid%albedo,grid%albbck,grid%GMT,grid%JULYR,grid%JULDAY, &
696	grid%levsiz, num_ozmixm, num_aerosolc, grid%paerlev, &
697	grid%tmn,grid%xland,grid%znt,grid%z0,grid%ust,grid%mol,grid%pblh,grid%tke_pbl, &
698	grid%exch_h,grid%thc,grid%snowc,grid%mavail,grid%hfx,grid%qfx,grid%rainbl, &
699	grid%tslb,grid%zs,grid%dzs,config_flags%num_soil_layers,grid%warm_rain, &
700	grid%adv_moist_cond, &
701	grid%apr_gr,grid%apr_w,grid%apr_mc,grid%apr_st,grid%apr_as, &
702	grid%apr_capma,grid%apr_capme,grid%apr_capmi, &
703	grid%xice,grid%xicem,grid%vegfra,grid%snow,grid%canwat,grid%smstav, &
704	grid%smstot, grid%sfcrunoff,grid%udrunoff,grid%grdflx,grid%acsnow, &
705	grid%acsnom,grid%ivgtyp,grid%isltyp, grid%sfcevp,grid%smois, &
706	grid%sh2o, grid%snowh, grid%smfr3d, &
707	grid%snoalb, &
708	grid%DX,grid%DY,grid%f_ice_phy,grid%f_rain_phy,grid%f_rimef_phy, &
709	grid%mp_restart_state,grid%tbpvs_state,grid%tbpvs0_state,&
710	allowed_to_read, grid%moved, start_of_simulation, &
711	grid%LAGDAY, &
712	ids, ide, jds, jde, kds, kde, &
713	ims, ime, jms, jme, kms, kme, &
714	grid%i_start(ij), grid%i_end(ij), grid%j_start(ij), grid%j_end(ij), kts, kte, &
715	config_flags%num_urban_layers, & !multi-layer urban
716	ozmixm,grid%pin, & ! Optional
717	grid%m_ps_1,grid%m_ps_2,grid%m_hybi,aerosolc_1,aerosolc_2,& ! Optional
718	grid%rundgdten,grid%rvndgdten,grid%rthndgdten, & ! Optional
719	grid%rphndgdten,grid%rqvndgdten,grid%rmundgdten, & ! Optional
720	grid%FGDT,grid%stepfg, & ! Optional
721	grid%cugd_tten,grid%cugd_ttens,grid%cugd_qvten, & ! Optional
722	grid%cugd_qvtens,grid%cugd_qcten, & ! Optional
723	grid%DZR, grid%DZB, grid%DZG, & !Optional urban
724	grid%TR_URB2D,grid%TB_URB2D,grid%TG_URB2D,grid%TC_URB2D, & !Optional urban
725	grid%QC_URB2D, grid%XXXR_URB2D,grid%XXXB_URB2D, & !Optional urban
726	grid%XXXG_URB2D, grid%XXXC_URB2D, & !Optional urban
727	grid%TRL_URB3D, grid%TBL_URB3D, grid%TGL_URB3D, & !Optional urban
728	grid%SH_URB2D, grid%LH_URB2D, grid%G_URB2D, grid%RN_URB2D, & !Optional urban
729	grid%TS_URB2D, grid%FRC_URB2D, grid%UTYPE_URB2D, & !Optional urban
730	grid%TRB_URB4D,grid%TW1_URB4D,grid%TW2_URB4D,grid%TGB_URB4D,grid%TLEV_URB3D, & !multi-layer urban
731	grid%QLEV_URB3D,grid%TW1LEV_URB3D,grid%TW2LEV_URB3D, & !multi-layer urban
732	grid%TGLEV_URB3D,grid%TFLEV_URB3D,grid%SF_AC_URB3D, & !multi-layer urban
733	grid%LF_AC_URB3D,grid%CM_AC_URB3D,grid%SFVENT_URB3D,grid%LFVENT_URB3D, & !multi-layer urban
734	grid%SFWIN1_URB3D,grid%SFWIN2_URB3D, & !multi-layer urban
735	grid%SFW1_URB3D,grid%SFW2_URB3D,grid%SFR_URB3D,grid%SFG_URB3D, & !multi-layer urban
736	grid%A_U_BEP,grid%A_V_BEP,grid%A_T_BEP,grid%A_Q_BEP, & !multi-layer urban
737	grid%A_E_BEP,grid%B_U_BEP,grid%B_V_BEP,grid%B_T_BEP, & !multi-layer urban
738	grid%B_Q_BEP,grid%B_E_BEP,grid%DLG_BEP, & !multi-layer urban
739	grid%DL_U_BEP,grid%SF_BEP,grid%VL_BEP, & !multi-layer urban
740	grid%TML,grid%T0ML,grid%HML,grid%H0ML,grid%HUML,grid%HVML,grid%TMOML, & !Optional oml
741	grid%itimestep, grid%fdob, &
742	t00, p00, a, & ! for obs_nudge base state
743	grid%TYR, grid%TYRA, grid%TDLY, grid%TLAG, grid%NYEAR, grid%NDAY,grid%tmn_update, &
744	grid%achfx, grid%aclhf, grid%acgrdflx &
745	,grid%te_temf,grid%cf3d_temf,grid%wm_temf & ! WA
746	)
747	ENDIF !tgs
748
749	ENDDO
750
751	CALL wrf_debug ( 100 , 'start_domain_em: After call to phy_init' )
752
753	IF (config_flags%do_avgflx_em .EQ. 1) THEN
754	WRITE ( message , FMT = '("start_em: initializing avgflx on domain ",I3)' ) &
755	& grid%id
756	CALL wrf_message(trim(message))
757	grid%avgflx_count = 0
758	f_flux = config_flags%do_avgflx_cugd .EQ. 1 ! WA 9/24/10
759	DO ij = 1, grid%num_tiles
760	call wrf_debug(200,'In start_em, before zero_avgflx call')
761	if (.not. grid%restart) call zero_avgflx(grid%avgflx_rum,grid%avgflx_rvm,grid%avgflx_wwm, &
762	& ids, ide, jds, jde, kds, kde, &
763	& ims, ime, jms, jme, kms, kme, &
764	& grid%i_start(ij), grid%i_end(ij), grid%j_start(ij), grid%j_end(ij), kts, kte, f_flux, &
765	& grid%avgflx_cfu1,grid%avgflx_cfd1,grid%avgflx_dfu1, &
766	& grid%avgflx_efu1,grid%avgflx_dfd1,grid%avgflx_efd1 )
767	call wrf_debug(200,'In start_em, after zero_avgflx call')
768	ENDDO
769	ENDIF
770
771	#ifdef MCELIO
772	grid%LU_MASK = 0.
773	WHERE ( grid%lu_index .EQ. 16 ) grid%LU_MASK = 1.
774	#endif
775
776	END IF
777
778	#if 0
779	#include "CYCLE_TEST.inc"
780	#endif
781
782	!
783	!
784
785	! set physical boundary conditions for all initialized variables
786
787	!-----------------------------------------------------------------------
788	! Stencils for patch communications (WCS, 29 June 2001)
789	! Note: the size of this halo exchange reflects the
790	! fact that we are carrying the uncoupled variables
791	! as state variables in the mass coordinate model, as
792	! opposed to the coupled variables as in the height
793	! coordinate model.
794	!
795	! * * * * *
796	! * * * * * * * * *
797	! * + * * + * * * + * *
798	! * * * * * * * * *
799	! * * * * *
800	!
801	!j grid%u_1 x
802	!j grid%u_2 x
803	!j grid%v_1 x
804	!j grid%v_2 x
805	!j grid%w_1 x
806	!j grid%w_2 x
807	!j grid%t_1 x
808	!j grid%t_2 x
809	!j grid%ph_1 x
810	!j grid%ph_2 x
811	!
812	!j grid%t_init x
813	!
814	!j grid%phb x
815	!j grid%ph0 x
816	!j grid%php x
817	!j grid%pb x
818	!j grid%al x
819	!j grid%alt x
820	!j grid%alb x
821	!
822	! the following are 2D (xy) variables
823	!
824	!j grid%mu_1 x
825	!j grid%mu_2 x
826	!j grid%mub x
827	!j grid%mu0 x
828	!j grid%ht x
829	!j grid%msftx x
830	!j grid%msfty x
831	!j grid%msfux x
832	!j grid%msfuy x
833	!j grid%msfvx x
834	!j grid%msfvy x
835	!j grid%sina x
836	!j grid%cosa x
837	!j grid%e x
838	!j grid%f x
839	!
840	! 4D variables
841	!
842	! moist x
843	! chem x
844	!scalar x
845
846	!--------------------------------------------------------------
847
848	#ifdef DM_PARALLEL
849	# include "HALO_EM_INIT_1.inc"
850	# include "HALO_EM_INIT_2.inc"
851	# include "HALO_EM_INIT_3.inc"
852	# include "HALO_EM_INIT_4.inc"
853	# include "HALO_EM_INIT_5.inc"
854	# include "PERIOD_BDY_EM_INIT.inc"
855	# include "PERIOD_BDY_EM_MOIST.inc"
856	# include "PERIOD_BDY_EM_TKE.inc"
857	# include "PERIOD_BDY_EM_SCALAR.inc"
858	# include "PERIOD_BDY_EM_CHEM.inc"
859	#endif
860
861
862	CALL set_physical_bc3d( grid%u_1 , 'U' , config_flags , &
863	ids , ide , jds , jde , kds , kde , &
864	ims , ime , jms , jme , kms , kme , &
865	its , ite , jts , jte , kts , kte , &
866	its , ite , jts , jte , kts , kte )
867	CALL set_physical_bc3d( grid%u_2 , 'U' , config_flags , &
868	ids , ide , jds , jde , kds , kde , &
869	ims , ime , jms , jme , kms , kme , &
870	its , ite , jts , jte , kts , kte , &
871	its , ite , jts , jte , kts , kte )
872
873	CALL set_physical_bc3d( grid%v_1 , 'V' , config_flags , &
874	ids , ide , jds , jde , kds , kde , &
875	ims , ime , jms , jme , kms , kme , &
876	its , ite , jts , jte , kts , kte , &
877	its , ite , jts , jte , kts , kte )
878	CALL set_physical_bc3d( grid%v_2 , 'V' , config_flags , &
879	ids , ide , jds , jde , kds , kde , &
880	ims , ime , jms , jme , kms , kme , &
881	its , ite , jts , jte , kts , kte , &
882	its , ite , jts , jte , kts , kte )
883
884	! set kinematic condition for w
885
886	CALL set_physical_bc2d( grid%ht , 'r' , config_flags , &
887	ids , ide , jds , jde , &
888	ims , ime , jms , jme , &
889	its , ite , jts , jte , &
890	its , ite , jts , jte )
891
892	! Set the w at the surface. If this is the start of a forecast, or if this is a cycled run
893	! and the start of that forecast, we define the w field. However, if this a restart, then
894	! we leave w alone. Initial value is V dot grad(topo) at surface, then smoothly decaying
895	! above that.
896
897	IF ( ( start_of_simulation .OR. config_flags%cycling ) .AND. ( .NOT. config_flags%restart ) ) THEN
898	fill_w_flag = .true.
899	CALL set_w_surface( config_flags, grid%znw, fill_w_flag, &
900	grid%w_1, grid%ht, grid%u_1, grid%v_1, grid%cf1, &
901	grid%cf2, grid%cf3, grid%rdx, grid%rdy, grid%msftx, grid%msfty, &
902	ids, ide, jds, jde, kds, kde, &
903	ims, ime, jms, jme, kms, kme, &
904	its, ite, jts, jte, kts, kte )
905	CALL set_w_surface( config_flags, grid%znw, fill_w_flag, &
906	grid%w_2, grid%ht, grid%u_2, grid%v_2, grid%cf1, &
907	grid%cf2, grid%cf3, grid%rdx, grid%rdy, grid%msftx, grid%msfty, &
908	ids, ide, jds, jde, kds, kde, &
909	ims, ime, jms, jme, kms, kme, &
910	its, ite, jts, jte, kts, kte )
911
912	! finished setting kinematic condition for w at the surface
913
914	! set up slope-radiation constant arrays based on topography
915	DO j = jts,min(jte,jde-1)
916	DO i = its, min(ite,ide-1)
917	im1 = max(i-1,ids)
918	ip1 = min(i+1,ide-1)
919	jm1 = max(j-1,jds)
920	jp1 = min(j+1,jde-1)
921	grid%toposlpx(i,j)=(grid%ht(ip1,j)-grid%ht(im1,j))grid%msftx(i,j)grid%rdx/(ip1-im1)
922	grid%toposlpy(i,j)=(grid%ht(i,jp1)-grid%ht(i,jm1))grid%msfty(i,j)grid%rdy/(jp1-jm1)
923	hx = grid%toposlpx(i,j)
924	hy = grid%toposlpy(i,j)
925	pi = 4.*atan(1.)
926	grid%slope(i,j) = atan((hx2+hy2)**.5)
927	if (grid%slope(i,j).lt.1.e-4) then
928	grid%slope(i,j) = 0.
929	grid%slp_azi(i,j) = 0.
930	else
931	grid%slp_azi(i,j) = atan2(hx,hy)+pi
932
933	! Rotate slope azimuth to lat-lon grid
934	if (grid%cosa(i,j).ge.0) then
935	grid%slp_azi(i,j) = grid%slp_azi(i,j) - asin(grid%sina(i,j))
936	else
937	grid%slp_azi(i,j) = grid%slp_azi(i,j) - (pi - asin(grid%sina(i,j)))
938	endif
939	endif
940	ENDDO
941	ENDDO
942	END IF
943
944	CALL set_physical_bc3d( grid%w_1 , 'W' , config_flags , &
945	ids , ide , jds , jde , kds , kde , &
946	ims , ime , jms , jme , kms , kme , &
947	its , ite , jts , jte , kts , kte , &
948	its , ite , jts , jte , kts , kte )
949	CALL set_physical_bc3d( grid%w_2 , 'W' , config_flags , &
950	ids , ide , jds , jde , kds , kde , &
951	ims , ime , jms , jme , kms , kme , &
952	its , ite , jts , jte , kts , kte , &
953	its , ite , jts , jte , kts , kte )
954
955	CALL set_physical_bc3d( grid%ph_1 , 'W' , config_flags , &
956	ids , ide , jds , jde , kds , kde , &
957	ims , ime , jms , jme , kms , kme , &
958	its , ite , jts , jte , kts , kte , &
959	its , ite , jts , jte , kts , kte )
960
961	CALL set_physical_bc3d( grid%ph_2 , 'W' , config_flags , &
962	ids , ide , jds , jde , kds , kde , &
963	ims , ime , jms , jme , kms , kme , &
964	its , ite , jts , jte , kts , kte , &
965	its , ite , jts , jte , kts , kte )
966
967	CALL set_physical_bc3d( grid%t_1 , 't' , config_flags , &
968	ids , ide , jds , jde , kds , kde , &
969	ims , ime , jms , jme , kms , kme , &
970	its , ite , jts , jte , kts , kte , &
971	its , ite , jts , jte , kts , kte )
972
973	CALL set_physical_bc3d( grid%t_2 , 't' , config_flags , &
974	ids , ide , jds , jde , kds , kde , &
975	ims , ime , jms , jme , kms , kme , &
976	its , ite , jts , jte , kts , kte , &
977	its , ite , jts , jte , kts , kte )
978
979	CALL set_physical_bc2d( grid%mu_1, 't' , config_flags , &
980	ids , ide , jds , jde , &
981	ims , ime , jms , jme , &
982	its , ite , jts , jte , &
983	its , ite , jts , jte )
984	CALL set_physical_bc2d( grid%mu_2, 't' , config_flags , &
985	ids , ide , jds , jde , &
986	ims , ime , jms , jme , &
987	its , ite , jts , jte , &
988	its , ite , jts , jte )
989	CALL set_physical_bc2d( grid%mub , 't' , config_flags , &
990	ids , ide , jds , jde , &
991	ims , ime , jms , jme , &
992	its , ite , jts , jte , &
993	its , ite , jts , jte )
994	! CALL set_physical_bc2d( grid%mu0 , 't' , config_flags , &
995	! ids , ide , jds , jde , &
996	! ims , ime , jms , jme , &
997	! its , ite , jts , jte , &
998	! its , ite , jts , jte )
999
1000
1001	CALL set_physical_bc3d( grid%phb , 'W' , config_flags , &
1002	ids , ide , jds , jde , kds , kde , &
1003	ims , ime , jms , jme , kms , kme , &
1004	its , ite , jts , jte , kts , kte , &
1005	its , ite , jts , jte , kts , kte )
1006	CALL set_physical_bc3d( grid%ph0 , 'W' , config_flags , &
1007	ids , ide , jds , jde , kds , kde , &
1008	ims , ime , jms , jme , kms , kme , &
1009	its , ite , jts , jte , kts , kte , &
1010	its , ite , jts , jte , kts , kte )
1011	CALL set_physical_bc3d( grid%php , 'W' , config_flags , &
1012	ids , ide , jds , jde , kds , kde , &
1013	ims , ime , jms , jme , kms , kme , &
1014	its , ite , jts , jte , kts , kte , &
1015	its , ite , jts , jte , kts , kte )
1016
1017	CALL set_physical_bc3d( grid%pb , 't' , config_flags , &
1018	ids , ide , jds , jde , kds , kde , &
1019	ims , ime , jms , jme , kms , kme , &
1020	its , ite , jts , jte , kts , kte , &
1021	its , ite , jts , jte , kts , kte )
1022	CALL set_physical_bc3d( grid%al , 't' , config_flags , &
1023	ids , ide , jds , jde , kds , kde , &
1024	ims , ime , jms , jme , kms , kme , &
1025	its , ite , jts , jte , kts , kte , &
1026	its , ite , jts , jte , kts , kte )
1027	CALL set_physical_bc3d( grid%alt , 't' , config_flags , &
1028	ids , ide , jds , jde , kds , kde , &
1029	ims , ime , jms , jme , kms , kme , &
1030	its , ite , jts , jte , kts , kte , &
1031	its , ite , jts , jte , kts , kte )
1032	CALL set_physical_bc3d( grid%alb , 't' , config_flags , &
1033	ids , ide , jds , jde , kds , kde , &
1034	ims , ime , jms , jme , kms , kme , &
1035	its , ite , jts , jte , kts , kte , &
1036	its , ite , jts , jte , kts , kte )
1037	CALL set_physical_bc3d(grid%t_init, 't' , config_flags , &
1038	ids , ide , jds , jde , kds , kde , &
1039	ims , ime , jms , jme , kms , kme , &
1040	its , ite , jts , jte , kts , kte , &
1041	its , ite , jts , jte , kts , kte )
1042	CALL set_physical_bc3d(grid%tke_2, 't' , config_flags , &
1043	ids , ide , jds , jde , kds , kde , &
1044	ims , ime , jms , jme , kms , kme , &
1045	its , ite , jts , jte , kts , kte , &
1046	its , ite , jts , jte , kts , kte )
1047
1048	IF (num_moist > 0) THEN
1049
1050	! use of (:,:,:,loop) not efficient on DEC, but (ims,kms,jms,loop) not portable to SGI/Cray
1051
1052	loop_3d_m : DO loop = 1 , num_moist
1053	CALL set_physical_bc3d( moist(:,:,:,loop) , 'r' , config_flags , &
1054	ids , ide , jds , jde , kds , kde , &
1055	ims , ime , jms , jme , kms , kme , &
1056	its , ite , jts , jte , kts , kte , &
1057	its , ite , jts , jte , kts , kte )
1058	END DO loop_3d_m
1059
1060	ENDIF
1061
1062	IF (num_scalar > 0) THEN
1063
1064	! use of (:,:,:,loop) not efficient on DEC, but (ims,kms,jms,loop) not portable to SGI/Cray
1065
1066	loop_3d_s : DO loop = 1 , num_scalar
1067	CALL set_physical_bc3d( scalar(:,:,:,loop) , 'r' , config_flags , &
1068	ids , ide , jds , jde , kds , kde , &
1069	ims , ime , jms , jme , kms , kme , &
1070	its , ite , jts , jte , kts , kte , &
1071	its , ite , jts , jte , kts , kte )
1072	END DO loop_3d_s
1073
1074	ENDIF
1075
1076
1077	#ifdef WRF_CHEM
1078	if(config_flags%tracer_opt > 0 )then
1079	call initialize_tracer (tracer,config_flags%chem_in_opt, &
1080	config_flags%tracer_opt,num_tracer, &
1081	ids,ide, jds,jde, kds,kde, & ! domain dims
1082	ims,ime, jms,jme, kms,kme, & ! memory dims
1083	ips,ipe, jps,jpe, kps,kpe, & ! patch dims
1084	its,ite, jts,jte, kts,kte )
1085
1086	endif
1087	!
1088	! we do this here, so we only have one chem_init routine for either core....
1089	!
1090	do j=jts,min(jte,jde-1)
1091	do i=its,min(ite,ide-1)
1092	do k=kts,kte
1093	z_at_w(i,k,j)=(grid%ph_2(i,k,j)+grid%phb(i,k,j))/g
1094	enddo
1095	do k=kts,min(kte,kde-1)
1096	tempfac=(grid%t_1(i,k,j) + t0)((grid%p(i,k,j) + grid%pb(i,k,j))/p1000mb)*rcp
1097	convfac(i,k,j) = (grid%p(i,k,j)+grid%pb(i,k,j))/rgasuniv/tempfac
1098	enddo
1099	enddo
1100	enddo
1101
1102	CALL chem_init (grid%id,chem,emis_ant,scalar,grid%dt,grid%bioemdt,grid%photdt, &
1103	grid%chemdt, &
1104	grid%stepbioe,grid%stepphot,grid%stepchem,grid%stepfirepl, &
1105	grid%plumerisefire_frq,z_at_w,grid%xlat,grid%xlong,g, &
1106	grid%aerwrf,config_flags,grid, &
1107	grid%alt,grid%t_1,grid%p,convfac,grid%ttday,grid%tcosz, &
1108	grid%julday,grid%gmt,&
1109	grid%gd_cloud, grid%gd_cloud2,grid%raincv_a,grid%raincv_b, &
1110	grid%gd_cloud_a, grid%gd_cloud2_a, &
1111	grid%gd_cloud_b, grid%gd_cloud2_b, &
1112	grid%tauaer1,grid%tauaer2,grid%tauaer3,grid%tauaer4, &
1113	grid%gaer1,grid%gaer2,grid%gaer3,grid%gaer4, &
1114	grid%waer1,grid%waer2,grid%waer3,grid%waer4, &
1115	grid%l2aer,grid%l3aer,grid%l4aer,grid%l5aer,grid%l6aer,grid%l7aer, &
1116	grid%extaerlw1,grid%extaerlw2,grid%extaerlw3,grid%extaerlw4, &
1117	grid%extaerlw5,grid%extaerlw6,grid%extaerlw7,grid%extaerlw8, &
1118	grid%extaerlw9,grid%extaerlw10,grid%extaerlw11,grid%extaerlw12, &
1119	grid%extaerlw13,grid%extaerlw14,grid%extaerlw15,grid%extaerlw16, &
1120	grid%tauaerlw1,grid%tauaerlw2,grid%tauaerlw3,grid%tauaerlw4, &
1121	grid%tauaerlw5,grid%tauaerlw6,grid%tauaerlw7,grid%tauaerlw8, &
1122	grid%tauaerlw9,grid%tauaerlw10,grid%tauaerlw11,grid%tauaerlw12, &
1123	grid%tauaerlw13,grid%tauaerlw14,grid%tauaerlw15,grid%tauaerlw16, &
1124	grid%pm2_5_dry,grid%pm2_5_water,grid%pm2_5_dry_ec, &
1125	grid%last_chem_time_year,grid%last_chem_time_month, &
1126	grid%last_chem_time_day,grid%last_chem_time_hour, &
1127	grid%last_chem_time_minute,grid%last_chem_time_second, &
1128	grid%chem_in_opt,grid%kemit, &
1129	ids , ide , jds , jde , kds , kde , &
1130	ims , ime , jms , jme , kms , kme , &
1131	its , ite , jts , jte , kts , kte )
1132
1133	!
1134	! calculate initial pm
1135	!
1136	! print *,'calculating initial pm'
1137	select case (config_flags%chem_opt)
1138	case (GOCART_SIMPLE,GOCARTRACM_KPP,GOCARTRADM2,GOCARTRADM2_KPP)
1139	call sum_pm_gocart ( &
1140	grid%alt, chem, grid%pm2_5_dry, grid%pm2_5_dry_ec,grid%pm10,&
1141	ids,ide, jds,jde, kds,kde, &
1142	ims,ime, jms,jme, kms,kme, &
1143	its,ite, jts,jte, kts,kte-1 )
1144	case (RADM2SORG,RADM2SORG_AQ,RADM2SORG_KPP,RACMSORG_AQ,RACMSORG_KPP)
1145	call sum_pm_sorgam ( &
1146	grid%alt, chem, grid%h2oaj, grid%h2oai, &
1147	grid%pm2_5_dry, grid%pm2_5_water, grid%pm2_5_dry_ec, grid%pm10, &
1148	grid%dust_opt,ids,ide, jds,jde, kds,kde, &
1149	ims,ime, jms,jme, kms,kme, &
1150	its,ite, jts,jte, kts,kte-1 )
1151
1152	case (CBMZ_MOSAIC_4BIN,CBMZ_MOSAIC_8BIN,CBMZ_MOSAIC_4BIN_AQ,CBMZ_MOSAIC_8BIN_AQ, &
1153	CBMZ_MOSAIC_DMS_4BIN, CBMZ_MOSAIC_DMS_8BIN, CBMZ_MOSAIC_DMS_4BIN_AQ, &
1154	CBMZ_MOSAIC_DMS_8BIN_AQ)
1155	call sum_pm_mosaic ( &
1156	grid%alt, chem, &
1157	grid%pm2_5_dry, grid%pm2_5_water, grid%pm2_5_dry_ec, grid%pm10, &
1158	ids,ide, jds,jde, kds,kde, &
1159	ims,ime, jms,jme, kms,kme, &
1160	its,ite, jts,jte, kts,kte-1 )
1161
1162	case default
1163	do j=jts,min(jte,jde-1)
1164	do k=kts,min(kte,kde-1)
1165	do i=its,min(ite,ide-1)
1166	grid%pm2_5_dry(i,k,j) = 0.
1167	grid%pm2_5_water(i,k,j) = 0.
1168	grid%pm2_5_dry_ec(i,k,j) = 0.
1169	grid%pm10(i,k,j) = 0.
1170	enddo
1171	enddo
1172	enddo
1173	end select
1174
1175	! initialize advective tendency diagnostics for chem
1176	if ( grid%itimestep .eq. 0 .and. config_flags%chemdiag .eq. USECHEMDIAG ) then
1177	grid%advh_ct(:,:,:,:) = 0.
1178	grid%advz_ct(:,:,:,:) = 0.
1179	endif
1180	#endif
1181
1182	! initialize advective tendency diagnostics for non-chem
1183	if ( grid%itimestep .eq. 0 .and. config_flags%tenddiag .eq. USETENDDIAG ) then
1184	grid%advh_t(:,:,:,:) = 0.
1185	grid%advz_t(:,:,:,:) = 0.
1186	endif
1187
1188	IF (num_chem >= PARAM_FIRST_SCALAR ) THEN
1189	! use of (:,:,:,loop) not efficient on DEC, but (ims,kms,jms,loop) not portable to SGI/Cray
1190
1191	loop_3d_c : DO loop = PARAM_FIRST_SCALAR , num_chem
1192	CALL set_physical_bc3d( chem(:,:,:,loop) , 'r' , config_flags , &
1193	ids , ide , jds , jde , kds , kde , &
1194	ims , ime , jms , jme , kms , kme , &
1195	its , ite , jts , jte , kts , kte , &
1196	its , ite , jts , jte , kts , kte )
1197	END DO loop_3d_c
1198
1199	ENDIF
1200
1201	CALL set_physical_bc2d( grid%msftx , 'r' , config_flags , &
1202	ids , ide , jds , jde , &
1203	ims , ime , jms , jme , &
1204	its , ite , jts , jte , &
1205	its , ite , jts , jte )
1206	CALL set_physical_bc2d( grid%msfty , 'r' , config_flags , &
1207	ids , ide , jds , jde , &
1208	ims , ime , jms , jme , &
1209	its , ite , jts , jte , &
1210	its , ite , jts , jte )
1211	CALL set_physical_bc2d( grid%msfux , 'x' , config_flags , &
1212	ids , ide , jds , jde , &
1213	ims , ime , jms , jme , &
1214	its , ite , jts , jte , &
1215	its , ite , jts , jte )
1216	CALL set_physical_bc2d( grid%msfuy , 'x' , config_flags , &
1217	ids , ide , jds , jde , &
1218	ims , ime , jms , jme , &
1219	its , ite , jts , jte , &
1220	its , ite , jts , jte )
1221	CALL set_physical_bc2d( grid%msfvx , 'y' , config_flags , &
1222	ids , ide , jds , jde , &
1223	ims , ime , jms , jme , &
1224	its , ite , jts , jte , &
1225	its , ite , jts , jte )
1226	CALL set_physical_bc2d( grid%msfvy , 'y' , config_flags , &
1227	ids , ide , jds , jde , &
1228	ims , ime , jms , jme , &
1229	its , ite , jts , jte , &
1230	its , ite , jts , jte )
1231	CALL set_physical_bc2d( grid%sina , 'r' , config_flags , &
1232	ids , ide , jds , jde , &
1233	ims , ime , jms , jme , &
1234	its , ite , jts , jte , &
1235	its , ite , jts , jte )
1236	CALL set_physical_bc2d( grid%cosa , 'r' , config_flags , &
1237	ids , ide , jds , jde , &
1238	ims , ime , jms , jme , &
1239	its , ite , jts , jte , &
1240	its , ite , jts , jte )
1241	CALL set_physical_bc2d( grid%e , 'r' , config_flags , &
1242	ids , ide , jds , jde , &
1243	ims , ime , jms , jme , &
1244	its , ite , jts , jte , &
1245	its , ite , jts , jte )
1246	CALL set_physical_bc2d( grid%f , 'r' , config_flags , &
1247	ids , ide , jds , jde , &
1248	ims , ime , jms , jme , &
1249	its , ite , jts , jte , &
1250	its , ite , jts , jte )
1251
1252	#ifndef WRF_CHEM
1253	DEALLOCATE(CLDFRA_OLD)
1254	#endif
1255	#ifdef DM_PARALLEL
1256	# include "HALO_EM_INIT_1.inc"
1257	# include "HALO_EM_INIT_2.inc"
1258	# include "HALO_EM_INIT_3.inc"
1259	# include "HALO_EM_INIT_4.inc"
1260	# include "HALO_EM_INIT_5.inc"
1261	# include "PERIOD_BDY_EM_INIT.inc"
1262	# include "PERIOD_BDY_EM_MOIST.inc"
1263	# include "PERIOD_BDY_EM_TKE.inc"
1264	# include "PERIOD_BDY_EM_SCALAR.inc"
1265	# include "PERIOD_BDY_EM_CHEM.inc"
1266	#endif
1267
1268	! FIRE
1269	if(config_flags%ifire.eq.2)then
1270
1271	call sfire_driver_em_init ( grid , config_flags &
1272	,ids,ide, kds,kde, jds,jde &
1273	,ims,ime, kms,kme, jms,jme &
1274	,ips,ipe, kps,kpe, jps,jpe )
1275
1276	CALL wrf_debug ( 100 , 'start_domain_em: After call to sfire_driver_em_init' )
1277	endif
1278
1279
1280	CALL wrf_debug ( 100 , 'start_domain_em: Returning' )
1281
1282	RETURN
1283
1284	END SUBROUTINE start_domain_em
1285

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format