source: trunk/WRF.COMMON/WRFV2/dyn_em/module_initialize_quarter_ss.F_ref @ 3567

Last change on this file since 3567 was 31, checked in by aslmd, 14 years ago

LMD_MM_MARS: teste compilation, scripts additionnels et initialisation nouvelle physique

File size: 26.0 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
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                :: dyn_opt
51   INTEGER :: idum1, idum2
52
53   CALL nl_get_dyn_opt( 1, dyn_opt )
54   
55   CALL set_scalar_indices_from_config ( head_grid%id , idum1, idum2 )
56
57   IF (      dyn_opt .eq. 1 &
58        .or. dyn_opt .eq. 2 &
59        .or. dyn_opt .eq. 3 &
60                                       ) THEN
61     CALL init_domain_rk( grid &
62!
63#include <em_actual_new_args.inc>
64!
65                        )
66
67   ELSE
68     WRITE(0,*)' init_domain: unknown or unimplemented dyn_opt = ',dyn_opt
69     CALL wrf_error_fatal ( ' init_domain: unknown or unimplemented dyn_opt ' )
70   ENDIF
71
72   END SUBROUTINE init_domain
73
74!-------------------------------------------------------------------
75
76   SUBROUTINE init_domain_rk ( grid &
77!
78# include <em_dummy_new_args.inc>
79!
80)
81   IMPLICIT NONE
82
83   !  Input data.
84   TYPE (domain), POINTER :: grid
85
86# include <em_dummy_new_decl.inc>
87
88   TYPE (grid_config_rec_type)              :: config_flags
89
90   !  Local data
91   INTEGER                             ::                       &
92                                  ids, ide, jds, jde, kds, kde, &
93                                  ims, ime, jms, jme, kms, kme, &
94                                  its, ite, jts, jte, kts, kte, &
95                                  i, j, k
96
97   ! Local data
98
99   INTEGER, PARAMETER :: nl_max = 1000
100   REAL, DIMENSION(nl_max) :: zk, p_in, theta, rho, u, v, qv, pd_in
101   INTEGER :: nl_in
102
103
104   INTEGER :: icm,jcm, ii, im1, jj, jm1, loop, error, fid, nxc, nyc
105   REAL    :: u_mean,v_mean, f0, p_surf, p_level, qvf, z_at_v, z_at_u
106   REAL    :: z_scale, xrad, yrad, zrad, rad, delt, cof1, cof2
107!   REAL, EXTERNAL :: interp_0
108   REAL    :: hm
109   REAL    :: pi
110
111!  stuff from original initialization that has been dropped from the Registry
112   REAL    :: vnu, xnu, xnus, dinit0, cbh, p0_temp, t0_temp, zd, zt
113   REAL    :: qvf1, qvf2, pd_surf
114   INTEGER :: it
115   real :: thtmp, ptmp, temp(3)
116
117   LOGICAL :: moisture_init
118   LOGICAL :: stretch_grid, dry_sounding
119
120  INTEGER :: xs , xe , ys , ye
121  REAL :: mtn_ht
122   LOGICAL, EXTERNAL :: wrf_dm_on_monitor
123
124!!MARS
125 REAL :: lon_input, lat_input, alt_input, tsurf_input
126!!MARS
127
128
129#ifdef DM_PARALLEL
130#    include <em_data_calls.inc>
131#endif
132
133
134   SELECT CASE ( model_data_order )
135         CASE ( DATA_ORDER_ZXY )
136   kds = grid%sd31 ; kde = grid%ed31 ;
137   ids = grid%sd32 ; ide = grid%ed32 ;
138   jds = grid%sd33 ; jde = grid%ed33 ;
139
140   kms = grid%sm31 ; kme = grid%em31 ;
141   ims = grid%sm32 ; ime = grid%em32 ;
142   jms = grid%sm33 ; jme = grid%em33 ;
143
144   kts = grid%sp31 ; kte = grid%ep31 ;   ! note that tile is entire patch
145   its = grid%sp32 ; ite = grid%ep32 ;   ! note that tile is entire patch
146   jts = grid%sp33 ; jte = grid%ep33 ;   ! note that tile is entire patch
147         CASE ( DATA_ORDER_XYZ )
148   ids = grid%sd31 ; ide = grid%ed31 ;
149   jds = grid%sd32 ; jde = grid%ed32 ;
150   kds = grid%sd33 ; kde = grid%ed33 ;
151
152   ims = grid%sm31 ; ime = grid%em31 ;
153   jms = grid%sm32 ; jme = grid%em32 ;
154   kms = grid%sm33 ; kme = grid%em33 ;
155
156   its = grid%sp31 ; ite = grid%ep31 ;   ! note that tile is entire patch
157   jts = grid%sp32 ; jte = grid%ep32 ;   ! note that tile is entire patch
158   kts = grid%sp33 ; kte = grid%ep33 ;   ! note that tile is entire patch
159         CASE ( DATA_ORDER_XZY )
160   ids = grid%sd31 ; ide = grid%ed31 ;
161   kds = grid%sd32 ; kde = grid%ed32 ;
162   jds = grid%sd33 ; jde = grid%ed33 ;
163
164   ims = grid%sm31 ; ime = grid%em31 ;
165   kms = grid%sm32 ; kme = grid%em32 ;
166   jms = grid%sm33 ; jme = grid%em33 ;
167
168   its = grid%sp31 ; ite = grid%ep31 ;   ! note that tile is entire patch
169   kts = grid%sp32 ; kte = grid%ep32 ;   ! note that tile is entire patch
170   jts = grid%sp33 ; jte = grid%ep33 ;   ! note that tile is entire patch
171
172   END SELECT
173
174
175   stretch_grid = .true.
176!   stretch_grid = .false.
177
178   delt = 3.
179!   delt = 10.
180!   z_scale = .50
181   z_scale = .40
182
183!!!MARS
184!!!MARS
185!  open(unit=16,file='input_vert',form='formatted',status='old')
186!  rewind(16)
187!  read(16,*) delt, z_scale
188!  write(6,*) 'delt, z_scale are  ', delt, z_scale
189!  close(16)
190!!!MARS
191!!!MARS
192
193   pi = 2.*asin(1.0)
194   write(6,*) ' pi is ',pi
195   nxc = (ide-ids)/2
196   nyc = (jde-jds)/2
197
198   CALL model_to_grid_config_rec ( grid%id , model_config_rec , config_flags )
199
200! here we check to see if the boundary conditions are set properly
201
202   CALL boundary_condition_check( config_flags, bdyzone, error, grid%id )
203
204   moisture_init = .true.
205
206    grid%itimestep=0
207
208#ifdef DM_PARALLEL
209   CALL wrf_dm_bcast_bytes( icm , IWORDSIZE )
210   CALL wrf_dm_bcast_bytes( jcm , IWORDSIZE )
211#endif
212
213    CALL nl_set_mminlu(1, '    ')
214    CALL nl_set_iswater(1,0)
215    CALL nl_set_cen_lat(1,40.)
216    CALL nl_set_cen_lon(1,-105.)
217    CALL nl_set_truelat1(1,0.)
218    CALL nl_set_truelat2(1,0.)
219    CALL nl_set_moad_cen_lat (1,0.)
220    CALL nl_set_stand_lon (1,0.)
221    CALL nl_set_map_proj(1,0)
222
223
224!  here we initialize data we currently is not initialized
225!  in the input data
226
227    DO j = jts, jte
228      DO i = its, ite
229         grid%msft(i,j)     = 1.
230         grid%msfu(i,j)     = 1.
231         grid%msfv(i,j)     = 1.
232         grid%sina(i,j)     = 0.
233         grid%cosa(i,j)     = 1.
234         grid%e(i,j)        = 0.
235         grid%f(i,j)        = 0.         !! MARS: put coriolis here if needed
236
237      END DO
238   END DO
239
240    DO j = jts, jte
241    DO k = kts, kte
242      DO i = its, ite
243         grid%em_ww(i,k,j)     = 0.
244      END DO
245   END DO
246   END DO
247
248   grid%step_number = 0
249
250! set up the grid
251
252   IF (stretch_grid) THEN ! exponential stretch for eta (nearly constant dz)
253     DO k=1, kde
254      grid%em_znw(k) = (exp(-(k-1)/float(kde-1)/z_scale) - exp(-1./z_scale))/ &
255                                (1.-exp(-1./z_scale))
256     ENDDO
257   ELSE
258     DO k=1, kde
259      grid%em_znw(k) = 1. - float(k-1)/float(kde-1)
260     ENDDO
261   ENDIF
262
263!!!MARS
264!!!MARS
265!  open(unit=12,file='levels',form='formatted',status='old')
266!  rewind(12)
267!  DO k=1, kde
268!  read(12,*) grid%em_znw(k)
269!  write(6,*) 'read level ', k,grid%em_znw(k)
270!  ENDDO
271!  close(12)
272!!!MARS
273!!!MARS
274   
275
276   DO k=1, kde-1
277    grid%em_dnw(k) = grid%em_znw(k+1) - grid%em_znw(k)
278    grid%em_rdnw(k) = 1./grid%em_dnw(k)
279    grid%em_znu(k) = 0.5*(grid%em_znw(k+1)+grid%em_znw(k))
280   ENDDO
281   DO k=2, kde-1
282    grid%em_dn(k) = 0.5*(grid%em_dnw(k)+grid%em_dnw(k-1))
283    grid%em_rdn(k) = 1./grid%em_dn(k)
284    grid%em_fnp(k) = .5* grid%em_dnw(k  )/grid%em_dn(k)
285    grid%em_fnm(k) = .5* grid%em_dnw(k-1)/grid%em_dn(k)
286   ENDDO
287
288   cof1 = (2.*grid%em_dn(2)+grid%em_dn(3))/(grid%em_dn(2)+grid%em_dn(3))*grid%em_dnw(1)/grid%em_dn(2)
289   cof2 =     grid%em_dn(2)        /(grid%em_dn(2)+grid%em_dn(3))*grid%em_dnw(1)/grid%em_dn(3)
290   grid%cf1  = grid%em_fnp(2) + cof1
291   grid%cf2  = grid%em_fnm(2) - cof1 - cof2
292   grid%cf3  = cof2       
293
294   grid%cfn  = (.5*grid%em_dnw(kde-1)+grid%em_dn(kde-1))/grid%em_dn(kde-1)
295   grid%cfn1 = -.5*grid%em_dnw(kde-1)/grid%em_dn(kde-1)
296   grid%rdx = 1./config_flags%dx
297   grid%rdy = 1./config_flags%dy
298
299!  get the sounding from the ascii sounding file, first get dry sounding and
300!  calculate base state
301
302  dry_sounding = .true.
303  IF ( wrf_dm_on_monitor() ) THEN
304  write(6,*) ' getting dry sounding for base state '
305
306  CALL get_sounding( zk, p_in, pd_in, theta, rho, u, v, qv, dry_sounding, nl_max, nl_in )
307  ENDIF
308  CALL wrf_dm_bcast_real( zk , nl_max )
309  CALL wrf_dm_bcast_real( p_in , nl_max )
310  CALL wrf_dm_bcast_real( pd_in , nl_max )
311  CALL wrf_dm_bcast_real( theta , nl_max )
312  CALL wrf_dm_bcast_real( rho , nl_max )
313  CALL wrf_dm_bcast_real( u , nl_max )
314  CALL wrf_dm_bcast_real( v , nl_max )
315  CALL wrf_dm_bcast_real( qv , nl_max )
316  CALL wrf_dm_bcast_integer ( nl_in , 1 )
317
318  write(6,*) ' returned from reading sounding, nl_in is ',nl_in
319
320!  find ptop for the desired ztop (ztop is input from the namelist),
321!  and find surface pressure
322
323  grid%p_top = interp_0( p_in, zk, config_flags%ztop, nl_in )
324
325!!MARS
326!!MARS
327  open(unit=14,file='input_coord',form='formatted',status='old')
328  rewind(14)
329  read(14,*) lon_input
330  read(14,*) lat_input
331  close(14)
332  write(6,*) ' lon is ',lon_input
333  write(6,*) ' lat is ',lat_input
334!!MARS
335!!MARS
336
337!!MARS
338!!MARS
339  open(unit=18,file='input_more',form='formatted',status='old')
340  rewind(18)
341  read(18,*) alt_input, tsurf_input
342  close(18)
343  write(6,*) ' alt is ',alt_input
344  write(6,*) ' tsurf is ',tsurf_input
345!!MARS
346!!MARS
347
348  DO j=jts,jte
349  DO i=its,ite
350!!MARS
351    grid%ht(i,j) = alt_input
352    grid%tsk(i,j) = tsurf_input
353!!MARS
354    grid%xlat(i,j) = lat_input
355    grid%xlong(i,j) = lon_input!+float(i)*config_flags%dx/59000.
356    grid%mars_emiss(i,j)=0.95
357    grid%mars_cice(i,j)=0.
358    grid%slpx(i,j) = 0.
359    grid%slpy(i,j) = 0.
360!!MARS
361  ENDDO
362  ENDDO
363
364  xs=ide/2 -3
365  xs=ids   -3
366  xe=xs + 6
367  ys=jde/2 -3
368  ye=ys + 6
369  mtn_ht = 500
370#ifdef MTN
371  DO j=max(ys,jds),min(ye,jde-1)
372  DO i=max(xs,ids),min(xe,ide-1)
373     grid%ht(i,j) = mtn_ht * 0.25 * &
374               ( 1. + COS ( 2*pi/(xe-xs) * ( i-xs ) + pi ) ) * &
375               ( 1. + COS ( 2*pi/(ye-ys) * ( j-ys ) + pi ) )
376  ENDDO
377  ENDDO
378#endif
379#ifdef EW_RIDGE
380  DO j=max(ys,jds),min(ye,jde-1)
381  DO i=ids,ide
382     grid%ht(i,j) = mtn_ht * 0.50 * &
383               ( 1. + COS ( 2*pi/(ye-ys) * ( j-ys ) + pi ) )
384  ENDDO
385  ENDDO
386#endif
387#ifdef NS_RIDGE
388  DO j=jds,jde
389  DO i=max(xs,ids),min(xe,ide-1)
390     grid%ht(i,j) = mtn_ht * 0.50 * &
391               ( 1. + COS ( 2*pi/(xe-xs) * ( i-xs ) + pi ) )
392  ENDDO
393  ENDDO
394#endif
395  DO j=jts,jte
396  DO i=its,ite
397    grid%em_phb(i,1,j) = g * grid%ht(i,j)
398    grid%em_ph0(i,1,j) = g * grid%ht(i,j)
399  ENDDO
400  ENDDO
401
402  DO J = jts, jte
403  DO I = its, ite
404
405    p_surf = interp_0( p_in, zk, grid%em_phb(i,1,j)/g, nl_in )
406    grid%em_mub(i,j) = p_surf-grid%p_top
407
408!  this is dry hydrostatic sounding (base state), so given grid%em_p (coordinate),
409!  interp theta (from interp) and compute 1/rho from eqn. of state
410
411    DO K = 1, kte-1
412      p_level = grid%em_znu(k)*(p_surf - grid%p_top) + grid%p_top
413      grid%em_pb(i,k,j) = p_level
414      grid%em_t_init(i,k,j) = interp_0( theta, p_in, p_level, nl_in ) - t0
415      grid%em_alb(i,k,j) = (r_d/p1000mb)*(grid%em_t_init(i,k,j)+t0)*(grid%em_pb(i,k,j)/p1000mb)**cvpm
416    ENDDO
417
418!  calc hydrostatic balance (alternatively we could interp the geopotential from the
419!  sounding, but this assures that the base state is in exact hydrostatic balance with
420!  respect to the model eqns.
421
422    DO k  = 2,kte
423      grid%em_phb(i,k,j) = grid%em_phb(i,k-1,j) - grid%em_dnw(k-1)*grid%em_mub(i,j)*grid%em_alb(i,k-1,j)
424    ENDDO
425
426  ENDDO
427  ENDDO
428
429  IF ( wrf_dm_on_monitor() ) THEN
430    write(6,*) ' ptop is ',grid%p_top
431    write(6,*) ' base state grid%em_mub(1,1), p_surf is ',grid%em_mub(1,1),grid%em_mub(1,1)+grid%p_top
432  ENDIF
433
434!  calculate full state for each column - this includes moisture.
435
436!!!!!MARS MARS
437!  write(6,*) ' getting moist sounding for full state '
438!  dry_sounding = .false.
439  dry_sounding = .true.
440  CALL get_sounding( zk, p_in, pd_in, theta, rho, u, v, qv, dry_sounding, nl_max, nl_in )
441
442  DO J = jts, min(jde-1,jte)
443  DO I = its, min(ide-1,ite)
444
445!  At this point grid%p_top is already set. find the DRY mass in the column
446!  by interpolating the DRY pressure. 
447
448   pd_surf = interp_0( pd_in, zk, grid%em_phb(i,1,j)/g, nl_in )
449
450!  compute the perturbation mass and the full mass
451
452    grid%em_mu_1(i,j) = pd_surf-grid%p_top - grid%em_mub(i,j)
453    grid%em_mu_2(i,j) = grid%em_mu_1(i,j)
454    grid%em_mu0(i,j) = grid%em_mu_1(i,j) + grid%em_mub(i,j)
455
456! given the dry pressure and coordinate system, interp the potential
457! temperature and qv
458
459    do k=1,kde-1
460
461      p_level = grid%em_znu(k)*(pd_surf - grid%p_top) + grid%p_top
462
463      moist(i,k,j,P_QV) = interp_0( qv, pd_in, p_level, nl_in )
464      grid%em_t_1(i,k,j)          = interp_0( theta, pd_in, p_level, nl_in ) - t0
465      grid%em_t_2(i,k,j)          = grid%em_t_1(i,k,j)
466     
467
468    enddo
469
470!  integrate the hydrostatic equation (from the RHS of the bigstep
471!  vertical momentum equation) down from the top to get grid%em_p.
472!  first from the top of the model to the top pressure
473
474    k = kte-1  ! top level
475
476    qvf1 = 0.5*(moist(i,k,j,P_QV)+moist(i,k,j,P_QV))
477    qvf2 = 1./(1.+qvf1)
478    qvf1 = qvf1*qvf2
479
480!    grid%em_p(i,k,j) = - 0.5*grid%em_mu_1(i,j)/grid%em_rdnw(k)
481    grid%em_p(i,k,j) = - 0.5*(grid%em_mu_1(i,j)+qvf1*grid%em_mub(i,j))/grid%em_rdnw(k)/qvf2
482    qvf = 1. + rvovrd*moist(i,k,j,P_QV)
483    grid%em_alt(i,k,j) = (r_d/p1000mb)*(grid%em_t_1(i,k,j)+t0)*qvf* &
484                (((grid%em_p(i,k,j)+grid%em_pb(i,k,j))/p1000mb)**cvpm)
485    grid%em_al(i,k,j) = grid%em_alt(i,k,j) - grid%em_alb(i,k,j)
486
487!  down the column
488
489    do k=kte-2,1,-1
490      qvf1 = 0.5*(moist(i,k,j,P_QV)+moist(i,k+1,j,P_QV))
491      qvf2 = 1./(1.+qvf1)
492      qvf1 = qvf1*qvf2
493      grid%em_p(i,k,j) = grid%em_p(i,k+1,j) - (grid%em_mu_1(i,j) + qvf1*grid%em_mub(i,j))/qvf2/grid%em_rdn(k+1)
494      qvf = 1. + rvovrd*moist(i,k,j,P_QV)
495      grid%em_alt(i,k,j) = (r_d/p1000mb)*(grid%em_t_1(i,k,j)+t0)*qvf* &
496                  (((grid%em_p(i,k,j)+grid%em_pb(i,k,j))/p1000mb)**cvpm)
497      grid%em_al(i,k,j) = grid%em_alt(i,k,j) - grid%em_alb(i,k,j)
498    enddo
499
500!  this is the hydrostatic equation used in the model after the
501!  small timesteps.  In the model, grid%em_al (inverse density)
502!  is computed from the geopotential.
503
504
505    grid%em_ph_1(i,1,j) = 0.
506    DO k  = 2,kte
507      grid%em_ph_1(i,k,j) = grid%em_ph_1(i,k-1,j) - (1./grid%em_rdnw(k-1))*(       &
508                   (grid%em_mub(i,j)+grid%em_mu_1(i,j))*grid%em_al(i,k-1,j)+ &
509                    grid%em_mu_1(i,j)*grid%em_alb(i,k-1,j)  )
510                                                   
511      grid%em_ph_2(i,k,j) = grid%em_ph_1(i,k,j)
512      grid%em_ph0(i,k,j) = grid%em_ph_1(i,k,j) + grid%em_phb(i,k,j)
513    ENDDO
514
515    IF ( wrf_dm_on_monitor() ) THEN
516    if((i==2) .and. (j==2)) then
517     write(6,*) ' grid%em_ph_1 calc ',grid%em_ph_1(2,1,2),grid%em_ph_1(2,2,2),&
518                              grid%em_mu_1(2,2)+grid%em_mub(2,2),grid%em_mu_1(2,2), &
519                              grid%em_alb(2,1,2),grid%em_al(1,2,1),grid%em_rdnw(1)
520    endif
521    ENDIF
522
523  ENDDO
524  ENDDO
525
526!#if 0
527
528!  thermal perturbation to kick off convection
529
530  write(6,*) ' nxc, nyc for perturbation ',nxc,nyc
531  write(6,*) ' delt for perturbation ',delt
532
533  DO J = jts, min(jde-1,jte)
534    yrad = config_flags%dy*float(j-nyc)/10000.
535!   yrad = 0.
536    DO I = its, min(ide-1,ite)
537      xrad = config_flags%dx*float(i-nxc)/10000.
538!     xrad = 0.
539      DO K = 1, kte-1
540
541!  put in preturbation theta (bubble) and recalc density.  note,
542!  the mass in the column is not changing, so when theta changes,
543!  we recompute density and geopotential
544
545        zrad = 0.5*(grid%em_ph_1(i,k,j)+grid%em_ph_1(i,k+1,j)  &
546                   +grid%em_phb(i,k,j)+grid%em_phb(i,k+1,j))/g
547        zrad = (zrad-1500.)/1500.
548        RAD=SQRT(xrad*xrad+yrad*yrad+zrad*zrad)
549        IF(RAD <= 1.) THEN
550           grid%em_t_1(i,k,j)=grid%em_t_1(i,k,j)+delt*COS(.5*PI*RAD)**2
551           grid%em_t_2(i,k,j)=grid%em_t_1(i,k,j)
552           qvf = 1. + rvovrd*moist(i,k,j,P_QV)
553           grid%em_alt(i,k,j) = (r_d/p1000mb)*(grid%em_t_1(i,k,j)+t0)*qvf* &
554                        (((grid%em_p(i,k,j)+grid%em_pb(i,k,j))/p1000mb)**cvpm)
555           grid%em_al(i,k,j) = grid%em_alt(i,k,j) - grid%em_alb(i,k,j)
556        ENDIF
557      ENDDO
558
559!  rebalance hydrostatically
560
561      DO k  = 2,kte
562        grid%em_ph_1(i,k,j) = grid%em_ph_1(i,k-1,j) - (1./grid%em_rdnw(k-1))*(       &
563                     (grid%em_mub(i,j)+grid%em_mu_1(i,j))*grid%em_al(i,k-1,j)+ &
564                      grid%em_mu_1(i,j)*grid%em_alb(i,k-1,j)  )
565                                                   
566        grid%em_ph_2(i,k,j) = grid%em_ph_1(i,k,j)
567        grid%em_ph0(i,k,j) = grid%em_ph_1(i,k,j) + grid%em_phb(i,k,j)
568      ENDDO
569
570    ENDDO
571  ENDDO
572
573!#endif
574
575   IF ( wrf_dm_on_monitor() ) THEN
576   write(6,*) ' grid%em_mu_1 from comp ', grid%em_mu_1(1,1)
577   write(6,*) ' full state sounding from comp, ph, grid%em_p, grid%em_al, grid%em_t_1, qv '
578   do k=1,kde-1
579     write(6,'(i3,1x,5(1x,1pe10.3))') k, grid%em_ph_1(1,k,1)+grid%em_phb(1,k,1), &
580                                      grid%em_p(1,k,1)+grid%em_pb(1,k,1), grid%em_alt(1,k,1), &
581                                      grid%em_t_1(1,k,1)+t0, moist(1,k,1,P_QV)
582   enddo
583
584   write(6,*) ' pert state sounding from comp, grid%em_ph_1, pp, alp, grid%em_t_1, qv '
585   do k=1,kde-1
586     write(6,'(i3,1x,5(1x,1pe10.3))') k, grid%em_ph_1(1,k,1), &
587                                      grid%em_p(1,k,1), grid%em_al(1,k,1), &
588                                      grid%em_t_1(1,k,1), moist(1,k,1,P_QV)
589   enddo
590   ENDIF
591
592! interp v
593
594  DO J = jts, jte
595  DO I = its, min(ide-1,ite)
596
597    IF (j == jds) THEN
598      z_at_v = grid%em_phb(i,1,j)/g
599    ELSE IF (j == jde) THEN
600      z_at_v = grid%em_phb(i,1,j-1)/g
601    ELSE
602      z_at_v = 0.5*(grid%em_phb(i,1,j)+grid%em_phb(i,1,j-1))/g
603    END IF
604    p_surf = interp_0( p_in, zk, z_at_v, nl_in )
605
606    DO K = 1, kte-1
607      p_level = grid%em_znu(k)*(p_surf - grid%p_top) + grid%p_top
608      grid%em_v_1(i,k,j) = interp_0( v, p_in, p_level, nl_in )
609      grid%em_v_2(i,k,j) = grid%em_v_1(i,k,j)
610    ENDDO
611
612  ENDDO
613  ENDDO
614
615! interp u
616
617  DO J = jts, min(jde-1,jte)
618  DO I = its, ite
619
620    IF (i == ids) THEN
621      z_at_u = grid%em_phb(i,1,j)/g
622    ELSE IF (i == ide) THEN
623      z_at_u = grid%em_phb(i-1,1,j)/g
624    ELSE
625      z_at_u = 0.5*(grid%em_phb(i,1,j)+grid%em_phb(i-1,1,j))/g
626    END IF
627
628    p_surf = interp_0( p_in, zk, z_at_u, nl_in )
629
630    DO K = 1, kte-1
631      p_level = grid%em_znu(k)*(p_surf - grid%p_top) + grid%p_top
632      grid%em_u_1(i,k,j) = interp_0( u, p_in, p_level, nl_in )
633      grid%em_u_2(i,k,j) = grid%em_u_1(i,k,j)
634    ENDDO
635
636  ENDDO
637  ENDDO
638
639!  set w
640
641  DO J = jts, min(jde-1,jte)
642  DO K = kts, kte
643  DO I = its, min(ide-1,ite)
644    grid%em_w_1(i,k,j) = 0.
645    grid%em_w_2(i,k,j) = 0.
646  ENDDO
647  ENDDO
648  ENDDO
649
650!  set a few more things
651
652  DO J = jts, min(jde-1,jte)
653  DO K = kts, kte-1
654  DO I = its, min(ide-1,ite)
655    grid%h_diabatic(i,k,j) = 0.
656  ENDDO
657  ENDDO
658  ENDDO
659
660  IF ( wrf_dm_on_monitor() ) THEN
661  DO k=1,kte-1
662    grid%em_t_base(k) = grid%em_t_1(1,k,1)
663    grid%qv_base(k) = moist(1,k,1,P_QV)
664    grid%u_base(k) = grid%em_u_1(1,k,1)
665    grid%v_base(k) = grid%em_v_1(1,k,1)
666    grid%z_base(k) = 0.5*(grid%em_phb(1,k,1)+grid%em_phb(1,k+1,1)+grid%em_ph_1(1,k,1)+grid%em_ph_1(1,k+1,1))/g
667  ENDDO
668  ENDIF
669  CALL wrf_dm_bcast_real( grid%em_t_base , kte )
670  CALL wrf_dm_bcast_real( grid%qv_base , kte )
671  CALL wrf_dm_bcast_real( grid%u_base , kte )
672  CALL wrf_dm_bcast_real( grid%v_base , kte )
673  CALL wrf_dm_bcast_real( grid%z_base , kte )
674
675  DO J = jts, min(jde-1,jte)
676  DO I = its, min(ide-1,ite)
677     thtmp   = grid%em_t_2(i,1,j)+t0
678     ptmp    = grid%em_p(i,1,j)+grid%em_pb(i,1,j)
679     temp(1) = thtmp * (ptmp/p1000mb)**rcp
680     thtmp   = grid%em_t_2(i,2,j)+t0
681     ptmp    = grid%em_p(i,2,j)+grid%em_pb(i,2,j)
682     temp(2) = thtmp * (ptmp/p1000mb)**rcp
683     thtmp   = grid%em_t_2(i,3,j)+t0
684     ptmp    = grid%em_p(i,3,j)+grid%em_pb(i,3,j)
685     temp(3) = thtmp * (ptmp/p1000mb)**rcp
686
687!!MARS
688!     grid%tsk(I,J)=grid%cf1*temp(1)+grid%cf2*temp(2)+grid%cf3*temp(3)
689     grid%tmn(I,J)=grid%tsk(I,J)-0.5
690!!!MARS
691!!TODO: passer la valeur a partir des donnees
692!grid%mars_tsoil(I,:,J)=grid%tsk(I,J)
693!!!MARS
694  ENDDO
695  ENDDO
696
697 END SUBROUTINE init_domain_rk
698
699   SUBROUTINE init_module_initialize
700   END SUBROUTINE init_module_initialize
701
702!---------------------------------------------------------------------
703
704!  test driver for get_sounding
705!
706!      implicit none
707!      integer n
708!      parameter(n = 1000)
709!      real zk(n),p(n),theta(n),rho(n),u(n),v(n),qv(n),pd(n)
710!      logical dry
711!      integer nl,k
712!
713!      dry = .false.
714!      dry = .true.
715!      call get_sounding( zk, p, pd, theta, rho, u, v, qv, dry, n, nl )
716!      write(6,*) ' input levels ',nl
717!      write(6,*) ' sounding '
718!      write(6,*) '  k  height(m)  press (Pa) pd(Pa) theta (K) den(kg/m^3)  u(m/s)     v(m/s)    qv(g/g) '
719!      do k=1,nl
720!        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)
721!      enddo
722!      end
723!
724!---------------------------------------------------------------------------
725
726      subroutine get_sounding( zk, p, p_dry, theta, rho, &
727                               u, v, qv, dry, nl_max, nl_in )
728      implicit none
729
730      integer nl_max, nl_in
731      real zk(nl_max), p(nl_max), theta(nl_max), rho(nl_max), &
732           u(nl_max), v(nl_max), qv(nl_max), p_dry(nl_max)
733      logical dry
734
735      integer n
736      parameter(n=1000)
737      logical debug
738      parameter( debug = .true.)
739
740! input sounding data
741
742      real p_surf, th_surf, qv_surf
743      real pi_surf, pi(n)
744      real h_input(n), th_input(n), qv_input(n), u_input(n), v_input(n)
745
746! diagnostics
747
748      real rho_surf, p_input(n), rho_input(n)
749      real pm_input(n)  !  this are for full moist sounding
750
751! local data
752
753      real p1000mb,cv,cp,r,cvpm,g
754!      parameter (p1000mb = 1.e+05, r = 287, cp = 1003., cv = cp-r, cvpm = -cv/cp, g=9.81 )
755      parameter (p1000mb = 610., r = 192., cp = 844.6, cv = cp-r, cvpm = -cv/cp, g=3.72)
756      integer k, it, nl
757      real qvf, qvf1, dz
758
759!  first, read the sounding
760
761      call read_sounding( p_surf, th_surf, qv_surf, &
762                          h_input, th_input, qv_input, u_input, v_input,n, nl, debug )
763
764      if(dry) then
765       do k=1,nl
766         qv_input(k) = 0.
767       enddo
768      endif
769
770      if(debug) write(6,*) ' number of input levels = ',nl
771
772        nl_in = nl
773        if(nl_in .gt. nl_max ) then
774          write(6,*) ' too many levels for input arrays ',nl_in,nl_max
775          call wrf_error_fatal ( ' too many levels for input arrays ' )
776        end if
777
778!  compute diagnostics,
779!  first, convert qv(g/kg) to qv(g/g)
780
781      do k=1,nl
782        qv_input(k) = 0.001*qv_input(k)
783      enddo
784
785      p_surf = 100.*p_surf  ! convert to pascals
786      qvf = 1. + rvovrd*qv_input(1)
787      rho_surf = 1./((r/p1000mb)*th_surf*qvf*((p_surf/p1000mb)**cvpm))
788      pi_surf = (p_surf/p1000mb)**(r/cp)
789
790      if(debug) then
791        write(6,*) ' surface density is ',rho_surf
792        write(6,*) ' surface pi is      ',pi_surf
793      end if
794
795
796!  integrate moist sounding hydrostatically, starting from the
797!  specified surface pressure
798!  -> first, integrate from surface to lowest level
799
800          qvf = 1. + rvovrd*qv_input(1)
801          qvf1 = 1. + qv_input(1)
802          rho_input(1) = rho_surf
803          dz = h_input(1)
804          do it=1,10
805!            pm_input(1) = p_surf &
806!                    - 0.5*dz*(rho_surf+rho_input(1))*g*qvf1
807!!!MARS MARS MARS
808            pm_input(1) = p_surf
809            rho_input(1) = 1./((r/p1000mb)*th_input(1)*qvf*((pm_input(1)/p1000mb)**cvpm))
810          enddo
811
812! integrate up the column
813
814          do k=2,nl
815            rho_input(k) = rho_input(k-1)
816            dz = h_input(k)-h_input(k-1)
817            qvf1 = 0.5*(2.+(qv_input(k-1)+qv_input(k)))
818            qvf = 1. + rvovrd*qv_input(k)   ! qv is in g/kg here
819 
820            do it=1,10
821              pm_input(k) = pm_input(k-1) &
822                      - 0.5*dz*(rho_input(k)+rho_input(k-1))*g*qvf1
823              rho_input(k) = 1./((r/p1000mb)*th_input(k)*qvf*((pm_input(k)/p1000mb)**cvpm))
824            enddo
825          enddo
826
827!  we have the moist sounding
828
829!  next, compute the dry sounding using p at the highest level from the
830!  moist sounding and integrating down.
831
832        p_input(nl) = pm_input(nl)
833
834          do k=nl-1,1,-1
835            dz = h_input(k+1)-h_input(k)
836            p_input(k) = p_input(k+1) + 0.5*dz*(rho_input(k)+rho_input(k+1))*g
837          enddo
838
839
840        do k=1,nl
841
842          zk(k) = h_input(k)
843          p(k) = pm_input(k)
844          p_dry(k) = p_input(k)
845          theta(k) = th_input(k)
846          rho(k) = rho_input(k)
847          u(k) = u_input(k)
848          v(k) = v_input(k)
849          qv(k) = qv_input(k)
850
851        enddo
852
853     if(debug) then
854      write(6,*) ' sounding '
855      write(6,*) '  k  height(m)  press (Pa) pd(Pa) theta (K) den(kg/m^3)  u(m/s)     v(m/s)    qv(g/g) '
856      do k=1,nl
857        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)
858      enddo
859
860     end if
861
862      end subroutine get_sounding
863
864!-------------------------------------------------------
865
866      subroutine read_sounding( ps,ts,qvs,h,th,qv,u,v,n,nl,debug )
867      implicit none
868      integer n,nl
869      real ps,ts,qvs,h(n),th(n),qv(n),u(n),v(n)
870      logical end_of_file
871      logical debug
872
873      integer k
874
875      open(unit=10,file='input_sounding',form='formatted',status='old')
876      rewind(10)
877      read(10,*) ps, ts, qvs
878      if(debug) then
879        write(6,*) ' input sounding surface parameters '
880        write(6,*) ' surface pressure (mb) ',ps
881        write(6,*) ' surface pot. temp (K) ',ts
882        write(6,*) ' surface mixing ratio (g/kg) ',qvs
883      end if
884
885      end_of_file = .false.
886      k = 0
887
888      do while (.not. end_of_file)
889
890        read(10,*,end=100) h(k+1), th(k+1), qv(k+1), u(k+1), v(k+1)
891        k = k+1
892        if(debug) write(6,'(1x,i3,5(1x,e10.3))') k, h(k), th(k), qv(k), u(k), v(k)
893        go to 110
894 100    end_of_file = .true.
895 110    continue
896      enddo
897
898      nl = k
899
900      close(unit=10,status = 'keep')
901
902      end subroutine read_sounding
903
904END MODULE module_initialize
Note: See TracBrowser for help on using the repository browser.