Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

tc_em.F @ 1577

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

Line
1	! Create an initial data set for the WRF model based on real data. This
2	! program is specifically set up for the Eulerian, mass-based coordinate.
3	PROGRAM tc_data
4	USE module_machine
5	USE module_domain, ONLY : domain, alloc_and_configure_domain, &
6	domain_clock_set, head_grid, program_name, domain_clockprint, &
7	set_current_grid_ptr
8	USE module_io_domain
9	USE module_initialize_real, ONLY : wrfu_initialize
10	USE module_driver_constants
11	USE module_configure, ONLY : grid_config_rec_type, model_config_rec, &
12	initial_config, get_config_as_buffer, set_config_as_buffer
13	USE module_timing
14	USE module_state_description, ONLY : realonly
15	USE module_symbols_util, ONLY: wrfu_cal_gregorian
16	USE module_utility, ONLY : WRFU_finalize
17
18	IMPLICIT NONE
19
20
21	REAL :: time , bdyfrq
22
23	INTEGER :: loop , levels_to_process , debug_level
24
25
26	TYPE(domain) , POINTER :: null_domain
27	TYPE(domain) , POINTER :: grid , another_grid
28	TYPE(domain) , POINTER :: grid_ptr , grid_ptr2
29	TYPE (grid_config_rec_type) :: config_flags
30	INTEGER :: number_at_same_level
31
32	INTEGER :: max_dom, domain_id , grid_id , parent_id , parent_id1 , id
33	INTEGER :: e_we , e_sn , i_parent_start , j_parent_start
34	INTEGER :: idum1, idum2
35	#ifdef DM_PARALLEL
36	INTEGER :: nbytes
37	INTEGER, PARAMETER :: configbuflen = 4* CONFIG_BUF_LEN
38	INTEGER :: configbuf( configbuflen )
39	LOGICAL , EXTERNAL :: wrf_dm_on_monitor
40	#endif
41	LOGICAL found_the_id
42
43	INTEGER :: ids , ide , jds , jde , kds , kde
44	INTEGER :: ims , ime , jms , jme , kms , kme
45	INTEGER :: ips , ipe , jps , jpe , kps , kpe
46	INTEGER :: ijds , ijde , spec_bdy_width
47	INTEGER :: i , j , k , idts, rc
48	INTEGER :: sibling_count , parent_id_hold , dom_loop
49
50	CHARACTER (LEN=80) :: message
51
52	INTEGER :: start_year , start_month , start_day , start_hour , start_minute , start_second
53	INTEGER :: end_year , end_month , end_day , end_hour , end_minute , end_second
54	INTEGER :: interval_seconds , real_data_init_type
55	INTEGER :: time_loop_max , time_loop, bogus_id, storm
56	real::t1,t2
57	real :: latc_loc(max_bogus),lonc_loc(max_bogus),vmax(max_bogus),rmax(max_bogus)
58	real :: rankine_lid
59	INTERFACE
60	SUBROUTINE Setup_Timekeeping( grid )
61	USE module_domain, ONLY : domain
62	TYPE(domain), POINTER :: grid
63	END SUBROUTINE Setup_Timekeeping
64	END INTERFACE
65
66	#include "version_decl"
67
68	! Define the name of this program (program_name defined in module_domain)
69
70	program_name = "TC_EM " // TRIM(release_version) // " PREPROCESSOR"
71
72	! The TC bogus algorithm assumes that the user defines a central point, and then
73	! allows the program to remove a typhoon based on a distance in km. This is
74	! implemented on a single processor only.
75
76	#ifdef DM_PARALLEL
77	IF ( .NOT. wrf_dm_on_monitor() ) THEN
78	CALL wrf_error_fatal( 'TC bogus must run with a single processor only, re-run with num procs set to 1' )
79	END IF
80	#endif
81
82	#ifdef DM_PARALLEL
83	CALL disable_quilting
84	#endif
85
86	! Initialize the modules used by the WRF system. Many of the CALLs made from the
87	! init_modules routine are NO-OPs. Typical initializations are: the size of a
88	! REAL, setting the file handles to a pre-use value, defining moisture and
89	! chemistry indices, etc.
90
91	CALL wrf_debug ( 100 , 'real_em: calling init_modules ' )
92	CALL init_modules(1) ! Phase 1 returns after MPI_INIT() (if it is called)
93	#ifdef NO_LEAP_CALENDAR
94	CALL WRFU_Initialize( defaultCalendar=WRFU_CAL_NOLEAP, rc=rc )
95	#else
96	CALL WRFU_Initialize( defaultCalendar=WRFU_CAL_GREGORIAN, rc=rc )
97	#endif
98	CALL init_modules(2) ! Phase 2 resumes after MPI_INIT() (if it is called)
99
100	! The configuration switches mostly come from the NAMELIST input.
101
102	#ifdef DM_PARALLEL
103	IF ( wrf_dm_on_monitor() ) THEN
104	CALL initial_config
105	END IF
106	CALL get_config_as_buffer( configbuf, configbuflen, nbytes )
107	CALL wrf_dm_bcast_bytes( configbuf, nbytes )
108	CALL set_config_as_buffer( configbuf, configbuflen )
109	! CALL wrf_dm_initialize
110	#else
111	CALL initial_config
112	#endif
113
114
115	CALL nl_get_debug_level ( 1, debug_level )
116	CALL set_wrf_debug_level ( debug_level )
117
118	CALL wrf_message ( program_name )
119
120	! There are variables in the Registry that are only required for the real
121	! program, fields that come from the WPS package. We define the run-time
122	! flag that says to allocate space for these input-from-WPS-only arrays.
123
124	CALL nl_set_use_wps_input ( 1 , REALONLY )
125
126	! Allocate the space for the mother of all domains.
127
128	NULLIFY( null_domain )
129	CALL wrf_debug ( 100 , 'real_em: calling alloc_and_configure_domain ' )
130	CALL alloc_and_configure_domain ( domain_id = 1 , &
131	grid = head_grid , &
132	parent = null_domain , &
133	kid = -1 )
134
135	grid => head_grid
136	CALL nl_get_max_dom ( 1 , max_dom )
137
138	IF ( model_config_rec%interval_seconds .LE. 0 ) THEN
139	CALL wrf_error_fatal( 'namelist value for interval_seconds must be > 0')
140	END IF
141
142	all_domains : DO domain_id = 1 , max_dom
143
144	IF ( ( model_config_rec%input_from_file(domain_id) ) .OR. &
145	( domain_id .EQ. 1 ) ) THEN
146
147	CALL Setup_Timekeeping ( grid )
148	CALL set_current_grid_ptr( grid )
149	CALL domain_clockprint ( 150, grid, &
150	'DEBUG real: clock after Setup_Timekeeping,' )
151	CALL domain_clock_set( grid, &
152	time_step_seconds=model_config_rec%interval_seconds )
153	CALL domain_clockprint ( 150, grid, &
154	'DEBUG real: clock after timeStep set,' )
155
156
157	CALL wrf_debug ( 100 , 'tc_em: calling set_scalar_indices_from_config ' )
158	CALL set_scalar_indices_from_config ( grid%id , idum1, idum2 )
159
160	!This is goofy but we need to loop through the number of storms to get
161	!the namelist variables for the tc_bogus. But then we need to
162	!call model_to_grid_config_rec with the grid%id = to 1 in order to
163	!reset to the correct information.
164	CALL wrf_debug ( 100 , 'tc_em: calling model_to_grid_config_rec ' )
165	lonc_loc(:) = -999.
166	latc_loc(:) = -999.
167	vmax(:) = -999.
168	rmax(:) = -999.
169	CALL model_to_grid_config_rec ( grid%id , model_config_rec , config_flags )
170	lonc_loc(1) = config_flags%lonc_loc
171	latc_loc(1) = config_flags%latc_loc
172	vmax(1) = config_flags%vmax_meters_per_second
173	rmax(1) = config_flags%rmax
174	rankine_lid = config_flags%rankine_lid
175	do storm = 2,config_flags%num_storm
176	bogus_id = storm
177	CALL model_to_grid_config_rec ( bogus_id , model_config_rec , config_flags )
178	lonc_loc(storm) = config_flags%lonc_loc
179	latc_loc(storm) = config_flags%latc_loc
180	vmax(storm) = config_flags%vmax_meters_per_second
181	rmax(storm) = config_flags%rmax
182	! print *,"in loop ",storm,lonc_loc(storm),latc_loc(storm),vmax(storm),rmax(storm)
183	end do
184	CALL model_to_grid_config_rec ( grid%id , model_config_rec , config_flags )
185
186	! Initialize the WRF IO: open files, init file handles, etc.
187
188	CALL wrf_debug ( 100 , 'tc_em: calling init_wrfio' )
189	CALL init_wrfio
190
191	! Some of the configuration values may have been modified from the initial READ
192	! of the NAMELIST, so we re-broadcast the configuration records.
193
194	#ifdef DM_PARALLEL
195	CALL wrf_debug ( 100 , 'tc_em: re-broadcast the configuration records' )
196	CALL get_config_as_buffer( configbuf, configbuflen, nbytes )
197	CALL wrf_dm_bcast_bytes( configbuf, nbytes )
198	CALL set_config_as_buffer( configbuf, configbuflen )
199	#endif
200
201	! No looping in this layer.
202
203	CALL wrf_debug ( 100 , 'calling tc_med_sidata_input' )
204	CALL tc_med_sidata_input ( grid , config_flags, latc_loc, lonc_loc, &
205	vmax,rmax,rankine_lid)
206	CALL wrf_debug ( 100 , 'backfrom tc_med_sidata_input' )
207
208	ELSE
209	CYCLE all_domains
210	END IF
211
212	END DO all_domains
213
214	CALL set_current_grid_ptr( head_grid )
215
216	! We are done.
217
218	CALL wrf_debug ( 0 , 'tc_em: SUCCESS COMPLETE TC BOGUS' )
219
220	CALL wrf_shutdown
221
222	CALL WRFU_Finalize( rc=rc )
223
224
225	END PROGRAM tc_data
226
227
228	!-----------------------------------------------------------------
229	SUBROUTINE tc_med_sidata_input ( grid , config_flags, latc_loc, lonc_loc, &
230	vmax, rmax,rankine_lid)
231	! Driver layer
232	USE module_domain
233	USE module_io_domain
234	! Model layer
235	USE module_configure
236	USE module_bc_time_utilities
237	USE module_optional_input
238
239	USE module_date_time
240	USE module_utility
241
242	IMPLICIT NONE
243
244
245	! Interface
246	INTERFACE
247	SUBROUTINE start_domain ( grid , allowed_to_read ) ! comes from module_start in appropriate dyn_ directory
248	USE module_domain
249	TYPE (domain) grid
250	LOGICAL, INTENT(IN) :: allowed_to_read
251	END SUBROUTINE start_domain
252	END INTERFACE
253
254	! Arguments
255	TYPE(domain) :: grid
256	TYPE (grid_config_rec_type) :: config_flags
257	! Local
258	INTEGER :: time_step_begin_restart
259	INTEGER :: idsi , ierr , myproc, internal_time_loop,iflag
260	! Declarations for the netcdf routines.
261	INTEGER ::nf_inq
262	!
263	CHARACTER (LEN=80) :: si_inpname
264	CHARACTER (LEN=80) :: message
265
266	CHARACTER(LEN=19) :: start_date_char , end_date_char , current_date_char , next_date_char
267	CHARACTER(LEN=8) :: flag_name
268
269	INTEGER :: time_loop_max , loop, rc,icnt,itmp
270	INTEGER :: julyr , julday ,metndims, metnvars, metngatts, nunlimdimid,rcode
271	REAL :: gmt
272	real :: t1,t2,t3,t4
273	real :: latc_loc(max_bogus), lonc_loc(max_bogus)
274	real :: vmax(max_bogus),rmax(max_bogus),rankine_lid
275
276	grid%input_from_file = .true.
277	grid%input_from_file = .false.
278
279	CALL tc_compute_si_start ( model_config_rec%start_year (grid%id) , &
280	model_config_rec%start_month (grid%id) , &
281	model_config_rec%start_day (grid%id) , &
282	model_config_rec%start_hour (grid%id) , &
283	model_config_rec%start_minute(grid%id) , &
284	model_config_rec%start_second(grid%id) , &
285	model_config_rec%interval_seconds , &
286	model_config_rec%real_data_init_type , &
287	start_date_char)
288
289	end_date_char = start_date_char
290	IF ( end_date_char .LT. start_date_char ) THEN
291	CALL wrf_error_fatal( 'Ending date in namelist ' // end_date_char // ' prior to beginning date ' // start_date_char )
292	END IF
293	print *,"the start date char ",start_date_char
294	print *,"the end date char ",end_date_char
295
296	time_loop_max = 1
297	! Override stop time with value computed above.
298	CALL domain_clock_set( grid, stop_timestr=end_date_char )
299
300	! TBH: for now, turn off stop time and let it run data-driven
301	CALL WRFU_ClockStopTimeDisable( grid%domain_clock, rc=rc )
302	CALL wrf_check_error( WRFU_SUCCESS, rc, &
303	'WRFU_ClockStopTimeDisable(grid%domain_clock) FAILED', &
304	__FILE__ , &
305	__LINE__ )
306	CALL domain_clockprint ( 150, grid, &
307	'DEBUG med_sidata_input: clock after stopTime set,' )
308
309	! Here we define the initial time to process, for later use by the code.
310
311	current_date_char = start_date_char
312	start_date = start_date_char // '.0000'
313	current_date = start_date
314
315	CALL nl_set_bdyfrq ( grid%id , REAL(model_config_rec%interval_seconds) )
316
317
318	CALL cpu_time ( t1 )
319	DO loop = 1 , time_loop_max
320
321	internal_time_loop = loop
322	IF ( ( grid%id .GT. 1 ) .AND. ( loop .GT. 1 ) .AND. &
323	( model_config_rec%grid_fdda(grid%id) .EQ. 0 ) .AND. &
324	( model_config_rec%sst_update .EQ. 0 ) ) EXIT
325
326	print *,' '
327	print *,'-----------------------------------------------------------------------------'
328	print *,' '
329	print '(A,I2,A,A,A,I4,A,I4)' , &
330	' Domain ',grid%id,': Current date being processed: ',current_date, ', which is loop #',loop,' out of ',time_loop_max
331
332	! After current_date has been set, fill in the julgmt stuff.
333
334	CALL geth_julgmt ( config_flags%julyr , config_flags%julday , config_flags%gmt )
335
336	print *,'configflags%julyr, %julday, %gmt:',config_flags%julyr, config_flags%julday, config_flags%gmt
337	! Now that the specific Julian info is available, save these in the model config record.
338
339	CALL nl_set_gmt (grid%id, config_flags%gmt)
340	CALL nl_set_julyr (grid%id, config_flags%julyr)
341	CALL nl_set_julday (grid%id, config_flags%julday)
342
343	! Open the input file for tc stuff. Either the "new" one or the "old" one. The "new" one could have
344	! a suffix for the type of the data format. Check to see if either is around.
345
346	CALL cpu_time ( t3 )
347	WRITE ( wrf_err_message , FMT='(A,A)' )'med_sidata_input: calling open_r_dataset for ', &
348	TRIM(config_flags%auxinput1_inname)
349	CALL wrf_debug ( 100 , wrf_err_message )
350	IF ( config_flags%auxinput1_inname(1:8) .NE. 'wrf_real' ) THEN
351	CALL construct_filename4a( si_inpname , config_flags%auxinput1_inname , grid%id , 2 , &
352	current_date_char , config_flags%io_form_auxinput1 )
353	ELSE
354	CALL construct_filename2a( si_inpname , config_flags%auxinput1_inname , grid%id , 2 , &
355	current_date_char )
356	END IF
357	CALL open_r_dataset ( idsi, TRIM(si_inpname) , grid , config_flags , "DATASET=AUXINPUT1", ierr )
358	IF ( ierr .NE. 0 ) THEN
359	CALL wrf_error_fatal( 'error opening ' // TRIM(si_inpname) // &
360	' for input; bad date in namelist or file not in directory' )
361	END IF
362
363	! Input data.
364
365	CALL wrf_debug ( 100 , 'med_sidata_input: calling input_auxinput1' )
366	CALL input_auxinput1 ( idsi , grid , config_flags , ierr )
367	WRITE ( wrf_err_message , FMT='(A,I10,A)' ) 'Timing for input ',NINT(t4-t3) ,' s.'
368	CALL wrf_debug( 0, wrf_err_message )
369
370	! Possible optional SI input. This sets flags used by init_domain.
371
372	CALL cpu_time ( t3 )
373	CALL wrf_debug ( 100 , 'med_sidata_input: calling init_module_optional_input' )
374	CALL init_module_optional_input ( grid , config_flags )
375	CALL wrf_debug ( 100 , 'med_sidata_input: calling optional_input' )
376	CALL optional_input ( grid , idsi , config_flags )
377
378	!Here we check the flags yet again. The flags are checked in optional_input but
379	!the grid% flags are not set.
380	flag_name(1:8) = 'SM000010'
381	CALL wrf_get_dom_ti_integer ( idsi, 'FLAG_' // flag_name, itmp, 1, icnt, ierr )
382	IF ( ierr .EQ. 0 ) THEN
383	grid%flag_sm000010 = 1
384	end if
385
386	flag_name(1:8) = 'SM010040'
387	CALL wrf_get_dom_ti_integer ( idsi, 'FLAG_' // flag_name, itmp, 1, icnt, ierr )
388	IF ( ierr .EQ. 0 ) THEN
389	grid%flag_sm010040 = 1
390	end if
391
392	flag_name(1:8) = 'SM040100'
393	CALL wrf_get_dom_ti_integer ( idsi, 'FLAG_' // flag_name, itmp, 1, icnt, ierr )
394	IF ( ierr .EQ. 0 ) THEN
395	grid%flag_sm040100 = itmp
396	end if
397
398
399	flag_name(1:8) = 'SM100200'
400	CALL wrf_get_dom_ti_integer ( idsi, 'FLAG_' // flag_name, itmp, 1, icnt, ierr )
401	IF ( ierr .EQ. 0 ) THEN
402	grid%flag_sm100200 = itmp
403	end if
404
405	! flag_name(1:8) = 'SM010200'
406	! CALL wrf_get_dom_ti_integer ( idsi, 'FLAG_' // flag_name, itmp, 1, icnt, ierr )
407	! IF ( ierr .EQ. 0 ) THEN
408	! config_flags%flag_sm010200 = itmp
409	! print *,"found the flag_sm010200 "
410	! end if
411
412	!Now the soil temperature flags
413	flag_name(1:8) = 'ST000010'
414	CALL wrf_get_dom_ti_integer ( idsi, 'FLAG_' // flag_name, itmp, 1, icnt, ierr )
415	IF ( ierr .EQ. 0 ) THEN
416	grid%flag_st000010 = 1
417	END IF
418
419
420	flag_name(1:8) = 'ST010040'
421	CALL wrf_get_dom_ti_integer ( idsi, 'FLAG_' // flag_name, itmp, 1, icnt, ierr )
422	IF ( ierr .EQ. 0 ) THEN
423	grid%flag_st010040 = 1
424	END IF
425
426	flag_name(1:8) = 'ST040100'
427	CALL wrf_get_dom_ti_integer ( idsi, 'FLAG_' // flag_name, itmp, 1, icnt, ierr )
428	IF ( ierr .EQ. 0 ) THEN
429	grid%flag_st040100 = 1
430	END IF
431
432
433	flag_name(1:8) = 'ST100200'
434	CALL wrf_get_dom_ti_integer ( idsi, 'FLAG_' // flag_name, itmp, 1, icnt, ierr )
435	IF ( ierr .EQ. 0 ) THEN
436	grid%flag_st100200 = 1
437	END IF
438
439
440
441	CALL close_dataset ( idsi , config_flags , "DATASET=AUXINPUT1" )
442	CALL cpu_time ( t4 )
443
444	! Possible optional SI input. This sets flags used by init_domain.
445	! We need to call the optional input routines to get the flags that
446	! are in the metgrid output file so they can be put in the tc bogus
447	! output file for real to read.
448	CALL cpu_time ( t3 )
449	already_been_here = .FALSE.
450	CALL model_to_grid_config_rec ( grid%id , model_config_rec , config_flags )
451
452
453	CALL cpu_time ( t3 )
454
455	CALL assemble_output ( grid , config_flags , loop , time_loop_max, current_date_char, &
456	latc_loc, lonc_loc, vmax, rmax, rankine_lid,si_inpname)
457	CALL cpu_time ( t4 )
458	WRITE ( wrf_err_message , FMT='(A,I10,A)' ) 'Timing for output ',NINT(t4-t3) ,' s.'
459	CALL wrf_debug( 0, wrf_err_message )
460	CALL cpu_time ( t2 )
461	WRITE ( wrf_err_message , FMT='(A,I4,A,I10,A)' ) 'Timing for loop # ',loop,' = ',NINT(t2-t1) ,' s.'
462	CALL wrf_debug( 0, wrf_err_message )
463
464	CALL cpu_time ( t1 )
465	END DO
466
467	END SUBROUTINE tc_med_sidata_input
468
469
470	!-------------------------------------------------------------------------------------
471	SUBROUTINE tc_compute_si_start( &
472	start_year , start_month , start_day , start_hour , start_minute , start_second , &
473	interval_seconds , real_data_init_type , &
474	start_date_char)
475
476	USE module_date_time
477
478	IMPLICIT NONE
479
480	INTEGER :: start_year , start_month , start_day , start_hour , start_minute , start_second
481	INTEGER :: end_year , end_month , end_day , end_hour , end_minute , end_second
482	INTEGER :: interval_seconds , real_data_init_type
483	INTEGER :: time_loop_max , time_loop
484
485	CHARACTER(LEN=19) :: current_date_char , start_date_char , end_date_char , next_date_char
486
487	#ifdef PLANET
488	WRITE ( start_date_char , FMT = '(I4.4,"-",I5.5,"_",I2.2,":",I2.2,":",I2.2)' ) &
489	start_year,start_day,start_hour,start_minute,start_second
490	#else
491	WRITE ( start_date_char , FMT = '(I4.4,"-",I2.2,"-",I2.2,"_",I2.2,":",I2.2,":",I2.2)' ) &
492	start_year,start_month,start_day,start_hour,start_minute,start_second
493	#endif
494
495
496	END SUBROUTINE tc_compute_si_start
497
498	!-----------------------------------------------------------------------
499	SUBROUTINE assemble_output ( grid , config_flags , loop , time_loop_max,current_date_char, &
500	latc_loc, lonc_loc,vmax,rmax,rankine_lid,si_inpname)
501
502	USE module_big_step_utilities_em
503	USE module_domain
504	USE module_io_domain
505	USE module_configure
506	USE module_date_time
507	USE module_bc
508	IMPLICIT NONE
509
510	TYPE(domain) :: grid
511	TYPE (grid_config_rec_type) :: config_flags
512
513	INTEGER , INTENT(IN) :: loop , time_loop_max
514
515	!These values are in the name list and are avaiable from
516	!from the config_flags.
517	real :: vmax(max_bogus),vmax_ratio,rankine_lid
518	real :: rmax(max_bogus),stand_lon,cen_lat,ptop_in_pa
519	real :: latc_loc(max_bogus),lonc_loc(max_bogus)
520
521	INTEGER :: ijds , ijde , spec_bdy_width
522	INTEGER :: i , j , k , idts,map_proj,remove_only,storms
523
524	INTEGER :: id1 , interval_seconds , ierr, rc, sst_update, grid_fdda
525	INTEGER , SAVE :: id, id2, id4
526	CHARACTER (LEN=80) :: tcoutname , bdyname,si_inpname
527	CHARACTER(LEN= 4) :: loop_char
528	CHARACTER(LEN=19) :: current_date_char
529
530	character *19 :: temp19
531	character *24 :: temp24 , temp24b
532
533	real::t1,t2,truelat1,truelat2
534
535
536	! Boundary width, scalar value.
537
538	spec_bdy_width = model_config_rec%spec_bdy_width
539	interval_seconds = model_config_rec%interval_seconds
540	sst_update = model_config_rec%sst_update
541	grid_fdda = model_config_rec%grid_fdda(grid%id)
542	truelat1 = config_flags%truelat1
543	truelat2 = config_flags%truelat2
544
545	stand_lon = config_flags%stand_lon
546	cen_lat = config_flags%cen_lat
547	map_proj = config_flags%map_proj
548
549	vmax_ratio = config_flags%vmax_ratio
550	ptop_in_pa = config_flags%p_top_requested
551	remove_only = 0
552	if(config_flags%remove_storm) then
553	remove_only = 1
554	end if
555
556	storms = config_flags%num_storm
557	print *,"number of storms ",config_flags%num_storm
558	call tc_bogus(cen_lat,stand_lon,map_proj,truelat1,truelat2, &
559	grid%dx,grid%e_we,grid%e_sn,grid%num_metgrid_levels,ptop_in_pa, &
560	rankine_lid,latc_loc,lonc_loc,vmax,vmax_ratio,rmax,remove_only, &
561	storms,grid)
562
563
564
565	! Open the tc bogused output file. cd
566	CALL construct_filename4a( tcoutname , config_flags%auxinput1_outname , grid%id , 2 , &
567	current_date_char , config_flags%io_form_auxinput1 )
568
569	print *,"outfile name from construct filename ",tcoutname
570	CALL open_w_dataset ( id1, TRIM(tcoutname) , grid , config_flags ,output_auxinput1,"DATASET=AUXINPUT1",ierr )
571	IF ( ierr .NE. 0 ) THEN
572	CALL wrf_error_fatal( 'tc_em: error opening tc bogus file for writing' )
573	END IF
574	CALL output_auxinput1( id1, grid , config_flags , ierr )
575	CALL close_dataset ( id1 , config_flags , "DATASET=AUXINPUT1" )
576
577
578	END SUBROUTINE assemble_output
579
580	!----------------------------------------------------------------------------------------------
581
582	SUBROUTINE tc_bogus(centerlat,stdlon,nproj,truelat1,truelat2,dsm,ew,ns,nz,ptop_in_pa, &
583	rankine_lid,latc_loc,lonc_loc,vmax,vmax_ratio,rmax,remove_only, &
584	storms,grid)
585
586	!!Original Author Dave Gill. Modified by Sherrie Fredrick
587	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
588	!These are read in from the netcdf file.
589	!centerlat The center latitude from the global attributes in the netcdf file.
590	!stdlon The center longitude from the global attributes in the netcdf file.
591	!nproj The map projection from the global attributes in the netcdf file.
592	!dsm The spacing in meters from the global attributes in the netcdf file.
593	!ew The west_east_stag from the dimensions in the netcdf file..
594	!ns The south_north_stag from the dimensions in the netcdf file. .
595	!nz The number of metgrid levels from the dimensions in the netcdf file.
596
597	!ptop_in_pa This is part of the namelist.input file under the &domains section.
598
599	!These values are part of the namelist.input file under the &tc section specifically
600	!for the tc bogus code.
601	!NOTES: There can be up to five bogus storms. The variable max_bogus is set in
602	!the WRF subroutine called module_driver_constants.F in the ./WRFV3/frame directory.
603
604	!latc_loc The center latitude of the bogus strorm. This is an array dimensioned max_bogus.
605
606	!lonc_loc The center longitude of the bogus strorm. This is an array dimensioned max_bogus.
607
608	!vmax The max vortex in meters/second it comes from the namelist entry.
609	! This is an array dimensioned max_bogus.
610
611	!vmax_ratio This comes from the namelist entry.
612
613	!rmax The maximum radius this comes from the namelist entry.
614	! This is an array dimensioned max_bogus
615
616	!remove_only If this is set to true in the namelist.input file a value of 0.1
617	! is automatically assigned to vmax.
618
619	!rankine_lid This comes from the namelist entry. It can be used to determine
620	! what model levels the bogus storm affects.
621
622	!storms The number of bogus storms.
623
624	!grid This is a Fortran structure which holds all of the field data values
625	! for the netcdf that was read in.
626	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
627
628
629	!module_llxy resides in the share directory.
630	USE module_llxy
631	!This is for the large structure (grid)
632	USE module_domain
633
634
635
636	IMPLICIT NONE
637	TYPE(domain) :: grid
638	integer ew,ns,nz
639	integer nproj
640	integer storms,nstrm
641	real :: centerlat,stdlon,conef,truelat1,truelat2,dsm,dx,rankine_lid
642	real :: latc_loc(max_bogus),lonc_loc(max_bogus),vmax(max_bogus),vmax_ratio,rmax(max_bogus)
643
644	real :: press(ew-1,nz,ns-1),rhmx(nz), vwgt(nz),old_slp(ew-1,ns-1)
645	real, dimension(:,:,:) , allocatable :: u11,v11,t11,rh11,phi11
646	real, dimension(:,:,:) , allocatable :: u1 , v1 , t1 , rh1 , phi1
647	real, dimension(ew-1,ns-1) :: lond,terrain,cor,pslx
648
649
650	!The map scale factors.
651	real, dimension(ew,ns-1) :: msfu !The mapscale factor for the ew wind staggered grid
652	real, dimension(ew-1,ns) :: msfv !The mapscale factor for the ns wind staggered grid
653	real, dimension(ew-1,ns-1) :: msfm !The mapscale factor for the unstaggered grid.
654
655	CHARACTER*2 jproj
656	LOGICAL :: l_tcbogus
657
658
659	real :: r_search,r_vor,beta,devps,humidity_max
660	real :: devpc,const,r_vor2,cnst,alphar,epsilon,vormx , rad , sum_q
661	real :: avg_q ,q_old,ror,q_new,dph,dphx0
662	real :: rh_max,min_RH_value,ps
663	integer :: vert_variation
664	integer :: i,k,j,kx,remove_only
665	integer :: k00,kfrm ,kto ,k85,n_iter,ew_mvc,ns_mvc,nct,itr
666	integer :: strmci(nz), strmcj(nz)
667	real :: disx,disy,alpha,degran,pie,rovcp,cp
668	REAL :: rho,pprm,phip0,x0,y0,vmx,xico,xjco,xicn,xjcn,p85,xlo,rconst,ew_gcntr,ns_gcntr
669	real :: ptop_in_pa,themax,themin
670	real :: latinc,loninc
671	real :: rtemp,colat0,colat
672	REAL :: q1(ew-1,nz,ns-1), psi1(ew-1,nz,ns-1)
673
674	! This is the entire map projection enchilada.
675	TYPE(proj_info) :: proj
676
677
678
679	REAL :: lat1 , lon1
680	! These values are read in from the data set.
681	real :: knowni,knownj
682
683	! TC bogus
684	REAL utcr(ew,nz,ns-1), vtcr(ew-1,nz,ns)
685	REAL utcp(ew,nz,ns-1), vtcp(ew-1,nz,ns)
686	REAL psitc(ew-1,nz,ns-1), psiv(nz)
687	REAL vortc(ew-1,nz,ns-1), vorv(nz)
688	REAL tptc(ew-1,nz,ns-1)
689	REAL phiptc(ew-1,nz,ns-1)
690
691	! Work arrays
692	REAL uuwork(nz), vvwork(nz), temp2(ew,ns)
693	REAL vort(ew-1,nz,ns-1), div(ew-1,nz,ns-1)
694	REAL vortsv(ew-1,nz,ns-1)
695	REAL theta(ew-1,nz,ns-1), t_reduce(ew-1,nz,ns-1)
696	REAL ug(ew,nz,ns-1), vg(ew-1,nz,ns), vorg(ew-1,nz,ns-1)
697	REAL delpx(ew-1,ns-1)
698
699	!subroutines for relaxation
700	REAL outold(ew-1,ns-1)
701	REAL rd(ew-1,ns-1), ff(ew-1,ns-1)
702	REAL tmp1(ew-1,ns-1), tmp2(ew-1,ns-1)
703
704	! Background fields.
705	REAL , DIMENSION (ew-1,nz,ns-1) :: t0, t00, rh0, q0, phi0, psi0, chi
706
707	! Perturbations
708	REAL , DIMENSION (ew-1,nz,ns-1) :: psipos, tpos, psi ,phipos, phip
709
710	! Final fields.
711	REAL u2(ew,nz,ns-1), v2(ew-1,nz,ns)
712	REAL t2(ew-1,nz,ns-1),z2(ew-1,nz,ns-1)
713	REAL phi2(ew-1,nz,ns-1),rh2(ew-1,nz,ns-1)
714
715	print *,"the dimensions: north-south = ",ns," east-west =",ew
716	IF (nproj .EQ. 1) THEN
717	jproj = 'LC'
718	print *,"Lambert Conformal projection"
719	ELSE IF (nproj .EQ. 2) THEN
720	jproj = 'ST'
721	ELSE IF (nproj .EQ. 3) THEN
722	jproj = 'ME'
723	print *,"A mercator projection"
724	END IF
725
726
727	knowni = 1.
728	knownj = 1.
729	pie = 3.141592653589793
730	degran = pie/180.
731	rconst = 287.04
732	min_RH_value = 5.0
733	cp = 1004.0
734	rovcp = rconst/cp
735
736	r_search = 400000.0
737	r_vor = 300000.0
738	r_vor2 = r_vor * 4
739	beta = 0.5
740	devpc= 40.0
741	vert_variation = 1
742	humidity_max = 95.0
743	alphar = 1.8
744	latinc = -999.
745	loninc = -999.
746
747	if(remove_only .eq. 1) then
748	do nstrm=1,storms
749	vmax(nstrm) = 0.1
750	end do
751	end if
752
753	! Set up initializations for map projection so that the lat/lon
754	! of the tropical storm can be put into model (i,j) space. This needs to be done once per
755	! map projection definition. Since this is the domain that we are "GOING TO", it is a once
756	! per regridder requirement. If the user somehow ends up calling this routine for several
757	! time periods, there is no problemos, just a bit of overhead with redundant calls.
758
759	dx = dsm
760	lat1 = grid%xlat_gc(1,1)
761	lon1 = grid%xlong_gc(1,1)
762	IF( jproj .EQ. 'ME' )THEN
763	IF ( lon1 .LT. -180. ) lon1 = lon1 + 360.
764	IF ( lon1 .GT. 180. ) lon1 = lon1 - 360.
765	IF ( stdlon .LT. -180. ) stdlon = stdlon + 360.
766	IF ( stdlon .GT. 180. ) stdlon = stdlon - 360.
767	CALL map_set ( proj_merc, proj, lat1, lon1, lat1, lon1, knowni, knownj, dx, &
768	latinc,loninc,stdlon , truelat1 , truelat2)
769	conef = 0.
770	ELSE IF ( jproj .EQ. 'LC' ) THEN
771	if((truelat1 .eq. 0.0) .and. (truelat2 .eq. 0.0)) then
772	print *,"Truelat1 and Truelat2 are both 0"
773	stop
774	end if
775	CALL map_set (proj_lc,proj, lat1, lon1, lat1, lon1, knowni, knownj, dx, &
776	latinc,loninc,stdlon , truelat1 , truelat2)
777	conef = proj%cone
778	ELSE IF ( jproj .EQ. 'ST' ) THEN
779	conef = 1.
780	CALL map_set ( proj_ps,proj,lat1, lon1, lat1, lon1, knowni, knownj, dx, &
781	latinc,loninc,stdlon , truelat1 , truelat2)
782	END IF
783
784	! Load the pressure array.
785	kx = nz
786	do j = 1,ns-1
787	do k = 1,nz
788	do i = 1,ew-1
789	press(i,k,j) = grid%p_gc(i,k,j)*0.01
790	end do
791	end do
792	end do
793
794
795	! Initialize the vertical profiles for humidity and weighting.
796	!The ptop variable will be read in from the namelist
797	IF ( ( ptop_in_pa .EQ. 40000. ) .OR. ( ptop_in_pa .EQ. 60000. ) ) THEN
798	PRINT '(A)','Hold on pardner, your value for PTOP is gonna cause problems for the TC bogus option.'
799	PRINT '(A)','Make it higher up than 400 mb.'
800	STOP 'ptop_woes_for_tc_bogus'
801	END IF
802
803	IF ( vert_variation .EQ. 1 ) THEN
804	DO k=1,kx
805	IF ( press(1,k,1) .GT. 400. ) THEN
806	rhmx(k) = humidity_max
807	ELSE
808	rhmx(k) = humidity_max * MAX( 0.1 , (press(1,k,1) - ptop_in_pa/100.)/(400.-ptop_in_pa/100.) )
809	END IF
810
811	IF ( press(1,k,1) .GT. 600. ) THEN
812	vwgt(k) = 1.0
813	ELSE IF ( press(1,k,1) .LE. 100. ) THEN
814	vwgt(k) = 0.0001
815	ELSE
816	vwgt(k) = MAX ( 0.0001 , (press(1,k,1)-ptop_in_pa/100.)/(600.-ptop_in_pa/100.) )
817	END IF
818	END DO
819
820	ELSE IF ( vert_variation .EQ. 2 ) THEN
821	IF ( kx .eq. 24 ) THEN
822	rhmx = (/ 95., 95., 95., 95., 95., 95., 95., 95., &
823	95., 95., 95., 95., 95., 90., 85., 80., 75., &
824	70., 66., 60., 39., 10., 10., 10./)
825	vwgt = (/ 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 0.9850, &
826	0.9680, 0.9500, 0.9290, 0.9060, 0.8810, 0.8500, 0.7580, 0.6500, 0.5100, &
827	0.3500, 0.2120, 0.0500, 0.0270, 0.0001, 0.0001, 0.0001/)
828	ELSE
829	PRINT '(A)','Number of vertical levels assumed to be 24 for AFWA TC bogus option'
830	STOP 'AFWA_TC_BOGUS_LEVEL_ERROR'
831	END IF
832	END IF
833
834	!Remember that ns = the north south staggered. This is one more than the ns mass point grid.
835	! ew = the east west staggered. This is one more than the ew mass point grid.
836
837
838	!Put the U and V into the new arrays.
839	!Remember that the WRF ordering is ew,vert level,ns
840	!Vorticity and Divergence calculatins are done on
841	!the staggered grids so the winds are not destaggered
842	allocate(u11 (1:ew, 1:nz, 1:ns-1))
843	allocate(u1 (1:ew, 1:nz, 1:ns-1))
844	allocate(v11 (1:ew-1, 1:nz, 1:ns))
845	allocate(v1 (1:ew-1, 1:nz, 1:ns))
846	do j = 1,ns-1
847	do k = 1,nz
848	do i = 1,ew
849	u11(i,k,j) = grid%u_gc(i,k,j)
850	u1(i,k,j) = grid%u_gc(i,k,j)
851	msfu(i,j) = grid%msfu(i,j) !map scale factor on the U staggered grid
852	end do
853	end do
854	end do
855
856
857	do j = 1,ns
858	do k = 1,nz
859	do i = 1,ew-1
860	v11(i,k,j) = grid%v_gc(i,k,j)
861	v1(i,k,j) = grid%v_gc(i,k,j)
862	msfv(i,j) = grid%msfv(i,j) !map scale factor on the V staggered grid
863	end do
864	end do
865	end do
866
867
868	!Put the temperature, relative humidity and height fields
869	!into arrays. Save the initial fields also.
870	!These arrays are on the WRF mass points
871	allocate(t11 (1:ew-1, 1:nz, 1:ns-1))
872	allocate(t1 (1:ew-1, 1:nz, 1:ns-1))
873	allocate(rh11 (1:ew-1, 1:nz, 1:ns-1))
874	allocate(rh1 (1:ew-1, 1:nz, 1:ns-1))
875	allocate(phi11(1:ew-1, 1:nz, 1:ns-1))
876	allocate(phi1 (1:ew-1, 1:nz, 1:ns-1))
877	do j = 1,ns-1
878	do k = 1,nz
879	do i = 1,ew-1
880	t11(i,k,j) = grid%t_gc(i,k,j)
881	t1(i,k,j) = grid%t_gc(i,k,j)
882	rh11(i,k,j) = grid%rh_gc(i,k,j)
883	rh1(i,k,j) = grid%rh_gc(i,k,j)
884	msfm(i,j) = grid%msft(i,j)
885	if(k .eq. 1)then
886	phi11(i,k,j) = grid%ht_gc(i,j)
887	phi1(i,k,j) = grid%ht_gc(i,j) * 9.81
888	else
889	phi11(i,k,j) = grid%ght_gc(i,k,j)
890	phi1(i,k,j) = grid%ght_gc(i,k,j) * 9.81
891	end if
892	end do
893	end do
894	end do
895
896	!The two D fields
897	!The terrain soil height is from ght at level 1
898	do j = 1,ns-1
899	do i = 1,ew-1
900	pslx(i,j) = grid%pslv_gc(i,j) * 0.01
901	cor(i,j) = grid%f(i,j) !coreolous
902	lond(i,j) = grid%xlong_gc(i,j)
903	terrain(i,j) = grid%ht_gc(i,j)
904	old_slp(i,j) = grid%pslv_gc(i,j)
905	end do
906	end do
907
908
909
910	! Loop over the number of storms to process.
911
912	l_tcbogus = .FALSE.
913	all_storms : DO nstrm=1,storms
914
915
916	!Make sure the user has defined the rmax variable
917	if(rmax(nstrm) .eq. -999.) then
918	print *,"Please enter a value for rmax in the namelist"
919	stop
920	end if
921
922
923	k00 = 2
924	kfrm = k00
925	p85 = 850.
926
927	kto = kfrm
928	DO k=kfrm+1,kx
929	IF ( press(1,k,1) .GE. p85 ) THEN
930	kto = kto + 1
931	END IF
932	END DO
933	k85 = kto
934
935
936	! Parameters for max wind
937	rho = 1.2
938	pprm = devpc*100.
939	phip0= pprm/rho
940
941
942	!latc_loc and lonc_loc come in from the namelist.
943	!These x0 and y0 points are relative to the mass points.
944	CALL latlon_to_ij ( proj , latc_loc(nstrm) , lonc_loc(nstrm) , x0 , y0 )
945	IF ( ( x0 .LT. 1. ) .OR. ( x0 .GT. REAL(ew-1) ) .OR. &
946	( y0 .LT. 1. ) .OR. ( y0 .GT. REAL(ns-1) ) ) THEN
947	PRINT '(A,I3,A,A,A)',' Storm position is outside the computational domain.'
948	PRINT '(A,2F6.2,A)' ,' Storm postion: (x,y) = ',x0,y0,'.'
949	stop
950	END IF
951
952	l_tcbogus = .TRUE.
953	! Bogus vortex specifications, vmax (m/s); rmax (m);
954	vmx = vmax(nstrm) * vmax_ratio
955
956	IF ( latc_loc(nstrm) .LT. 0. ) THEN
957	vmx = -vmx
958	END IF
959
960	IF ( vmax(nstrm) .LE. 0. ) THEN
961	vmx = SQRT( 2.(1-beta)ABS(phip0) )
962	END IF
963
964	ew_gcntr = x0 !ew center grid location
965	ns_gcntr = y0 !ns center grid location
966	!For right now we are adding 0.5 to the grid location this
967	!makes the output of the wrf tc_bogus scheme analogous to the
968	!ouput of the MM5 tc_bogus scheme.
969	ew_gcntr = x0 + 0.5
970	ns_gcntr = y0 + 0.5
971
972	n_iter = 1
973
974	! Start computing.
975
976	PRINT '(/,A,I3,A,A,A)' ,'---> TC: Processing storm number= ',nstrm
977	PRINT '(A,F6.2,A,F7.2,A)' ,' Storm center lat= ',latc_loc(nstrm),' lon= ',lonc_loc(nstrm),'.'
978	PRINT '(A,2F6.2,A)' ,' Storm center grid position (x,y)= ',ew_gcntr,ns_gcntr,'.'
979	PRINT '(A,F5.2,F9.2,A)' ,' Storm max wind (m/s) and max radius (m)= ',vmx,rmax(nstrm),'.'
980	PRINT '(A,F5.2,A)' ,' Estimated central press dev (mb)= ',devpc,'.'
981
982
983	! Initialize storm center to (1,1)
984
985	DO k=1,kx
986	strmci(k) = 1
987	strmcj(k) = 1
988	END DO
989
990	! Define complete field of bogus storm
991	!Note dx is spacing in meters.
992	!The output arrays from the rankine subroutine vvwork,uuwork,psiv and vorv
993	!are defined on the WRF mass points.
994	utcp(:,:,:) = 0.0
995	vtcp(:,:,:) = 0.0
996	print *,"nstrm ",rmax(nstrm),ew_gcntr,ns_gcntr
997	DO j=1,ns-1
998	DO i=1,ew-1
999	disx = REAL(i) - ew_gcntr
1000	disy = REAL(j) - ns_gcntr
1001	CALL rankine(disx,disy,dx,kx,vwgt,rmax(nstrm),vmx,uuwork,vvwork,psiv,vorv)
1002	DO k=1,kx
1003	utcp(i,k,j) = uuwork(k)
1004	vtcp(i,k,j) = vvwork(k)
1005	psitc(i,k,j) = psiv(k)
1006	vortc(i,k,j) = vorv(k)
1007	END DO
1008	END DO
1009	END DO
1010	call stagger_rankine_winds(utcp,vtcp,ew,ns,nz)
1011
1012
1013	utcr(:,:,:) = 0.0
1014	vtcr(:,:,:) = 0.0
1015	! dave Rotate wind to map proj, on the correct staggering
1016	DO j=1,ns-1
1017	DO i=2,ew-1
1018	xlo = stdlon-grid%xlong_u(i,j)
1019	IF ( xlo .GT. 180.)xlo = xlo-360.
1020	IF ( xlo .LT.-180.)xlo = xlo+360.
1021
1022	alpha = xloconefdegran*SIGN(1.,centerlat)
1023	DO k=1,kx
1024	utcr(i,k,j) = (vtcp(i-1,k,j)+vtcp(i,k,j)+vtcp(i,k,j+1)+vtcp(i-1,k,j+1))/4 SIN(alpha)+utcp(i,k,j)COS(alpha)
1025	if(utcr(i,k,j) .gt. 300.) then
1026	print *,i,k,j,"a very bad value of utcr"
1027	stop
1028	end if
1029	END DO
1030	END DO
1031	END DO
1032
1033
1034	DO j=2,ns-1
1035	DO i=1,ew-1
1036	xlo = stdlon-grid%xlong_v(i,j)
1037	IF ( xlo .GT. 180.)xlo = xlo-360.
1038	IF ( xlo .LT.-180.)xlo = xlo+360.
1039
1040	alpha = xloconefdegran*SIGN(1.,centerlat)
1041	DO k=1,kx
1042	vtcr(i,k,j) = vtcp(i,k,j)COS(alpha)-(utcp(i,k,j-1)+utcp(i+1,k,j-1)+utcp(i+1,k,j)+utcp(i,k,j))/4SIN(alpha)
1043	if(vtcr(i,k,j) .gt. 300.) then
1044	print *,i,k,j,"a very bad value of vtcr"
1045	stop
1046	end if
1047	END DO
1048	END DO
1049	END DO
1050
1051
1052	!Fill in UTCR's along the left and right side.
1053	do j = 1,ns-1
1054	utcr(1,:,j) = utcr(2,:,j)
1055	utcr(ew,:,j) = utcr(ew-1,:,j)
1056	end do
1057
1058	!Fill in V's along the bottom and top.
1059	do i = 1,ew-1
1060	vtcr(i,:,1) = vtcr(i,:,2)
1061	vtcr(i,:,ns) = vtcr(i,:,ns-1)
1062	end do
1063
1064
1065	! Compute vorticity of FG. This is the vorticity of the original winds
1066	! on the staggered grid. The vorticity and divergence are defined at
1067	! the mass points when done.
1068	CALL vor(u1,v1,msfu,msfv,msfm,ew,ns,kx,dx,vort)
1069
1070
1071	! Compute divergence of FG
1072	CALL diverg(u1,v1,msfu,msfv,msfm,ew,ns,kx,dx,div)
1073
1074
1075	! Compute mixing ratio of FG
1076	CALL mxratprs(rh1,t1,press*100.,ew,ns,kx,q1,min_RH_value)
1077	q1(:,1,:) = q1(:,2,:)
1078
1079
1080	! Compute initial streamfunction - PSI1
1081	vortsv = vort
1082	q0 = q1
1083
1084
1085	! Solve for streamfunction.
1086	DO k=1,kx
1087	DO j=1,ns-1
1088	DO i=1,ew-1
1089	ff(i,j) = vort(i,k,j)
1090	tmp1(i,j)= 0.0
1091	END DO
1092	END DO
1093	epsilon = 1.E-2
1094	CALL relax(tmp1,ff,rd,ew,ns,dx,epsilon,alphar)
1095	DO j=1,ns-1
1096	DO i=1,ew-1
1097	psi1(i,k,j) = tmp1(i,j)
1098	END DO
1099	END DO
1100	END DO
1101
1102
1103	DO k=1,kx !start of the k loop
1104	IF ( latc_loc(nstrm) .GE. 0. ) THEN
1105	vormx = -1.e10
1106	ELSE
1107	vormx = 1.e10
1108	END IF
1109
1110	ew_mvc = 1
1111	ns_mvc = 1
1112
1113	DO j=1,ns-1
1114	DO i=1,ew-1
1115	rad = SQRT((REAL(i)-ew_gcntr)2.+(REAL(j)-ns_gcntr)2.)*dx
1116	IF ( rad .LE. r_search ) THEN
1117	IF ( latc_loc(nstrm) .GE. 0. ) THEN
1118	IF ( vortsv(i,k,j) .GT. vormx ) THEN
1119	vormx = vortsv(i,k,j)
1120	ew_mvc = i
1121	ns_mvc = j
1122	END IF
1123	ELSE IF (latc_loc(nstrm) .LT. 0. ) THEN
1124	IF ( vortsv(i,k,j) .LT. vormx ) THEN
1125	vormx = vortsv(i,k,j)
1126	ew_mvc = i
1127	ns_mvc = j
1128	END IF
1129	END IF
1130	END IF
1131	END DO
1132	END DO
1133
1134	strmci(k) = ew_mvc
1135	strmcj(k) = ns_mvc
1136
1137	DO j=1,ns-1
1138	DO i=1,ew-1
1139	rad = SQRT(REAL((i-ew_mvc)2.+(j-ns_mvc)2.))*dx
1140	IF ( rad .GT. r_vor ) THEN
1141	vort(i,k,j) = 0.
1142	div(i,k,j) = 0.
1143	END IF
1144	END DO
1145	END DO
1146
1147	DO itr=1,n_iter
1148	sum_q = 0.
1149	nct = 0
1150	DO j=1,ns-1
1151	DO i=1,ew-1
1152	rad = SQRT(REAL(i-ew_mvc)2.+REAL(j-ns_mvc)2.)*dx
1153	IF ( (rad .LT. r_vor2).AND.(rad .GE. 0.8*r_vor2) ) THEN
1154	sum_q = sum_q + q0(i,k,j)
1155	nct = nct + 1
1156	END IF
1157	END DO
1158	END DO
1159	avg_q = sum_q/MAX(REAL(nct),1.)
1160
1161	DO j=1,ns-1
1162	DO i=1,ew-1
1163	q_old = q0(i,k,j)
1164	rad = SQRT(REAL(i-ew_mvc)2.+REAL(j-ns_mvc)2.)*dx
1165	IF ( rad .LT. r_vor2 ) THEN
1166	ror = rad/r_vor2
1167	q_new = ((1.-ror)avg_q) + (rorq_old)
1168	q0(i,k,j) = q_new
1169	END IF
1170	END DO
1171	END DO
1172	END DO !end of itr loop
1173	END DO !of the k loop
1174
1175
1176	! Compute divergent wind (chi) at the mass points
1177	DO k=1,kx
1178	DO j=1,ns-1
1179	DO i=1,ew-1
1180	ff(i,j) = div(i,k,j)
1181	tmp1(i,j)= 0.0
1182	END DO
1183	END DO
1184
1185	epsilon = 1.e-2
1186	CALL relax(tmp1,ff,rd,ew,ns,dx,epsilon,alphar)
1187	DO j=1,ns-1
1188	DO i=1,ew-1
1189	chi(i,k,j) = tmp1(i,j)
1190	END DO
1191	END DO
1192	END DO !of the k loop for divergent winds
1193
1194
1195
1196	! Compute background streamfunction (PSI0) and perturbation field (PSI)
1197	! print *,"perturbation field (PSI) relax three"
1198	DO k=1,kx
1199	DO j=1,ns-1
1200	DO i=1,ew-1
1201	ff(i,j)=vort(i,k,j)
1202	tmp1(i,j)=0.0
1203	END DO
1204	END DO
1205	epsilon = 1.e-2
1206	CALL relax(tmp1,ff,rd,ew,ns,dx,epsilon,alphar)
1207	DO j=1,ns-1
1208	DO i=1,ew-1
1209	psi(i,k,j)=tmp1(i,j)
1210	END DO
1211	END DO
1212	END DO
1213
1214
1215	!We can now calculate the final wind fields.
1216	call final_ew_velocity(u2,u1,chi,psi,utcr,dx,ew,ns,nz)
1217	call final_ns_velocity(v2,v1,chi,psi,vtcr,dx,ew,ns,nz)
1218
1219	DO k=1,kx
1220	DO j=1,ns-1
1221	DO i=1,ew-1
1222	psi0(i,k,j) = psi1(i,k,j)-psi(i,k,j)
1223	END DO
1224	END DO
1225	END DO
1226
1227	DO k=k00,kx
1228	DO j=1,ns-1
1229	DO i=1,ew-1
1230	psipos(i,k,j)=psi(i,k,j)
1231	END DO
1232	END DO
1233	END DO
1234
1235
1236	! Geostrophic vorticity.
1237	!We calculate the ug and vg on the wrf U and V staggered grids
1238	!since this is where the vorticity subroutine expects them.
1239
1240	CALL geowind(phi1,ew,ns,kx,dx,ug,vg)
1241	CALL vor(ug,vg,msfu,msfv,msfm,ew,ns,kx,dx,vorg)
1242
1243	DO k=1,kx
1244	ew_mvc = strmci(k)
1245	ns_mvc = strmcj(k)
1246
1247	DO j=1,ns-1
1248	DO i=1,ew-1
1249	rad = SQRT(REAL(i-ew_mvc)2.+REAL(j-ns_mvc)2.)*dx
1250	IF ( rad .GT. r_vor ) THEN
1251	vorg(i,k,j) = 0.
1252	END IF
1253	END DO
1254	END DO
1255	END DO
1256
1257	DO k=k00,kx
1258	DO j=1,ns-1
1259	DO i=1,ew-1
1260	ff(i,j) = vorg(i,k,j)
1261	tmp1(i,j)= 0.0
1262	END DO
1263	END DO
1264	epsilon = 1.e-3
1265	CALL relax(tmp1,ff,rd,ew,ns,dx,epsilon,alphar)
1266	DO j=1,ns-1
1267	DO i=1,ew-1
1268	phip(i,k,j) = tmp1(i,j)
1269	END DO
1270	END DO
1271	END DO
1272
1273
1274	! Background geopotential.
1275	DO k=k00,kx
1276	DO j=1,ns-1
1277	DO i=1,ew-1
1278	phi0(i,k,j) = phi1(i,k,j) - phip(i,k,j)
1279	END DO
1280	END DO
1281	END DO
1282
1283
1284	! Background temperature
1285	DO k=k00,kx
1286	DO j=1,ns-1
1287	DO i=1,ew-1
1288	IF( k .EQ. 2 ) THEN
1289	tpos(i,k,j) = (-1./rconst)*(phip(i,k+1,j)-phip(i,k,j ))/LOG(press(i,k+1,j)/press(i,k,j))
1290	ELSE IF ( k .EQ. kx ) THEN
1291	tpos(i,k,j) = (-1./rconst)*(phip(i,k ,j)-phip(i,k-1,j))/LOG(press(i,k,j )/press(i,k-1,j))
1292	ELSE
1293	tpos(i,k,j) = (-1./rconst)*(phip(i,k+1,j)-phip(i,k-1,j))/LOG(press(i,k+1,j)/press(i,k-1,j))
1294	END IF
1295	t0(i,k,j) = t1(i,k,j)-tpos(i,k,j)
1296	t00(i,k,j) = t0(i,k,j)
1297	if(t0(i,k,j) .gt. 400) then
1298	print *,"interesting temperature ",t0(i,k,j)," at ",i,j,k
1299	stop
1300	end if
1301	END DO
1302	END DO
1303	END DO
1304
1305	! New RH.
1306	CALL qvtorh (q0,t0,press*100.,k00,ew,ns,kx,rh0,min_RH_value)
1307	call final_RH(rh2,rh0,rhmx,strmci,strmcj,rmax(nstrm),ew,ns,nz,k00,dx,ew_gcntr,ns_gcntr,r_vor2)
1308
1309
1310
1311	! adjust T0
1312	DO k=k00,kx
1313	DO j=1,ns-1
1314	DO i=1,ew-1
1315	theta(i,k,j) = t1(i,k,j)(1000./press(i,k,j))*rovcp
1316	END DO
1317	END DO
1318	END DO
1319
1320
1321	ew_mvc = strmci(k00)
1322	ns_mvc = strmcj(k00)
1323	DO k=kfrm,kto
1324	DO j=1,ns-1
1325	DO i=1,ew-1
1326	rad = SQRT(REAL(i-ew_mvc)2.+REAL(j-ns_mvc)2.)*dx
1327	IF ( rad .LT. r_vor2 ) THEN
1328	t_reduce(i,k,j) = theta(i,k85,j)-0.03*(press(i,k,j)-press(i,k85,j))
1329	t0(i,k,j) = t00(i,k,j)(rad/r_vor2) + (((press(i,k,j)/1000.)rovcp)t_reduce(i,k,j))*(1.-(rad/r_vor2))
1330	END IF
1331	END DO
1332	END DO
1333	END DO
1334
1335	! Geopotential perturbation
1336	DO k=1,kx
1337	DO j=1,ns-1
1338	DO i=1,ew-1
1339	tmp1(i,j)=psitc(i,k,j)
1340	END DO
1341	END DO
1342	CALL balance(cor,tmp1,ew,ns,dx,outold)
1343	DO j=1,ns-1
1344	DO i=1,ew-1
1345	ff(i,j)=outold(i,j)
1346	tmp1(i,j)=0.0
1347	END DO
1348	END DO
1349	epsilon = 1.e-3
1350	CALL relax (tmp1,ff,rd,ew,ns,dx,epsilon,alphar)
1351	DO j=1,ns-1
1352	DO i=1,ew-1
1353	phiptc(i,k,j) = tmp1(i,j)
1354	END DO
1355	END DO
1356	END DO
1357
1358
1359	! New geopotential field.
1360	DO j=1,ns-1
1361	DO k=1,kx
1362	DO i=1,ew-1
1363	phi2(i,k,j) = phi0(i,k,j) + phiptc(i,k,j)
1364	END DO
1365	END DO
1366	END DO
1367
1368
1369	! New temperature field.
1370	DO j=1,ns-1
1371	DO k=k00,kx
1372	DO i=1,ew-1
1373	IF( k .EQ. 2 ) THEN
1374	tptc(i,k,j)=(-1./rconst)*(phiptc(i,k+1,j)-phiptc(i,k,j ))/LOG(press(i,k+1,j)/press(i,k,j))
1375	ELSE IF ( k .EQ. kx ) THEN
1376	tptc(i,k,j)=(-1./rconst)*(phiptc(i,k,j )-phiptc(i,k-1,j))/LOG(press(i,k,j)/press(i,k-1,j))
1377	ELSE
1378	tptc(i,k,j)=(-1./rconst)*(phiptc(i,k+1,j)-phiptc(i,k-1,j))/LOG(press(i,k+1,j)/press(i,k-1,j))
1379	END IF
1380	t2(i,k,j) = t0(i,k,j) + tptc(i,k,j)
1381	if(t2(i,k,j) .gt. 400) then
1382	print *,"interesting temperature "
1383	print *,t2(i,k,j),i,k,j,tptc(i,k,j)
1384	stop
1385	end if
1386	END DO
1387	END DO
1388	END DO
1389
1390
1391	! Sea level pressure change.
1392	DO j=1,ns-1
1393	DO i=1,ew-1
1394	dph = phi2(i,k00,j)-phi1(i,k00,j)
1395	delpx(i,j) = rhodph0.01
1396	END DO
1397	END DO
1398
1399
1400	! New SLP.
1401	! print *,"new slp",nstrm
1402	DO j=1,ns-1
1403	DO i=1,ew-1
1404	pslx(i,j) = pslx(i,j)+delpx(i,j)
1405	grid%pslv_gc(i,j) = pslx(i,j) * 100.
1406	! print *,pslx(i,j)
1407	END DO
1408	END DO
1409
1410	! Set new geopotential at surface to terrain elevation.
1411	DO j=1,ns-1
1412	DO i=1,ew-1
1413	z2(i,1,j) = terrain(i,j)
1414	END DO
1415	END DO
1416
1417	! Geopotential back to height.
1418
1419	DO j=1,ns-1
1420	DO k=k00,kx
1421	DO i=1,ew-1
1422	z2(i,k,j) = phi2(i,k,j)/9.81
1423	END DO
1424	END DO
1425	END DO
1426
1427
1428	! New surface temperature, assuming same theta as from 1000 mb.
1429	! print *,"new surface temperature"
1430	DO j=1,ns-1
1431	DO i=1,ew-1
1432	ps = pslx(i,j)
1433	t2(i,1,j) = t2(i,k00,j)((ps/1000.)*rovcp)
1434	if(t2(i,1,j) .gt. 400) then
1435	print *,"Interesting surface temperature"
1436	print *,t2(i,1,j),t2(i,k00,j),ps,i,j
1437	stop
1438	end if
1439	END DO
1440	END DO
1441
1442
1443	! Set surface RH to the value from 1000 mb.
1444	DO j=1,ns-1
1445	DO i=1,ew-1
1446	rh2(i,1,j) = rh2(i,k00,j)
1447	END DO
1448	END DO
1449
1450	! Modification of tropical storm complete.
1451	PRINT '(A,I3,A)' ,' Bogus storm number ',nstrm,' completed.'
1452
1453	do j = 1,ns-1
1454	do k = 1,nz
1455	do i = 1,ew
1456	u1(i,k,j) = u2(i,k,j)
1457	grid%u_gc(i,k,j) = u2(i,k,j)
1458	end do
1459	end do
1460	end do
1461
1462	do j = 1,ns
1463	do k = 1,nz
1464	do i = 1,ew-1
1465	v1(i,k,j) = v2(i,k,j)
1466	grid%v_gc(i,k,j) = v2(i,k,j)
1467	end do
1468	end do
1469	end do
1470
1471	do j = 1,ns-1
1472	do k = 1,nz
1473	do i = 1,ew-1
1474	t1(i,k,j) = t2(i,k,j)
1475	grid%t_gc(i,k,j) = t2(i,k,j)
1476	rh1(i,k,j) = rh2(i,k,j)
1477	grid%rh_gc(i,k,j) = rh2(i,k,j)
1478	phi1(i,k,j) = phi2(i,k,j)
1479	grid%ght_gc(i,k,j) = z2(i,k,j)
1480	END DO
1481	END DO
1482	END DO
1483
1484
1485	END DO all_storms
1486	deallocate(u11)
1487	deallocate(v11)
1488	deallocate(t11)
1489	deallocate(rh11)
1490	deallocate(phi11)
1491	deallocate(u1)
1492	deallocate(v1)
1493	deallocate(t1)
1494	deallocate(rh1)
1495	deallocate(phi1)
1496
1497	do j = 1,ns-1
1498	do i = 1,ew-1
1499	if(grid%ht_gc(i,j) .gt. 1) then
1500	grid%p_gc(i,1,j) = grid%p_gc(i,1,j) + (pslx(i,j) * 100. - old_slp(i,j))
1501	grid%psfc(i,j) = grid%psfc(i,j) + (pslx(i,j) * 100. - old_slp(i,j))
1502	else
1503	grid%p_gc(i,1,j) = pslx(i,j) * 100.
1504	grid%psfc(i,j) = pslx(i,j) * 100.
1505	end if
1506	end do
1507	end do
1508
1509	END SUBROUTINE tc_bogus
1510
1511	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1512	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1513
1514	SUBROUTINE rankine(dx,dy,ds,nlvl,vwgt,rmax,vmax,uu,vv,psi,vor)
1515
1516	! Define analytical bogus vortex
1517
1518	IMPLICIT NONE
1519
1520	INTEGER nlvl
1521	REAL , DIMENSION(nlvl) :: uu, vv, psi, vor
1522	REAL , DIMENSION(nlvl) :: vwgt
1523	REAL :: dx,dy,ds,rmax,vmax
1524
1525	REAL , PARAMETER :: alpha1= 1.
1526	REAL , PARAMETER :: alpha2= -0.75
1527	real :: pi
1528
1529
1530	INTEGER :: k
1531	REAL :: vr , ang , rr , term1 , bb , term2 , alpha
1532
1533
1534	pi = 3.141592653589793
1535	! Wind component
1536
1537	DO k=1,nlvl
1538	rr = SQRT(dx2+dy2)*ds
1539	IF ( rr .LT. rmax ) THEN
1540	alpha = 1.
1541	ELSE IF ( rr .GE. rmax ) THEN
1542	alpha = alpha2
1543	END IF
1544	vr = vmax * (rr/rmax)**(alpha)
1545	IF ( dx.GE.0. ) THEN
1546	ang = (pi/2.) - ATAN2(dy,MAX(dx,1.e-6))
1547	uu(k) = vwgt(k)(-vrCOS(ang))
1548	vv(k) = vwgt(k)( vrSIN(ang))
1549	ELSE IF ( dx.LT.0. ) THEN
1550	ang = ((3.*pi)/2.) + ATAN2(dy,dx)
1551	uu(k) = vwgt(k)(-vrCOS(ang))
1552	vv(k) = vwgt(k)(-vrSIN(ang))
1553	END IF
1554	END DO
1555
1556	! psi
1557
1558	DO k=1,nlvl
1559	rr = SQRT(dx2+dy2)*ds
1560	IF ( rr .LT. rmax ) THEN
1561	psi(k) = vwgt(k) * (vmaxrrrr)/(2.*rmax)
1562	ELSE IF ( rr .GE. rmax ) THEN
1563	IF (alpha1.EQ.1.0 .AND. alpha2.eq.-1.0) THEN
1564	psi(k) = vwgt(k) * vmaxrmax(0.5+LOG(rr/rmax))
1565	ELSE IF (alpha1.EQ.1.0 .AND. alpha2.NE.-1.0) THEN
1566	term1 = vmax/(rmax*alpha1)(rmax**(alpha1+1)/(alpha1+1))
1567	bb = (rr(alpha2+1)/(alpha2+1))-(rmax(alpha2+1))/(alpha2+1)
1568	term2 = vmax/(rmax*alpha2)bb
1569	psi(k) = vwgt(k) * (term1 + term2)
1570	END IF
1571	END IF
1572	END DO
1573
1574	! vort
1575
1576	DO k=1,nlvl
1577	rr = SQRT(dx2+dy2)*ds
1578	IF ( rr .LT. rmax ) THEN
1579	vor(k) = vwgt(k) * (2.*vmax)/rmax
1580	ELSE IF ( rr .GE. rmax ) THEN
1581	vor(k) = vwgt(k) * ( (vmax/rmax*alpha2)(rr*(alpha2-1.))(1.+alpha2) )
1582	END IF
1583	END DO
1584
1585	END SUBROUTINE rankine
1586
1587	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1588	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1589
1590	SUBROUTINE vor(uin,vin,msfu,msfv,msfm,ew,ns,nz,ds,vort)
1591
1592	!Here we assume that the U and V's are still on the WRF staggered grid.
1593	!The vorticity is then calculated at the mass points on the WRF grid.
1594
1595
1596	IMPLICIT NONE
1597
1598	INTEGER :: jp1,jm1,ip1,im1,i,j,k
1599	INTEGER :: ns, ew, nz, k1
1600
1601	REAL , DIMENSION(ew,nz,ns-1) :: uin !u values on unstaggered U grid
1602	REAL , DIMENSION(ew-1,nz,ns) :: vin !v values on unstaggered V grid
1603	REAL , DIMENSION(ew-1,nz,ns-1) :: vort !vort is defined on the mass points
1604
1605	REAL , DIMENSION(ew,ns-1) :: msfu !map scale factors on U staggered grid
1606	REAL , DIMENSION(ew-1,ns) :: msfv !map scale factors on V staggered grid
1607	REAL , DIMENSION(ew-1,ns-1) :: msfm !map scale factors on unstaggered grid
1608
1609	real :: u(ew,ns-1),v(ew-1,ns)
1610
1611
1612	REAL :: ds
1613
1614	REAL :: dsx,dsy , u1 , u2 , u3 , u4 , v1 , v2 , v3 , v4
1615	real :: dudy,dvdx,mm
1616
1617
1618	vort(:,:,:) = -999.
1619	do k = 1,nz
1620
1621	do j = 1,ns-1
1622	do i = 1,ew
1623	u(i,j) = uin(i,k,j)
1624	end do
1625	end do
1626
1627
1628	do j = 1,ns
1629	do i = 1,ew-1
1630	v(i,j) = vin(i,k,j)
1631	end do
1632	end do
1633
1634	!Our indicies are from 2 to ns-2 and ew-2. This is because out
1635	!map scale factors are not defined for the entire grid.
1636	do j = 2,ns-2
1637	do i = 2,ew-2
1638	mm = msfm(i,j) * msfm(i,j)
1639	u1 = u(i ,j-1)/msfu(i ,j-1)
1640	u2 = u(i+1,j-1)/msfu(i+1,j-1)
1641	u3 = u(i+1,j+1)/msfu(i+1,j+1)
1642	u4 = u(i ,j+1)/msfu(i ,j+1)
1643	dudy = mm * (u4 + u3 -(u1 + u2)) /(4*ds)
1644
1645	v1 = v(i-1,j )/msfv(i-1,j)
1646	v2 = v(i+1,j )/msfv(i+1,j)
1647	v3 = v(i-1 ,j+1)/msfv(i-1,j+1)
1648	v4 = v(i+1,j+1)/msfv(i+1,j+1)
1649	dvdx = mm * (v4 + v2 - (v1 + v3))/(4*ds)
1650
1651	vort(i,k,j) = dvdx - dudy
1652	end do
1653	end do
1654	!Our vorticity array goes out to ew-1 and ns-1 which is the
1655	!mass point grid dimensions.
1656	do i = 2,ew-2
1657	vort(i,k,1) = vort(i,k,2) !bottom not corners
1658	vort(i,k,ns-1) = vort(i,k,ns-2) !top not corners
1659	end do
1660
1661	do j = 1,ns-1
1662	vort(ew-1,k,j) = vort(ew-2,k,j) !right side including corners
1663	vort(1,k,j) = vort(2,k,j) !left side including corners
1664	end do
1665
1666	end do ! this is the k loop end
1667
1668	END SUBROUTINE
1669
1670	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1671	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1672
1673	SUBROUTINE diverg(uin,vin,msfu,msfv,msfm,ew,ns,nz,ds,div)
1674
1675	! Computes divergence on unstaggered grid. The divergence is calculated
1676	! at the mass points on the WRF grid.
1677	! div = m*m (du/dx + dv/dy)
1678
1679	IMPLICIT NONE
1680
1681	INTEGER :: jp1,jm1,ip1,im1,i,j,k
1682	INTEGER :: ns, ew, nz, k1
1683
1684	REAL , DIMENSION(ew,nz,ns-1) :: uin !u values on unstaggered U grid
1685	REAL , DIMENSION(ew-1,nz,ns) :: vin !v values on unstaggered V grid
1686	REAL , DIMENSION(ew-1,nz,ns-1) :: div !divergence is calculate on the mass points
1687	REAL , DIMENSION(ew,ns-1) :: msfu !map scale factors on U staggered grid
1688	REAL , DIMENSION(ew-1,ns) :: msfv !map scale factors on V staggered grid
1689	REAL , DIMENSION(ew-1,ns-1) :: msfm !map scale factors on unstaggered grid
1690
1691	real :: u(ew,ns-1),v(ew-1,ns)
1692
1693
1694	REAL :: ds
1695
1696	REAL :: dsr,u1,u2,v1,v2
1697	real :: dudx,dvdy,mm,arg1,arg2
1698
1699	dsr = 1/ds
1700	do k = 1,nz
1701
1702	do j = 1,ns-1
1703	do i = 1,ew
1704	u(i,j) = uin(i,k,j)
1705	end do
1706	end do
1707
1708
1709	do j = 1,ns
1710	do i = 1,ew-1
1711	v(i,j) = vin(i,k,j)
1712	end do
1713	end do
1714	!Our indicies are from 2 to ns-2 and ew-2. This is because out
1715	!map scale factors are not defined for the entire grid.
1716	do j = 2,ns-2
1717	do i = 2,ew-2
1718	mm = msfm(i,j) * msfm(i,j)
1719	u1 = u(i+1,j)/msfu(i+1,j)
1720	u2 = u(i ,j)/msfu(i ,j)
1721
1722	v1 = v(i,j+1)/msfv(i,j+1)
1723	v2 = v(i,j) /msfv(i,j)
1724
1725	div(i,k,j) = mm * (u1 - u2 + v1 - v2) * dsr
1726	end do
1727	end do
1728
1729	!Our divergence array is defined on the mass points.
1730	do i = 2,ew-2
1731	div(i,k,1) = div(i,k,2) !bottom not corners
1732	div(i,k,ns-1) = div(i,k,ns-2) !top not corners
1733	end do
1734
1735	do j = 1,ns-1
1736	div(ew-1,k,j) = div(ew-2,k,j) !right side including corners
1737	div(1,k,j) = div(2,k,j) !left side including corners
1738	end do
1739
1740	end do !end for the k loop
1741
1742	END SUBROUTINE diverg
1743
1744	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1745	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1746
1747	SUBROUTINE mxratprs (rh, t, ppa, ew, ns, nz, q, min_RH_value)
1748
1749
1750	IMPLICIT NONE
1751
1752	INTEGER :: i , ew , j , ns , k , nz
1753
1754
1755	REAL :: min_RH_value
1756	REAL :: ppa(ew-1,nz,ns-1)
1757	REAL :: p( ew-1,nz,ns-1 )
1758	REAL :: q (ew-1,nz,ns-1),rh(ew-1,nz,ns-1),t(ew-1,nz,ns-1)
1759
1760	REAL :: es
1761	REAL :: qs
1762	REAL :: cp = 1004.0
1763	REAL :: svp1,svp2,svp3
1764	REAL :: celkel
1765	REAL :: eps
1766
1767
1768	! This function is designed to compute (q) from basic variables
1769	! p (mb), t(K) and rh(0-100%) to give (q) in (kg/kg).
1770
1771
1772	p = ppa * 0.01
1773
1774	DO j = 1, ns - 1
1775	DO k = 1, nz
1776	DO i = 1, ew - 1
1777	rh(i,k,j) = MIN ( MAX ( rh(i,k,j) ,min_RH_value ) , 100. )
1778	END DO
1779	END DO
1780	END DO
1781
1782	svp3 = 29.65
1783	svp1 = 0.6112
1784	svp2 = 17.67
1785	celkel = 273.15
1786	eps = 0.622
1787
1788	DO j = 1, ns-1
1789	DO k = 1, nz
1790	DO i = 1,ew-1
1791	es = svp1 * 10. * EXP(svp2 * (t(i,k,j) - celkel ) / (t(i,k,j) - svp3 ))
1792	qs = eps * es / (p(i,k,j) - es)
1793	q(i,k,j) = MAX(0.01 * rh(i,k,j) * qs,0.0)
1794	END DO
1795	END DO
1796	END DO
1797
1798	END SUBROUTINE mxratprs
1799
1800	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1801	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1802	SUBROUTINE mass2_Ustag(field,dim1,dim2,dim3)
1803
1804	IMPLICIT NONE
1805
1806	INTEGER :: dim1 , dim2 , dim3
1807	REAL , DIMENSION(dim1,dim2,dim3) :: field,dummy
1808
1809	dummy = 0.0
1810	dummy(:,2:dim2-1,:) = ( field(:,1:dim2-2,:) + &
1811	field(:,2:dim2-1,:) ) * 0.5
1812	dummy(:,1,:) = field(:,1,:)
1813	dummy(:,dim2,:) = field(:,dim2-1,:)
1814
1815	field = dummy
1816
1817	END SUBROUTINE mass2_Ustag
1818
1819	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1820	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1821	SUBROUTINE mass2_Vstag(field,dim1,dim2,dim3)
1822
1823	IMPLICIT NONE
1824
1825	INTEGER :: dim1 , dim2 , dim3
1826	REAL , DIMENSION(dim1,dim2,dim3) :: field,dummy
1827
1828	dummy = 0.0
1829	dummy(2:dim1-1,:,:) = ( field(1:dim1-2,:,:) + &
1830	field(2:dim1-1,:,:) ) * 0.5
1831	dummy(1,:,:) = field(1,:,:)
1832	dummy(dim1,:,:) = field(dim1-1,:,:)
1833
1834	field = dummy
1835
1836	END SUBROUTINE mass2_Vstag
1837
1838
1839	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1840	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1841
1842	SUBROUTINE relax (chi, ff, rd, ew, ns, ds, smallres, alpha)
1843
1844	IMPLICIT NONE
1845
1846	INTEGER, PARAMETER :: mm = 20000
1847
1848	INTEGER :: i
1849	INTEGER :: ie
1850	INTEGER :: ew !ew direction
1851	INTEGER :: iter
1852	INTEGER :: j
1853	INTEGER :: je
1854	INTEGER :: jm
1855	INTEGER :: ns !ns direction
1856	INTEGER :: mi
1857
1858	REAL :: alpha
1859	REAL :: alphaov4
1860	REAL :: chi(ew-1,ns-1)
1861	REAL :: chimx(ns-1)
1862	REAL :: ds
1863	REAL :: epx
1864	REAL :: fac
1865	REAL :: ff(ew-1,ns-1)
1866	REAL :: rd(ew-1,ns-1)
1867	REAL :: rdmax(ns-1)
1868	REAL :: smallres
1869
1870	LOGICAL :: converged = .FALSE.
1871
1872	fac = ds * ds
1873	alphaov4 = alpha * 0.25
1874
1875	ie=ew-2
1876	je=ns-2
1877
1878	DO j = 1, ns-1
1879	DO i = 1, ew-1
1880	ff(i,j) = fac * ff(i,j)
1881	rd(i,j) = 0.0
1882	END DO
1883	END DO
1884
1885	iter_loop : DO iter = 1, mm
1886	mi = iter
1887	chimx = 0.0
1888
1889
1890	DO j = 2, ns-1
1891	DO i = 2, ew-1
1892	chimx(j) = MAX(ABS(chi(i,j)),chimx(j))
1893	END DO
1894	END DO
1895
1896	epx = MAXVAL(chimx) * SMALLRES * 4.0 / alpha
1897
1898	DO j = 2, ns-2
1899	DO i = 2, ew-2
1900	rd(i,j) = chi(i,j+1) + chi(i,j-1) + chi(i+1,j) + chi(i-1,j) - 4.0 * chi(i,j) - ff(i,j)
1901	chi(i,j) = chi(i,j) + rd(i,j) * alphaov4
1902	END DO
1903	END DO
1904
1905	rdmax = 0.0
1906
1907	DO j = 2, ns-2
1908	DO i = 2, ew-2
1909	rdmax(j) = MAX(ABS(rd(i,j)),rdmax(j))
1910	END DO
1911	END DO
1912
1913
1914	IF (MAXVAL(rdmax) .lt. epx) THEN
1915	converged = .TRUE.
1916	EXIT iter_loop
1917	END IF
1918
1919	END DO iter_loop
1920
1921	IF (converged ) THEN
1922	! PRINT '(A,I5,A)','Relaxation converged in ',mi,' iterations.'
1923	ELSE
1924	PRINT '(A,I5,A)','Relaxation did not converge in',mm,' iterations.'
1925	STOP 'no_converge'
1926	END IF
1927
1928
1929	do i = 2,ew-2
1930	chi(i,ns-1) = chi(i,ns-2) !top not including corners
1931	chi(i,1) = chi(i,2) !bottom not including corners
1932	end do
1933
1934	do j = 2,ns-2
1935	chi(ew-1,j) = chi(ew-2,j) !right side not including corners
1936	chi(1,j) = chi(2,j) !left side not including corners
1937	end do
1938
1939	!Fill in the corners
1940	chi(1,1) = chi(2,1)
1941	chi(ew-1,1) = chi(ew-2,1)
1942	chi(1,ns-1) = chi(2,ns-1)
1943	chi(ew-1,ns-1) = chi(ew-2,ns-1)
1944
1945
1946
1947	END SUBROUTINE relax
1948	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1949	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1950	SUBROUTINE geowind(height,ew,ns,nz,ds,ug,vg)
1951
1952	IMPLICIT NONE
1953
1954	! input height geopotential on wrf mass grid points
1955	! ns wrf staggered V dimension n-s
1956	! ew wrf staggered U dimension e-w
1957	! nz number of vertical levels
1958	!
1959	! output ug u component of geo wind on wrf staggered V points
1960	! vg v component of geo wind on wrf staggered U points
1961
1962	INTEGER :: ew , ns , nz
1963	REAL :: ds
1964	REAL , DIMENSION(ew-1,nz,ns-1) :: height
1965	REAL , DIMENSION(ew,nz,ns-1) :: ug
1966	REAL , DIMENSION(ew-1,nz,ns) :: vg
1967
1968	REAL :: ds2r , h1 , h2 , h3 , h4, ds4r
1969	INTEGER :: i , j , k
1970
1971	ds4r=1./(4.*ds)
1972
1973	! The height field comes in on the WRF mass points.
1974
1975
1976
1977	! ug is the derivative of height in the ns direction ug = -dheight/dy
1978	ug(:,:,:) = -999.
1979	do j=2,ns-2
1980	do k=1,nz
1981	do i=2,ew-1
1982	h1 = height(i,k,j+1)
1983	h2 = height(i-1,k,j+1)
1984	h3 = height(i ,k,j-1)
1985	h4 = height(i-1,k,j-1)
1986	ug(i,k,j) = -( (h1 + h2) - ( h3 + h4) ) * ds4r
1987	end do
1988	end do
1989	end do
1990
1991	do i = 2,ew-1
1992	ug(i,:,1) = ug(i,:,2) !bottom not including corner points
1993	ug(i,:,ns-1) = ug(i,:,ns-2) !top not including corner points
1994	end do
1995
1996	do j = 2,ns-2
1997	ug(1,:,j) = ug(2,:,j) !left side
1998	ug(ew,:,j) = ug(ew-1,:,j) !right side
1999	end do
2000
2001	ug(1,:,1) = ug(2,:,1) !Lower left hand corner
2002	ug(1,:,ns-1) = ug(2,:,ns-1) !upper left hand corner
2003	ug(ew,:,1) = ug(ew-1,:,1) !Lower right hand corner
2004	ug(ew,:,ns-1) = ug(ew-1,:,ns-1) !upper right hand corner
2005
2006
2007	! ug is the derivative of height in the ns direction vg = dheight/dx
2008	vg(:,:,:) = -999.
2009	DO j=2,ns-1
2010	DO k=1,nz
2011	DO i=2,ew-2
2012	h1 = height(i+1,k,j)
2013	h2 = height(i-1,k,j)
2014	h3 = height(i+1,k,j-1)
2015	h4 = height(i-1,k,j-1)
2016	vg(i,k,j) = ( (h1 + h3) - ( h2 + h4) ) * ds4r
2017	end do
2018	end do
2019	end do
2020
2021	do i = 2,ew-2
2022	vg(i,:,1) = vg(i,:,2) !bottom not including corner points
2023	vg(i,:,ns) = vg(i,:,ns-1) !top not including corner points
2024	end do
2025
2026	do j = 2,ns-1
2027	vg(1,:,j) = vg(2,:,j) !left side not including corner points
2028	vg(ew-1,:,j) = vg(ew-2,:,j) !right side not including corner points
2029	end do
2030
2031	vg(1,:,1) = vg(2,:,1) !Lower left hand corner
2032	vg(1,:,ns) = vg(2,:,ns) !upper left hand corner
2033	vg(ew-1,:,1) = vg(ew-2,:,1) !Lower right hand corner
2034	vg(ew-1,:,ns) = vg(ew-2,:,ns) !upper right hand corner
2035
2036
2037	END SUBROUTINE geowind
2038	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2039	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2040
2041	SUBROUTINE balance (f,psi,ew,ns,ds,out)
2042
2043	! Calculates the forcing terms in balance equation
2044
2045	IMPLICIT NONE
2046
2047	! f coriolis force
2048	! psi stream function
2049	! ew, ns grid points in east west, north south direction, respectively
2050	! ds grid distance
2051	! out output array
2052
2053	INTEGER :: ew , ns,nslast,ewlast,ifill
2054	REAL , DIMENSION(ew-1,ns-1) :: f,psi,out
2055	REAL :: ds
2056
2057	REAL :: psixx , psiyy , psiy , psix, psixy
2058	REAL :: dssq , ds2 , dssq4,arg1,arg2,arg3,arg4
2059
2060	INTEGER :: i , j
2061
2062	dssq = ds * ds
2063	ds2 = ds * 2.
2064	dssq4 = ds * ds * 4.
2065
2066	!The forcing terms are calculated on the WRF mass points.
2067	out(:,:) = -999.0
2068	DO j=2,ns-2
2069	DO i=2,ew-2
2070	psiyy = ( psi(i,j+1) + psi(i,j-1) - 2.*psi(i,j) ) / dssq
2071	psixx = ( psi(i+1,j) + psi(i-1,j) - 2.*psi(i,j) ) / dssq
2072	psiy = ( psi(i,j+1) - psi(i,j-1) ) / ds2
2073	psix = ( psi(i+1,j) - psi(i-1,j) ) / ds2
2074	psixy = ( psi(i+1,j+1)+psi(i-1,j-1)-psi(i-1,j+1)-psi(i+1,j-1)) / dssq4
2075
2076	arg1 = f(i,j)* (psixx+psiyy)
2077	arg2 = ( ( f(i,j+1) - f(i,j-1)) / ds2 ) * psiy
2078	arg3 = ( ( f(i+1,j) - f(i-1,j)) / ds2 ) * psix
2079	arg4 = 2 (psixypsixy-psixx*psiyy)
2080	out(i,j)= arg1 + arg2 + arg3 - arg4
2081	END DO
2082	END DO
2083
2084	do i = 2,ew-2
2085	out(i,ns-1) = out(i,ns-2) !top not including corners
2086	out(i,1) = out(i,2) !bottom not including corners
2087	end do
2088
2089	do j = 2,ns-2
2090	out(ew-1,j) = out(ew-2,j) !right side not including corners
2091	out(1,j) = out(2,j) !left side not including corners
2092	end do
2093
2094	!Fill in the corners
2095	out(1,1) = out(2,1)
2096	out(ew-1,1) = out(ew-2,1)
2097	out(1,ns-1) = out(2,ns-1)
2098	out(ew-1,ns-1) = out(ew-2,ns-1)
2099
2100	END SUBROUTINE balance
2101
2102	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2103	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2104
2105	SUBROUTINE qvtorh ( q , t , p , k00, ew , ns , nz , rh, min_RH_value )
2106
2107	IMPLICIT NONE
2108
2109	INTEGER , INTENT(IN) :: ew , ns , nz , k00
2110	REAL , INTENT(IN) , DIMENSION(ew-1,nz,ns-1) :: q ,t, p
2111	REAL , INTENT(OUT) , DIMENSION(ew-1,nz,ns-1) :: rh
2112
2113	real min_RH_value
2114
2115	! Local variables.
2116
2117	INTEGER :: i , j , k,fill
2118	REAL :: es
2119	REAL :: qs
2120	REAL :: cp = 1004.0
2121	REAL :: svp1,svp2,svp3
2122	REAL :: celkel
2123	REAL :: eps
2124	svp3 = 29.65
2125	svp1 = 0.6112
2126	svp2 = 17.67
2127	celkel = 273.15
2128	eps = 0.622
2129
2130	DO j = 1 , ns - 1
2131	DO k = k00 , nz
2132	DO i = 1 , ew -1
2133	es = svp1 * 10. * EXP(svp2 * (t(i,k,j) - celkel ) / (t(i,k,j) - svp3 ))
2134	qs = epses/(0.01p(i,k,j) - es)
2135	rh(i,k,j) = MIN ( 100. , MAX ( 100.*q(i,k,j)/qs , min_RH_value ) )
2136	END DO
2137	END DO
2138	END DO
2139
2140	END SUBROUTINE qvtorh
2141
2142	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2143	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2144
2145	SUBROUTINE stagger_rankine_winds(utcp,vtcp,ew,ns,nz)
2146
2147
2148	!utcp and vtcp are the output winds of the rankine subroutine
2149	!The winds are calculated on the mass points of the WRF grid
2150	!so they need to be staggered out to the WRF staggering.
2151	!The utcp is placed on the staggered ew wind grid.
2152	!The vtcp is placed on the staggered ns wind grid.
2153
2154	!ew is the full grid dimension in the ew direction.
2155	!ns is the full grid dimension in the ns direction.
2156
2157	!nz is the number of vertical levels.
2158
2159	INTEGER :: ew, ns, nz, i,k,j
2160	REAL utcp(ew,nz,ns-1), vtcp(ew-1,nz,ns)
2161
2162	!----------------------------------------------------
2163	!Stagger UTCP
2164	DO j=1,ns-1
2165	DO i=2,ew-1
2166	DO k=1,nz
2167	utcp(i,k,j) = ( utcp(i-1,k,j) + utcp(i,k,j) ) /2
2168	end do
2169	end do
2170	end do
2171
2172	!Fill in U's along the left and right side.
2173	do j = 1,ns
2174	utcp(1,:,j) = utcp(2,:,j)
2175	utcp(ew,:,j) = utcp(ew-1,:,j)
2176	end do
2177
2178
2179	!Stagger VTCP
2180	DO j=2,ns-1
2181	DO i=1,ew-1
2182	DO k=1,nz
2183	vtcp(i,k,j) = ( vtcp(i,k,j+1) + vtcp(i,k,j-1) ) /2
2184	end do
2185	end do
2186	end do
2187
2188	!Fill in V's along the bottom and bottom.
2189	do i = 1,ew
2190	vtcp(i,:,1) = vtcp(i,:,2)
2191	vtcp(i,:,ns) = vtcp(i,:,ns-1)
2192	end do
2193
2194
2195	END SUBROUTINE stagger_rankine_winds
2196
2197	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2198	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2199
2200	subroutine final_ew_velocity(u2,u1,chi,psi,utcr,dx,ew,ns,nz)
2201
2202
2203	integer :: ew,ns,nz,i,j,k
2204	real :: u1(ew,nz,ns-1),utcr(ew,nz,ns-1)
2205	real :: psi(ew-1,nz,ns-1),chi(ew-1,nz,ns-1)
2206	! input arrays:
2207	! u1 is the original wind coming in from the metgrid file.
2208	! utcr is the rankine winds rotated to the map projection put on u WRF staggered grid points.
2209
2210	! psi comes in on the WRF mass points. psi is the perturbation field
2211	! calculated from the relaxation of the vorticity.
2212
2213	! chi is the relaxation of the divergent winds on WRF mass points.
2214
2215
2216	! ew is the grid dimension of the WRF ew staggered wind
2217	! ns is the grid dimension of the WRF ns staggered wind
2218	! nz is the number of vertical levels
2219	! dx is the grid spacing
2220	!-------------------------------------------------------------------------------------------
2221
2222	real :: u2(ew,nz,ns-1)
2223	! output array u2 is the new wind in the ew direction. Is is on WRF staggering.
2224	!-------------------------------------------------------------------------------------------
2225
2226	real upos(ew,nz,ns-1),u0(ew,nz,ns-1),uchi(ew,nz,ns-1)
2227	! upos is the derivative of psi in the ns direction u = -dpsi/dy
2228	! u0 is the background ew velocity
2229	! uchi is the array chi on the u staggered grid.
2230
2231	real :: dx,arg1,arg2
2232
2233	!-------------------------------------------------------------
2234	!Take the derivative of chi in the ew direction.
2235	uchi(:,:,:) = -999.
2236	DO k=1,nz !start of k loop
2237	DO j=1,ns-1
2238	DO i=2,ew-1
2239	uchi(i,k,j) = ( chi(i,k,j) - chi(i-1,k,j) )/dx
2240	END DO
2241	END DO
2242
2243	do j = 1,ns-1
2244	uchi(1,k,j) = uchi(2,k,j) !fill in the left side
2245	uchi(ew,k,j) = uchi(ew-1,k,j) !fill in the right side
2246	end do
2247	end do !k loop
2248
2249	!-----------------------------------------------------------------------------------------
2250	! Take the derivative of psi in the ns direction.
2251	upos = - dpsi/dy
2252	upos(:,:,:) = -999.
2253	DO k=1,nz
2254
2255	DO j=2,ns-2
2256	DO i=2,ew-1
2257	arg1 = psi(i,k,j+1) + psi(i-1,k,j+1)
2258	arg2 = psi(i,k,j-1) + psi(i-1,k,j-1)
2259	upos(i,k,j) = -( arg1 - arg2 )/(4.*dx)
2260	END DO
2261	END DO
2262
2263	do i = 2,ew-1
2264	upos(i,k,1) = upos(i,k,2) !bottom not including corner points
2265	upos(i,k,ns-1) = upos(i,k,ns-2) !top not including corner points
2266	end do
2267
2268	do j = 1,ns-2
2269	upos(1,k,j) = upos(2,k,j) !left side not including corners
2270	upos(ew,k,j) = upos(ew-1,k,j) !right side not including corners
2271	end do
2272
2273
2274	upos(1,k,1) = upos(2,k,1) !Lower left hand corner
2275	upos(1,k,ns-1) = upos(2,k,ns-1) !upper left hand corner
2276	upos(ew,k,1) = upos(ew-1,k,1) !Lower right hand corner
2277	upos(ew,k,ns-1) = upos(ew-1,k,ns-1) !upper right hand corner
2278
2279	end do !k loop for dspi/dy
2280
2281
2282
2283	!-----------------------------------------------------------------------------------------
2284
2285	! Background u field.
2286	! Subtract the first quess wind field from the original winds.
2287	do j=1,ns-1
2288	do k=1,nz
2289	do i=1,ew
2290	u0(i,k,j) = u1(i,k,j)-(upos(i,k,j)+uchi(i,k,j))
2291	end do
2292	end do
2293	end do
2294
2295
2296	! Calculate the final velocity
2297	do j=1,ns-1
2298	do k=1,nz
2299	do i=1,ew
2300	u2(i,k,j) = u0(i,k,j)+utcr(i,k,j)
2301	end do
2302	end do
2303	end do
2304
2305	end subroutine final_ew_velocity
2306
2307	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2308	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2309
2310	subroutine final_ns_velocity(v2,v1,chi,psi,vtcr,dx,ew,ns,nz)
2311
2312
2313	integer :: ew,ns,nz,i,j,k
2314	real :: v1(ew-1,nz,ns),vtcr(ew-1,nz,ns)
2315	real :: psi(ew-1,nz,ns-1),chi(ew-1,nz,ns-1)
2316	! input arrays:
2317	! v1 is the original wind coming in from the metgrid file.
2318	! vtcr is the is the rankine winds rotated to the map projection put on v WRF staggered grid points.
2319
2320	! psi comes on the WRF mass points. psi is the perturbation field
2321	! calculated from the relaxation of the vorticity.
2322
2323	! chi is the relaxation of the divergent winds on WRF mass points.
2324
2325	! ew is the grid dimension of the WRF ew staggered wind
2326	! ns is the grid dimension of the WRF ns staggered wind
2327	! nz is the number of vertical levels
2328
2329
2330	real :: v2(ew-1,nz,ns)
2331	! output array v2 is the new wind in the ns direction. Is is on WRF staggering.
2332
2333
2334	real vpos(ew-1,nz,ns),v0(ew-1,nz,ns),vchi(ew-1,nz,ns)
2335	! vpos is the derivative of psi in the ew direction v = dpsi/dx
2336	! v0 is the background ns velocity
2337	! vchi is the relaxation of the divergent wind put on v WRF staggered grid points.
2338
2339	real :: dx,arg1,arg2
2340
2341
2342	!-----------------------------------------------------------------------------------------
2343	vchi(:,:,:) = -999.0
2344	!The derivative dchi in the ns direction.
2345	do k = 1,nz
2346	DO j=2,ns-1
2347	DO i=1,ew-1
2348	vchi(i,k,j) = ( chi(i,k,j) - chi(i,k,j-1))/dx
2349	END DO
2350	END DO
2351
2352	do i = 1,ew-1
2353	vchi(i,k,1) = vchi(i,k,2)
2354	vchi(i,k,ns) = vchi(i,k,ns-1)
2355	end do
2356
2357	end do !end of k loop
2358
2359	!-----------------------------------------------------------------------------------------
2360	!Take the derivative of psi in the ew direction.
2361	vpos(:,:,:) = -999.
2362
2363	DO k=1,nz
2364	DO j=2,ns-1
2365	DO i=2,ew-2
2366	arg1 = psi(i+1,k,j) + psi(i+1,k,j-1)
2367	arg2 = psi(i-1,k,j) + psi(i-1,k,j-1)
2368	vpos(i,k,j) = ( arg1 - arg2 )/(4.*dx)
2369	END DO
2370	END DO
2371
2372	do i = 2,ew-2
2373	vpos(i,k,1) = vpos(i,k,2) !bottom not including corner points
2374	vpos(i,k,ns) = vpos(i,k,ns-1) !top not including corner points
2375	end do
2376
2377	do j = 1,ns
2378	vpos(1,k,j) = vpos(2,k,j) !left side not including corner points
2379	vpos(ew-1,k,j) = vpos(ew-2,k,j) !right side not including corner points
2380	end do
2381
2382
2383	vpos(1,k,1) = vpos(2,k,1) !Lower left hand corner
2384	vpos(1,k,ns) = vpos(2,k,ns) !upper left hand corner
2385	vpos(ew-1,k,1) = vpos(ew-2,k,1) !Lower right hand corner
2386	vpos(ew-1,k,ns) = vpos(ew-2,k,ns) !upper right hand corner
2387
2388	END DO!k loop for dspi/dx
2389
2390
2391	do j=1,ns
2392	do k=1,nz
2393	do i=1,ew-1
2394	v0(i,k,j) = v1(i,k,j)-(vpos(i,k,j)+vchi(i,k,j))
2395	if( v0(i,k,j) .gt. 100.) then
2396	print *,vchi(i,k,j),i,k,j
2397	stop
2398	end if
2399	end do
2400	end do
2401	end do
2402
2403
2404	! Calculate the final velocity
2405	do j=1,ns
2406	do k=1,nz
2407	do i=1,ew-1
2408	v2(i,k,j) = v0(i,k,j)+vtcr(i,k,j)
2409	end do
2410	end do
2411	end do
2412
2413	end subroutine final_ns_velocity
2414	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2415	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2416	subroutine final_RH(rh2,rh0,rhmx,strmci,strmcj,rmax_nstrm,ew,ns,nz,k00, &
2417	dx,ew_gcntr,ns_gcntr,r_vor2)
2418
2419
2420
2421	integer :: ew,ns,nz
2422	real :: rh2(ew-1,nz,ns-1) !The final output relative humidity.
2423	real :: rh0(ew-1,nz,ns-1) !First quess rh read from the metgrid input file.
2424	real :: rhmx(nz)
2425	real :: ew_gcntr !ew grid center as returned from the map projection routines.
2426	real :: ns_gcntr !ns grid center as returned from the map projection routines.
2427	real :: dx !grid spacing
2428	real :: rmax_nstrm
2429
2430
2431	!Local real variables
2432	real :: sum_rh,avg_rh,rh_min,rhbkg,rhbog,r_ratio
2433	real :: rad
2434	real :: rhtc(ew-1,nz,ns-1)
2435
2436	integer :: nct,k00,i,j,k,ew_mvc,ns_mvc
2437	integer :: strmci(nz), strmcj(nz)
2438
2439
2440	!-----------------------------------------------------------------------
2441	DO k=k00,nz
2442	rh_max= rhmx(k)
2443	ew_mvc = strmci(k)
2444	ns_mvc = strmcj(k)
2445
2446
2447	sum_rh = 0.
2448	nct = 0
2449	DO j=1,ns-1
2450	DO i=1,ew-1
2451	rad = SQRT(REAL(i-ew_mvc)2.+REAL(j-ns_mvc)2.)*dx
2452	IF ( (rad .LT. r_vor2).AND.(rad .GE. 0.8*r_vor2) ) THEN
2453	sum_rh = sum_rh + rh0(i,k,j)
2454	nct = nct + 1
2455	END IF
2456	END DO
2457	END DO
2458	avg_rh = sum_rh/MAX(REAL(nct),1.)
2459
2460	DO j=1,ns-1
2461	DO i=1,ew-1
2462	rh_min = avg_rh
2463	rad = SQRT((REAL(i)-ew_gcntr)2.+(REAL(j)-ns_gcntr)2.)*dx
2464	IF ( rad .LE. rmax_nstrm ) THEN
2465	rhtc(i,k,j) = rh_max
2466	ELSE
2467	rhtc(i,k,j) = (rmax_nstrm/rad)rh_max+(1.-(rmax_nstrm/rad))rh_min
2468	END IF
2469	END DO
2470	END DO
2471	END DO
2472
2473
2474	! New RH.
2475	DO j=1,ns-1
2476	DO k=k00,nz
2477	DO i=1,ew-1
2478	rhbkg = rh0(i,k,j)
2479	rhbog = rhtc(i,k,j)
2480	rad = SQRT((REAL(i)-ew_mvc)2.+(REAL(j)-ns_mvc)2.)*dx
2481	IF ( (rad.GT.rmax_nstrm) .AND. (rad.LE.r_vor2) ) THEN
2482	r_ratio = (rad-rmax_nstrm)/(r_vor2-rmax_nstrm)
2483	rh2(i,k,j) = ((1.-r_ratio)rhbog) + (r_ratiorhbkg)
2484	ELSE IF (rad .LE. rmax_nstrm ) THEN
2485	rh2(i,k,j) = rhbog
2486	ELSE
2487	rh2(i,k,j) = rhbkg
2488	END IF
2489
2490	END DO
2491	END DO
2492	END DO
2493
2494
2495
2496	end subroutine final_RH
2497
2498	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2499	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

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

Download in other formats:

Original Format