da_setup_firstguess_wrf.inc
References to this file elsewhere.
1 subroutine da_setup_firstguess_wrf(xbx, grid)
2
3 !---------------------------------------------------------------------------
4 ! Purpose: Define/allocate components of WRF model state.
5 !---------------------------------------------------------------------------
6
7 implicit none
8
9 type (xbx_type), intent(out) :: xbx ! Header & non-gridded vars.
10
11 type (domain), intent(inout) :: grid
12
13 integer :: map_util_project
14 real :: x, y, xxc, yyc, xxx, yyy, lat_cen, lon_cen
15
16 real :: conv
17 real :: buf(2)
18
19 character(len=24) :: xb_date, an_date
20 integer :: len, seconds, i_grid, j_grid, m_expand
21
22 if (trace_use) call da_trace_entry("da_setup_firstguess_wrf")
23
24 !-----------------------------------------------------------------------
25 ! [0.0] check the xb_date for 3DVAR
26 !-----------------------------------------------------------------------
27
28 if (.NOT. var4d) then
29 write(unit=xb_date,fmt='(i4.4,2("-",i2.2),"_",i2.2,2(":",i2.2),".0000")') &
30 grid%start_year, grid%start_month, grid%start_day, &
31 grid%start_hour, grid%start_minute,grid%start_second
32
33 len = len_trim(ANALYSIS_DATE)
34
35 write(unit=an_date(1:len), fmt='(a)') trim(ANALYSIS_DATE)
36
37 seconds = int(da_diff_seconds(an_date, xb_date))
38
39 if (seconds > ANALYSIS_ACCU) then
40 write(unit=message(1),fmt='(A,A,A,A)') &
41 "xb_date=",xb_date," an_date=", an_date
42 write(unit=message(2),fmt='(A,I6,A,I6)') &
43 "diff=",seconds," ANALYSIS_ACCU=",ANALYSIS_ACCU
44 message(3)="=======> Wrong xb time found???"
45 call da_warning(__FILE__,__LINE__,message(1:3))
46 end if
47 end if
48
49 !------------------------------------------------------------------------
50 ! [1.0] Read original WRF format first guess:
51 !------------------------------------------------------------------------
52
53 conv = 180.0 / pi
54
55 !------------------------------------------------------------------------
56 ! [2.0] Copy header info:
57 !------------------------------------------------------------------------
58
59 if ((grid%xp%its == grid%xp%ids) .and. (grid%xp%jts == grid%xp%jds)) then
60 buf(1) = grid%xlat(grid%xp%its, grid%xp%jts)
61 buf(2) = grid%xlong(grid%xp%its, grid%xp%jts)
62 end if
63
64 call wrf_dm_bcast_real(buf, 2)
65 start_lat=buf(1)
66 start_lon=buf(2)
67
68 !------------------------------------------------------------------------
69 ! Setup map utility
70 !------------------------------------------------------------------------
71
72 call nl_get_map_proj (grid%id , grid%map_proj)
73 call nl_get_truelat1 (grid%id , grid%truelat1)
74 call nl_get_truelat2 (grid%id , grid%truelat2)
75 call nl_get_dx (grid%id , grid%dx)
76 call nl_get_cen_lat (grid%id , grid%cen_lat)
77 call nl_get_cen_lon (grid%id , grid%cen_lon)
78 call nl_get_moad_cen_lat (grid%id , grid%moad_cen_lat)
79 call nl_get_stand_lon (grid%id , grid%stand_lon)
80
81 phic = grid%moad_cen_lat
82 xlonc = grid%stand_lon
83
84 truelat1_3dv = grid%truelat1
85 truelat2_3dv = grid%truelat2
86 pole = 90.0
87 dsm = 0.001 * grid%dx
88
89 map_util_project = grid%map_proj
90
91 if (print_detail_map) then
92 write(unit=stdout, fmt='(a, i6)') &
93 'map_proj =', grid%map_proj
94
95 write(unit=stdout, fmt='(a, e16.6)') &
96 'cen_lat =', grid%cen_lat, &
97 'cen_lon =', grid%cen_lon, &
98 'truelat1 =', grid%truelat1, &
99 'truelat2 =', grid%truelat2, &
100 'start_lat =', start_lat, &
101 'start_lon =', start_lon, &
102 'dsm =', dsm
103 end if
104
105 ! Set map projection in WRFSI world.
106 map_util_project = PROJ_LC
107
108 if (grid%map_proj == 0) then
109 map_util_project = PROJ_LATLON
110 else if (grid%map_proj == 1) then
111 map_util_project = PROJ_LC
112 else if (grid%map_proj == 2) then
113 map_util_project = PROJ_PS
114 else if (grid%map_proj == 3) then
115 map_util_project = PROJ_MERC
116 end if
117
118 call da_map_set(map_util_project,grid%cen_lat,grid%cen_lon, &
119 real(grid%xp%ide-grid%xp%ids+2)/2., real(grid%xp%jde-grid%xp%jds+2)/2., &
120 grid%dx,grid%stand_lon,grid%truelat1,grid%truelat2,map_info)
121
122 ! Need to set map projection in WRF world.
123 map_projection = grid%map_proj
124
125 cone_factor = map_info%cone
126
127 if (print_detail_map) then
128
129 !----------------------------------------------------------------------
130 ! Check the ll_to_ij:
131 !----------------------------------------------------------------------
132
133 message(1)="Check the map_set correctness::::::::::::::::::::::::"
134
135 ! Domain center:
136 call da_latlon_to_ij(map_info, grid%cen_lat, grid%cen_lon, start_x, start_y)
137 write(unit=message(2),fmt='("Center: latc,lonc,x,y, Xc, Yc:",6f10.3)') &
138 grid%cen_lat, grid%cen_lon, start_x, start_y, &
139 real(grid%xp%ide-grid%xp%ids+2)/2., real(grid%xp%jde-grid%xp%jds+2)/2.
140
141 start_x = real(grid%xp%ide-grid%xp%ids+2)/2
142 start_y = real(grid%xp%jde-grid%xp%jds+2)/2
143 lat_cen = -999.9
144 lon_cen = -999.9
145 call da_ij_to_latlon(map_info, start_x, start_y, lat_cen, lon_cen)
146 write(unit=message(3), &
147 fmt='("Center: X, Y, latc, lonc, phic, xlonc:",6f10.3)') &
148 start_x, start_y, lat_cen, lon_cen, &
149 grid%cen_lat, grid%cen_lon
150 call da_message(message(1:3))
151 end if
152
153 ! Setup the domain definition for use of the GRAPH:
154
155 coarse_ds = 0.001 * grid%dx
156 coarse_ix = grid%e_we - grid%s_we + 1
157 coarse_jy = grid%e_sn - grid%s_sn + 1
158 start_x = 1.0
159 start_y = 1.0
160
161 xxc = real(coarse_ix)/2.0
162 yyc = real(coarse_jy)/2.0
163
164 ! Warn if grid centres differ by more than .1%
165 if (abs(grid%cen_lon-grid%stand_lon)/grid%cen_lon > 0.001 .or. &
166 abs(grid%cen_lat-grid%moad_cen_lat)/grid%cen_lat > 0.001) then
167
168 write(unit=message(1),fmt='(a,2e20.12)')'grid%cen_lon , grid%cen_lat :', &
169 grid%cen_lon, grid%cen_lat
170 write(unit=message(2),fmt='(a,2e20.12)')'grid%stand_lon, grid%moad_cen_lat:', &
171 grid%stand_lon, grid%moad_cen_lat
172 write(unit=message(3),fmt='(a)') &
173 '## Compute the coarse_ix, coarse_jy, start_x, and start_y: ##'
174
175 call da_latlon_to_ij(map_info, phic , xlonc , xxx, yyy)
176
177 i_grid = nint(xxx-xxc)
178 j_grid = nint(yyy-yyc)
179 m_expand = 16
180
181 coarse_ix = coarse_ix + 2*abs(i_grid) + m_expand
182 coarse_jy = coarse_jy + 2*abs(j_grid) + m_expand
183
184 start_x = m_expand/2 + 1
185 if (i_grid < 0) start_x = start_x - 2*i_grid
186 start_y = m_expand/2 + 1
187 if (j_grid < 0) start_y = start_y - 2*j_grid
188
189 write(unit=message(4),fmt='(a,2i5,2x,a,2f6.1)') &
190 'Dimension of MOAD: ix, jy:',coarse_ix, coarse_jy, &
191 'parent_x, parent_y:', start_x, start_y
192 call da_warning(__FILE__,__LINE__,message(1:4))
193 end if
194
195 call da_set_map_para ! set up the map background parameters
196
197 ! call da_llxy(PHIC,XLONC,X,Y)
198 call da_llxy(grid%cen_lat,grid%cen_lon,x,y)
199
200 !--------------------------------------------------------------------------
201 ! [3.0] Interpolate WRF C-grid winds to p points of WRFVAR grid (interpolate
202 ! u to west, v to south?
203 !---------------------------------------------------------------------------
204
205 grid%xb % mix = grid%xp%ide - grid%xp%ids + 1
206 grid%xb % mjy = grid%xp%jde - grid% xp%jds + 1
207 grid%xb % mkz = grid%xp%kde - grid%xp%kds + 1
208
209 grid%xb % ds = 0.001 * grid%dx
210
211 mix = grid%xb % mix
212 mjy = grid%xb % mjy
213 mkz = grid%xb % mkz
214
215 call da_transfer_wrftoxb(xbx, grid)
216
217 if (trace_use) call da_trace_exit("da_setup_firstguess_wrf")
218
219 end subroutine da_setup_firstguess_wrf
220
221