source: lmdz_wrf/trunk/WRFV3/dyn_em/module_initialize_les.F @ 1419

Last change on this file since 1419 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: 24.1 KB
Line 
1!IDEAL:MODEL_LAYER:INITIALIZATION
2!
3
4!  This MODULE holds the routines which are used to perform various initializations
5!  for the individual domains. 
6
7!  This MODULE CONTAINS the following routines:
8
9!  initialize_field_test - 1. Set different fields to different constant
10!                             values.  This is only a test.  If the correct
11!                             domain is not found (based upon the "id")
12!                             then a fatal error is issued.               
13
14!-----------------------------------------------------------------------
15
16MODULE module_initialize_ideal
17
18   USE module_domain
19   USE module_io_domain
20   USE module_state_description
21   USE module_model_constants
22   USE module_bc
23   USE module_timing
24   USE module_configure
25   USE module_init_utilities
26#ifdef DM_PARALLEL
27   USE module_dm
28#endif
29
30
31CONTAINS
32
33
34!-------------------------------------------------------------------
35! this is a wrapper for the solver-specific init_domain routines.
36! Also dereferences the grid variables and passes them down as arguments.
37! This is crucial, since the lower level routines may do message passing
38! and this will get fouled up on machines that insist on passing down
39! copies of assumed-shape arrays (by passing down as arguments, the
40! data are treated as assumed-size -- ie. f77 -- arrays and the copying
41! business is avoided).  Fie on the F90 designers.  Fie and a pox.
42
43   SUBROUTINE init_domain ( grid )
44
45   IMPLICIT NONE
46
47   !  Input data.
48   TYPE (domain), POINTER :: grid
49   !  Local data.
50   INTEGER :: idum1, idum2
51
52   CALL set_scalar_indices_from_config ( head_grid%id , idum1, idum2 )
53
54     CALL init_domain_rk( grid &
55!
56#include <actual_new_args.inc>
57!
58                        )
59
60   END SUBROUTINE init_domain
61
62!-------------------------------------------------------------------
63
64   SUBROUTINE init_domain_rk ( grid &
65!
66# include <dummy_new_args.inc>
67!
68)
69   IMPLICIT NONE
70
71   !  Input data.
72   TYPE (domain), POINTER :: grid
73
74# include <dummy_new_decl.inc>
75
76   TYPE (grid_config_rec_type)              :: config_flags
77
78   !  Local data
79   INTEGER                             ::                       &
80                                  ids, ide, jds, jde, kds, kde, &
81                                  ims, ime, jms, jme, kms, kme, &
82                                  its, ite, jts, jte, kts, kte, &
83                                  i, j, k
84
85   ! Local data
86
87   INTEGER, PARAMETER :: nl_max = 1000
88   REAL, DIMENSION(nl_max) :: zk, p_in, theta, rho, u, v, qv, pd_in
89   INTEGER :: nl_in
90
91
92   INTEGER :: icm,jcm, ii, im1, jj, jm1, loop, error, fid, nxc, nyc
93   REAL    :: u_mean,v_mean, f0, p_surf, p_level, qvf, z_at_v, z_at_u
94   REAL    :: z_scale, xrad, yrad, zrad, rad, delt, cof1, cof2
95!   REAL, EXTERNAL :: interp_0
96   REAL    :: hm
97   REAL    :: pi
98
99!  stuff from original initialization that has been dropped from the Registry
100   REAL    :: vnu, xnu, xnus, dinit0, cbh, p0_temp, t0_temp, zd, zt
101   REAL    :: qvf1, qvf2, pd_surf
102   INTEGER :: it
103   real :: thtmp, ptmp, temp(3)
104
105   LOGICAL :: moisture_init
106   LOGICAL :: stretch_grid, dry_sounding
107
108  INTEGER :: xs , xe , ys , ye
109  REAL :: mtn_ht
110   LOGICAL, EXTERNAL :: wrf_dm_on_monitor
111!  For LES, add randx
112   real :: randx
113
114#ifdef DM_PARALLEL
115#    include <data_calls.inc>
116#endif
117
118
119   SELECT CASE ( model_data_order )
120         CASE ( DATA_ORDER_ZXY )
121   kds = grid%sd31 ; kde = grid%ed31 ;
122   ids = grid%sd32 ; ide = grid%ed32 ;
123   jds = grid%sd33 ; jde = grid%ed33 ;
124
125   kms = grid%sm31 ; kme = grid%em31 ;
126   ims = grid%sm32 ; ime = grid%em32 ;
127   jms = grid%sm33 ; jme = grid%em33 ;
128
129   kts = grid%sp31 ; kte = grid%ep31 ;   ! note that tile is entire patch
130   its = grid%sp32 ; ite = grid%ep32 ;   ! note that tile is entire patch
131   jts = grid%sp33 ; jte = grid%ep33 ;   ! note that tile is entire patch
132         CASE ( DATA_ORDER_XYZ )
133   ids = grid%sd31 ; ide = grid%ed31 ;
134   jds = grid%sd32 ; jde = grid%ed32 ;
135   kds = grid%sd33 ; kde = grid%ed33 ;
136
137   ims = grid%sm31 ; ime = grid%em31 ;
138   jms = grid%sm32 ; jme = grid%em32 ;
139   kms = grid%sm33 ; kme = grid%em33 ;
140
141   its = grid%sp31 ; ite = grid%ep31 ;   ! note that tile is entire patch
142   jts = grid%sp32 ; jte = grid%ep32 ;   ! note that tile is entire patch
143   kts = grid%sp33 ; kte = grid%ep33 ;   ! note that tile is entire patch
144         CASE ( DATA_ORDER_XZY )
145   ids = grid%sd31 ; ide = grid%ed31 ;
146   kds = grid%sd32 ; kde = grid%ed32 ;
147   jds = grid%sd33 ; jde = grid%ed33 ;
148
149   ims = grid%sm31 ; ime = grid%em31 ;
150   kms = grid%sm32 ; kme = grid%em32 ;
151   jms = grid%sm33 ; jme = grid%em33 ;
152
153   its = grid%sp31 ; ite = grid%ep31 ;   ! note that tile is entire patch
154   kts = grid%sp32 ; kte = grid%ep32 ;   ! note that tile is entire patch
155   jts = grid%sp33 ; jte = grid%ep33 ;   ! note that tile is entire patch
156
157   END SELECT
158
159
160!  stretch_grid = .true.
161!  FOR LES, set stretch to false
162   stretch_grid = .false.
163   delt = 3.
164!   z_scale = .50
165   z_scale = .40
166   pi = 2.*asin(1.0)
167   write(6,*) ' pi is ',pi
168   nxc = (ide-ids)/2
169   nyc = (jde-jds)/2
170
171   CALL model_to_grid_config_rec ( grid%id , model_config_rec , config_flags )
172
173! here we check to see if the boundary conditions are set properly
174
175   CALL boundary_condition_check( config_flags, bdyzone, error, grid%id )
176
177   moisture_init = .true.
178
179    grid%itimestep=0
180
181#ifdef DM_PARALLEL
182   CALL wrf_dm_bcast_bytes( icm , IWORDSIZE )
183   CALL wrf_dm_bcast_bytes( jcm , IWORDSIZE )
184#endif
185
186    CALL nl_set_mminlu(1, '    ')
187    CALL nl_set_iswater(1,0)
188    CALL nl_set_cen_lat(1,40.)
189    CALL nl_set_cen_lon(1,-105.)
190    CALL nl_set_truelat1(1,0.)
191    CALL nl_set_truelat2(1,0.)
192    CALL nl_set_moad_cen_lat (1,0.)
193    CALL nl_set_stand_lon (1,0.)
194    CALL nl_set_pole_lon (1,0.)
195    CALL nl_set_pole_lat (1,90.)
196    CALL nl_set_map_proj(1,0)
197
198
199!  here we initialize data we currently is not initialized
200!  in the input data
201
202    DO j = jts, jte
203      DO i = its, ite
204         grid%msftx(i,j)    = 1.
205         grid%msfty(i,j)    = 1.
206         grid%msfux(i,j)    = 1.
207         grid%msfuy(i,j)    = 1.
208         grid%msfvx(i,j)    = 1.
209         grid%msfvx_inv(i,j)= 1.
210         grid%msfvy(i,j)    = 1.
211         grid%sina(i,j)     = 0.
212         grid%cosa(i,j)     = 1.
213         grid%e(i,j)        = 0.
214!  for LES, include Coriolis force
215         grid%f(i,j)        = 1.e-4
216
217      END DO
218   END DO
219
220    DO j = jts, jte
221    DO k = kts, kte
222      DO i = its, ite
223         grid%ww(i,k,j)     = 0.
224      END DO
225   END DO
226   END DO
227
228   grid%step_number = 0
229
230! set up the grid
231
232   IF (stretch_grid) THEN ! exponential stretch for eta (nearly constant dz)
233     DO k=1, kde
234      grid%znw(k) = (exp(-(k-1)/float(kde-1)/z_scale) - exp(-1./z_scale))/ &
235                                (1.-exp(-1./z_scale))
236     ENDDO
237   ELSE
238     DO k=1, kde
239      grid%znw(k) = 1. - float(k-1)/float(kde-1)
240     ENDDO
241   ENDIF
242
243   DO k=1, kde-1
244    grid%dnw(k) = grid%znw(k+1) - grid%znw(k)
245    grid%rdnw(k) = 1./grid%dnw(k)
246    grid%znu(k) = 0.5*(grid%znw(k+1)+grid%znw(k))
247   ENDDO
248   DO k=2, kde-1
249    grid%dn(k) = 0.5*(grid%dnw(k)+grid%dnw(k-1))
250    grid%rdn(k) = 1./grid%dn(k)
251    grid%fnp(k) = .5* grid%dnw(k  )/grid%dn(k)
252    grid%fnm(k) = .5* grid%dnw(k-1)/grid%dn(k)
253   ENDDO
254
255   cof1 = (2.*grid%dn(2)+grid%dn(3))/(grid%dn(2)+grid%dn(3))*grid%dnw(1)/grid%dn(2)
256   cof2 =     grid%dn(2)        /(grid%dn(2)+grid%dn(3))*grid%dnw(1)/grid%dn(3)
257   grid%cf1  = grid%fnp(2) + cof1
258   grid%cf2  = grid%fnm(2) - cof1 - cof2
259   grid%cf3  = cof2       
260
261   grid%cfn  = (.5*grid%dnw(kde-1)+grid%dn(kde-1))/grid%dn(kde-1)
262   grid%cfn1 = -.5*grid%dnw(kde-1)/grid%dn(kde-1)
263   grid%rdx = 1./config_flags%dx
264   grid%rdy = 1./config_flags%dy
265
266!  get the sounding from the ascii sounding file, first get dry sounding and
267!  calculate base state
268
269  dry_sounding = .true.
270  IF ( wrf_dm_on_monitor() ) THEN
271  write(6,*) ' getting dry sounding for base state '
272
273  CALL get_sounding( zk, p_in, pd_in, theta, rho, u, v, qv, dry_sounding, nl_max, nl_in )
274  ENDIF
275  CALL wrf_dm_bcast_real( zk , nl_max )
276  CALL wrf_dm_bcast_real( p_in , nl_max )
277  CALL wrf_dm_bcast_real( pd_in , nl_max )
278  CALL wrf_dm_bcast_real( theta , nl_max )
279  CALL wrf_dm_bcast_real( rho , nl_max )
280  CALL wrf_dm_bcast_real( u , nl_max )
281  CALL wrf_dm_bcast_real( v , nl_max )
282  CALL wrf_dm_bcast_real( qv , nl_max )
283  CALL wrf_dm_bcast_integer ( nl_in , 1 )
284
285  write(6,*) ' returned from reading sounding, nl_in is ',nl_in
286
287!  find ptop for the desired ztop (ztop is input from the namelist),
288!  and find surface pressure
289
290  grid%p_top = interp_0( p_in, zk, config_flags%ztop, nl_in )
291
292  DO j=jts,jte
293  DO i=its,ite
294    grid%ht(i,j) = 0.
295  ENDDO
296  ENDDO
297
298  xs=ide/2 -3
299  xs=ids   -3
300  xe=xs + 6
301  ys=jde/2 -3
302  ye=ys + 6
303  mtn_ht = 500
304#ifdef MTN
305  DO j=max(ys,jds),min(ye,jde-1)
306  DO i=max(xs,ids),min(xe,ide-1)
307     grid%ht(i,j) = mtn_ht * 0.25 * &
308               ( 1. + COS ( 2*pi/(xe-xs) * ( i-xs ) + pi ) ) * &
309               ( 1. + COS ( 2*pi/(ye-ys) * ( j-ys ) + pi ) )
310  ENDDO
311  ENDDO
312#endif
313#ifdef EW_RIDGE
314  DO j=max(ys,jds),min(ye,jde-1)
315  DO i=ids,ide
316     grid%ht(i,j) = mtn_ht * 0.50 * &
317               ( 1. + COS ( 2*pi/(ye-ys) * ( j-ys ) + pi ) )
318  ENDDO
319  ENDDO
320#endif
321#ifdef NS_RIDGE
322  DO j=jds,jde
323  DO i=max(xs,ids),min(xe,ide-1)
324     grid%ht(i,j) = mtn_ht * 0.50 * &
325               ( 1. + COS ( 2*pi/(xe-xs) * ( i-xs ) + pi ) )
326  ENDDO
327  ENDDO
328#endif
329  DO j=jts,jte
330  DO i=its,ite
331    grid%phb(i,1,j) = g * grid%ht(i,j)
332    grid%ph0(i,1,j) = g * grid%ht(i,j)
333  ENDDO
334  ENDDO
335
336  DO J = jts, jte
337  DO I = its, ite
338
339    p_surf = interp_0( p_in, zk, grid%phb(i,1,j)/g, nl_in )
340    grid%mub(i,j) = p_surf-grid%p_top
341
342!  this is dry hydrostatic sounding (base state), so given grid%p (coordinate),
343!  interp theta (from interp) and compute 1/rho from eqn. of state
344
345    DO K = 1, kte-1
346      p_level = grid%znu(k)*(p_surf - grid%p_top) + grid%p_top
347      grid%pb(i,k,j) = p_level
348      grid%t_init(i,k,j) = interp_0( theta, p_in, p_level, nl_in ) - t0
349      grid%alb(i,k,j) = (r_d/p1000mb)*(grid%t_init(i,k,j)+t0)*(grid%pb(i,k,j)/p1000mb)**cvpm
350    ENDDO
351
352!  calc hydrostatic balance (alternatively we could interp the geopotential from the
353!  sounding, but this assures that the base state is in exact hydrostatic balance with
354!  respect to the model eqns.
355
356    DO k  = 2,kte
357      grid%phb(i,k,j) = grid%phb(i,k-1,j) - grid%dnw(k-1)*grid%mub(i,j)*grid%alb(i,k-1,j)
358    ENDDO
359
360  ENDDO
361  ENDDO
362
363  IF ( wrf_dm_on_monitor() ) THEN
364    write(6,*) ' ptop is ',grid%p_top
365    write(6,*) ' base state grid%mub(1,1), p_surf is ',grid%mub(1,1),grid%mub(1,1)+grid%p_top
366  ENDIF
367
368!  calculate full state for each column - this includes moisture.
369
370  write(6,*) ' getting moist sounding for full state '
371  dry_sounding = .false.
372  CALL get_sounding( zk, p_in, pd_in, theta, rho, u, v, qv, dry_sounding, nl_max, nl_in )
373
374  DO J = jts, min(jde-1,jte)
375  DO I = its, min(ide-1,ite)
376
377!  At this point grid%p_top is already set. find the DRY mass in the column
378!  by interpolating the DRY pressure. 
379
380   pd_surf = interp_0( pd_in, zk, grid%phb(i,1,j)/g, nl_in )
381
382!  compute the perturbation mass and the full mass
383
384    grid%mu_1(i,j) = pd_surf-grid%p_top - grid%mub(i,j)
385    grid%mu_2(i,j) = grid%mu_1(i,j)
386    grid%mu0(i,j) = grid%mu_1(i,j) + grid%mub(i,j)
387
388! given the dry pressure and coordinate system, interp the potential
389! temperature and qv
390
391    do k=1,kde-1
392
393      p_level = grid%znu(k)*(pd_surf - grid%p_top) + grid%p_top
394
395      moist(i,k,j,P_QV) = interp_0( qv, pd_in, p_level, nl_in )
396      grid%t_1(i,k,j)          = interp_0( theta, pd_in, p_level, nl_in ) - t0
397      grid%t_2(i,k,j)          = grid%t_1(i,k,j)
398     
399
400    enddo
401
402!  integrate the hydrostatic equation (from the RHS of the bigstep
403!  vertical momentum equation) down from the top to get grid%p.
404!  first from the top of the model to the top pressure
405
406    k = kte-1  ! top level
407
408    qvf1 = 0.5*(moist(i,k,j,P_QV)+moist(i,k,j,P_QV))
409    qvf2 = 1./(1.+qvf1)
410    qvf1 = qvf1*qvf2
411
412!    grid%p(i,k,j) = - 0.5*grid%mu_1(i,j)/grid%rdnw(k)
413    grid%p(i,k,j) = - 0.5*(grid%mu_1(i,j)+qvf1*grid%mub(i,j))/grid%rdnw(k)/qvf2
414    qvf = 1. + rvovrd*moist(i,k,j,P_QV)
415    grid%alt(i,k,j) = (r_d/p1000mb)*(grid%t_1(i,k,j)+t0)*qvf* &
416                (((grid%p(i,k,j)+grid%pb(i,k,j))/p1000mb)**cvpm)
417    grid%al(i,k,j) = grid%alt(i,k,j) - grid%alb(i,k,j)
418
419!  down the column
420
421    do k=kte-2,1,-1
422      qvf1 = 0.5*(moist(i,k,j,P_QV)+moist(i,k+1,j,P_QV))
423      qvf2 = 1./(1.+qvf1)
424      qvf1 = qvf1*qvf2
425      grid%p(i,k,j) = grid%p(i,k+1,j) - (grid%mu_1(i,j) + qvf1*grid%mub(i,j))/qvf2/grid%rdn(k+1)
426      qvf = 1. + rvovrd*moist(i,k,j,P_QV)
427      grid%alt(i,k,j) = (r_d/p1000mb)*(grid%t_1(i,k,j)+t0)*qvf* &
428                  (((grid%p(i,k,j)+grid%pb(i,k,j))/p1000mb)**cvpm)
429      grid%al(i,k,j) = grid%alt(i,k,j) - grid%alb(i,k,j)
430    enddo
431
432!  this is the hydrostatic equation used in the model after the
433!  small timesteps.  In the model, grid%al (inverse density)
434!  is computed from the geopotential.
435
436
437    grid%ph_1(i,1,j) = 0.
438    DO k  = 2,kte
439      grid%ph_1(i,k,j) = grid%ph_1(i,k-1,j) - (1./grid%rdnw(k-1))*(       &
440                   (grid%mub(i,j)+grid%mu_1(i,j))*grid%al(i,k-1,j)+ &
441                    grid%mu_1(i,j)*grid%alb(i,k-1,j)  )
442                                                   
443      grid%ph_2(i,k,j) = grid%ph_1(i,k,j)
444      grid%ph0(i,k,j) = grid%ph_1(i,k,j) + grid%phb(i,k,j)
445    ENDDO
446
447    IF ( wrf_dm_on_monitor() ) THEN
448    if((i==2) .and. (j==2)) then
449     write(6,*) ' grid%ph_1 calc ',grid%ph_1(2,1,2),grid%ph_1(2,2,2),&
450                              grid%mu_1(2,2)+grid%mub(2,2),grid%mu_1(2,2), &
451                              grid%alb(2,1,2),grid%al(1,2,1),grid%rdnw(1)
452    endif
453    ENDIF
454
455  ENDDO
456  ENDDO
457
458!#if 0
459
460!  thermal perturbation to kick off convection
461
462  write(6,*) ' nxc, nyc for perturbation ',nxc,nyc
463  write(6,*) ' delt for perturbation ',delt
464
465! For LES, change the initial random perturbations
466! For 2D test, call randx outside I-loop
467! For 3D runs, call randx inside both I-J loops
468
469  DO J = jts, min(jde-1,jte)
470!   yrad = config_flags%dy*float(j-nyc)/10000.
471    yrad = 0.
472    DO I = its, min(ide-1,ite)
473!     xrad = config_flags%dx*float(i-nxc)/10000.
474      xrad = 0.
475      call random_number (randx)
476      randx = randx - 0.5
477!     DO K = 1, kte-1
478      DO K = 1, 4
479
480!  No bubbles for LES!
481!  put in preturbation theta (bubble) and recalc density.  note,
482!  the mass in the column is not changing, so when theta changes,
483!  we recompute density and geopotential
484
485!       zrad = 0.5*(grid%ph_1(i,k,j)+grid%ph_1(i,k+1,j)  &
486!                  +grid%phb(i,k,j)+grid%phb(i,k+1,j))/g
487!       zrad = (zrad-1500.)/1500.
488        zrad = 0.
489        RAD=SQRT(xrad*xrad+yrad*yrad+zrad*zrad)
490        IF(RAD <= 1.) THEN
491!          grid%t_1(i,k,j)=grid%t_1(i,k,j)+delt*COS(.5*PI*RAD)**2
492           grid%t_1(i,k,j)=grid%t_1(i,k,j)+ 0.1 *randx
493           grid%t_2(i,k,j)=grid%t_1(i,k,j)
494           qvf = 1. + rvovrd*moist(i,k,j,P_QV)
495           grid%alt(i,k,j) = (r_d/p1000mb)*(grid%t_1(i,k,j)+t0)*qvf* &
496                        (((grid%p(i,k,j)+grid%pb(i,k,j))/p1000mb)**cvpm)
497           grid%al(i,k,j) = grid%alt(i,k,j) - grid%alb(i,k,j)
498        ENDIF
499      ENDDO
500
501!  rebalance hydrostatically
502
503      DO k  = 2,kte
504        grid%ph_1(i,k,j) = grid%ph_1(i,k-1,j) - (1./grid%rdnw(k-1))*(       &
505                     (grid%mub(i,j)+grid%mu_1(i,j))*grid%al(i,k-1,j)+ &
506                      grid%mu_1(i,j)*grid%alb(i,k-1,j)  )
507                                                   
508        grid%ph_2(i,k,j) = grid%ph_1(i,k,j)
509        grid%ph0(i,k,j) = grid%ph_1(i,k,j) + grid%phb(i,k,j)
510      ENDDO
511
512    ENDDO
513  ENDDO
514
515!#endif
516
517   IF ( wrf_dm_on_monitor() ) THEN
518   write(6,*) ' grid%mu_1 from comp ', grid%mu_1(1,1)
519   write(6,*) ' full state sounding from comp, ph, grid%p, grid%al, grid%t_1, qv '
520   do k=1,kde-1
521     write(6,'(i3,1x,5(1x,1pe10.3))') k, grid%ph_1(1,k,1)+grid%phb(1,k,1), &
522                                      grid%p(1,k,1)+grid%pb(1,k,1), grid%alt(1,k,1), &
523                                      grid%t_1(1,k,1)+t0, moist(1,k,1,P_QV)
524   enddo
525
526   write(6,*) ' pert state sounding from comp, grid%ph_1, pp, alp, grid%t_1, qv '
527   do k=1,kde-1
528     write(6,'(i3,1x,5(1x,1pe10.3))') k, grid%ph_1(1,k,1), &
529                                      grid%p(1,k,1), grid%al(1,k,1), &
530                                      grid%t_1(1,k,1), moist(1,k,1,P_QV)
531   enddo
532   ENDIF
533
534! interp v
535
536  DO J = jts, jte
537  DO I = its, min(ide-1,ite)
538
539    IF (j == jds) THEN
540      z_at_v = grid%phb(i,1,j)/g
541    ELSE IF (j == jde) THEN
542      z_at_v = grid%phb(i,1,j-1)/g
543    ELSE
544      z_at_v = 0.5*(grid%phb(i,1,j)+grid%phb(i,1,j-1))/g
545    END IF
546    p_surf = interp_0( p_in, zk, z_at_v, nl_in )
547
548    DO K = 1, kte-1
549      p_level = grid%znu(k)*(p_surf - grid%p_top) + grid%p_top
550      grid%v_1(i,k,j) = interp_0( v, p_in, p_level, nl_in )
551      grid%v_2(i,k,j) = grid%v_1(i,k,j)
552    ENDDO
553
554  ENDDO
555  ENDDO
556
557! interp u
558
559  DO J = jts, min(jde-1,jte)
560  DO I = its, ite
561
562    IF (i == ids) THEN
563      z_at_u = grid%phb(i,1,j)/g
564    ELSE IF (i == ide) THEN
565      z_at_u = grid%phb(i-1,1,j)/g
566    ELSE
567      z_at_u = 0.5*(grid%phb(i,1,j)+grid%phb(i-1,1,j))/g
568    END IF
569
570    p_surf = interp_0( p_in, zk, z_at_u, nl_in )
571
572    DO K = 1, kte-1
573      p_level = grid%znu(k)*(p_surf - grid%p_top) + grid%p_top
574      grid%u_1(i,k,j) = interp_0( u, p_in, p_level, nl_in )
575      grid%u_2(i,k,j) = grid%u_1(i,k,j)
576    ENDDO
577
578  ENDDO
579  ENDDO
580
581!  set w
582
583  DO J = jts, min(jde-1,jte)
584  DO K = kts, kte
585  DO I = its, min(ide-1,ite)
586    grid%w_1(i,k,j) = 0.
587    grid%w_2(i,k,j) = 0.
588  ENDDO
589  ENDDO
590  ENDDO
591
592!  set a few more things
593
594  DO J = jts, min(jde-1,jte)
595  DO K = kts, kte-1
596  DO I = its, min(ide-1,ite)
597    grid%h_diabatic(i,k,j) = 0.
598  ENDDO
599  ENDDO
600  ENDDO
601
602  IF ( wrf_dm_on_monitor() ) THEN
603  DO k=1,kte-1
604    grid%t_base(k) = grid%t_1(1,k,1)
605    grid%qv_base(k) = moist(1,k,1,P_QV)
606    grid%u_base(k) = grid%u_1(1,k,1)
607    grid%v_base(k) = grid%v_1(1,k,1)
608    grid%z_base(k) = 0.5*(grid%phb(1,k,1)+grid%phb(1,k+1,1)+grid%ph_1(1,k,1)+grid%ph_1(1,k+1,1))/g
609  ENDDO
610  ENDIF
611  CALL wrf_dm_bcast_real( grid%t_base , kte )
612  CALL wrf_dm_bcast_real( grid%qv_base , kte )
613  CALL wrf_dm_bcast_real( grid%u_base , kte )
614  CALL wrf_dm_bcast_real( grid%v_base , kte )
615  CALL wrf_dm_bcast_real( grid%z_base , kte )
616
617  DO J = jts, min(jde-1,jte)
618  DO I = its, min(ide-1,ite)
619     thtmp   = grid%t_2(i,1,j)+t0
620     ptmp    = grid%p(i,1,j)+grid%pb(i,1,j)
621     temp(1) = thtmp * (ptmp/p1000mb)**rcp
622     thtmp   = grid%t_2(i,2,j)+t0
623     ptmp    = grid%p(i,2,j)+grid%pb(i,2,j)
624     temp(2) = thtmp * (ptmp/p1000mb)**rcp
625     thtmp   = grid%t_2(i,3,j)+t0
626     ptmp    = grid%p(i,3,j)+grid%pb(i,3,j)
627     temp(3) = thtmp * (ptmp/p1000mb)**rcp
628
629!    For LES-CBL, add 5 degrees to the surface temperature!
630!
631!    grid%tsk(I,J)=grid%cf1*temp(1)+grid%cf2*temp(2)+grid%cf3*temp(3)
632     grid%tsk(I,J)=grid%cf1*temp(1)+grid%cf2*temp(2)+grid%cf3*temp(3)+5.
633     grid%tmn(I,J)=grid%tsk(I,J)-0.5
634  ENDDO
635  ENDDO
636
637 END SUBROUTINE init_domain_rk
638
639   SUBROUTINE init_module_initialize
640   END SUBROUTINE init_module_initialize
641
642!---------------------------------------------------------------------
643
644!  test driver for get_sounding
645!
646!      implicit none
647!      integer n
648!      parameter(n = 1000)
649!      real zk(n),p(n),theta(n),rho(n),u(n),v(n),qv(n),pd(n)
650!      logical dry
651!      integer nl,k
652!
653!      dry = .false.
654!      dry = .true.
655!      call get_sounding( zk, p, pd, theta, rho, u, v, qv, dry, n, nl )
656!      write(6,*) ' input levels ',nl
657!      write(6,*) ' sounding '
658!      write(6,*) '  k  height(m)  press (Pa) pd(Pa) theta (K) den(kg/m^3)  u(m/s)     v(m/s)    qv(g/g) '
659!      do k=1,nl
660!        write(6,'(1x,i3,8(1x,1pe10.3))') k, zk(k), p(k), pd(k), theta(k), rho(k), u(k), v(k), qv(k)
661!      enddo
662!      end
663!
664!---------------------------------------------------------------------------
665
666      subroutine get_sounding( zk, p, p_dry, theta, rho, &
667                               u, v, qv, dry, nl_max, nl_in )
668      implicit none
669
670      integer nl_max, nl_in
671      real zk(nl_max), p(nl_max), theta(nl_max), rho(nl_max), &
672           u(nl_max), v(nl_max), qv(nl_max), p_dry(nl_max)
673      logical dry
674
675      integer n
676      parameter(n=1000)
677      logical debug
678      parameter( debug = .true.)
679
680! input sounding data
681
682      real p_surf, th_surf, qv_surf
683      real pi_surf, pi(n)
684      real h_input(n), th_input(n), qv_input(n), u_input(n), v_input(n)
685
686! diagnostics
687
688      real rho_surf, p_input(n), rho_input(n)
689      real pm_input(n)  !  this are for full moist sounding
690
691! local data
692
693      real r
694      parameter (r = r_d)
695      integer k, it, nl
696      real qvf, qvf1, dz
697
698!  first, read the sounding
699
700      call read_sounding( p_surf, th_surf, qv_surf, &
701                          h_input, th_input, qv_input, u_input, v_input,n, nl, debug )
702
703      if(dry) then
704       do k=1,nl
705         qv_input(k) = 0.
706       enddo
707      endif
708
709      if(debug) write(6,*) ' number of input levels = ',nl
710
711        nl_in = nl
712        if(nl_in .gt. nl_max ) then
713          write(6,*) ' too many levels for input arrays ',nl_in,nl_max
714          call wrf_error_fatal ( ' too many levels for input arrays ' )
715        end if
716
717!  compute diagnostics,
718!  first, convert qv(g/kg) to qv(g/g)
719
720      do k=1,nl
721        qv_input(k) = 0.001*qv_input(k)
722      enddo
723
724      p_surf = 100.*p_surf  ! convert to pascals
725      qvf = 1. + rvovrd*qv_input(1)
726      rho_surf = 1./((r/p1000mb)*th_surf*qvf*((p_surf/p1000mb)**cvpm))
727      pi_surf = (p_surf/p1000mb)**(r/cp)
728
729      if(debug) then
730        write(6,*) ' surface density is ',rho_surf
731        write(6,*) ' surface pi is      ',pi_surf
732      end if
733
734
735!  integrate moist sounding hydrostatically, starting from the
736!  specified surface pressure
737!  -> first, integrate from surface to lowest level
738
739          qvf = 1. + rvovrd*qv_input(1)
740          qvf1 = 1. + qv_input(1)
741          rho_input(1) = rho_surf
742          dz = h_input(1)
743          do it=1,10
744            pm_input(1) = p_surf &
745                    - 0.5*dz*(rho_surf+rho_input(1))*g*qvf1
746            rho_input(1) = 1./((r/p1000mb)*th_input(1)*qvf*((pm_input(1)/p1000mb)**cvpm))
747          enddo
748
749! integrate up the column
750
751          do k=2,nl
752            rho_input(k) = rho_input(k-1)
753            dz = h_input(k)-h_input(k-1)
754            qvf1 = 0.5*(2.+(qv_input(k-1)+qv_input(k)))
755            qvf = 1. + rvovrd*qv_input(k)   ! qv is in g/kg here
756 
757            do it=1,10
758              pm_input(k) = pm_input(k-1) &
759                      - 0.5*dz*(rho_input(k)+rho_input(k-1))*g*qvf1
760              rho_input(k) = 1./((r/p1000mb)*th_input(k)*qvf*((pm_input(k)/p1000mb)**cvpm))
761            enddo
762          enddo
763
764!  we have the moist sounding
765
766!  next, compute the dry sounding using p at the highest level from the
767!  moist sounding and integrating down.
768
769        p_input(nl) = pm_input(nl)
770
771          do k=nl-1,1,-1
772            dz = h_input(k+1)-h_input(k)
773            p_input(k) = p_input(k+1) + 0.5*dz*(rho_input(k)+rho_input(k+1))*g
774          enddo
775
776
777        do k=1,nl
778
779          zk(k) = h_input(k)
780          p(k) = pm_input(k)
781          p_dry(k) = p_input(k)
782          theta(k) = th_input(k)
783          rho(k) = rho_input(k)
784          u(k) = u_input(k)
785          v(k) = v_input(k)
786          qv(k) = qv_input(k)
787
788        enddo
789
790     if(debug) then
791      write(6,*) ' sounding '
792      write(6,*) '  k  height(m)  press (Pa) pd(Pa) theta (K) den(kg/m^3)  u(m/s)     v(m/s)    qv(g/g) '
793      do k=1,nl
794        write(6,'(1x,i3,8(1x,1pe10.3))') k, zk(k), p(k), p_dry(k), theta(k), rho(k), u(k), v(k), qv(k)
795      enddo
796
797     end if
798
799      end subroutine get_sounding
800
801!-------------------------------------------------------
802
803      subroutine read_sounding( ps,ts,qvs,h,th,qv,u,v,n,nl,debug )
804      implicit none
805      integer n,nl
806      real ps,ts,qvs,h(n),th(n),qv(n),u(n),v(n)
807      logical end_of_file
808      logical debug
809
810      integer k
811
812      open(unit=10,file='input_sounding',form='formatted',status='old')
813      rewind(10)
814      read(10,*) ps, ts, qvs
815      if(debug) then
816        write(6,*) ' input sounding surface parameters '
817        write(6,*) ' surface pressure (mb) ',ps
818        write(6,*) ' surface pot. temp (K) ',ts
819        write(6,*) ' surface mixing ratio (g/kg) ',qvs
820      end if
821
822      end_of_file = .false.
823      k = 0
824
825      do while (.not. end_of_file)
826
827        read(10,*,end=100) h(k+1), th(k+1), qv(k+1), u(k+1), v(k+1)
828        k = k+1
829        if(debug) write(6,'(1x,i3,5(1x,e10.3))') k, h(k), th(k), qv(k), u(k), v(k)
830        go to 110
831 100    end_of_file = .true.
832 110    continue
833      enddo
834
835      nl = k
836
837      close(unit=10,status = 'keep')
838
839      end subroutine read_sounding
840
841END MODULE module_initialize_ideal
Note: See TracBrowser for help on using the repository browser.