Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_initialize_fire.F @ 1

Last change on this file since 1 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: 37.0 KB

Line
1	!IDEAL:MODEL_LAYER:INITIALIZATION
2	!
3
4	! This MODULE holds the routines which are used to perform various initializations
5	! for the individual domains.
6
7	! This MODULE CONTAINS the following routines:
8
9	! initialize_field_test - 1. Set different fields to different constant
10	! values. This is only a test. If the correct
11	! domain is not found (based upon the "id")
12	! then a fatal error is issued.
13
14	!-----------------------------------------------------------------------
15
16	MODULE module_initialize_ideal
17
18	USE module_domain
19	USE module_io_domain
20	USE module_state_description
21	USE module_model_constants
22	USE module_bc
23	USE module_timing
24	USE module_configure
25	USE module_init_utilities
26	USE module_soil_pre !AK/ak for full surface initialization
27	#ifdef DM_PARALLEL
28	USE module_dm
29	#endif
30	USE module_fr_sfire_util, ONLY: continue_at_boundary,crash,read_array_2d_real, &
31	read_array_2d_integer,interpolate_2d,set_ideal_coord
32
33	CONTAINS
34
35
36	!-------------------------------------------------------------------
37	! this is a wrapper for the solver-specific init_domain routines.
38	! Also dereferences the grid variables and passes them down as arguments.
39	! This is crucial, since the lower level routines may do message passing
40	! and this will get fouled up on machines that insist on passing down
41	! copies of assumed-shape arrays (by passing down as arguments, the
42	! data are treated as assumed-size -- ie. f77 -- arrays and the copying
43	! business is avoided). Fie on the F90 designers. Fie and a pox.
44
45	SUBROUTINE init_domain ( grid )
46
47	IMPLICIT NONE
48
49	! Input data.
50	TYPE (domain), POINTER :: grid
51	! Local data.
52	INTEGER :: idum1, idum2
53
54	CALL set_scalar_indices_from_config ( head_grid%id , idum1, idum2 )
55
56	CALL init_domain_rk( grid &
57	!
58	#include <actual_new_args.inc>
59	!
60	)
61
62	END SUBROUTINE init_domain
63
64	!-------------------------------------------------------------------
65
66	SUBROUTINE init_domain_rk ( grid &
67	!
68	# include <dummy_new_args.inc>
69	!
70	)
71	IMPLICIT NONE
72
73	! Input data.
74	TYPE (domain), POINTER :: grid
75
76	# include <dummy_new_decl.inc>
77
78	TYPE (grid_config_rec_type) :: config_flags
79
80	LOGICAL, EXTERNAL :: wrf_dm_on_monitor
81
82	! Local data
83	INTEGER :: &
84	ids, ide, jds, jde, kds, kde, &
85	ims, ime, jms, jme, kms, kme, &
86	its, ite, jts, jte, kts, kte, &
87	i, j, k
88
89	INTEGER, PARAMETER :: nl_max = 1000
90	REAL, DIMENSION(nl_max) :: zk, p_in, theta, rho, u, v, qv, pd_in
91	INTEGER :: nl_in
92
93
94	INTEGER :: icm,jcm, ii, im1, jj, jm1, loop, error, fid, nxc, nyc
95	REAL :: u_mean,v_mean, f0, p_surf, p_level, qvf, z_at_v, z_at_u
96	REAL :: z_scale, xrad, yrad, zrad, rad, delt, cof1, cof2
97	REAL :: x_rad, y_rad, z_rad, hght_pert !Ak/ak
98	character (len=256) :: mminlu2 !AK/ak land use scheme (USGS)
99	! REAL, EXTERNAL :: interp_0
100	REAL :: hm
101	REAL :: pi
102
103	! stuff from original initialization that has been dropped from the Registry
104	REAL :: vnu, xnu, xnus, dinit0, cbh, p0_temp, t0_temp, zd, zt
105	REAL :: qvf1, qvf2, pd_surf
106	INTEGER :: it
107	real :: thtmp, ptmp, temp(3)
108
109	LOGICAL :: moisture_init
110	LOGICAL :: stretch_grd, dry_sounding
111	LOGICAL :: stretch_hyp, sfc_init !AK/ak switches for hyperbolic grid streching and surface initialization
112
113	INTEGER :: xs , xe , ys , ye
114	INTEGER :: mtn_type
115	INTEGER :: & ! fire mesh sizes
116	iots,iote,jots,jote, & ! tile dims out
117	ifds,ifde, kfds,kfde, jfds,jfde, &
118	ifms,ifme, kfms,kfme, jfms,jfme, &
119	ifts,ifte, kfts,kfte, jfts,jfte
120
121	REAL :: mtn_ht, mtn_max, mtn_x, mtn_y, mtn_z, grad_max
122	REAL :: mtn_axs, mtn_ays, mtn_axe, mtn_aye
123	REAL :: mtn_fxs, mtn_fys, mtn_fxe, mtn_fye
124	REAL :: mtn_xs, mtn_ys, mtn_xe, mtn_ye
125	REAL :: fdx,fdy ! fire mesh step
126	INTEGER:: ir,jr ! refinement factors
127
128	logical have_fire_ht,have_fire_grad,have_atm_grad
129
130	!*** executable
131
132	SELECT CASE ( model_data_order )
133	CASE ( DATA_ORDER_ZXY )
134	kds = grid%sd31 ; kde = grid%ed31 ;
135	ids = grid%sd32 ; ide = grid%ed32 ;
136	jds = grid%sd33 ; jde = grid%ed33 ;
137
138	kms = grid%sm31 ; kme = grid%em31 ;
139	ims = grid%sm32 ; ime = grid%em32 ;
140	jms = grid%sm33 ; jme = grid%em33 ;
141
142	kts = grid%sp31 ; kte = grid%ep31 ; ! note that tile is entire patch
143	its = grid%sp32 ; ite = grid%ep32 ; ! note that tile is entire patch
144	jts = grid%sp33 ; jte = grid%ep33 ; ! note that tile is entire patch
145	CASE ( DATA_ORDER_XYZ )
146	ids = grid%sd31 ; ide = grid%ed31 ;
147	jds = grid%sd32 ; jde = grid%ed32 ;
148	kds = grid%sd33 ; kde = grid%ed33 ;
149
150	ims = grid%sm31 ; ime = grid%em31 ;
151	jms = grid%sm32 ; jme = grid%em32 ;
152	kms = grid%sm33 ; kme = grid%em33 ;
153
154	its = grid%sp31 ; ite = grid%ep31 ; ! note that tile is entire patch
155	jts = grid%sp32 ; jte = grid%ep32 ; ! note that tile is entire patch
156	kts = grid%sp33 ; kte = grid%ep33 ; ! note that tile is entire patch
157	CASE ( DATA_ORDER_XZY )
158	ids = grid%sd31 ; ide = grid%ed31 ;
159	kds = grid%sd32 ; kde = grid%ed32 ;
160	jds = grid%sd33 ; jde = grid%ed33 ;
161
162	ims = grid%sm31 ; ime = grid%em31 ;
163	kms = grid%sm32 ; kme = grid%em32 ;
164	jms = grid%sm33 ; jme = grid%em33 ;
165
166	its = grid%sp31 ; ite = grid%ep31 ; ! note that tile is entire patch
167	kts = grid%sp32 ; kte = grid%ep32 ; ! note that tile is entire patch
168	jts = grid%sp33 ; jte = grid%ep33 ; ! note that tile is entire patch
169
170	END SELECT
171
172	! z_scale = .40
173	pi = 2.*asin(1.0)
174	write(6,*) ' pi is ',pi
175	nxc = (ide-ids)/2
176	nyc = (jde-jds)/2
177
178	CALL model_to_grid_config_rec ( grid%id , model_config_rec , config_flags )
179
180	! here we check to see if the boundary conditions are set properly
181
182	CALL boundary_condition_check( config_flags, bdyzone, error, grid%id )
183
184	delt = config_flags%delt_perturbation
185	x_rad = config_flags%xrad_perturbation
186	y_rad = config_flags%yrad_perturbation
187	z_rad = config_flags%zrad_perturbation
188	hght_pert = config_flags%hght_perturbation
189
190	stretch_grd = config_flags%stretch_grd
191	stretch_hyp = config_flags%stretch_hyp
192	z_scale = config_flags%z_grd_scale
193	sfc_init = config_flags%sfc_full_init
194
195	moisture_init = .true. !AK/ak
196
197	grid%itimestep=0
198
199	#ifdef DM_PARALLEL
200	CALL wrf_dm_bcast_bytes( icm , IWORDSIZE )
201	CALL wrf_dm_bcast_bytes( jcm , IWORDSIZE )
202	#endif
203
204	!AK/ak land use initialization (USGS)
205	IF (sfc_init) THEN
206	mminlu2=' '
207	mminlu2(1:4) = 'USGS' !Ak/ak
208	CALL nl_set_mminlu(1, mminlu2) !Ak/ak
209	CALL nl_get_iswater(1,grid%iswater) ! Ak/ak
210	ENDIF
211
212	CALL nl_set_iswater(1,0)
213	CALL nl_set_cen_lat(1,40.)
214	CALL nl_set_cen_lon(1,-105.)
215	CALL nl_set_truelat1(1,0.)
216	CALL nl_set_truelat2(1,0.)
217	CALL nl_set_moad_cen_lat (1,0.)
218	CALL nl_set_stand_lon (1,0.)
219	CALL nl_set_pole_lon (1,0.)
220	CALL nl_set_pole_lat (1,90.)
221	CALL nl_set_map_proj(1,0)
222
223	! here we initialize data we currently is not initialized
224	! in the input data
225
226	DO j = jts, jte
227	DO i = its, ite
228	grid%msftx(i,j) = 1.
229	grid%msfty(i,j) = 1.
230	grid%msfux(i,j) = 1.
231	grid%msfuy(i,j) = 1.
232	grid%msfvx(i,j) = 1.
233	grid%msfvx_inv(i,j)= 1.
234	grid%msfvy(i,j) = 1.
235	grid%sina(i,j) = 0.
236	grid%cosa(i,j) = 1.
237	grid%e(i,j) = 0.
238	grid%f(i,j) = 0.
239	END DO
240	END DO
241
242	!AK/ak surface initialization latitude, longitude, landuse index from from LANDUSE.TBL skin temperature and soil temperature
243	write(6,) '************************************'
244	IF (sfc_init) THEN
245	DO j = jts, jte
246	DO i = its, ite
247	grid%xlat(i,j) = config_flags%fire_lat_init !Ak/sk (35)
248	grid%xlong(i,j) = config_flags%fire_lon_init !Ak/ak (-111)
249	grid%xland(i,j) = 1. !Ak/ak
250	grid%lu_index(i,j) = config_flags%sfc_lu_index !AK/ak land use index (28)
251	grid%tsk(i,j) = config_flags%sfc_tsk !AK/ak surface skin temperature [K] (285)
252	grid%tmn(i,j) = config_flags%sfc_tmn !AK/ak soil temperature at lower boundary [K] (285)
253	END DO
254	END DO
255	! read land use data from files, overwriting the constant
256	if(config_flags%fire_read_lu) &
257	call read_array_2d_real('input_lu',grid%lu_index,ids,ide,jds,jde,ims,ime,jms,jme)
258	if(config_flags%fire_read_tsk) &
259	call read_array_2d_real ('input_tsk',grid%tsk, ids,ide,jds,jde,ims,ime,jms,jme)
260	if(config_flags%fire_read_tmn) &
261	call read_array_2d_real ('input_tmn',grid%tmn, ids,ide,jds,jde,ims,ime,jms,jme)
262
263	! for Noah LSM, additional variables need to be initializedi !AK/ak \|
264
265	other_masked_fields : SELECT CASE ( model_config_rec%sf_surface_physics(grid%id) )
266
267	CASE (SLABSCHEME)
268	write(6,*) ' SLAB surface scheme activated'
269
270	CASE (LSMSCHEME)
271	write(6,*) ' Noah unified LSM scheme activated with:'
272	write(6,*) ' vegetation fraction=',config_flags%sfc_vegfra
273	write(6,*) ' canopy water=',config_flags%sfc_canwat
274	write(6,*) ' dominant veg. type=',config_flags%sfc_ivgtyp
275	write(6,*) ' dominant soil type=',config_flags%sfc_isltyp
276
277	DO j = jts , MIN(jde-1,jte)
278	DO i = its , MIN(ide-1,ite)
279	grid%vegfra(i,j) = config_flags%sfc_vegfra !0.5
280	grid%canwat(i,j) = config_flags%sfc_canwat !0.
281	grid%ivgtyp(i,j) = config_flags%sfc_ivgtyp !18
282	grid%isltyp(i,j) = config_flags%sfc_isltyp !7
283	grid%xice(i,j) = 0.
284	grid%snow(i,j) = 0.
285	END DO
286	END DO
287
288	CASE (RUCLSMSCHEME)
289	write(6,*) ' RUS surface scheme activated'
290	END SELECT other_masked_fields !AK/ak \|
291
292	ENDIF
293
294	DO j = jts, jte
295	DO k = kts, kte
296	DO i = its, ite
297	grid%ww(i,k,j) = 0.
298	END DO
299	END DO
300	END DO
301
302	grid%step_number = 0
303
304	IF (sfc_init) THEN
305
306	write(6,*) ' full surface initialization activated '
307	! write(6,*) ' land use index =', config_flags%sfc_lu_index
308	! write(6,*) ' skin temperature=',grid%tsk(10,10),&
309	! '[K] soil temperature=', grid%tmn(10,10),'[K]'
310	! Process the soil; note that there are some things hard-wired into share/module_soil_pre.F
311	CALL process_soil_ideal(grid%xland,grid%xice,grid%vegfra,grid%snow,grid%canwat, &
312	grid%ivgtyp,grid%isltyp,grid%tslb,grid%smois, &
313	grid%tsk,grid%tmn,grid%zs,grid%dzs,model_config_rec%num_soil_layers, &
314	model_config_rec%sf_surface_physics(grid%id), &
315	ids,ide, jds,jde, kds,kde,&
316	ims,ime, jms,jme, kms,kme,&
317	its,ite, jts,jte, kts,kte )
318
319
320	ELSE
321	write(6,*) 'full surface initialization is turned off!! '
322	ENDIF !end of surface initialization
323
324	! set up the grid
325	write(6,) '************************************'
326
327	IF (stretch_grd) THEN ! exponential or hyperbolic tangential stretch for eta
328
329	IF (stretch_hyp) THEN ! hyperbolic tangential stretch (more levels at the surface)
330	write(6,*) ' hyperbolic tangential stretching activated with z_scale =',z_scale
331	DO k=1, kde
332	grid%znw(k) = -1.* (tanh(z_scale*(float(k-1) / float(kde-1) -1.)))/ &
333	(tanh(z_scale))
334	ENDDO
335	ELSE ! exponential stretch for eta (nearly constant dz)
336	write(6,*) ' exponential grid stretching activated with z_scale =',z_scale
337	DO k=1, kde
338	grid%znw(k) = (exp(-(k-1)/float(kde-1)/z_scale) - exp(-1./z_scale))/ &
339	(1.-exp(-1./z_scale))
340	ENDDO
341	ENDIF
342	ELSE
343	write(6,*) ' no grid stretching'
344	DO k=1, kde
345	grid%znw(k) = 1. - float(k-1)/float(kde-1)
346	ENDDO
347	ENDIF
348	write(6,) '************************************'
349
350	DO k=1, kde-1
351	grid%dnw(k) = grid%znw(k+1) - grid%znw(k)
352	grid%rdnw(k) = 1./grid%dnw(k)
353	grid%znu(k) = 0.5*(grid%znw(k+1)+grid%znw(k))
354	ENDDO
355	DO k=2, kde-1
356	grid%dn(k) = 0.5*(grid%dnw(k)+grid%dnw(k-1))
357	grid%rdn(k) = 1./grid%dn(k)
358	grid%fnp(k) = .5* grid%dnw(k )/grid%dn(k)
359	grid%fnm(k) = .5* grid%dnw(k-1)/grid%dn(k)
360	ENDDO
361
362	cof1 = (2.grid%dn(2)+grid%dn(3))/(grid%dn(2)+grid%dn(3))grid%dnw(1)/grid%dn(2)
363	cof2 = grid%dn(2) /(grid%dn(2)+grid%dn(3))*grid%dnw(1)/grid%dn(3)
364	grid%cf1 = grid%fnp(2) + cof1
365	grid%cf2 = grid%fnm(2) - cof1 - cof2
366	grid%cf3 = cof2
367
368	grid%cfn = (.5*grid%dnw(kde-1)+grid%dn(kde-1))/grid%dn(kde-1)
369	grid%cfn1 = -.5*grid%dnw(kde-1)/grid%dn(kde-1)
370	grid%rdx = 1./config_flags%dx
371	grid%rdy = 1./config_flags%dy
372
373	! get the sounding from the ascii sounding file, first get dry sounding and
374	! calculate base state
375
376	dry_sounding = .true.
377	IF ( wrf_dm_on_monitor() ) THEN
378	write(6,*) ' getting dry sounding for base state '
379
380	CALL get_sounding( zk, p_in, pd_in, theta, rho, u, v, qv, dry_sounding, nl_max, nl_in )
381	ENDIF
382	CALL wrf_dm_bcast_real( zk , nl_max )
383	CALL wrf_dm_bcast_real( p_in , nl_max )
384	CALL wrf_dm_bcast_real( pd_in , nl_max )
385	CALL wrf_dm_bcast_real( theta , nl_max )
386	CALL wrf_dm_bcast_real( rho , nl_max )
387	CALL wrf_dm_bcast_real( u , nl_max )
388	CALL wrf_dm_bcast_real( v , nl_max )
389	CALL wrf_dm_bcast_real( qv , nl_max )
390	CALL wrf_dm_bcast_integer ( nl_in , 1 )
391
392	write(6,*) ' returned from reading sounding, nl_in is ',nl_in
393
394	! find ptop for the desired ztop (ztop is input from the namelist),
395	! and find surface pressure
396
397	grid%p_top = interp_0( p_in, zk, config_flags%ztop, nl_in )
398
399	! get fire mesh dimensions
400	CALL get_ijk_from_subgrid ( grid , &
401	ifds,ifde, jfds,jfde,kfds,kfde, &
402	ifms,ifme, jfms,jfme,kfms,kfme, &
403	ifts,ifte, jfts,jfte,kfts,kfte)
404
405	write (6,)' ***** SFIRE ideal initialization ******'
406
407	! fire grid step size
408	fdx = grid%dx/grid%sr_x
409	fdy = grid%dy/grid%sr_y
410	! refinement ratios
411	ir = grid%sr_x
412	jr = grid%sr_y
413
414	write (6,*)' atm mesh step ',grid%dx,grid%dy
415	write (6,*)' fire mesh step ',fdx,fdy
416	write (6,*)' refinement ratio ',grid%sr_x,grid%sr_y
417	write (6,*)' atm domain bounds ',ids,ide, jds,jde,kds,kde
418	write (6,*)' atm memory bounds ',ims,ime, jms,jme,kms,kme
419	write (6,*)' atm tile bounds ',its,ite, jts,jte,kts,kte
420	write (6,*)' fire domain bounds ',ifds,ifde, jfds,jfde,kfds,kfde
421	write (6,*)' fire memory bounds ',ifms,ifme, jfms,jfme,kfms,kfme
422	write (6,*)' fire tile bounds ',ifts,ifte, jfts,jfte,kfts,kfte
423	write (6,*)' Note that atm mesh and fire mesh are cell-centered'
424
425	! set ideal coordinates
426	call set_ideal_coord( fdx,fdy, &
427	ifds,ifde,jfds,jfde, &
428	ifms,ifme,jfms,jfme, &
429	ifts,ifte,jfts,jfte, &
430	grid%fxlong,grid%fxlat )
431	call set_ideal_coord( grid%dx,grid%dy, &
432	ids,ide,jds,jde, &
433	ims,ime,jms,jme, &
434	its,ite,jts,jte, &
435	grid%xlong,grid%xlat )
436
437	! set terrain height
438
439	DO j=jts,jte
440	DO i=its,ite
441	grid%ht(i,j) = 0.
442	ENDDO
443	ENDDO
444
445	if(config_flags%fire_fuel_read.eq.2) &
446	call read_array_2d_real('input_fc',grid%nfuel_cat,ifds,ifde,jfds,jfde,ifms,ifme,jfms,jfme)
447
448
449	have_fire_grad=.false.
450	have_atm_grad=.false.
451	have_fire_ht=.false.
452
453	!******* set terrain height
454
455	! copy params from the namelist
456	mtn_type = config_flags%fire_mountain_type
457	mtn_xs = config_flags%fire_mountain_start_x
458	mtn_ys = config_flags%fire_mountain_start_y
459	mtn_xe = config_flags%fire_mountain_end_x
460	mtn_ye = config_flags%fire_mountain_end_y
461	mtn_ht = config_flags%fire_mountain_height
462
463	IF(mtn_type .ne. 0)THEN
464
465	! idealized mountain
466
467	! atmospheric grid coordinates of the mountain
468	mtn_axs = mtn_xs/grid%dx + ids - 0.5
469	mtn_axe = mtn_xe/grid%dx + ids - 0.5
470	mtn_ays = mtn_ys/grid%dy + jds - 0.5
471	mtn_aye = mtn_ye/grid%dy + jds - 0.5
472
473	! fire grid coordinates of the mountain
474	mtn_fxs = mtn_xs/fdx + ifds - 0.5
475	mtn_fxe = mtn_xe/fdx + ifds - 0.5
476	mtn_fys = mtn_ys/fdy + jfds - 0.5
477	mtn_fye = mtn_ye/fdy + jfds - 0.5
478
479	write(6,*)' Mountain height ',mtn_ht,' type',mtn_type
480	write(6,*)' Mountain (m) LL=(0,0) ',mtn_xs,':',mtn_xe,' ',mtn_ys,':',mtn_ye
481	write(6,*)' Mountain on atm grid ',mtn_axs,':',mtn_axe,' ',mtn_ays,':',mtn_aye
482	write(6,*)' Mountain on fire grid ',mtn_fxs,':',mtn_fxe,' ',mtn_fys,':',mtn_fye
483
484	mtn_max = 0.
485	DO j=jts,jte
486	DO i=its,ite
487	mtn_x = pi + 2pi max(0. , min( (i - mtn_axs)/(mtn_axe - mtn_axs), 1. ))
488	mtn_y = pi + 2pi max(0. , min( (j - mtn_ays)/(mtn_aye - mtn_ays), 1. ))
489	SELECT CASE(mtn_type)
490	CASE (1) ! circ/elliptic mountain
491	mtn_z = mtn_ht * 0.25 * (1. + COS(mtn_x))*(1. + COS(mtn_y))
492	CASE (2) ! EW ridge
493	mtn_z = mtn_ht * 0.5 * (1. + COS(mtn_y))
494	CASE (3) ! NS ridge
495	mtn_z = mtn_ht * 0.5 * (1. + COS(mtn_x))
496	CASE DEFAULT
497	call wrf_error_fatal ( ' bad fire_mountain_type ' )
498	END SELECT
499	mtn_max = max(mtn_max, mtn_z)
500	grid%ht(i,j) = mtn_z
501	ENDDO
502	ENDDO
503
504	write(6, *)' Atm tile ',its,':',ite,' ',jts,':',jte,' max terrain height ',mtn_max
505
506	DO j=jfts,jfte
507	DO i=ifts,ifte
508	mtn_x = pi + 2pi max(0. , min( (i - mtn_fxs)/(mtn_fxe - mtn_fxs), 1. ))
509	mtn_y = pi + 2pi max(0. , min( (j - mtn_fys)/(mtn_fye - mtn_fys), 1. ))
510	SELECT CASE(mtn_type)
511	CASE (1) ! circ/elliptic mountain
512	mtn_z = mtn_ht * 0.25 * (1. + COS(mtn_x))*(1. + COS(mtn_y))
513	CASE (2) ! EW ridge
514	mtn_z = mtn_ht * 0.5 * (1. + COS(mtn_y))
515	CASE (3) ! NS ridge
516	mtn_z = mtn_ht * 0.5 * (1. + COS(mtn_x))
517	CASE DEFAULT
518	call wrf_error_fatal ( ' bad fire_mountain_type ' )
519	END SELECT
520	grid%zsf(i,j) = mtn_z
521	ENDDO
522	ENDDO
523
524	have_fire_ht=.true.
525
526	ELSE ! mtn_type
527
528	if(config_flags%fire_read_atm_ht)then !
529	call read_array_2d_real('input_ht',grid%ht,ids,ide,jds,jde,ims,ime,jms,jme)
530	! no flag - we always have the terrain height on atm mesh, zero if not set
531	endif
532
533	if(config_flags%fire_read_fire_ht)then !
534	call read_array_2d_real('input_zsf',grid%zsf,ifds,ifde,jfds,jfde,ifms,ifme,jfms,jfme)
535	have_fire_ht=.true.
536	endif
537
538	if(config_flags%fire_read_atm_grad)then !
539	call crash('Reading terrain gradient on atm mesh from file not supported.')
540	have_atm_grad=.true.
541	endif
542
543	if(config_flags%fire_read_fire_grad)then !
544	call read_array_2d_real('input_dzdxf',grid%dzdxf,ifds,ifde,jfds,jfde,ifms,ifme,jfms,jfme)
545	call read_array_2d_real('input_dzdyf',grid%dzdyf,ifds,ifde,jfds,jfde,ifms,ifme,jfms,jfme)
546	have_fire_grad=.true.
547	endif
548	ENDIF ! mtn_type
549
550	if(have_fire_ht)then
551	write(6, *)'Fine-resolution terrain height on the fire mesh used.'
552	else
553	write(6,*)'Interpolating the terrain height from the atm mesh to the fire mesh'
554	call interpolate_2d( &
555	ims,ime,jms,jme, & ! memory dims atm grid tile
556	its,ite,jts,jte, & ! where atm grid values set
557	ifms,ifme,jfms,jfme, & ! array dims fire grid
558	ifts,ifte,jfts,jfte, & ! dimensions fire grid tile
559	ir,jr, & ! refinement ratio
560	real(ids),real(jds),ifds+(ir-1)0.5,jfds+(jr-1)0.5, & ! line up by lower left corner of domain
561	grid%ht, & ! atm grid arrays in
562	grid%zsf) ! fire grid arrays out
563	have_fire_ht=.true.
564	endif
565
566
567	if(have_fire_grad)then
568	write(6, *)'Fine-resolution terrain gradient on the fire mesh used.'
569	else
570
571	write(6,*)'Computing the terrain gradient from fire mesh height'
572	if(.not.have_fire_ht)then
573	write(6,*)'WARNING: Fire mesh terrain height not given, setting to zero'
574	do j=jfts,jfte
575	do i=ifts,ifte
576	grid%zsf(i,j) = 0.
577	enddo
578	enddo
579	endif
580
581	! extend the terrain height one beyond the domain
582	call continue_at_boundary(1,1,0., & ! do x direction or y direction
583	ifms,ifme,jfms,jfme, & ! memory dims
584	ifds,ifde,jfds,jfde, & ! domain dims
585	ifds,ifde,jfds,jfde, & ! patch dims = domain, not parallel!
586	ifts,ifte,jfts,jfte, & ! tile dims
587	iots,iote,jots,jote, & ! tile dims out
588	grid%zsf) ! array
589
590	! compute the terrain gradient by differencing
591	do j=jfts,jfte
592	do i=ifts,ifte
593	grid%dzdxf(i,j) = (grid%zsf(i+1,j)-grid%zsf(i-1,j))/(2.*fdx)
594	grid%dzdyf(i,j) = (grid%zsf(i,j+1)-grid%zsf(i,j-1))/(2.*fdy)
595	enddo
596	enddo
597	have_fire_grad=.true.
598	endif ! have_fire_grad
599
600	if(.not.have_fire_grad)call crash('Fire mesh terrain gradient not set')
601
602	mtn_max = 0.
603	DO j=jts,jte
604	DO i=its,ite
605	mtn_max = max(mtn_max, grid%ht(i,j))
606	ENDDO
607	ENDDO
608	write(6, *)' Max terrain height on the atmosphere mesh ',mtn_max
609
610	mtn_max = 0.
611	grad_max =0.
612	DO j=jfts,jfte
613	DO i=ifts,ifte
614	mtn_max = max(mtn_max, grid%zsf(i,j))
615	grad_max = max( grad_max, sqrt(grid%dzdxf(i,j)2+grid%dzdyf(i,j)2) )
616	ENDDO
617	ENDDO
618	write(6, *)' Max terrain height on the fire mesh ',mtn_max
619	write(6, *)' Max terrain gradient on the fire mesh ',grad_max
620
621	! the rest of initialization dependent on the atmosphere grid terrain height set
622
623	DO j=jts,jte
624	DO i=its,ite
625	grid%phb(i,1,j) = g * grid%ht(i,j)
626	grid%ph0(i,1,j) = g * grid%ht(i,j)
627	ENDDO
628	ENDDO
629
630	DO J = jts, jte
631	DO I = its, ite
632
633	p_surf = interp_0( p_in, zk, grid%phb(i,1,j)/g, nl_in )
634	grid%mub(i,j) = p_surf-grid%p_top
635
636	! this is dry hydrostatic sounding (base state), so given grid%p (coordinate),
637	! interp theta (from interp) and compute 1/rho from eqn. of state
638
639	DO K = 1, kte-1
640	p_level = grid%znu(k)*(p_surf - grid%p_top) + grid%p_top
641	grid%pb(i,k,j) = p_level
642	grid%t_init(i,k,j) = interp_0( theta, p_in, p_level, nl_in ) - t0
643	grid%alb(i,k,j) = (r_d/p1000mb)(grid%t_init(i,k,j)+t0)(grid%pb(i,k,j)/p1000mb)**cvpm
644	ENDDO
645
646	! calc hydrostatic balance (alternatively we could interp the geopotential from the
647	! sounding, but this assures that the base state is in exact hydrostatic balance with
648	! respect to the model eqns.
649
650	DO k = 2,kte
651	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)
652	ENDDO
653
654	ENDDO
655	ENDDO
656
657	IF ( wrf_dm_on_monitor() ) THEN
658	write(6,*) ' ptop is ',grid%p_top
659	write(6,*) ' base state grid%mub(1,1), p_surf is ',grid%mub(1,1),grid%mub(1,1)+grid%p_top
660	ENDIF
661
662	! calculate full state for each column - this includes moisture.
663
664	write(6,*) ' getting moist sounding for full state '
665	dry_sounding = .false.
666	CALL get_sounding( zk, p_in, pd_in, theta, rho, u, v, qv, dry_sounding, nl_max, nl_in )
667
668	DO J = jts, min(jde-1,jte)
669	DO I = its, min(ide-1,ite)
670
671	! At this point grid%p_top is already set. find the DRY mass in the column
672	! by interpolating the DRY pressure.
673
674	pd_surf = interp_0( pd_in, zk, grid%phb(i,1,j)/g, nl_in )
675
676	! compute the perturbation mass and the full mass
677
678	grid%mu_1(i,j) = pd_surf-grid%p_top - grid%mub(i,j)
679	grid%mu_2(i,j) = grid%mu_1(i,j)
680	grid%mu0(i,j) = grid%mu_1(i,j) + grid%mub(i,j)
681
682	! given the dry pressure and coordinate system, interp the potential
683	! temperature and qv
684
685	do k=1,kde-1
686
687	p_level = grid%znu(k)*(pd_surf - grid%p_top) + grid%p_top
688
689	moist(i,k,j,P_QV) = interp_0( qv, pd_in, p_level, nl_in )
690	grid%t_1(i,k,j) = interp_0( theta, pd_in, p_level, nl_in ) - t0
691	grid%t_2(i,k,j) = grid%t_1(i,k,j)
692
693
694	enddo
695
696	! integrate the hydrostatic equation (from the RHS of the bigstep
697	! vertical momentum equation) down from the top to get grid%p.
698	! first from the top of the model to the top pressure
699
700	k = kte-1 ! top level
701
702	qvf1 = 0.5*(moist(i,k,j,P_QV)+moist(i,k,j,P_QV))
703	qvf2 = 1./(1.+qvf1)
704	qvf1 = qvf1*qvf2
705
706	! grid%p(i,k,j) = - 0.5*grid%mu_1(i,j)/grid%rdnw(k)
707	grid%p(i,k,j) = - 0.5(grid%mu_1(i,j)+qvf1grid%mub(i,j))/grid%rdnw(k)/qvf2
708	qvf = 1. + rvovrd*moist(i,k,j,P_QV)
709	grid%alt(i,k,j) = (r_d/p1000mb)(grid%t_1(i,k,j)+t0)qvf* &
710	(((grid%p(i,k,j)+grid%pb(i,k,j))/p1000mb)**cvpm)
711	grid%al(i,k,j) = grid%alt(i,k,j) - grid%alb(i,k,j)
712
713	! down the column
714
715	do k=kte-2,1,-1
716	qvf1 = 0.5*(moist(i,k,j,P_QV)+moist(i,k+1,j,P_QV))
717	qvf2 = 1./(1.+qvf1)
718	qvf1 = qvf1*qvf2
719	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)
720	qvf = 1. + rvovrd*moist(i,k,j,P_QV)
721	grid%alt(i,k,j) = (r_d/p1000mb)(grid%t_1(i,k,j)+t0)qvf* &
722	(((grid%p(i,k,j)+grid%pb(i,k,j))/p1000mb)**cvpm)
723	grid%al(i,k,j) = grid%alt(i,k,j) - grid%alb(i,k,j)
724	enddo
725
726	! this is the hydrostatic equation used in the model after the
727	! small timesteps. In the model, grid%al (inverse density)
728	! is computed from the geopotential.
729
730
731	grid%ph_1(i,1,j) = 0.
732	DO k = 2,kte
733	grid%ph_1(i,k,j) = grid%ph_1(i,k-1,j) - (1./grid%rdnw(k-1))*( &
734	(grid%mub(i,j)+grid%mu_1(i,j))*grid%al(i,k-1,j)+ &
735	grid%mu_1(i,j)*grid%alb(i,k-1,j) )
736
737	grid%ph_2(i,k,j) = grid%ph_1(i,k,j)
738	grid%ph0(i,k,j) = grid%ph_1(i,k,j) + grid%phb(i,k,j)
739	ENDDO
740
741	IF ( wrf_dm_on_monitor() ) THEN
742	if((i==2) .and. (j==2)) then
743	write(6,*) ' grid%ph_1 calc ',grid%ph_1(2,1,2),grid%ph_1(2,2,2),&
744	grid%mu_1(2,2)+grid%mub(2,2),grid%mu_1(2,2), &
745	grid%alb(2,1,2),grid%al(1,2,1),grid%rdnw(1)
746	endif
747	ENDIF
748
749	ENDDO
750	ENDDO
751
752	! checking if the perturbation (bubble) is to be applied
753	IF ((delt/=0.) .and. (x_rad > 0.) &
754	.and. (y_rad > 0.) &
755	.and. (z_rad > 0.)) THEN
756	! thermal perturbation to kick off convection
757
758	write(6,*) ' nxc, nyc for perturbation ',nxc,nyc
759	write(6,'(A23,f18.16)') ' delt for perturbation ',delt
760	write(6,'(A30,f18.12)') ' x radius of the perturbation ' ,x_rad
761	write(6,'(A30,f18.12)') ' y radius of the perturbation ' ,y_rad
762	write(6,'(A30,f18.12)') ' z radius of the perturbation ' ,z_rad
763	write(6,'(A30,f18.12)') ' height of the perturbation ' ,hght_pert
764
765	DO J = jts, min(jde-1,jte)
766	yrad = config_flags%dy*float(j-nyc)/y_rad
767	! yrad = 0.
768	DO I = its, min(ide-1,ite)
769	xrad = config_flags%dx*float(i-nxc)/x_rad
770	! xrad = 0.
771	DO K = 1, kte-1
772
773	! put in preturbation theta (bubble) and recalc density. note,
774	! the mass in the column is not changing, so when theta changes,
775	! we recompute density and geopotential
776
777	zrad = 0.5*(grid%ph_1(i,k,j)+grid%ph_1(i,k+1,j) &
778	+grid%phb(i,k,j)+grid%phb(i,k+1,j))/g
779	zrad = (zrad-hght_pert)/z_rad
780	RAD=SQRT(xradxrad+yradyrad+zrad*zrad)
781	IF(RAD <= 1.) THEN
782	grid%t_1(i,k,j)=grid%t_1(i,k,j)+deltCOS(.5PIRAD)*2
783	grid%t_2(i,k,j)=grid%t_1(i,k,j)
784	qvf = 1. + rvovrd*moist(i,k,j,P_QV)
785	grid%alt(i,k,j) = (r_d/p1000mb)(grid%t_1(i,k,j)+t0)qvf* &
786	(((grid%p(i,k,j)+grid%pb(i,k,j))/p1000mb)**cvpm)
787	grid%al(i,k,j) = grid%alt(i,k,j) - grid%alb(i,k,j)
788	ENDIF
789	ENDDO
790
791	! rebalance hydrostatically
792
793	DO k = 2,kte
794	grid%ph_1(i,k,j) = grid%ph_1(i,k-1,j) - (1./grid%rdnw(k-1))*( &
795	(grid%mub(i,j)+grid%mu_1(i,j))*grid%al(i,k-1,j)+ &
796	grid%mu_1(i,j)*grid%alb(i,k-1,j) )
797
798	grid%ph_2(i,k,j) = grid%ph_1(i,k,j)
799	grid%ph0(i,k,j) = grid%ph_1(i,k,j) + grid%phb(i,k,j)
800	ENDDO
801
802	ENDDO
803	ENDDO
804
805	!End of setting up the perturbation (bubble)
806	ENDIF
807
808	IF ( wrf_dm_on_monitor() ) THEN
809	write(6,*) ' grid%mu_1 from comp ', grid%mu_1(1,1)
810	write(6,*) ' full state sounding from comp, ph, grid%p, grid%al, grid%t_1, qv '
811	do k=1,kde-1
812	write(6,'(i3,1x,5(1x,1pe10.3))') k, grid%ph_1(1,k,1)+grid%phb(1,k,1), &
813	grid%p(1,k,1)+grid%pb(1,k,1), grid%alt(1,k,1), &
814	grid%t_1(1,k,1)+t0, moist(1,k,1,P_QV)
815	enddo
816
817	write(6,*) ' pert state sounding from comp, grid%ph_1, pp, alp, grid%t_1, qv '
818	do k=1,kde-1
819	write(6,'(i3,1x,5(1x,1pe10.3))') k, grid%ph_1(1,k,1), &
820	grid%p(1,k,1), grid%al(1,k,1), &
821	grid%t_1(1,k,1), moist(1,k,1,P_QV)
822	enddo
823	ENDIF
824
825	! interp v
826
827	DO J = jts, jte
828	DO I = its, min(ide-1,ite)
829
830	IF (j == jds) THEN
831	z_at_v = grid%phb(i,1,j)/g
832	ELSE IF (j == jde) THEN
833	z_at_v = grid%phb(i,1,j-1)/g
834	ELSE
835	z_at_v = 0.5*(grid%phb(i,1,j)+grid%phb(i,1,j-1))/g
836	END IF
837	p_surf = interp_0( p_in, zk, z_at_v, nl_in )
838
839	DO K = 1, kte-1
840	p_level = grid%znu(k)*(p_surf - grid%p_top) + grid%p_top
841	grid%v_1(i,k,j) = interp_0( v, p_in, p_level, nl_in )
842	grid%v_2(i,k,j) = grid%v_1(i,k,j)
843	ENDDO
844
845	ENDDO
846	ENDDO
847
848	! interp u
849
850	DO J = jts, min(jde-1,jte)
851	DO I = its, ite
852
853	IF (i == ids) THEN
854	z_at_u = grid%phb(i,1,j)/g
855	ELSE IF (i == ide) THEN
856	z_at_u = grid%phb(i-1,1,j)/g
857	ELSE
858	z_at_u = 0.5*(grid%phb(i,1,j)+grid%phb(i-1,1,j))/g
859	END IF
860
861	p_surf = interp_0( p_in, zk, z_at_u, nl_in )
862
863	DO K = 1, kte-1
864	p_level = grid%znu(k)*(p_surf - grid%p_top) + grid%p_top
865	grid%u_1(i,k,j) = interp_0( u, p_in, p_level, nl_in )
866	grid%u_2(i,k,j) = grid%u_1(i,k,j)
867	ENDDO
868
869	ENDDO
870	ENDDO
871
872	! set w
873
874	DO J = jts, min(jde-1,jte)
875	DO K = kts, kte
876	DO I = its, min(ide-1,ite)
877	grid%w_1(i,k,j) = 0.
878	grid%w_2(i,k,j) = 0.
879	ENDDO
880	ENDDO
881	ENDDO
882
883	! set a few more things
884
885	DO J = jts, min(jde-1,jte)
886	DO K = kts, kte-1
887	DO I = its, min(ide-1,ite)
888	grid%h_diabatic(i,k,j) = 0.
889	ENDDO
890	ENDDO
891	ENDDO
892
893	IF ( wrf_dm_on_monitor() ) THEN
894	DO k=1,kte-1
895	grid%t_base(k) = grid%t_1(1,k,1)
896	grid%qv_base(k) = moist(1,k,1,P_QV)
897	grid%u_base(k) = grid%u_1(1,k,1)
898	grid%v_base(k) = grid%v_1(1,k,1)
899	grid%z_base(k) = 0.5*(grid%phb(1,k,1)+grid%phb(1,k+1,1)+grid%ph_1(1,k,1)+grid%ph_1(1,k+1,1))/g
900	ENDDO
901	ENDIF
902	CALL wrf_dm_bcast_real( grid%t_base , kte )
903	CALL wrf_dm_bcast_real( grid%qv_base , kte )
904	CALL wrf_dm_bcast_real( grid%u_base , kte )
905	CALL wrf_dm_bcast_real( grid%v_base , kte )
906	CALL wrf_dm_bcast_real( grid%z_base , kte )
907
908	DO J = jts, min(jde-1,jte)
909	DO I = its, min(ide-1,ite)
910	thtmp = grid%t_2(i,1,j)+t0
911	ptmp = grid%p(i,1,j)+grid%pb(i,1,j)
912	temp(1) = thtmp * (ptmp/p1000mb)**rcp
913	thtmp = grid%t_2(i,2,j)+t0
914	ptmp = grid%p(i,2,j)+grid%pb(i,2,j)
915	temp(2) = thtmp * (ptmp/p1000mb)**rcp
916	thtmp = grid%t_2(i,3,j)+t0
917	ptmp = grid%p(i,3,j)+grid%pb(i,3,j)
918	temp(3) = thtmp * (ptmp/p1000mb)**rcp
919
920	! grid%tsk(I,J)=grid%cf1temp(1)+grid%cf2temp(2)+grid%cf3*temp(3) !AK/AK it is already declared via namelist.input if sfc_init=.true.
921	! grid%tmn(I,J)=grid%tsk(I,J)-0.5 !AK/AK it is already declared via namelist.input if sfc_init=.true.
922	ENDDO
923	ENDDO
924
925	IF (.NOT.sfc_init) THEN
926	write(6,*) ' setting tsk and tmn default'
927	DO J = jts, min(jde-1,jte)
928	DO I = its, min(ide-1,ite)
929	grid%tsk(I,J)=grid%cf1temp(1)+grid%cf2temp(2)+grid%cf3*temp(3)
930	grid%tmn(I,J)=grid%tsk(I,J)-0.5
931	ENDDO
932	ENDDO
933	ENDIF
934
935	END SUBROUTINE init_domain_rk
936
937	SUBROUTINE init_module_initialize
938	END SUBROUTINE init_module_initialize
939
940	!---------------------------------------------------------------------
941
942	! test driver for get_sounding
943	!
944	! implicit none
945	! integer n
946	! parameter(n = 1000)
947	! real zk(n),p(n),theta(n),rho(n),u(n),v(n),qv(n),pd(n)
948	! logical dry
949	! integer nl,k
950	!
951	! dry = .false.
952	! dry = .true.
953	! call get_sounding( zk, p, pd, theta, rho, u, v, qv, dry, n, nl )
954	! write(6,*) ' input levels ',nl
955	! write(6,*) ' sounding '
956	! write(6,*) ' k height(m) press (Pa) pd(Pa) theta (K) den(kg/m^3) u(m/s) v(m/s) qv(g/g) '
957	! do k=1,nl
958	! write(6,'(1x,i3,8(1x,1pe10.3))') k, zk(k), p(k), pd(k), theta(k), rho(k), u(k), v(k), qv(k)
959	! enddo
960	! end
961	!
962	!---------------------------------------------------------------------------
963
964	subroutine get_sounding( zk, p, p_dry, theta, rho, &
965	u, v, qv, dry, nl_max, nl_in )
966	implicit none
967
968	integer nl_max, nl_in
969	real zk(nl_max), p(nl_max), theta(nl_max), rho(nl_max), &
970	u(nl_max), v(nl_max), qv(nl_max), p_dry(nl_max)
971	logical dry
972
973	integer n
974	parameter(n=1000)
975	logical debug
976	parameter( debug = .true.)
977
978	! input sounding data
979
980	real p_surf, th_surf, qv_surf
981	real pi_surf, pi(n)
982	real h_input(n), th_input(n), qv_input(n), u_input(n), v_input(n)
983
984	! diagnostics
985
986	real rho_surf, p_input(n), rho_input(n)
987	real pm_input(n) ! this are for full moist sounding
988
989	! local data
990
991	real r
992	parameter (r = r_d)
993	integer k, it, nl
994	real qvf, qvf1, dz
995
996	! first, read the sounding
997
998	call read_sounding( p_surf, th_surf, qv_surf, &
999	h_input, th_input, qv_input, u_input, v_input,n, nl, debug )
1000
1001	if(dry) then
1002	do k=1,nl
1003	qv_input(k) = 0.
1004	enddo
1005	endif
1006
1007	if(debug) write(6,*) ' number of input levels = ',nl
1008
1009	nl_in = nl
1010	if(nl_in .gt. nl_max ) then
1011	write(6,*) ' too many levels for input arrays ',nl_in,nl_max
1012	call wrf_error_fatal ( ' too many levels for input arrays ' )
1013	end if
1014
1015	! compute diagnostics,
1016	! first, convert qv(g/kg) to qv(g/g)
1017
1018	do k=1,nl
1019	qv_input(k) = 0.001*qv_input(k)
1020	enddo
1021
1022	p_surf = 100.*p_surf ! convert to pascals
1023	qvf = 1. + rvovrd*qv_input(1)
1024	rho_surf = 1./((r/p1000mb)th_surfqvf((p_surf/p1000mb)*cvpm))
1025	pi_surf = (p_surf/p1000mb)**(r/cp)
1026
1027	if(debug) then
1028	write(6,*) ' surface density is ',rho_surf
1029	write(6,*) ' surface pi is ',pi_surf
1030	end if
1031
1032
1033	! integrate moist sounding hydrostatically, starting from the
1034	! specified surface pressure
1035	! -> first, integrate from surface to lowest level
1036
1037	qvf = 1. + rvovrd*qv_input(1)
1038	qvf1 = 1. + qv_input(1)
1039	rho_input(1) = rho_surf
1040	dz = h_input(1)
1041	do it=1,10
1042	pm_input(1) = p_surf &
1043	- 0.5dz(rho_surf+rho_input(1))gqvf1
1044	rho_input(1) = 1./((r/p1000mb)th_input(1)qvf((pm_input(1)/p1000mb)*cvpm))
1045	enddo
1046
1047	! integrate up the column
1048
1049	do k=2,nl
1050	rho_input(k) = rho_input(k-1)
1051	dz = h_input(k)-h_input(k-1)
1052	qvf1 = 0.5*(2.+(qv_input(k-1)+qv_input(k)))
1053	qvf = 1. + rvovrd*qv_input(k) ! qv is in g/kg here
1054
1055	do it=1,10
1056	pm_input(k) = pm_input(k-1) &
1057	- 0.5dz(rho_input(k)+rho_input(k-1))gqvf1
1058	rho_input(k) = 1./((r/p1000mb)th_input(k)qvf((pm_input(k)/p1000mb)*cvpm))
1059	enddo
1060	enddo
1061
1062	! we have the moist sounding
1063
1064	! next, compute the dry sounding using p at the highest level from the
1065	! moist sounding and integrating down.
1066
1067	p_input(nl) = pm_input(nl)
1068
1069	do k=nl-1,1,-1
1070	dz = h_input(k+1)-h_input(k)
1071	p_input(k) = p_input(k+1) + 0.5dz(rho_input(k)+rho_input(k+1))*g
1072	enddo
1073
1074
1075	do k=1,nl
1076
1077	zk(k) = h_input(k)
1078	p(k) = pm_input(k)
1079	p_dry(k) = p_input(k)
1080	theta(k) = th_input(k)
1081	rho(k) = rho_input(k)
1082	u(k) = u_input(k)
1083	v(k) = v_input(k)
1084	qv(k) = qv_input(k)
1085
1086	enddo
1087
1088	if(debug) then
1089	write(6,*) ' sounding '
1090	write(6,*) ' k height(m) press (Pa) pd(Pa) theta (K) den(kg/m^3) u(m/s) v(m/s) qv(g/g) '
1091	do k=1,nl
1092	write(6,'(1x,i3,8(1x,1pe10.3))') k, zk(k), p(k), p_dry(k), theta(k), rho(k), u(k), v(k), qv(k)
1093	enddo
1094
1095	end if
1096
1097	end subroutine get_sounding
1098
1099	!-------------------------------------------------------
1100
1101	subroutine read_sounding( ps,ts,qvs,h,th,qv,u,v,n,nl,debug )
1102	implicit none
1103	integer n,nl
1104	real ps,ts,qvs,h(n),th(n),qv(n),u(n),v(n)
1105	logical end_of_file
1106	logical debug
1107
1108	integer k
1109
1110	open(unit=10,file='input_sounding',form='formatted',status='old')
1111	rewind(10)
1112	read(10,*) ps, ts, qvs
1113	if(debug) then
1114	write(6,*) ' input sounding surface parameters '
1115	write(6,*) ' surface pressure (mb) ',ps
1116	write(6,*) ' surface pot. temp (K) ',ts
1117	write(6,*) ' surface mixing ratio (g/kg) ',qvs
1118	end if
1119
1120	end_of_file = .false.
1121	k = 0
1122
1123	do while (.not. end_of_file)
1124
1125	read(10,*,end=100) h(k+1), th(k+1), qv(k+1), u(k+1), v(k+1)
1126	k = k+1
1127	if(debug) write(6,'(1x,i3,5(1x,e10.3))') k, h(k), th(k), qv(k), u(k), v(k)
1128	go to 110
1129	100 end_of_file = .true.
1130	110 continue
1131	enddo
1132
1133	nl = k
1134
1135	close(unit=10,status = 'keep')
1136
1137	end subroutine read_sounding
1138
1139	END MODULE module_initialize_ideal

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

Download in other formats:

Original Format