source: lmdz_wrf/WRFV3/dyn_em/nest_init_utils.F @ 1

Last change on this file since 1 was 1, checked in by lfita, 10 years ago
  • -- --- Opening of the WRF+LMDZ coupling repository --- -- -

WRF: version v3.3
LMDZ: version v1818

More details in:

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