source: trunk/WRF.COMMON/WRFV2/share/mediation_nest_move.F @ 2756

Last change on this file since 2756 was 11, checked in by aslmd, 14 years ago

spiga@svn-planeto:ajoute le modele meso-echelle martien

File size: 34.1 KB
Line 
1
2SUBROUTINE med_nest_move ( parent, nest )
3  ! Driver layer
4   USE module_domain
5   USE module_timing
6   USE module_configure
7   USE module_io_domain
8   USE module_dm
9   TYPE(domain) , POINTER                     :: parent, nest, grid
10   INTEGER dx, dy       ! number of parent domain points to move
11#ifdef MOVE_NESTS
12  ! Local
13   CHARACTER*256 mess
14   INTEGER i, j, p, parent_grid_ratio, dyn_opt
15   INTEGER px, py       ! number and direction of nd points to move
16   INTEGER                         :: ids , ide , jds , jde , kds , kde , &
17                                      ims , ime , jms , jme , kms , kme , &
18                                      ips , ipe , jps , jpe , kps , kpe
19   INTEGER ierr, fid
20   LOGICAL input_from_hires
21   LOGICAL saved_restart_value
22   TYPE (grid_config_rec_type)   :: config_flags
23   LOGICAL, EXTERNAL :: wrf_dm_on_monitor
24   LOGICAL, EXTERNAL :: should_not_move
25
26   INTERFACE
27     SUBROUTINE med_interp_domain ( parent , nest )
28        USE module_domain
29        TYPE(domain) , POINTER                 :: parent , nest
30     END SUBROUTINE med_interp_domain
31     SUBROUTINE start_domain ( grid , allowed_to_move )
32        USE module_domain
33        TYPE(domain) :: grid
34        LOGICAL, INTENT(IN) :: allowed_to_move
35     END SUBROUTINE start_domain
36#if ( EM_CORE == 1 )
37     SUBROUTINE shift_domain_em ( grid, disp_x, disp_y  &
38!
39# include <em_dummy_new_args.inc>
40!
41                           )
42        USE module_domain
43        USE module_configure
44        USE module_timing
45        IMPLICIT NONE
46        INTEGER disp_x, disp_y
47        TYPE(domain) , POINTER                 :: grid
48# include <em_dummy_new_decl.inc>
49     END SUBROUTINE shift_domain_em
50#endif
51#if ( NMM_CORE == 1 )
52#endif
53     LOGICAL FUNCTION time_for_move ( parent , nest , dx , dy )
54        USE module_domain
55        USE module_utility
56        TYPE(domain) , POINTER    :: parent , nest
57        INTEGER, INTENT(OUT)      :: dx , dy
58     END FUNCTION time_for_move
59     SUBROUTINE  input_terrain_rsmas ( grid ,                  &
60                           ids , ide , jds , jde , kds , kde , &
61                           ims , ime , jms , jme , kms , kme , &
62                           ips , ipe , jps , jpe , kps , kpe )
63       USE module_domain
64       TYPE ( domain ) :: grid
65       INTEGER                           :: ids , ide , jds , jde , kds , kde , &
66                                            ims , ime , jms , jme , kms , kme , &
67                                            ips , ipe , jps , jpe , kps , kpe
68     END SUBROUTINE input_terrain_rsmas
69     SUBROUTINE med_nest_feedback ( parent , nest , config_flags )
70       USE module_domain
71       USE module_configure
72       TYPE (domain), POINTER ::  nest , parent
73       TYPE (grid_config_rec_type) config_flags
74     END SUBROUTINE med_nest_feedback
75     SUBROUTINE  blend_terrain ( ter_interpolated , ter_input , &
76                           ids , ide , jds , jde , kds , kde , &
77                           ims , ime , jms , jme , kms , kme , &
78                           ips , ipe , jps , jpe , kps , kpe )
79       INTEGER                           :: ids , ide , jds , jde , kds , kde , &
80                                            ims , ime , jms , jme , kms , kme , &
81                                            ips , ipe , jps , jpe , kps , kpe
82       REAL , DIMENSION(ims:ime,jms:jme) :: ter_interpolated
83       REAL , DIMENSION(ims:ime,jms:jme) :: ter_input
84     END SUBROUTINE blend_terrain
85     SUBROUTINE  store_terrain ( ter_interpolated , ter_input , &
86                           ids , ide , jds , jde , kds , kde , &
87                           ims , ime , jms , jme , kms , kme , &
88                           ips , ipe , jps , jpe , kps , kpe )
89       INTEGER                           :: ids , ide , jds , jde , kds , kde , &
90                                            ims , ime , jms , jme , kms , kme , &
91                                            ips , ipe , jps , jpe , kps , kpe
92       REAL , DIMENSION(ims:ime,jms:jme) :: ter_interpolated
93       REAL , DIMENSION(ims:ime,jms:jme) :: ter_input
94     END SUBROUTINE store_terrain
95   END INTERFACE
96
97  ! set grid pointer for code in deref_kludge (if used)
98   grid => nest
99
100   CALL nl_get_dyn_opt ( 1, dyn_opt )
101
102   IF ( should_not_move( nest%id ) ) THEN
103      CALL wrf_message( 'Nest movement is disabled because of namelist settings' )
104      RETURN
105   ENDIF
106
107! if the nest has stopped don't do all this
108   IF ( WRFU_ClockIsStopTime(nest%domain_clock ,rc=ierr) ) RETURN
109
110! mask should be defined in nest domain
111
112   check_for_move: IF ( time_for_move ( parent , nest , dx, dy ) ) THEN
113
114     IF ( (dx .gt. 1 .or. dx .lt. -1 ) .or.  &
115          (dy .gt. 1 .or. dy .lt. -1 ) ) THEN
116       WRITE(mess,*)' invalid move: dx, dy ', dx, dy
117       CALL wrf_error_fatal( mess )
118     ENDIF
119
120     WRITE(mess,*)' moving ',grid%id,dx,dy
121     CALL wrf_message(mess)
122
123     CALL get_ijk_from_grid (  grid ,                   &
124                               ids, ide, jds, jde, kds, kde,    &
125                               ims, ime, jms, jme, kms, kme,    &
126                               ips, ipe, jps, jpe, kps, kpe    )
127
128     CALL wrf_dm_move_nest ( parent, nest%intermediate_grid, dx, dy )
129
130     CALL adjust_domain_dims_for_move( nest%intermediate_grid , dx, dy )
131
132     CALL get_ijk_from_grid (  grid ,                   &
133                               ids, ide, jds, jde, kds, kde,    &
134                               ims, ime, jms, jme, kms, kme,    &
135                               ips, ipe, jps, jpe, kps, kpe    )
136
137     grid => nest
138
139#if ( EM_CORE == 1 )
140     IF ( dyn_opt .EQ. DYN_EM ) THEN
141       CALL shift_domain_em( grid, dx, dy  &
142!
143# include <em_actual_new_args.inc>
144!
145                           )
146     ENDIF
147#endif
148#if ( WRF_NMM_CORE == 1 )
149     IF ( dyn_opt .EQ. DYN_NMM ) THEN
150     ENDIF
151#endif
152
153     px = grid%parent_grid_ratio*dx
154     py = grid%parent_grid_ratio*dy
155
156     grid%i_parent_start = grid%i_parent_start + px / grid%parent_grid_ratio
157     CALL nl_set_i_parent_start( grid%id, grid%i_parent_start )
158     grid%j_parent_start = grid%j_parent_start + py / grid%parent_grid_ratio
159     CALL nl_set_j_parent_start( grid%id, grid%j_parent_start )
160
161     IF ( wrf_dm_on_monitor() ) THEN
162       write(mess,*)  &
163         'Grid ',grid%id,' New SW corner (in parent x and y):',grid%i_parent_start, grid%j_parent_start
164       CALL wrf_message(TRIM(mess))
165     ENDIF
166
167     CALL med_interp_domain( parent, nest )
168
169     CALL nl_get_input_from_hires( nest%id , input_from_hires )
170     IF ( input_from_hires ) THEN
171
172! store horizontally interpolated terrain in temp location
173       CALL  store_terrain ( nest%ht_fine , nest%ht , &
174                             ids , ide , jds , jde , 1   , 1   , &
175                             ims , ime , jms , jme , 1   , 1   , &
176                             ips , ipe , jps , jpe , 1   , 1   )
177       CALL  store_terrain ( nest%em_mub_fine , nest%em_mub , &
178                             ids , ide , jds , jde , 1   , 1   , &
179                             ims , ime , jms , jme , 1   , 1   , &
180                             ips , ipe , jps , jpe , 1   , 1   )
181       CALL  store_terrain ( nest%em_phb_fine , nest%em_phb , &
182                             ids , ide , jds , jde , kds , kde , &
183                             ims , ime , jms , jme , kms , kme , &
184                             ips , ipe , jps , jpe , kps , kpe )
185
186       CALL  input_terrain_rsmas ( nest,                               &
187                                   ids , ide , jds , jde , 1   , 1   , &
188                                   ims , ime , jms , jme , 1   , 1   , &
189                                   ips , ipe , jps , jpe , 1   , 1   )
190
191       CALL  blend_terrain ( nest%ht_fine , nest%ht , &
192                             ids , ide , jds , jde , 1   , 1   , &
193                             ims , ime , jms , jme , 1   , 1   , &
194                             ips , ipe , jps , jpe , 1   , 1   )
195       CALL  blend_terrain ( nest%em_mub_fine , nest%em_mub , &
196                             ids , ide , jds , jde , 1   , 1   , &
197                             ims , ime , jms , jme , 1   , 1   , &
198                             ips , ipe , jps , jpe , 1   , 1   )
199       CALL  blend_terrain ( nest%em_phb_fine , nest%em_phb , &
200                             ids , ide , jds , jde , kds , kde , &
201                             ims , ime , jms , jme , kms , kme , &
202                             ips , ipe , jps , jpe , kps , kpe )
203!
204       CALL model_to_grid_config_rec ( parent%id , model_config_rec , config_flags )
205
206       CALL med_nest_feedback ( parent , nest , config_flags )
207       parent%imask_nostag = 1
208       parent%imask_xstag = 1
209       parent%imask_ystag = 1
210       parent%imask_xystag = 1
211
212! start_domain will key off "restart". Even if this is a restart run
213! we don't want it to here. Save the value, set it to false, and restore afterwards
214       saved_restart_value = config_flags%restart
215       config_flags%restart = .FALSE.
216       grid%restart = .FALSE.
217       CALL nl_set_restart ( 1, .FALSE. )
218       CALL start_domain ( parent , .FALSE. )
219       config_flags%restart = saved_restart_value
220       grid%restart = saved_restart_value
221       CALL nl_set_restart ( 1,  saved_restart_value )
222
223     ENDIF
224
225!
226! masks associated with nest will have been set by shift_domain_em above
227     nest%moved = .true.
228! start_domain will key off "restart". Even if this is a restart run
229! we don't want it to here. Save the value, set it to false, and restore afterwards
230     saved_restart_value = config_flags%restart
231     config_flags%restart = .FALSE.
232     CALL nl_set_restart ( 1, .FALSE. )
233     grid%restart = .FALSE.
234     CALL start_domain ( nest , .FALSE. )
235     config_flags%restart = saved_restart_value
236     grid%restart = saved_restart_value
237     CALL nl_set_restart ( 1,  saved_restart_value )
238     nest%moved = .false.
239     
240!
241! copy time level 2 to time level 1 in new regions of multi-time level fields
242! this should be registry generated.
243!
244#if ( EM_CORE == 1 )
245     IF ( dyn_opt .EQ. DYN_EM ) THEN
246       do k = kms,kme
247         where ( nest%imask_xstag  .EQ. 1 ) nest%em_u_1(:,k,:)   = nest%em_u_2(:,k,:)
248         where ( nest%imask_ystag  .EQ. 1 ) nest%em_v_1(:,k,:)   = nest%em_v_2(:,k,:)
249         where ( nest%imask_nostag .EQ. 1 ) nest%em_t_1(:,k,:)   = nest%em_t_2(:,k,:)
250         where ( nest%imask_nostag .EQ. 1 ) nest%em_w_1(:,k,:)   = nest%em_w_2(:,k,:)
251         where ( nest%imask_nostag .EQ. 1 ) nest%em_ph_1(:,k,:)  = nest%em_ph_2(:,k,:)
252         where ( nest%imask_nostag .EQ. 1 ) nest%em_tp_1(:,k,:)  = nest%em_tp_2(:,k,:)
253         where ( nest%imask_nostag .EQ. 1 ) nest%em_tke_1(:,k,:) = nest%em_tke_2(:,k,:)
254       enddo
255       where ( nest%imask_nostag .EQ. 1 ) nest%em_mu_1(:,:)  = nest%em_mu_2(:,:)
256     ENDIF
257#endif
258#if ( WRF_NMM_CORE == 1 )
259     IF ( dyn_opt .EQ. DYN_NMM ) THEN
260     ENDIF
261#endif
262!
263   ENDIF check_for_move
264#endif
265END SUBROUTINE med_nest_move
266
267LOGICAL FUNCTION time_for_move2 ( parent , grid , move_cd_x, move_cd_y )
268  ! Driver layer
269   USE module_domain
270   USE module_configure
271   USE module_compute_geop
272   USE module_dm
273   USE module_utility
274   IMPLICIT NONE
275! Arguments
276   TYPE(domain) , POINTER    :: parent, grid
277   INTEGER, INTENT(OUT)      :: move_cd_x , move_cd_y
278#ifdef MOVE_NESTS
279! Local
280   INTEGER  num_moves, rc
281   INTEGER  move_interval , move_id
282   TYPE(WRFU_Time) :: ct, st
283   TYPE(WRFU_TimeInterval) :: ti
284   CHARACTER*256 mess, timestr
285   INTEGER     :: ids, ide, jds, jde, kds, kde, &
286                  ims, ime, jms, jme, kms, kme, &
287                  ips, ipe, jps, jpe, kps, kpe
288   INTEGER :: is, ie, js, je, ierr
289   REAL    :: ipbar, pbar, jpbar, fact
290   REAL    :: last_vc_i , last_vc_j
291
292   REAL, ALLOCATABLE, DIMENSION(:,:) :: height_l, height
293   REAL, ALLOCATABLE, DIMENSION(:,:) :: psfc, xlat, xlong, terrain
294   REAL :: minh, maxh
295   INTEGER :: mini, minj, maxi, maxj, i, j, pgr, irad
296   REAL :: disp_x, disp_y, lag, radius, center_i, center_j, dx
297   REAL :: dijsmooth, vmax, vmin, a, b
298   REAL :: dc_i, dc_j   ! domain center
299   REAL :: maxws, ws
300   REAL :: pmin
301   INTEGER imploc, jmploc
302
303   INTEGER :: fje, fjs, fie, fis, fimloc, fjmloc, imloc, jmloc, dyn_opt
304   INTEGER :: i_parent_start, j_parent_start
305   INTEGER :: max_vortex_speed, vortex_interval  ! meters per second and seconds
306   REAL    :: rsmooth = 100.  ! kilometers?
307
308   LOGICAL, EXTERNAL :: wrf_dm_on_monitor
309
310character*256 message, message2
311
312!#define MOVING_DIAGS
313# ifdef VORTEX_CENTER
314
315
316   CALL nl_get_dyn_opt ( 1, dyn_opt )
317   CALL nl_get_parent_grid_ratio ( grid%id , pgr )
318   CALL nl_get_i_parent_start    ( grid%id , i_parent_start )
319   CALL nl_get_j_parent_start    ( grid%id , j_parent_start )
320
321   CALL get_ijk_from_grid (  grid ,                        &
322                             ids, ide, jds, jde, kds, kde, &
323                             ims, ime, jms, jme, kms, kme, &
324                             ips, ipe, jps, jpe, kps, kpe  )
325
326! If the alarm is ringing, recompute the Vortex Center (VC); otherwise
327! use the previous position of VC.  VC is not recomputed ever step to
328! save on cost for global collection of height field and broadcast
329! of new center.
330
331#  ifdef MOVING_DIAGS
332write(0,*)'Check to see if COMPUTE_VORTEX_CENTER_ALARM is ringing? '
333#  endif
334   CALL nl_get_parent_grid_ratio ( grid%id , pgr )
335   CALL nl_get_dx ( grid%id , dx )
336
337   IF ( WRFU_AlarmIsRinging( grid%alarms( COMPUTE_VORTEX_CENTER_ALARM ), rc=rc ) ) THEN
338
339#  ifdef MOVING_DIAGS
340     write(0,*)'COMPUTE_VORTEX_CENTER_ALARM is ringing  '
341#  endif
342     CALL WRFU_AlarmRingerOff( grid%alarms( COMPUTE_VORTEX_CENTER_ALARM ), rc=rc )
343     CALL domain_clock_get( grid, current_timestr=timestr )
344
345     last_vc_i = grid%vc_i
346     last_vc_j = grid%vc_j
347
348     ALLOCATE ( height_l ( ims:ime , jms:jme ), STAT=ierr )
349     IF ( ierr .ne. 0 ) CALL wrf_error_fatal ('allocating height_l in time_for_move2')
350     IF ( wrf_dm_on_monitor() ) THEN
351       ALLOCATE ( height   ( ids:ide , jds:jde ), STAT=ierr )
352       IF ( ierr .ne. 0 ) CALL wrf_error_fatal ('allocating height in time_for_move2')
353       ALLOCATE ( psfc     ( ids:ide , jds:jde ), STAT=ierr )
354       IF ( ierr .ne. 0 ) CALL wrf_error_fatal ('allocating psfc in time_for_move2')
355       ALLOCATE ( xlat     ( ids:ide , jds:jde ), STAT=ierr )
356       IF ( ierr .ne. 0 ) CALL wrf_error_fatal ('allocating xlat in time_for_move2')
357       ALLOCATE ( xlong    ( ids:ide , jds:jde ), STAT=ierr )
358       IF ( ierr .ne. 0 ) CALL wrf_error_fatal ('allocating xlong in time_for_move2')
359       ALLOCATE ( terrain  ( ids:ide , jds:jde ), STAT=ierr )
360       IF ( ierr .ne. 0 ) CALL wrf_error_fatal ('allocating terrain in time_for_move2')
361     ELSE
362       ALLOCATE ( height   ( 1:1 , 1:1 ), STAT=ierr )
363       IF ( ierr .ne. 0 ) CALL wrf_error_fatal ('allocating height in time_for_move2')
364       ALLOCATE ( psfc     ( 1:1 , 1:1 ), STAT=ierr )
365       IF ( ierr .ne. 0 ) CALL wrf_error_fatal ('allocating psfc in time_for_move2')
366       ALLOCATE ( xlat     ( 1:1 , 1:1 ), STAT=ierr )
367       IF ( ierr .ne. 0 ) CALL wrf_error_fatal ('allocating xlat in time_for_move2')
368       ALLOCATE ( xlong    ( 1:1 , 1:1 ), STAT=ierr )
369       IF ( ierr .ne. 0 ) CALL wrf_error_fatal ('allocating xlong in time_for_move2')
370       ALLOCATE ( terrain  ( 1:1 , 1:1 ), STAT=ierr )
371       IF ( ierr .ne. 0 ) CALL wrf_error_fatal ('allocating terrain in time_for_move2')
372     ENDIF
373
374#  if (EM_CORE == 1)
375     IF ( dyn_opt .EQ. DYN_EM ) THEN
376       CALL compute_500mb_height ( grid%em_ph_2 , grid%em_phb, grid%em_p, grid%em_pb, height_l , &
377                                   ids, ide, jds, jde, kds, kde, &
378                                   ims, ime, jms, jme, kms, kme, &
379                                   ips, ipe, jps, jpe, kps, kpe  )
380     ENDIF
381#  endif
382#  if (WRF_NMM_CORE == 1)
383     IF ( dyn_opt .EQ. DYN_NMM ) THEN
384     ENDIF
385#  endif
386
387     CALL wrf_patch_to_global_real ( height_l , height , grid%domdesc, "z", "xy", &
388                                     ids, ide-1 , jds , jde-1 , 1 , 1 , &
389                                     ims, ime   , jms , jme   , 1 , 1 , &
390                                     ips, ipe   , jps , jpe   , 1 , 1   )
391     CALL wrf_patch_to_global_real ( grid%psfc , psfc , grid%domdesc, "z", "xy", &
392                                     ids, ide-1 , jds , jde-1 , 1 , 1 , &
393                                     ims, ime   , jms , jme   , 1 , 1 , &
394                                     ips, ipe   , jps , jpe   , 1 , 1   )
395     CALL wrf_patch_to_global_real ( grid%xlat , xlat , grid%domdesc, "z", "xy", &
396                                     ids, ide-1 , jds , jde-1 , 1 , 1 , &
397                                     ims, ime   , jms , jme   , 1 , 1 , &
398                                     ips, ipe   , jps , jpe   , 1 , 1   )
399     CALL wrf_patch_to_global_real ( grid%xlong , xlong , grid%domdesc, "z", "xy", &
400                                     ids, ide-1 , jds , jde-1 , 1 , 1 , &
401                                     ims, ime   , jms , jme   , 1 , 1 , &
402                                     ips, ipe   , jps , jpe   , 1 , 1   )
403     CALL wrf_patch_to_global_real ( grid%ht , terrain , grid%domdesc, "z", "xy", &
404                                     ids, ide-1 , jds , jde-1 , 1 , 1 , &
405                                     ims, ime   , jms , jme   , 1 , 1 , &
406                                     ips, ipe   , jps , jpe   , 1 , 1   )
407
408! calculate max wind speed
409     maxws = 0.
410     do j = jps, jpe
411       do i = ips, ipe
412         ws = grid%u10(i,j) * grid%u10(i,j) + grid%v10(i,j) * grid%v10(i,j)
413         if ( ws > maxws ) maxws = ws
414       enddo
415     enddo
416     maxws = sqrt ( maxws )
417     maxws = wrf_dm_max_real ( maxws )
418
419     monitor_only : IF ( wrf_dm_on_monitor() ) THEN
420
421!
422! This vortex center finding code adapted from the Hurricane version of MM5,
423! Courtesy:
424!
425!   Shuyi Chen et al., Rosenstiel School of Marine and Atmos. Sci., U. Miami.
426!   Spring, 2005
427!
428! Get the first guess vortex center about which we do our search
429! as mini and minh; minimum value is minh
430!
431
432       CALL nl_get_vortex_interval( grid%id , vortex_interval )
433       CALL nl_get_max_vortex_speed( grid%id , max_vortex_speed )
434
435       IF ( grid%vc_i < 0. .AND. grid%vc_j < 0. ) THEN
436          ! first time through
437          is = ids
438          ie = ide-1
439          js = jds
440          je = jde-1
441       ELSE
442          ! limit the search to an area around the vortex
443          ! that is limited by max_vortex_speed (default 40) meters per second from
444          ! previous location over vortex_interval (default 15 mins)
445
446          is = max( grid%vc_i - 60 * vortex_interval * max_vortex_speed / dx , 1.0 * ids )
447          js = max( grid%vc_j - 60 * vortex_interval * max_vortex_speed / dx , 1.0 * jds )
448          ie = min( grid%vc_i + 60 * vortex_interval * max_vortex_speed / dx , 1.0 * (ide-1) )
449          je = min( grid%vc_j + 60 * vortex_interval * max_vortex_speed / dx , 1.0 * (jde-1) )
450
451       ENDIF
452
453#  ifdef MOVING_DIAGS
454write(0,*)'search set around last position '
455write(0,*)'   is, ids-1,  ie,  ide-1 ', is, ids-1, ie, ide-1
456write(0,*)'   js, jds-1,  je,  jde-1 ', js, jds-1, je, jde-1
457#  endif
458
459       imploc = -1
460       jmploc = -1
461
462       ! find minimum psfc
463       pmin = 99999999.0     ! make this very large to be sure we find a minumum
464       DO j = js, je
465       DO i = is, ie
466         IF ( psfc(i,j) .LT. pmin ) THEN
467           pmin = psfc(i,j)
468           imploc = i
469           jmploc = j
470         ENDIF
471       ENDDO
472       ENDDO
473
474       IF ( imploc .EQ. -1 .OR. jmploc .EQ. -1 ) THEN  ! if we fail to find a min there is something seriously wrong
475         WRITE(mess,*)'i,j,is,ie,js,je,imploc,jmploc ',i,j,is,ie,js,je,imploc,jmploc
476         CALL wrf_message(mess)
477         CALL wrf_error_fatal('time_for_move2: Method failure searching for minimum psfc.')
478       ENDIF
479
480       imloc = -1
481       jmloc = -1
482       maxi = -1
483       maxj = -1
484
485       ! find local min, max
486       vmin =  99999999.0
487       vmax = -99999999.0
488       DO j = js, je
489       DO i = is, ie
490         IF ( height(i,j) .LT. vmin ) THEN
491           vmin = height(i,j)
492           imloc = i
493           jmloc = j
494         ENDIF
495         IF ( height(i,j) .GT. vmax ) THEN
496           vmax = height(i,j)
497           maxi = i
498           maxj = j
499         ENDIF
500       ENDDO
501       ENDDO
502
503       IF ( imloc .EQ. -1 .OR. jmloc .EQ. -1 .OR. maxi .EQ. -1 .OR. maxj .EQ. -1 ) THEN
504         WRITE(mess,*)'i,j,is,ie,js,je,imloc,jmloc,maxi,maxj ',i,j,is,ie,js,je,imloc,jmloc,maxi,maxj
505         CALL wrf_message(mess)
506         CALL wrf_error_fatal('time_for_move2: Method failure searching max/min of height.')
507       ENDIF
508
509       fimloc = imloc
510       fjmloc = jmloc
511
512       if ( grid%xi .EQ. -1. ) grid%xi = fimloc
513       if ( grid%xj .EQ. -1. ) grid%xj = fjmloc
514
515       dijsmooth = rsmooth / dx
516
517       fjs = max(fjmloc-dijsmooth,1.0)
518       fje = min(fjmloc+dijsmooth,jde-2.0)
519       fis = max(fimloc-dijsmooth,1.0)
520       fie = min(fimloc+dijsmooth,ide-2.0)
521       js = fjs
522       je = fje
523       is = fis
524       ie = fie
525
526       vmin =  1000000.0
527       vmax = -1000000.0
528       DO j = js, je
529       DO i = is, ie
530         IF ( height(i,j) .GT. vmax ) THEN
531           vmax = height(i,j)
532         ENDIF
533       ENDDO
534       ENDDO
535
536       pbar  = 0.0
537       ipbar = 0.0
538       jpbar = 0.0
539
540       do j=js,je
541          do i=is,ie
542             fact = vmax - height(i,j)
543             pbar  = pbar + fact
544             ipbar = ipbar + fact*(i-is)
545             jpbar = jpbar + fact*(j-js)
546          enddo
547       enddo
548
549      IF ( pbar .NE. 0. ) THEN
550
551!     Compute an adjusted, smoothed, vortex center location in cross
552!     point index space.
553!     Time average. A is coef for old information; B is new
554!     If pbar is zero then just skip this, leave xi and xj alone,
555!     result will be no movement.
556         a = 0.0
557         b = 1.0
558         grid%xi =  (a * grid%xi + b * (is + ipbar / pbar))  / ( a + b )
559         grid%xj =  (a * grid%xj + b * (js + jpbar / pbar))  / ( a + b )
560
561         grid%vc_i = grid%xi + .5
562         grid%vc_j = grid%xj + .5
563
564
565      ENDIF
566
567#  ifdef MOVING_DIAGS
568write(0,*)'computed grid%vc_i, grid%vc_j ',grid%vc_i, grid%vc_j
569i = grid%vc_i ; j = grid%vc_j ; height( i,j ) = height(i,j) * 1.2   !mark the center
570CALL domain_clock_get( grid, current_timestr=message2 )
571WRITE ( message , FMT = '(A," on domain ",I3)' ) TRIM(message2), grid%id
572#  endif
573
574!
575        i = INT(grid%xi+.5)
576        j = INT(grid%xj+.5)
577        write(mess,'("ATCF"," ",A19," ",f8.2," ",f8.2," ",f6.1," ",f6.1)')                &
578                                       timestr(1:19),                               &
579                                       xlat(i,j),                                   &
580                                       xlong(i,j),                                  &
581                                       0.01*pmin+0.1138*terrain(imploc,jmploc),     &
582                                       maxws*1.94
583        CALL wrf_message(TRIM(mess))
584                           
585
586
587     ENDIF monitor_only
588
589     DEALLOCATE ( psfc )
590     DEALLOCATE ( xlat )
591     DEALLOCATE ( xlong )
592     DEALLOCATE ( terrain )
593     DEALLOCATE ( height )
594     DEALLOCATE ( height_l )
595
596     CALL wrf_dm_bcast_real( grid%vc_i , 1 )
597     CALL wrf_dm_bcast_real( grid%vc_j , 1 )
598
599     CALL wrf_dm_bcast_real( pmin , 1 )
600     CALL wrf_dm_bcast_integer( imploc , 1 )
601     CALL wrf_dm_bcast_integer( jmploc , 1 )
602
603#  ifdef MOVING_DIAGS
604write(0,*)'after bcast : grid%vc_i, grid%vc_j ',grid%vc_i, grid%vc_j
605#  endif
606
607
608   ENDIF   ! COMPUTE_VORTEX_CENTER_ALARM ringing
609
610#  ifdef MOVING_DIAGS
611write(0,*)'After ENDIF : grid%vc_i, grid%vc_j ',grid%vc_i, grid%vc_j
612#  endif
613
614   dc_i = (ide-ids+1)/2.    ! domain center
615   dc_j = (jde-jds+1)/2.
616
617   disp_x = grid%vc_i - dc_i * 1.0
618   disp_y = grid%vc_j - dc_j * 1.0
619
620#if 0
621! This appears to be an old, redundant, and perhaps even misnamed parameter.
622! Remove it from the namelist and Registry and just hard code it to
623! the default of 6. JM 20050721
624   CALL nl_get_vortex_search_radius( 1, irad )
625#else
626   irad = 6
627#endif
628
629   radius = irad
630
631   if ( disp_x .GT. 0 ) disp_x = min( disp_x , radius )
632   if ( disp_y .GT. 0 ) disp_y = min( disp_y , radius )
633
634   if ( disp_x .LT. 0 ) disp_x = max( disp_x , -radius )
635   if ( disp_y .LT. 0 ) disp_y = max( disp_y , -radius )
636
637   move_cd_x = int ( disp_x  / pgr )
638   move_cd_y = int ( disp_y  / pgr )
639
640   IF ( move_cd_x .GT. 0 ) move_cd_x = min ( move_cd_x , 1 )
641   IF ( move_cd_y .GT. 0 ) move_cd_y = min ( move_cd_y , 1 )
642   IF ( move_cd_x .LT. 0 ) move_cd_x = max ( move_cd_x , -1 )
643   IF ( move_cd_y .LT. 0 ) move_cd_y = max ( move_cd_y , -1 )
644
645   CALL domain_clock_get( grid, current_timestr=timestr )
646
647   WRITE(mess,*)timestr(1:19),' vortex center (in nest x and y): ',grid%vc_i, grid%vc_j
648   CALL wrf_message(TRIM(mess))
649   WRITE(mess,*)timestr(1:19),' grid   center (in nest x and y): ',     dc_i,      dc_j
650   CALL wrf_message(TRIM(mess))
651   WRITE(mess,*)timestr(1:19),' disp          : ',   disp_x,    disp_y
652   CALL wrf_message(TRIM(mess))
653   WRITE(mess,*)timestr(1:19),' move (rel cd) : ',move_cd_x, move_cd_y
654   CALL wrf_message(TRIM(mess))
655
656   grid%vc_i = grid%vc_i - move_cd_x * pgr
657   grid%vc_j = grid%vc_j - move_cd_y * pgr
658
659#  ifdef MOVING_DIAGS
660write(0,*)' changing grid%vc_i,  move_cd_x * pgr ', grid%vc_i, move_cd_x * pgr, move_cd_x, pgr
661#  endif
662
663   IF ( ( move_cd_x .NE. 0 ) .OR. ( move_cd_y .NE. 0 ) ) THEN
664     time_for_move2 = .TRUE.
665   ELSE
666     time_for_move2 = .FALSE.
667   ENDIF
668
669# else
670! from namelist
671   move_cd_x = 0
672   move_cd_y = 0
673   time_for_move2 = .FALSE.
674   CALL domain_clock_get( grid, current_time=ct, start_time=st )
675   CALL nl_get_num_moves( 1, num_moves )
676   IF ( num_moves .GT. max_moves ) THEN
677     WRITE(mess,*)'time_for_moves2: num_moves (',num_moves,') .GT. max_moves (',max_moves,')'
678     CALL wrf_error_fatal( TRIM(mess) )
679   ENDIF
680   DO i = 1, num_moves
681     CALL nl_get_move_id( i, move_id )
682     IF ( move_id .EQ. grid%id ) THEN
683       CALL nl_get_move_interval( i, move_interval )
684       IF ( move_interval .LT. 999999999 ) THEN
685         CALL WRFU_TimeIntervalSet ( ti, M=move_interval, rc=rc )
686         IF ( ct .GE. st + ti ) THEN
687           CALL nl_get_move_cd_x ( i, move_cd_x )
688           CALL nl_get_move_cd_y ( i, move_cd_y )
689           CALL nl_set_move_interval ( i, 999999999 )
690           time_for_move2 = .TRUE.
691           EXIT
692         ENDIF
693       ENDIF
694     ENDIF
695   ENDDO
696# endif
697   RETURN
698#endif
699END FUNCTION time_for_move2
700
701LOGICAL FUNCTION time_for_move ( parent , grid , move_cd_x, move_cd_y )
702   USE module_domain
703   USE module_configure
704   USE module_dm
705USE module_timing
706   USE module_utility
707   IMPLICIT NONE
708! arguments
709   TYPE(domain) , POINTER    :: parent, grid, par, nst
710   INTEGER, INTENT(OUT)      :: move_cd_x , move_cd_y
711#ifdef MOVE_NESTS
712! local
713   INTEGER     :: corral_dist, kid
714   INTEGER     :: dw, de, ds, dn, pgr
715   INTEGER     :: would_move_x, would_move_y
716   INTEGER     :: cids, cide, cjds, cjde, ckds, ckde, &
717                  cims, cime, cjms, cjme, ckms, ckme, &
718                  cips, cipe, cjps, cjpe, ckps, ckpe, &
719                  nids, nide, njds, njde, nkds, nkde, &
720                  nims, nime, njms, njme, nkms, nkme, &
721                  nips, nipe, njps, njpe, nkps, nkpe
722! interface
723   INTERFACE
724     LOGICAL FUNCTION time_for_move2 ( parent , nest , dx , dy )
725        USE module_domain
726        USE module_utility
727        TYPE(domain) , POINTER    :: parent , nest
728        INTEGER, INTENT(OUT)      :: dx , dy
729     END FUNCTION time_for_move2
730   END INTERFACE
731! executable
732!
733! Simplifying assumption: domains in moving nest simulations have only
734! one parent and only one child.
735
736   IF   ( grid%num_nests .GT. 1 ) THEN
737     CALL wrf_error_fatal ( 'domains in moving nest simulations can have only 1 nest' )
738   ENDIF
739   kid = 1
740!
741! find out if this is the innermost nest (will not have kids)
742   IF   ( grid%num_nests .EQ. 0 ) THEN
743     ! code that executes on innermost nest
744     time_for_move = time_for_move2 ( parent , grid , move_cd_x, move_cd_y )
745
746     ! Make sure the parent can move before allowing the nest to approach
747     ! its boundary
748     par => grid%parents(1)%ptr
749     nst => grid
750
751     would_move_x = move_cd_x
752     would_move_y = move_cd_y
753
754     ! top of until loop
755100  CONTINUE
756       CALL nl_get_corral_dist ( nst%id , corral_dist )
757       CALL get_ijk_from_grid (  nst ,                               &
758                                 nids, nide, njds, njde, nkds, nkde, &
759                                 nims, nime, njms, njme, nkms, nkme, &
760                                 nips, nipe, njps, njpe, nkps, nkpe  )
761       CALL get_ijk_from_grid (  par ,                               &
762                                 cids, cide, cjds, cjde, ckds, ckde, &
763                                 cims, cime, cjms, cjme, ckms, ckme, &
764                                 cips, cipe, cjps, cjpe, ckps, ckpe  )
765       CALL nl_get_parent_grid_ratio ( nst%id , pgr )
766       ! perform measurements...
767       !  from western boundary
768       dw = nst%i_parent_start + would_move_x - cids
769       !  from southern boundary
770       ds = nst%j_parent_start + would_move_y - cjds
771       !  from eastern boundary
772       de = cide - ( nst%i_parent_start + (nide-nids+1)/pgr + would_move_x )
773       !  from northern boundary
774       dn = cjde - ( nst%j_parent_start + (njde-njds+1)/pgr + would_move_y )
775
776       ! would this generate a move on the parent?
777       would_move_x = 0
778       would_move_y = 0
779       if ( dw .LE. corral_dist ) would_move_x = would_move_x - 1
780       if ( de .LE. corral_dist ) would_move_x = would_move_x + 1
781       if ( ds .LE. corral_dist ) would_move_y = would_move_y - 1
782       if ( dn .LE. corral_dist ) would_move_y = would_move_y + 1
783
784     IF ( par%id .EQ. 1 ) THEN
785         IF ( would_move_x .NE. 0 .AND. move_cd_x .NE. 0 ) THEN
786           CALL wrf_message('MOAD can not move. Cancelling nest move in X')
787           if ( grid%num_nests .eq. 0 ) grid%vc_i = grid%vc_i + move_cd_x * pgr  ! cancel effect of move
788           move_cd_x = 0
789         ENDIF
790         IF ( would_move_y .NE. 0 .AND. move_cd_y .NE. 0 ) THEN
791           CALL wrf_message('MOAD can not move. Cancelling nest move in Y')
792           if ( grid%num_nests .eq. 0 ) grid%vc_j = grid%vc_j + move_cd_y * pgr  ! cancel effect of move
793           move_cd_y = 0
794         ENDIF
795     ELSE
796         nst => par
797         par => nst%parents(1)%ptr
798         GOTO 100
799     ENDIF
800
801! bottom of until loop
802     time_for_move = ( move_cd_x .NE. 0 ) .OR. ( move_cd_y .NE. 0 )
803
804   ELSE
805     ! code that executes on parent to see if parent needs to move
806     ! get closest number of cells we'll allow nest edge to approach parent bdy
807     CALL nl_get_corral_dist ( grid%nests(kid)%ptr%id , corral_dist )
808     ! get dims
809     CALL get_ijk_from_grid (  grid%nests(kid)%ptr ,               &
810                               nids, nide, njds, njde, nkds, nkde, &
811                               nims, nime, njms, njme, nkms, nkme, &
812                               nips, nipe, njps, njpe, nkps, nkpe  )
813     CALL get_ijk_from_grid (  grid ,                              &
814                               cids, cide, cjds, cjde, ckds, ckde, &
815                               cims, cime, cjms, cjme, ckms, ckme, &
816                               cips, cipe, cjps, cjpe, ckps, ckpe  )
817     CALL nl_get_parent_grid_ratio ( grid%nests(kid)%ptr%id , pgr )
818     ! perform measurements...
819     !  from western boundary
820     dw = grid%nests(kid)%ptr%i_parent_start - 1
821     !  from southern boundary
822     ds = grid%nests(kid)%ptr%j_parent_start - 1
823     !  from eastern boundary
824     de = cide - ( grid%nests(kid)%ptr%i_parent_start + (nide-nids+1)/pgr )
825     !  from northern boundary
826     dn = cjde - ( grid%nests(kid)%ptr%j_parent_start + (njde-njds+1)/pgr )
827
828     ! move this domain (the parent containing the moving nest)
829     ! in a direction that reestablishes the distance from
830     ! the boundary.
831     move_cd_x = 0
832     move_cd_y = 0
833     if ( dw .LE. corral_dist ) move_cd_x = move_cd_x - 1
834     if ( de .LE. corral_dist ) move_cd_x = move_cd_x + 1
835     if ( ds .LE. corral_dist ) move_cd_y = move_cd_y - 1
836     if ( dn .LE. corral_dist ) move_cd_y = move_cd_y + 1
837
838     time_for_move = ( move_cd_x .NE. 0 ) .OR. ( move_cd_y .NE. 0 )
839
840     IF ( time_for_move ) THEN
841       IF ( grid%id .EQ. 1 ) THEN
842
843         CALL wrf_message( 'DANGER: Nest has moved too close to boundary of outermost domain.' )
844         time_for_move = .FALSE.
845
846       ELSE
847         ! need to adjust the intermediate domain of the nest in relation to this
848         ! domain since we're moving
849
850         CALL wrf_dm_move_nest ( grid , grid%nests(kid)%ptr%intermediate_grid , -move_cd_x*pgr, -move_cd_y*pgr )
851         CALL adjust_domain_dims_for_move( grid%nests(kid)%ptr%intermediate_grid , -move_cd_x*pgr, -move_cd_y*pgr )
852         grid%nests(kid)%ptr%i_parent_start = grid%nests(kid)%ptr%i_parent_start - move_cd_x*pgr
853         CALL nl_set_i_parent_start( grid%nests(kid)%ptr%id, grid%nests(kid)%ptr%i_parent_start )
854         grid%nests(kid)%ptr%j_parent_start = grid%nests(kid)%ptr%j_parent_start - move_cd_y*pgr
855         CALL nl_set_j_parent_start( grid%nests(kid)%ptr%id, grid%nests(kid)%ptr%j_parent_start )
856
857       ENDIF
858     ENDIF
859
860   ENDIF
861
862   RETURN
863#endif
864END FUNCTION time_for_move
865
866! Put any tests for non-moving options or conditions in here
867LOGICAL FUNCTION should_not_move ( id )
868  USE module_state_description
869  USE module_configure
870  IMPLICIT NONE
871  INTEGER, INTENT(IN) :: id
872 ! Local
873  LOGICAL retval
874  INTEGER cu_physics, ra_sw_physics, ra_lw_physics, ucmcall, obs_nudge_opt
875
876  retval = .FALSE.
877! check for GD ensemble cumulus, which can not move
878  CALL nl_get_cu_physics( id , cu_physics )
879  IF ( cu_physics .EQ. GDSCHEME ) THEN
880    CALL wrf_message('Grell cumulus can not be specified with moving nests. Movement disabled.')
881    retval = .TRUE.
882  ENDIF
883! check for CAM radiation scheme , which can not move
884  CALL nl_get_ra_sw_physics( id , ra_sw_physics )
885  IF ( ra_sw_physics .EQ. CAMSWSCHEME ) THEN
886    CALL wrf_message('CAM SW radiation can not be specified with moving nests. Movement disabled.')
887    retval = .TRUE.
888  ENDIF
889  CALL nl_get_ra_lw_physics( id , ra_lw_physics )
890  IF ( ra_lw_physics .EQ. CAMLWSCHEME ) THEN
891    CALL wrf_message('CAM LW radiation can not be specified with moving nests. Movement disabled.')
892    retval = .TRUE.
893  ENDIF
894! check for urban canopy Noah LSM, which can not move
895  CALL nl_get_ucmcall( id , ucmcall )
896  IF ( ucmcall .EQ. 1 ) THEN
897    CALL wrf_message('UCM Noah LSM can not be specified with moving nests. Movement disabled.')
898    retval = .TRUE.
899  ENDIF
900! check for observation nudging, which can not move
901  CALL nl_get_obs_nudge_opt( id , obs_nudge_opt )
902  IF ( obs_nudge_opt .EQ. 1 ) THEN
903    CALL wrf_message('Observation nudging can not be specified with moving nests. Movement disabled.')
904    retval = .TRUE.
905  ENDIF
906  should_not_move = retval
907END FUNCTION
908
Note: See TracBrowser for help on using the repository browser.