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