Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

tc_em.F @ 13

Last change on this file since 13 was 1, checked in by lfita, 11 years ago

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

WRF: version v3.3
LMDZ: version v1818

More details in:

File size: 78.0 KB

Rev	Line
[1]	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