merge standard release WRF/WPS V3.0.1.1 into wrffire
[wrffire.git] / wrfv2_fire / dyn_em / nest_init_utils.F
blobae07eb2282353a9a5ab3ef2b62ad1d95ea829fab
1 SUBROUTINE init_domain_constants_em ( parent , nest )
2    USE module_domain, ONLY : domain
3    IMPLICIT NONE
4    TYPE(domain)  :: parent , nest
6    INTEGER iswater , map_proj, julyr, julday
7    REAL    truelat1 , truelat2 , gmt , moad_cen_lat , stand_lon
8    CHARACTER (LEN=4) :: char_junk
10 ! single-value constants
12    nest%p_top   = parent%p_top
13    nest%cfn     = parent%cfn
14    nest%cfn1    = parent%cfn1
15    nest%rdx     = 1./nest%dx
16    nest%rdy     = 1./nest%dy
17 !  nest%dts     = nest%dt/float(nest%time_step_sound)
18    nest%dtseps  = parent%dtseps  ! used in height model only?
19    nest%resm    = parent%resm    ! used in height model only?
20    nest%zetatop = parent%zetatop ! used in height model only?
21    nest%cf1     = parent%cf1
22    nest%cf2     = parent%cf2
23    nest%cf3     = parent%cf3
24    nest%gmt     = parent%gmt
25    nest%julyr   = parent%julyr
26    nest%julday  = parent%julday
28    CALL nl_get_mminlu ( 1,char_junk(1:4) )
29    CALL nl_get_iswater (1, iswater )
30    CALL nl_get_truelat1 ( 1 , truelat1 )
31    CALL nl_get_truelat2 ( 1 , truelat2 )
32    CALL nl_get_moad_cen_lat ( 1 , moad_cen_lat )
33    CALL nl_get_stand_lon ( 1 , stand_lon )
34    CALL nl_get_map_proj ( 1 , map_proj )
35    CALL nl_get_gmt ( 1 , gmt)
36    CALL nl_get_julyr ( 1 , julyr)
37    CALL nl_get_julday ( 1 , julday)
38    IF ( nest%id .NE. 1 ) THEN
39      CALL nl_set_gmt (nest%id, gmt)
40      CALL nl_set_julyr (nest%id, julyr)
41      CALL nl_set_julday (nest%id, julday)
42      CALL nl_set_iswater (nest%id, iswater )
43      CALL nl_set_truelat1 ( nest%id , truelat1 )
44      CALL nl_set_truelat2 ( nest%id , truelat2 )
45      CALL nl_set_moad_cen_lat ( nest%id , moad_cen_lat )
46      CALL nl_set_stand_lon ( nest%id , stand_lon )
47      CALL nl_set_map_proj ( nest%id , map_proj )
48    END IF
49    nest%gmt     = gmt 
50    nest%julday  = julday
51    nest%julyr   = julyr
52    nest%iswater = iswater
53    nest%truelat1= truelat1
54    nest%truelat2= truelat2
55    nest%moad_cen_lat= moad_cen_lat
56    nest%stand_lon= stand_lon
57    nest%map_proj= map_proj
59    nest%step_number  = parent%step_number
61 ! 1D constants (Z)
63    nest%fnm    = parent%fnm
64    nest%fnp    = parent%fnp
65    nest%rdnw   = parent%rdnw
66    nest%rdn    = parent%rdn
67    nest%dnw    = parent%dnw
68    nest%dn     = parent%dn
69    nest%znu    = parent%znu
70    nest%znw    = parent%znw
71    nest%t_base = parent%t_base
72    nest%u_base    = parent%u_base
73    nest%v_base    = parent%v_base
74    nest%qv_base   = parent%qv_base
75    nest%z_base    = parent%z_base
76    nest%dzs       = parent%dzs
77    nest%zs        = parent%zs
79 END SUBROUTINE init_domain_constants_em
81 SUBROUTINE blend_terrain ( ter_interpolated , ter_input , &
82                            ids , ide , jds , jde , kds , kde , & 
83                            ims , ime , jms , jme , kms , kme , & 
84                            ips , ipe , jps , jpe , kps , kpe )
86    USE module_configure
87    IMPLICIT NONE
89    INTEGER , INTENT(IN)                       :: 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,kms:kme,jms:jme) , INTENT(IN)    :: ter_interpolated
93    REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(INOUT) :: ter_input
95    REAL , DIMENSION(ims:ime,kms:kme,jms:jme) :: ter_temp
96    INTEGER :: i , j , k , spec_bdy_width
97    REAL    :: r_blend_zones
98    INTEGER blend_cell, blend_width
100    !  The fine grid elevation comes from the horizontally interpolated
101    !  parent elevation for the first spec_bdy_width row/columns, so we need
102    !  to get that value.  We blend the coarse and fine in the next blend_width
103    !  rows and columns.  After that, in the interior, it is 100% fine grid.
105    CALL nl_get_spec_bdy_width ( 1, spec_bdy_width) 
106    CALL nl_get_blend_width ( 1, blend_width)
108    !  Initialize temp values to the nest ter elevation.  This fills in the values
109    !  that will not be modified below.  
111    DO j = jps , MIN(jpe, jde-1)
112       DO k = kps , kpe
113          DO i = ips , MIN(ipe, ide-1)
114             ter_temp(i,k,j) = ter_input(i,k,j)
115          END DO 
116       END DO 
117    END DO 
119    !  To avoid some tricky indexing, we fill in the values inside out.  This allows
120    !  us to overwrite incorrect assignments.  There are replicated assignments, and
121    !  there is much unnecessary "IF test inside of a loop" stuff.  For a large
122    !  domain, this is only a patch; for a small domain, this is not a biggy.
124    r_blend_zones = 1./(blend_width+1)
125    DO j = jps , MIN(jpe, jde-1)
126       DO k = kps , kpe
127          DO i = ips , MIN(ipe, ide-1)
128             DO blend_cell = blend_width,1,-1
129                IF   ( ( i .EQ.       spec_bdy_width + blend_cell ) .OR.  ( j .EQ.       spec_bdy_width + blend_cell ) .OR. &
130                       ( i .EQ. ide - spec_bdy_width - blend_cell ) .OR.  ( j .EQ. jde - spec_bdy_width - blend_cell ) ) THEN
131                   ter_temp(i,k,j) = ( (blend_cell)*ter_input(i,k,j) + (blend_width+1-blend_cell)*ter_interpolated(i,k,j) ) &
132                                     * r_blend_zones
133                END IF
134             ENDDO
135             IF      ( ( i .LE.       spec_bdy_width     ) .OR.  ( j .LE.       spec_bdy_width     ) .OR. &
136                       ( i .GE. ide - spec_bdy_width     ) .OR.  ( j .GE. jde - spec_bdy_width     ) ) THEN
137                ter_temp(i,k,j) =      ter_interpolated(i,k,j)
138             END IF
139          END DO 
140       END DO 
141    END DO 
143    !  Set nest elevation with temp values.  All values not overwritten in the above
144    !  loops have been previously set in the initial assignment.
146    DO j = jps , MIN(jpe, jde-1)
147       DO k = kps , kpe
148          DO i = ips , MIN(ipe, ide-1)
149             ter_input(i,k,j) = ter_temp(i,k,j)
150          END DO 
151       END DO 
152    END DO 
154 END SUBROUTINE blend_terrain
156 SUBROUTINE store_terrain ( ter_interpolated , ter_input , &
157                            ids , ide , jds , jde , kds , kde , & 
158                            ims , ime , jms , jme , kms , kme , & 
159                            ips , ipe , jps , jpe , kps , kpe )
161    IMPLICIT NONE
163    INTEGER , INTENT(IN)                       :: ids , ide , jds , jde , kds , kde , & 
164                                                  ims , ime , jms , jme , kms , kme , & 
165                                                  ips , ipe , jps , jpe , kps , kpe
166    REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(OUT) :: ter_interpolated
167    REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(IN)  :: ter_input
169    INTEGER :: i , j , k
171    DO j = jps , MIN(jpe, jde-1)
172       DO k = kps , kpe
173          DO i = ips , MIN(ipe, ide-1)
174             ter_interpolated(i,k,j) = ter_input(i,k,j)
175          END DO 
176       END DO 
177    END DO 
179 END SUBROUTINE store_terrain
181 SUBROUTINE adjust_tempqv ( mub, save_mub, znw, p_top, &
182                            th, pp, qv,  &
183                            ids , ide , jds , jde , kds , kde , &
184                            ims , ime , jms , jme , kms , kme , &
185                            ips , ipe , jps , jpe , kps , kpe )
187    !USE module_configure
188    !USE module_domain
189    USE module_model_constants
191    !USE module_bc
192    !USE module_io_domain
193    !USE module_state_description
194    !USE module_timing
195    !USE module_soil_pre
196    IMPLICIT NONE
198    INTEGER , INTENT(IN)                       :: ids , ide , jds , jde , kds , kde , &
199                                                  ims , ime , jms , jme , kms , kme , &
200                                                  ips , ipe , jps , jpe , kps , kpe
201    REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN)    :: mub, save_mub
202    REAL , DIMENSION(kms:kme) , INTENT(IN)    :: znw
203    REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(INOUT) :: th, pp, qv
205    REAL , DIMENSION(ims:ime,kms:kme,jms:jme) :: p_old, p_new, rh
206    REAL :: es,dth,tc,e,dth1
207    INTEGER :: i , j , k
209    real p_top
212 ! p_old = full pressure before terrain blending; also compute initial RH
213 ! which is going to be conserved during terrain blending
214    DO j = jps , MIN(jpe, jde-1)
215       DO k = kps , kpe-1
216          DO i = ips , MIN(ipe, ide-1)
217             p_old(i,k,j) = 0.5*(znw(k+1)+znw(k))*save_mub(i,j) + p_top + pp(i,k,j)
218             tc = (th(i,k,j)+300.)*(p_old(i,k,j)/1.e5)**(2./7.) - 273.15
219             es = 610.78*exp(17.0809*tc/(234.175+tc))
220             e = qv(i,k,j)*p_old(i,k,j)/(0.622+qv(i,k,j))
221             rh(i,k,j) = e/es
222          END DO
223       END DO
224    END DO
226 ! p_new = full pressure after terrain blending; also compute temperature correction and convert RH back to QV
227    DO j = jps , MIN(jpe, jde-1)
228       DO k = kps , kpe-1
229          DO i = ips , MIN(ipe, ide-1)
230             p_new(i,k,j) = 0.5*(znw(k+1)+znw(k))*mub(i,j) + p_top + pp(i,k,j)
231 ! 2*(g/cp-6.5e-3)*R_dry/g = -191.86e-3
232             dth1 = -191.86e-3*(th(i,k,j)+300.)/(p_new(i,k,j)+p_old(i,k,j))*(p_new(i,k,j)-p_old(i,k,j))
233             dth = -191.86e-3*(th(i,k,j)+0.5*dth1+300.)/(p_new(i,k,j)+p_old(i,k,j))*(p_new(i,k,j)-p_old(i,k,j))
234             th(i,k,j) = th(i,k,j)+dth
235             tc = (th(i,k,j)+300.)*(p_new(i,k,j)/1.e5)**(2./7.) - 273.15
236             es = 610.78*exp(17.0809*tc/(234.175+tc))
237             e = rh(i,k,j)*es
238             qv(i,k,j) = 0.622*e/(p_new(i,k,j)-e)
239          END DO
240       END DO
241    END DO
244 END SUBROUTINE adjust_tempqv
246 SUBROUTINE input_terrain_rsmas ( grid ,                        &
247                            ids , ide , jds , jde , kds , kde , &
248                            ims , ime , jms , jme , kms , kme , &
249                            ips , ipe , jps , jpe , kps , kpe )
251    USE module_domain, ONLY : domain
252    IMPLICIT NONE
253    TYPE ( domain ) :: grid
255    INTEGER , INTENT(IN)                       :: ids , ide , jds , jde , kds , kde , &
256                                                  ims , ime , jms , jme , kms , kme , &
257                                                  ips , ipe , jps , jpe , kps , kpe
259    LOGICAL, EXTERNAL ::  wrf_dm_on_monitor
261    INTEGER :: i , j , k , myproc
262    INTEGER, DIMENSION(256) :: ipath  ! array for integer coded ascii for passing path down to get_terrain
263    CHARACTER*256 :: message, message2
264    CHARACTER*256 :: rsmas_data_path
266 #if DM_PARALLEL
267 ! Local globally sized arrays
268    REAL , DIMENSION(ids:ide,jds:jde) :: ht_g, xlat_g, xlon_g
269 #endif
271    CALL wrf_get_myproc ( myproc ) 
273 #if 0
274 CALL domain_clock_get ( grid, current_timestr=message2 )
275 WRITE ( message , FMT = '(A," HT before ",I3)' ) TRIM(message2), grid%id
276 write(30+myproc,*)ipe-ips+1,jpe-jps+1,trim(message)
277 do j = jps,jpe
278 do i = ips,ipe
279 write(30+myproc,*)grid%ht(i,j)
280 enddo
281 enddo
282 #endif
284    CALL nl_get_rsmas_data_path(1,rsmas_data_path)
285    do i = 1, LEN(TRIM(rsmas_data_path))
286       ipath(i) = ICHAR(rsmas_data_path(i:i))
287    enddo
289 #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
291    CALL wrf_patch_to_global_real ( grid%xlat , xlat_g , grid%domdesc, ' ' , 'xy' ,       &
292                                    ids, ide-1 , jds , jde-1 , 1 , 1 , &
293                                    ims, ime   , jms , jme   , 1 , 1 , &
294                                    ips, ipe   , jps , jpe   , 1 , 1   ) 
295    CALL wrf_patch_to_global_real ( grid%xlong , xlon_g , grid%domdesc, ' ' , 'xy' ,       &
296                                    ids, ide-1 , jds , jde-1 , 1 , 1 , &
297                                    ims, ime   , jms , jme   , 1 , 1 , &
298                                    ips, ipe   , jps , jpe   , 1 , 1   ) 
300    IF ( wrf_dm_on_monitor() ) THEN
301      CALL get_terrain ( grid%dx/1000., xlat_g(ids:ide,jds:jde), xlon_g(ids:ide,jds:jde), ht_g(ids:ide,jds:jde), &
302                         ide-ids+1,jde-jds+1,ide-ids+1,jde-jds+1, ipath, LEN(TRIM(rsmas_data_path)) )
303      WHERE ( ht_g(ids:ide,jds:jde) < -1000. ) ht_g(ids:ide,jds:jde) = 0.
304    ENDIF
306    CALL wrf_global_to_patch_real ( ht_g , grid%ht , grid%domdesc, ' ' , 'xy' ,         &
307                                    ids, ide-1 , jds , jde-1 , 1 , 1 , &
308                                    ims, ime   , jms , jme   , 1 , 1 , &
309                                    ips, ipe   , jps , jpe   , 1 , 1   ) 
310 #else
312    CALL get_terrain ( grid%dx/1000., grid%xlat(ids:ide,jds:jde), grid%xlong(ids:ide,jds:jde), grid%ht(ids:ide,jds:jde), &
313                        ide-ids+1,jde-jds+1,ide-ids+1,jde-jds+1, ipath, LEN(TRIM(rsmas_data_path)) )
314    WHERE ( grid%ht(ids:ide,jds:jde) < -1000. ) grid%ht(ids:ide,jds:jde) = 0.
316 #endif
318 #if 0
319 CALL domain_clock_get ( grid, current_timestr=message2 )
320 WRITE ( message , FMT = '(A," HT after ",I3)' ) TRIM(message2), grid%id
321 write(30+myproc,*)ipe-ips+1,jpe-jps+1,trim(message)
322 do j = jps,jpe
323 do i = ips,ipe
324 write(30+myproc,*)grid%ht(i,j)
325 enddo
326 enddo
327 #endif
328                        
329 END SUBROUTINE input_terrain_rsmas
331 SUBROUTINE update_after_feedback_em ( grid  &
333 #include "dummy_new_args.inc"
335                  )
337 ! perform core specific updates, exchanges after
338 ! model feedback  (called from med_feedback_domain) -John
341 ! Driver layer modules
342    USE module_domain, ONLY : domain, get_ijk_from_grid
343    USE module_configure
344    USE module_driver_constants
345    USE module_machine
346    USE module_tiles
347    USE module_dm
348    USE module_bc
349 ! Mediation layer modules
350 ! Registry generated module
351    USE module_state_description
353    IMPLICIT NONE
355    !  Subroutine interface block.
357    TYPE(domain) , TARGET         :: grid
359    !  Definitions of dummy arguments
360 #include <dummy_new_decl.inc>
362    INTEGER                         :: ids , ide , jds , jde , kds , kde , &
363                                       ims , ime , jms , jme , kms , kme , &
364                                       ips , ipe , jps , jpe , kps , kpe
366   CALL wrf_debug( 500, "entering update_after_feedback_em" )
368 !  Obtain dimension information stored in the grid data structure.
369   CALL get_ijk_from_grid (  grid ,                   &
370                             ids, ide, jds, jde, kds, kde,    &
371                             ims, ime, jms, jme, kms, kme,    &
372                             ips, ipe, jps, jpe, kps, kpe    )
374   CALL wrf_debug( 500, "before HALO_EM_FEEDBACK.inc in update_after_feedback_em" )
375 #ifdef DM_PARALLEL
376 #include "HALO_EM_FEEDBACK.inc"
377 #endif
378   CALL wrf_debug( 500, "leaving update_after_feedback_em" )
380 END SUBROUTINE update_after_feedback_em