Specifying lat/lon on idealized Sea Breeze

Anything related to the idealized simulations in WRF.

Specifying lat/lon on idealized Sea Breeze

Postby ekimalab » Thu Jun 16, 2011 3:34 am

Hi!

I'm trying to run an idealized case in WRF where I can input different winds, temperature and moisture values over a specific location in the Philippines.

I'm fairly new with WRF and can't seem to work my way around this. The em_seabreeze2d_x idealized case doen't have terrain/topography input and I'm not sure if the real cases allow false sounding input.

Hope you can help me. Thank you so much!

Mike
ekimalab
 
Posts: 4
Joined: Wed Jun 15, 2011 10:34 pm

Re: Specifying lat/lon on idealized Sea Breeze

Postby Paul » Sun Jun 26, 2011 11:10 am

Hi Mike,

To set your domain's central lat./long. coordinate you need to edit these lines of code in the source code:
CALL nl_set_cen_lat(3,23.92)
CALL nl_set_cen_lon(2,7.50)

If you wish to include your own domain then you need to edit the array 'grid%ht(i,j)' within the source code by including your own function that defines your terrain.

I do not think a real WRF run will allow you to include a sounding file, like it does for the idealised cases. The procedure is more complex for the real cases and requires you to use the WPS system.

Hope this helps,
Paul
Paul
 
Posts: 1
Joined: Sun Jun 26, 2011 10:58 am

Re: Specifying lat/lon on idealized Sea Breeze

Postby andresperezcba » Wed Jun 27, 2012 6:03 pm

Hi Mike and Paul, The solution that I found is a little bit more complex.

nl_set_cen_lat(3,23.92) and nl_set_cen_lon(2,7.50) only set the lat an the long attribute, but the arrays xlat and xlong remain filled with 0 vaules. In order to fill the xlat and xlong variables, I used the WPS's maps_utils functions. At the end I attach the module_init code modified for the quarter_ss ideal case. Now ARWpost read the lat and long values and display them correctly with a mercator projector. I need to to do more test to check that everething is right and the values are correctly but this was a begining.

I modified em_ideal.F , module_initialize_quarter_ss.F of the code and postamble_new (in ./arch) and make (in main) to compile the code. The CONFIGFLAG (added by me) defined in config WRF (through postamble_new) is used to compile the ideal cases with the lat long specifications, if is undefined you can compile the modules like before. You need to do a ./clean -a and rebuild the code to run properly.
I hope that helps you. The !************************ something ************ comments type are my comments to the code.


em_ideal.F
Code: Select all

!IDEAL:DRIVER_LAYER
!
! create an initial data set for the WRF model based on an ideal condition

!*****************************************************************************************************************
!****Modification added in order to account the lat-lon and map projection set in the ideal init modules *********
!************************if they have the posibility to do that***************************************************
!****************Modified 24/06/2012   by A. Pérez Hortal - aperez1@famaf.unc.edu.ar******************************
!*****************************************************************************************************************

PROGRAM ideal

   USE module_domain , ONLY : domain
   USE module_initialize_ideal
   USE module_configure , ONLY : grid_config_rec_type

   USE module_timing
   USE module_wrf_error
#ifdef WRF_CHEM
   USE module_input_chem_data
   USE module_input_chem_bioemiss
#endif

   IMPLICIT NONE
#ifdef WRF_CHEM
  ! interface
   INTERFACE
     ! mediation-supplied
     SUBROUTINE med_read_wrf_chem_bioemiss ( grid , config_flags)
       USE module_domain
       TYPE (domain) grid
       TYPE (grid_config_rec_type) config_flags
     END SUBROUTINE med_read_wrf_chem_bioemiss
   END INTERFACE
#endif

   REAL    :: time

   INTEGER :: loop , &
              levels_to_process


   TYPE(domain) , POINTER :: keep_grid, grid_ptr, null_domain, grid
   TYPE(domain)           :: dummy
   TYPE (grid_config_rec_type)              :: config_flags
   TYPE (WRFU_Time) startTime, stopTime, currentTime
   TYPE (WRFU_TimeInterval) stepTime

   INTEGER :: max_dom , domain_id , fid , oid , idum1 , idum2 , ierr
   INTEGER :: debug_level, rc
   LOGICAL :: input_from_file

   INTERFACE
     SUBROUTINE med_initialdata_output ( grid , config_flags )
       USE module_domain , ONLY : domain
       USE module_configure , ONLY : grid_config_rec_type
       TYPE (domain) , POINTER :: grid
       TYPE (grid_config_rec_type) , INTENT(IN)   :: config_flags
     END SUBROUTINE med_initialdata_output
   END INTERFACE

#include "version_decl"


#ifdef DM_PARALLEL
   INTEGER                 :: nbytes
   INTEGER, PARAMETER      :: configbuflen = 4* CONFIG_BUF_LEN
   INTEGER                 :: configbuf( configbuflen )
   LOGICAL , EXTERNAL      :: wrf_dm_on_monitor
#endif

   CHARACTER (LEN=80)     :: message

   !  Define the name of this program (program_name defined in module_domain)

   program_name = "IDEAL " // TRIM(release_version) // " PREPROCESSOR"

   !  Get the NAMELIST data for input.

   CALL init_modules(1)   ! Phase 1 returns after MPI_INIT() (if it is called)
#ifdef NO_LEAP_CALENDAR
   CALL WRFU_Initialize( defaultCalKind=WRFU_CAL_NOLEAP, rc=rc )
#else
   CALL WRFU_Initialize( defaultCalKind=WRFU_CAL_GREGORIAN, rc=rc )
#endif
   CALL init_modules(2)   ! Phase 2 resumes after MPI_INIT() (if it is called)

#ifdef DM_PARALLEL
   IF ( wrf_dm_on_monitor() ) THEN
     CALL initial_config
   ENDIF
   CALL get_config_as_buffer( configbuf, configbuflen, nbytes )
   CALL wrf_dm_bcast_bytes( configbuf, nbytes )
   CALL set_config_as_buffer( configbuf, configbuflen )
   CALL wrf_dm_initialize
#else
   CALL initial_config
#endif
   CALL nl_get_debug_level ( 1, debug_level )
   CALL set_wrf_debug_level ( debug_level )

   CALL  wrf_message ( program_name )


   ! allocated and configure the mother domain

   NULLIFY( null_domain )

   CALL alloc_and_configure_domain ( domain_id  = 1 ,                  &
                                     grid       = head_grid ,          &
                                     parent     = null_domain ,        &
                                     kid        = -1                   )

   grid => head_grid
   ! TBH:  Note that historically, IDEAL did not set up clocks.  These
   ! TBH:  are explicit replacements for old default initializations...  They
   ! TBH:  are needed to ensure that time manager calls do not fail due to
   ! TBH:  uninitialized clock.  Clean this up later... 
   CALL WRFU_TimeSet(startTime, YY=1, MM=1, DD=1, H=0, M=0, S=0, rc=rc)
   stopTime = startTime
   currentTime = startTime
   ! TBH:  Bogus time step value -- clock is never advanced... 
   CALL WRFU_TimeIntervalSet(stepTime, S=180, rc=rc)
   grid%domain_clock = WRFU_ClockCreate( TimeStep= stepTime,  &
                                         StartTime=startTime, &
                                         StopTime= stopTime,  &
                                         rc=rc )
   CALL wrf_check_error( WRFU_SUCCESS, rc, &
                         'grid%domain_clock = WRFU_ClockCreate() FAILED', &
                         __FILE__ , &
                         __LINE__  )
   CALL       wrf_debug ( 100 , 'wrf: calling model_to_grid_config_rec ' )

   CALL model_to_grid_config_rec ( head_grid%id , model_config_rec , config_flags )
   CALL       wrf_debug ( 100 , 'wrf: calling set_scalar_indices_from_config ' )
   CALL set_scalar_indices_from_config ( head_grid%id , idum1, idum2 )


#ifdef PLANET
   WRITE ( current_date , FMT = '(I4.4,"-",I5.5,"_",I2.2,":",I2.2,":",I2.2,".0000")' ) &
           config_flags%start_year, &
           config_flags%start_day, &
           config_flags%start_hour, &
           config_flags%start_minute, &
           config_flags%start_second
#else
   WRITE ( current_date , FMT = '(I4.4,"-",I2.2,"-",I2.2,"_",I2.2,":",I2.2,":",I2.2,".0000")' ) &
           config_flags%start_year, &
           config_flags%start_month, &
           config_flags%start_day, &
           config_flags%start_hour, &
           config_flags%start_minute, &
           config_flags%start_second
#endif
   CALL domain_clockprint ( 150, grid, &
          'DEBUG assemble_output:  clock before 1st currTime set,' )
   WRITE (wrf_err_message,*) &
        'DEBUG assemble_output:  before 1st currTime set, current_date = ',TRIM(current_date)
   CALL wrf_debug ( 150 , wrf_err_message )
   CALL domain_clock_set( grid, current_timestr=current_date(1:19) )
   CALL domain_clockprint ( 150, grid, &
          'DEBUG assemble_output:  clock after 1st currTime set,' )

   CALL       wrf_debug ( 100 , 'wrf: calling init_wrfio' )
   CALL init_wrfio

#ifdef DM_PARALLEL
   CALL get_config_as_buffer( configbuf, configbuflen, nbytes )
   CALL wrf_dm_bcast_bytes( configbuf, nbytes )
   CALL set_config_as_buffer( configbuf, configbuflen )
#endif
     
#ifdef WRF_CHEM
         IF( grid%chem_opt > 0 ) then
           ! Read the chemistry data from a previous wrf forecast (wrfout file)
           IF(grid%chem_in_opt == 1 ) THEN
              message = 'INITIALIZING CHEMISTRY WITH OLD SIMULATION'
              CALL  wrf_message ( message )

              CALL med_read_wrf_chem_input ( grid , config_flags)
              IF(grid%emiss_opt == ECPTEC .or. grid%emiss_opt == GOCART_ECPTEC   &
                                         .or. grid%biomass_burn_opt == BIOMASSB) THEN
                 message = 'READING EMISSIONS DATA OPT 3'
                 CALL  wrf_message ( message )
!                CALL med_read_bin_chem_emissopt3 ( grid , config_flags)
                 CALL med_read_wrf_chem_emissopt3 ( grid , config_flags)
              END IF

              IF(grid%bio_emiss_opt == 2 ) THEN
                 message = 'READING BEIS3.11 EMISSIONS DATA'
                 CALL  wrf_message ( message )
                 CALL med_read_wrf_chem_bioemiss ( grid , config_flags)
              else IF(grid%bio_emiss_opt == 3 ) THEN !shc
                 message = 'READING MEGAN 2 EMISSIONS DATA'
                 CALL  wrf_message ( message )
                 CALL med_read_wrf_chem_bioemiss ( grid , config_flags)
              END IF

              IF(grid%dust_opt == 1 .or. grid%dmsemis_opt == 1 .or. grid%chem_opt == 300) THEN !shc
                 message = 'READING GOCART BG AND/OR DUST and DMS REF FIELDS'
                 CALL  wrf_message ( message )
                 CALL med_read_wrf_chem_gocart_bg ( grid , config_flags)
              END IF

           ELSEIF(grid%chem_in_opt == 0)then
              ! Generate chemistry data from a idealized vertical profile
              message = 'STARTING WITH BACKGROUND CHEMISTRY '
              CALL  wrf_message ( message )

              CALL input_chem_profile ( grid )

              IF(grid%bio_emiss_opt == 2 ) THEN
                 message = 'READING BEIS3.11 EMISSIONS DATA'
                 CALL  wrf_message ( message )
                 CALL med_read_wrf_chem_bioemiss ( grid , config_flags)
              else IF(grid%bio_emiss_opt == 3 ) THEN !shc
                 message = 'READING MEGAN 2 EMISSIONS DATA'
                 CALL  wrf_message ( message )
                 CALL med_read_wrf_chem_bioemiss ( grid , config_flags)
              END IF
              IF(grid%emiss_opt == ECPTEC .or. grid%emiss_opt == GOCART_ECPTEC   &
                                         .or. grid%biomass_burn_opt == BIOMASSB) THEN
                 message = 'READING EMISSIONS DATA OPT 3'
                 CALL  wrf_message ( message )
!                CALL med_read_bin_chem_emissopt3 ( grid , config_flags)
                 CALL med_read_wrf_chem_emissopt3 ( grid , config_flags)
              END IF

              IF(grid%dust_opt == 1 .or. grid%dmsemis_opt == 1 .or. grid%chem_opt == 300) THEN !shc
                 message = 'READING GOCART BG AND/OR DUST and DMS REF FIELDS'
                 CALL  wrf_message ( message )
                 CALL med_read_wrf_chem_gocart_bg ( grid , config_flags)
              END IF

           ELSE
             message = 'RUNNING WITHOUT CHEMISTRY INITIALIZATION'
             CALL  wrf_message ( message )
           END IF
         END IF
#endif

   CALL med_initialdata_output( head_grid , config_flags )

   CALL       wrf_debug (   0 , 'wrf: SUCCESS COMPLETE IDEAL INIT' )
   CALL med_shutdown_io ( head_grid , config_flags )
   CALL wrf_shutdown

   CALL WRFU_Finalize( rc=rc )

END PROGRAM ideal

SUBROUTINE med_initialdata_output ( grid , config_flags )
  ! Driver layer
   USE module_domain
   USE module_io_domain
   USE module_initialize_ideal
  ! Model layer
   USE module_configure

   IMPLICIT NONE

  ! Arguments
   TYPE(domain)  , POINTER                    :: grid
   TYPE (grid_config_rec_type) , INTENT(IN)   :: config_flags
  ! Local
   INTEGER                :: time_step_begin_restart
   INTEGER                :: fid , ierr , id
   CHARACTER (LEN=80)      :: rstname
   CHARACTER (LEN=80)      :: message
   CHARACTER (LEN=80)      :: inpname , bdyname

   !  Initialize the mother domain.

   grid%input_from_file = .false.

!*********************************************************************************************************************
!***if config flag present add the config_flag argument to init_domain************************************************
#ifdef CONFIGFLAG 
   CALL init_domain ( grid , config_flags)
#else
   CALL init_domain ( grid )
#endif
!**********************************************************************************************************************

   CALL calc_current_date ( grid%id, 0.)

   CALL construct_filename1 ( inpname , 'wrfinput' , grid%id , 2 )
   CALL open_w_dataset ( id, TRIM(inpname) , grid , config_flags , output_input , "DATASET=INPUT", ierr )
   IF ( ierr .NE. 0 ) THEN
     WRITE (wrf_err_message,*)'ideal: error opening wrfinput for writing ',ierr
     CALL wrf_error_fatal( wrf_err_message )
   ENDIF
   CALL output_input ( id, grid , config_flags , ierr )
   CALL close_dataset ( id , config_flags, "DATASET=INPUT" )


   IF ( config_flags%specified ) THEN
 
     CALL construct_filename1 ( bdyname , 'wrfbdy' , grid%id , 2 )
     CALL open_w_dataset ( id, TRIM(bdyname) , grid , config_flags , output_boundary , "DATASET=BOUNDARY", ierr )
     IF ( ierr .NE. 0 ) THEN
       WRITE (wrf_err_message,*)'ideal: error opening wrfbdy for writing ',ierr
       CALL wrf_error_fatal( wrf_err_message )
     ENDIF
     CALL output_boundary ( id, grid , config_flags , ierr )
     CALL close_dataset ( id , config_flags , "DATASET=BOUNDARY" )
 
   ENDIF
   


 
   RETURN
END SUBROUTINE med_initialdata_output





module_initialize_quarter_ss.F

Code: Select all
!IDEAL:MODEL_LAYER:INITIALIZATION
!

!  This MODULE holds the routines which are used to perform various initializations
!  for the individual domains. 

!  This MODULE CONTAINS the following routines:

!  initialize_field_test - 1. Set different fields to different constant
!                             values.  This is only a test.  If the correct
!                             domain is not found (based upon the "id")
!                             then a fatal error is issued.               
!-----------------------------------------------------------------------

!******************************************************************************************************
!*      Initialization with latitude and longitud specification                            ************
!*    This modified module set the lat-long and map projections setting in config_flags    ************
!*                         and is returned to em_ideal                                     ************
!*        It also set the xlong and xlat variable of grid datatype utilizing the           ************
!*   map_utils module from WPS. WPS need to be compiled with the same compiler as WRF      ************
!*                                                                                         ************
!*               Modified 24/06/2012   by A. Pérez Hortal - aperez1@famaf.unc.edu.ar       ************
!******************************************************************************************************

MODULE module_initialize_ideal
!******Add map_utils module to do the ij to latlong transforms*****************************************
   USE map_utils
!******************************************************************************************************
   USE module_domain
   USE module_io_domain
   USE module_state_description
   USE module_model_constants
   USE module_bc
   USE module_timing
   USE module_configure
   USE module_init_utilities
#ifdef DM_PARALLEL
   USE module_dm
#endif


CONTAINS


!-------------------------------------------------------------------
! this is a wrapper for the solver-specific init_domain routines.
! Also dereferences the grid variables and passes them down as arguments.
! This is crucial, since the lower level routines may do message passing
! and this will get fouled up on machines that insist on passing down
! copies of assumed-shape arrays (by passing down as arguments, the
! data are treated as assumed-size -- ie. f77 -- arrays and the copying
! business is avoided).  Fie on the F90 designers.  Fie and a pox.



!**********************************************************************************************************
!**The variable config_flags had been added in order to pass the lat-lon and map projection definitions*****
!**The type of this var is TYPE (grid_config_rec_type)  and its defined in module_configure****************
!**This var was also added to init_domain sobroutine*******************************************************
!**The init_domain call in ideal config_flag
!**********************************************************************************************************
   SUBROUTINE init_domain ( grid , config_flags)

   IMPLICIT NONE

   !  Input data.
   TYPE (domain), POINTER :: grid
   !Output data
   TYPE (grid_config_rec_type)    :: config_flags
   

   !  Local data.
   INTEGER :: idum1, idum2

   CALL set_scalar_indices_from_config ( head_grid%id , idum1, idum2 )

   CALL init_domain_rk( grid, config_flags &
!
#include <actual_new_args.inc>
!
                        )

   

   END SUBROUTINE init_domain
!-------------------------------------------------------------------

!**The variable config_flags had been added in order to pass the lat-lon and map projection definitions*****
   SUBROUTINE init_domain_rk ( grid , config_flags &
!
# include <dummy_new_args.inc>
!
)
   IMPLICIT NONE

   !  Input data.
   TYPE (domain), POINTER       :: grid
   TYPE (grid_config_rec_type)  :: config_flags

# include <dummy_new_decl.inc>



   !  Local data
   TYPE(proj_info)              :: proj

   INTEGER                             ::                       &
                                  ids, ide, jds, jde, kds, kde, &
                                  ims, ime, jms, jme, kms, kme, &
                                  its, ite, jts, jte, kts, kte, &
                                  i, j, k

   ! Local data

   INTEGER, PARAMETER :: nl_max = 1000
   REAL, DIMENSION(nl_max) :: zk, p_in, theta, rho, u, v, qv, pd_in
   INTEGER :: nl_in


   INTEGER :: icm,jcm, ii, im1, jj, jm1, loop, error, fid, nxc, nyc
   REAL    :: u_mean,v_mean, f0, p_surf, p_level, qvf, z_at_v, z_at_u
   REAL    :: z_scale, xrad, yrad, zrad, rad, delt, cof1, cof2
!   REAL, EXTERNAL :: interp_0
   REAL    :: hm
   REAL    :: pi

!  stuff from original initialization that has been dropped from the Registry
   REAL    :: vnu, xnu, xnus, dinit0, cbh, p0_temp, t0_temp, zd, zt
   REAL    :: qvf1, qvf2, pd_surf
   INTEGER :: it
   real :: thtmp, ptmp, temp(3)

   LOGICAL :: moisture_init
   LOGICAL :: stretch_grid, dry_sounding

  INTEGER :: xs , xe , ys , ye
  REAL :: mtn_ht
   LOGICAL, EXTERNAL :: wrf_dm_on_monitor

   SELECT CASE ( model_data_order )
         CASE ( DATA_ORDER_ZXY )
   kds = grid%sd31 ; kde = grid%ed31 ;
   ids = grid%sd32 ; ide = grid%ed32 ;
   jds = grid%sd33 ; jde = grid%ed33 ;

   kms = grid%sm31 ; kme = grid%em31 ;
   ims = grid%sm32 ; ime = grid%em32 ;
   jms = grid%sm33 ; jme = grid%em33 ;

   kts = grid%sp31 ; kte = grid%ep31 ;   ! note that tile is entire patch
   its = grid%sp32 ; ite = grid%ep32 ;   ! note that tile is entire patch
   jts = grid%sp33 ; jte = grid%ep33 ;   ! note that tile is entire patch
         CASE ( DATA_ORDER_XYZ )
   ids = grid%sd31 ; ide = grid%ed31 ;
   jds = grid%sd32 ; jde = grid%ed32 ;
   kds = grid%sd33 ; kde = grid%ed33 ;

   ims = grid%sm31 ; ime = grid%em31 ;
   jms = grid%sm32 ; jme = grid%em32 ;
   kms = grid%sm33 ; kme = grid%em33 ;

   its = grid%sp31 ; ite = grid%ep31 ;   ! note that tile is entire patch
   jts = grid%sp32 ; jte = grid%ep32 ;   ! note that tile is entire patch
   kts = grid%sp33 ; kte = grid%ep33 ;   ! note that tile is entire patch
         CASE ( DATA_ORDER_XZY )
   ids = grid%sd31 ; ide = grid%ed31 ;
   kds = grid%sd32 ; kde = grid%ed32 ;
   jds = grid%sd33 ; jde = grid%ed33 ;

   ims = grid%sm31 ; ime = grid%em31 ;
   kms = grid%sm32 ; kme = grid%em32 ;
   jms = grid%sm33 ; jme = grid%em33 ;

   its = grid%sp31 ; ite = grid%ep31 ;   ! note that tile is entire patch
   kts = grid%sp32 ; kte = grid%ep32 ;   ! note that tile is entire patch
   jts = grid%sp33 ; jte = grid%ep33 ;   ! note that tile is entire patch

   END SELECT


   stretch_grid = .true.
   delt = 3.
!   z_scale = .50
   z_scale = .40
   pi = 2.*asin(1.0)
   write(6,*) ' pi is ',pi
   nxc = (ide-ids)/2
   nyc = (jde-jds)/2


!*********Set of lat long and map projection info*******************
    CALL nl_set_mminlu(grid%id, '    ')
    CALL nl_set_iswater(grid%id,0)
    CALL nl_set_cen_lat(grid%id,-34.)
    CALL nl_set_cen_lon(grid%id,-68.)
    CALL nl_set_truelat1(grid%id,-34.)
    CALL nl_set_stand_lon (grid%id,-60.)
    CALL nl_set_truelat2(grid%id,0.)
    CALL nl_set_moad_cen_lat (grid%id,-34.)
    CALL nl_set_pole_lon (grid%id,0.)
    CALL nl_set_pole_lat (grid%id,90.)
    CALL nl_set_map_proj(grid%id,3)

!****        Save the lat, long and map projections info    in model_config_rec  ***********************
   CALL model_to_grid_config_rec ( grid%id , model_config_rec , config_flags )
   
!****  Populate the proj data type if undefined values are present   **********************************
   CALL map_set( proj_code=config_flags%map_proj, proj=proj, lat1=config_flags%cen_lat, truelat1=config_flags%truelat1, &
                 lon1=config_flags%cen_lon, knowni=real(nxc), knownj=real(nyc), dx=grid%dx, dy=grid%dy)
!****  Set the xlat and xlong values for all grid points   *********************************************
   DO j = jts, jte
    DO i = its, ite
      CALL ij_to_latlon(proj,real(i) , real(j), grid%xlat(i,j), grid%xlong(i,j))
    END DO
  END DO



   !here we check to see if the boundary conditions are set properly
   CALL boundary_condition_check( config_flags, bdyzone, error, grid%id )


!*********************************************************************************************************


   moisture_init = .true.

   grid%itimestep=0

#ifdef DM_PARALLEL
   CALL wrf_dm_bcast_bytes( icm , IWORDSIZE )
   CALL wrf_dm_bcast_bytes( jcm , IWORDSIZE )
#endif


   
!  here we initialize data we currently is not initialized
!  in the input data

    DO j = jts, jte
      DO i = its, ite
         grid%msftx(i,j)    = 1.
         grid%msfty(i,j)    = 1.
         grid%msfux(i,j)    = 1.
         grid%msfuy(i,j)    = 1.
         grid%msfvx(i,j)    = 1.
         grid%msfvx_inv(i,j)= 1.
         grid%msfvy(i,j)    = 1.
         grid%sina(i,j)     = 0.
         grid%cosa(i,j)     = 1.
         grid%e(i,j)        = 0.
         grid%f(i,j)        = 0.
             
      END DO
   END DO

    DO j = jts, jte
    DO k = kts, kte
      DO i = its, ite
         grid%ww(i,k,j)     = 0.
      END DO
   END DO
   END DO

   grid%step_number = 0

! set up the grid

   IF (stretch_grid) THEN ! exponential stretch for eta (nearly constant dz)
     DO k=1, kde
      grid%znw(k) = (exp(-(k-1)/float(kde-1)/z_scale) - exp(-1./z_scale))/ &
                                (1.-exp(-1./z_scale))
     ENDDO
   ELSE
     DO k=1, kde
      grid%znw(k) = 1. - float(k-1)/float(kde-1)
     ENDDO
   ENDIF

   DO k=1, kde-1
    grid%dnw(k) = grid%znw(k+1) - grid%znw(k)
    grid%rdnw(k) = 1./grid%dnw(k)
    grid%znu(k) = 0.5*(grid%znw(k+1)+grid%znw(k))
   ENDDO
   DO k=2, kde-1
    grid%dn(k) = 0.5*(grid%dnw(k)+grid%dnw(k-1))
    grid%rdn(k) = 1./grid%dn(k)
    grid%fnp(k) = .5* grid%dnw(k  )/grid%dn(k)
    grid%fnm(k) = .5* grid%dnw(k-1)/grid%dn(k)
   ENDDO

   cof1 = (2.*grid%dn(2)+grid%dn(3))/(grid%dn(2)+grid%dn(3))*grid%dnw(1)/grid%dn(2)
   cof2 =     grid%dn(2)        /(grid%dn(2)+grid%dn(3))*grid%dnw(1)/grid%dn(3)
   grid%cf1  = grid%fnp(2) + cof1
   grid%cf2  = grid%fnm(2) - cof1 - cof2
   grid%cf3  = cof2       

   grid%cfn  = (.5*grid%dnw(kde-1)+grid%dn(kde-1))/grid%dn(kde-1)
   grid%cfn1 = -.5*grid%dnw(kde-1)/grid%dn(kde-1)
   grid%rdx = 1./config_flags%dx
   grid%rdy = 1./config_flags%dy

!  get the sounding from the ascii sounding file, first get dry sounding and
!  calculate base state

  dry_sounding = .true.
  IF ( wrf_dm_on_monitor() ) THEN
  write(6,*) ' getting dry sounding for base state '

  CALL get_sounding( zk, p_in, pd_in, theta, rho, u, v, qv, dry_sounding, nl_max, nl_in )
  ENDIF
  CALL wrf_dm_bcast_real( zk , nl_max )
  CALL wrf_dm_bcast_real( p_in , nl_max )
  CALL wrf_dm_bcast_real( pd_in , nl_max )
  CALL wrf_dm_bcast_real( theta , nl_max )
  CALL wrf_dm_bcast_real( rho , nl_max )
  CALL wrf_dm_bcast_real( u , nl_max )
  CALL wrf_dm_bcast_real( v , nl_max )
  CALL wrf_dm_bcast_real( qv , nl_max )
  CALL wrf_dm_bcast_integer ( nl_in , 1 )

  write(6,*) ' returned from reading sounding, nl_in is ',nl_in

!  find ptop for the desired ztop (ztop is input from the namelist),
!  and find surface pressure

  grid%p_top = interp_0( p_in, zk, config_flags%ztop, nl_in )

  DO j=jts,jte
  DO i=its,ite
    grid%ht(i,j) = 0.
  ENDDO
  ENDDO

  xs=ide/2 -3
  xs=ids   -3
  xe=xs + 6
  ys=jde/2 -3
  ye=ys + 6
  mtn_ht = 500
#ifdef MTN
  DO j=max(ys,jds),min(ye,jde-1)
  DO i=max(xs,ids),min(xe,ide-1)
     grid%ht(i,j) = mtn_ht * 0.25 * &
               ( 1. + COS ( 2*pi/(xe-xs) * ( i-xs ) + pi ) ) * &
               ( 1. + COS ( 2*pi/(ye-ys) * ( j-ys ) + pi ) )
  ENDDO
  ENDDO
#endif
#ifdef EW_RIDGE
  DO j=max(ys,jds),min(ye,jde-1)
  DO i=ids,ide
     grid%ht(i,j) = mtn_ht * 0.50 * &
               ( 1. + COS ( 2*pi/(ye-ys) * ( j-ys ) + pi ) )
  ENDDO
  ENDDO
#endif
#ifdef NS_RIDGE
  DO j=jds,jde
  DO i=max(xs,ids),min(xe,ide-1)
     grid%ht(i,j) = mtn_ht * 0.50 * &
               ( 1. + COS ( 2*pi/(xe-xs) * ( i-xs ) + pi ) )
  ENDDO
  ENDDO
#endif
  DO j=jts,jte
  DO i=its,ite
    grid%phb(i,1,j) = g * grid%ht(i,j)
    grid%ph0(i,1,j) = g * grid%ht(i,j)
  ENDDO
  ENDDO

  DO J = jts, jte
  DO I = its, ite

    p_surf = interp_0( p_in, zk, grid%phb(i,1,j)/g, nl_in )
    grid%mub(i,j) = p_surf-grid%p_top

!  this is dry hydrostatic sounding (base state), so given grid%p (coordinate),
!  interp theta (from interp) and compute 1/rho from eqn. of state

    DO K = 1, kte-1
      p_level = grid%znu(k)*(p_surf - grid%p_top) + grid%p_top
      grid%pb(i,k,j) = p_level
      grid%t_init(i,k,j) = interp_0( theta, p_in, p_level, nl_in ) - t0
      grid%alb(i,k,j) = (r_d/p1000mb)*(grid%t_init(i,k,j)+t0)*(grid%pb(i,k,j)/p1000mb)**cvpm
    ENDDO

!  calc hydrostatic balance (alternatively we could interp the geopotential from the
!  sounding, but this assures that the base state is in exact hydrostatic balance with
!  respect to the model eqns.

    DO k  = 2,kte
      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)
    ENDDO

  ENDDO
  ENDDO

  IF ( wrf_dm_on_monitor() ) THEN
    write(6,*) ' ptop is ',grid%p_top
    write(6,*) ' base state grid%mub(1,1), p_surf is ',grid%mub(1,1),grid%mub(1,1)+grid%p_top
  ENDIF

!  calculate full state for each column - this includes moisture.

  write(6,*) ' getting moist sounding for full state '
  dry_sounding = .false.
  CALL get_sounding( zk, p_in, pd_in, theta, rho, u, v, qv, dry_sounding, nl_max, nl_in )

  DO J = jts, min(jde-1,jte)
  DO I = its, min(ide-1,ite)

!  At this point grid%p_top is already set. find the DRY mass in the column
!  by interpolating the DRY pressure. 

   pd_surf = interp_0( pd_in, zk, grid%phb(i,1,j)/g, nl_in )

!  compute the perturbation mass and the full mass

    grid%mu_1(i,j) = pd_surf-grid%p_top - grid%mub(i,j)
    grid%mu_2(i,j) = grid%mu_1(i,j)
    grid%mu0(i,j) = grid%mu_1(i,j) + grid%mub(i,j)

! given the dry pressure and coordinate system, interp the potential
! temperature and qv

    do k=1,kde-1

      p_level = grid%znu(k)*(pd_surf - grid%p_top) + grid%p_top

      moist(i,k,j,P_QV) = interp_0( qv, pd_in, p_level, nl_in )
      grid%t_1(i,k,j)          = interp_0( theta, pd_in, p_level, nl_in ) - t0
      grid%t_2(i,k,j)          = grid%t_1(i,k,j)
     

    enddo

!  integrate the hydrostatic equation (from the RHS of the bigstep
!  vertical momentum equation) down from the top to get grid%p.
!  first from the top of the model to the top pressure

    k = kte-1  ! top level

    qvf1 = 0.5*(moist(i,k,j,P_QV)+moist(i,k,j,P_QV))
    qvf2 = 1./(1.+qvf1)
    qvf1 = qvf1*qvf2

!    grid%p(i,k,j) = - 0.5*grid%mu_1(i,j)/grid%rdnw(k)
    grid%p(i,k,j) = - 0.5*(grid%mu_1(i,j)+qvf1*grid%mub(i,j))/grid%rdnw(k)/qvf2
    qvf = 1. + rvovrd*moist(i,k,j,P_QV)
    grid%alt(i,k,j) = (r_d/p1000mb)*(grid%t_1(i,k,j)+t0)*qvf* &
                (((grid%p(i,k,j)+grid%pb(i,k,j))/p1000mb)**cvpm)
    grid%al(i,k,j) = grid%alt(i,k,j) - grid%alb(i,k,j)

!  down the column

    do k=kte-2,1,-1
      qvf1 = 0.5*(moist(i,k,j,P_QV)+moist(i,k+1,j,P_QV))
      qvf2 = 1./(1.+qvf1)
      qvf1 = qvf1*qvf2
      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)
      qvf = 1. + rvovrd*moist(i,k,j,P_QV)
      grid%alt(i,k,j) = (r_d/p1000mb)*(grid%t_1(i,k,j)+t0)*qvf* &
                  (((grid%p(i,k,j)+grid%pb(i,k,j))/p1000mb)**cvpm)
      grid%al(i,k,j) = grid%alt(i,k,j) - grid%alb(i,k,j)
    enddo

!  this is the hydrostatic equation used in the model after the
!  small timesteps.  In the model, grid%al (inverse density)
!  is computed from the geopotential.


    grid%ph_1(i,1,j) = 0.
    DO k  = 2,kte
      grid%ph_1(i,k,j) = grid%ph_1(i,k-1,j) - (1./grid%rdnw(k-1))*(       &
                   (grid%mub(i,j)+grid%mu_1(i,j))*grid%al(i,k-1,j)+ &
                    grid%mu_1(i,j)*grid%alb(i,k-1,j)  )
                                                   
      grid%ph_2(i,k,j) = grid%ph_1(i,k,j)
      grid%ph0(i,k,j) = grid%ph_1(i,k,j) + grid%phb(i,k,j)
    ENDDO

    IF ( wrf_dm_on_monitor() ) THEN
    if((i==2) .and. (j==2)) then
     write(6,*) ' grid%ph_1 calc ',grid%ph_1(2,1,2),grid%ph_1(2,2,2),&
                              grid%mu_1(2,2)+grid%mub(2,2),grid%mu_1(2,2), &
                              grid%alb(2,1,2),grid%al(1,2,1),grid%rdnw(1)
    endif
    ENDIF

  ENDDO
  ENDDO

!#if 0

!  thermal perturbation to kick off convection

  write(6,*) ' nxc, nyc for perturbation ',nxc,nyc
  write(6,*) ' delt for perturbation ',delt

  DO J = jts, min(jde-1,jte)
    yrad = config_flags%dy*float(j-nyc)/10000.
!   yrad = 0.
    DO I = its, min(ide-1,ite)
      xrad = config_flags%dx*float(i-nxc)/10000.
!     xrad = 0.
      DO K = 1, kte-1

!  put in preturbation theta (bubble) and recalc density.  note,
!  the mass in the column is not changing, so when theta changes,
!  we recompute density and geopotential

        zrad = 0.5*(grid%ph_1(i,k,j)+grid%ph_1(i,k+1,j)  &
                   +grid%phb(i,k,j)+grid%phb(i,k+1,j))/g
        zrad = (zrad-1500.)/1500.
        RAD=SQRT(xrad*xrad+yrad*yrad+zrad*zrad)
        IF(RAD <= 1.) THEN
           grid%t_1(i,k,j)=grid%t_1(i,k,j)+delt*COS(.5*PI*RAD)**2
           grid%t_2(i,k,j)=grid%t_1(i,k,j)
           qvf = 1. + rvovrd*moist(i,k,j,P_QV)
           grid%alt(i,k,j) = (r_d/p1000mb)*(grid%t_1(i,k,j)+t0)*qvf* &
                        (((grid%p(i,k,j)+grid%pb(i,k,j))/p1000mb)**cvpm)
           grid%al(i,k,j) = grid%alt(i,k,j) - grid%alb(i,k,j)
        ENDIF
      ENDDO

!  rebalance hydrostatically

      DO k  = 2,kte
        grid%ph_1(i,k,j) = grid%ph_1(i,k-1,j) - (1./grid%rdnw(k-1))*(       &
                     (grid%mub(i,j)+grid%mu_1(i,j))*grid%al(i,k-1,j)+ &
                      grid%mu_1(i,j)*grid%alb(i,k-1,j)  )
                                                   
        grid%ph_2(i,k,j) = grid%ph_1(i,k,j)
        grid%ph0(i,k,j) = grid%ph_1(i,k,j) + grid%phb(i,k,j)
      ENDDO

    ENDDO
  ENDDO

!#endif

   IF ( wrf_dm_on_monitor() ) THEN
   write(6,*) ' grid%mu_1 from comp ', grid%mu_1(1,1)
   write(6,*) ' full state sounding from comp, ph, grid%p, grid%al, grid%t_1, qv '
   do k=1,kde-1
     write(6,'(i3,1x,5(1x,1pe10.3))') k, grid%ph_1(1,k,1)+grid%phb(1,k,1), &
                                      grid%p(1,k,1)+grid%pb(1,k,1), grid%alt(1,k,1), &
                                      grid%t_1(1,k,1)+t0, moist(1,k,1,P_QV)
   enddo

   write(6,*) ' pert state sounding from comp, grid%ph_1, pp, alp, grid%t_1, qv '
   do k=1,kde-1
     write(6,'(i3,1x,5(1x,1pe10.3))') k, grid%ph_1(1,k,1), &
                                      grid%p(1,k,1), grid%al(1,k,1), &
                                      grid%t_1(1,k,1), moist(1,k,1,P_QV)
   enddo
   ENDIF

! interp v

  DO J = jts, jte
  DO I = its, min(ide-1,ite)

    IF (j == jds) THEN
      z_at_v = grid%phb(i,1,j)/g
    ELSE IF (j == jde) THEN
      z_at_v = grid%phb(i,1,j-1)/g
    ELSE
      z_at_v = 0.5*(grid%phb(i,1,j)+grid%phb(i,1,j-1))/g
    END IF
    p_surf = interp_0( p_in, zk, z_at_v, nl_in )

    DO K = 1, kte-1
      p_level = grid%znu(k)*(p_surf - grid%p_top) + grid%p_top
      grid%v_1(i,k,j) = interp_0( v, p_in, p_level, nl_in )
      grid%v_2(i,k,j) = grid%v_1(i,k,j)
    ENDDO

  ENDDO
  ENDDO

! interp u

  DO J = jts, min(jde-1,jte)
  DO I = its, ite

    IF (i == ids) THEN
      z_at_u = grid%phb(i,1,j)/g
    ELSE IF (i == ide) THEN
      z_at_u = grid%phb(i-1,1,j)/g
    ELSE
      z_at_u = 0.5*(grid%phb(i,1,j)+grid%phb(i-1,1,j))/g
    END IF

    p_surf = interp_0( p_in, zk, z_at_u, nl_in )

    DO K = 1, kte-1
      p_level = grid%znu(k)*(p_surf - grid%p_top) + grid%p_top
      grid%u_1(i,k,j) = interp_0( u, p_in, p_level, nl_in )
      grid%u_2(i,k,j) = grid%u_1(i,k,j)
    ENDDO

  ENDDO
  ENDDO

!  set w

  DO J = jts, min(jde-1,jte)
  DO K = kts, kte
  DO I = its, min(ide-1,ite)
    grid%w_1(i,k,j) = 0.
    grid%w_2(i,k,j) = 0.
  ENDDO
  ENDDO
  ENDDO

!  set a few more things

  DO J = jts, min(jde-1,jte)
  DO K = kts, kte-1
  DO I = its, min(ide-1,ite)
    grid%h_diabatic(i,k,j) = 0.
    !scalar(j,k,i, P_XAGI) = 0.0;
  ENDDO
  ENDDO
  ENDDO

  DO K = kts, kte-1
     !scalar( int( (min(ide-1,ite)-its)/2), K , int( (min(jde-1,jte)-jts)/2) , P_XAGI) = 10.;
  END DO


  IF ( wrf_dm_on_monitor() ) THEN
  DO k=1,kte-1
    grid%t_base(k) = grid%t_1(1,k,1)
    grid%qv_base(k) = moist(1,k,1,P_QV)
    grid%u_base(k) = grid%u_1(1,k,1)
    grid%v_base(k) = grid%v_1(1,k,1)
    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
  ENDDO
  ENDIF
  CALL wrf_dm_bcast_real( grid%t_base , kte )
  CALL wrf_dm_bcast_real( grid%qv_base , kte )
  CALL wrf_dm_bcast_real( grid%u_base , kte )
  CALL wrf_dm_bcast_real( grid%v_base , kte )
  CALL wrf_dm_bcast_real( grid%z_base , kte )

  DO J = jts, min(jde-1,jte)
  DO I = its, min(ide-1,ite)
     thtmp   = grid%t_2(i,1,j)+t0
     ptmp    = grid%p(i,1,j)+grid%pb(i,1,j)
     temp(1) = thtmp * (ptmp/p1000mb)**rcp
     thtmp   = grid%t_2(i,2,j)+t0
     ptmp    = grid%p(i,2,j)+grid%pb(i,2,j)
     temp(2) = thtmp * (ptmp/p1000mb)**rcp
     thtmp   = grid%t_2(i,3,j)+t0
     ptmp    = grid%p(i,3,j)+grid%pb(i,3,j)
     temp(3) = thtmp * (ptmp/p1000mb)**rcp

     grid%tsk(I,J)=grid%cf1*temp(1)+grid%cf2*temp(2)+grid%cf3*temp(3)
     grid%tmn(I,J)=grid%tsk(I,J)-0.5
  ENDDO
  ENDDO

     !scalar(  jts, kts , its  , P_XAGI ) = 10.0;
     !scalar(  min(jde-1,jte),  kte-1, min(ide-1,ite)  , P_XAGI ) = 10.0;
 

 END SUBROUTINE init_domain_rk

   SUBROUTINE init_module_initialize
   END SUBROUTINE init_module_initialize

!---------------------------------------------------------------------

!  test driver for get_sounding
!
!      implicit none
!      integer n
!      parameter(n = 1000)
!      real zk(n),p(n),theta(n),rho(n),u(n),v(n),qv(n),pd(n)
!      logical dry
!      integer nl,k
!
!      dry = .false.
!      dry = .true.
!      call get_sounding( zk, p, pd, theta, rho, u, v, qv, dry, n, nl )
!      write(6,*) ' input levels ',nl
!      write(6,*) ' sounding '
!      write(6,*) '  k  height(m)  press (Pa) pd(Pa) theta (K) den(kg/m^3)  u(m/s)     v(m/s)    qv(g/g) '
!      do k=1,nl
!        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)
!      enddo
!      end
!
!---------------------------------------------------------------------------

      subroutine get_sounding( zk, p, p_dry, theta, rho, &
                               u, v, qv, dry, nl_max, nl_in )
      implicit none

      integer nl_max, nl_in
      real zk(nl_max), p(nl_max), theta(nl_max), rho(nl_max), &
           u(nl_max), v(nl_max), qv(nl_max), p_dry(nl_max)
      logical dry

      integer n
      parameter(n=1000)
      logical debug
      parameter( debug = .true.)

! input sounding data

      real p_surf, th_surf, qv_surf
      real pi_surf, pi(n)
      real h_input(n), th_input(n), qv_input(n), u_input(n), v_input(n)

! diagnostics

      real rho_surf, p_input(n), rho_input(n)
      real pm_input(n)  !  this are for full moist sounding

! local data

      real r
      parameter (r = r_d)
      integer k, it, nl
      real qvf, qvf1, dz

!  first, read the sounding

      call read_sounding( p_surf, th_surf, qv_surf, &
                          h_input, th_input, qv_input, u_input, v_input,n, nl, debug )

      if(dry) then
       do k=1,nl
         qv_input(k) = 0.
       enddo
      endif

      if(debug) write(6,*) ' number of input levels = ',nl

        nl_in = nl
        if(nl_in .gt. nl_max ) then
          write(6,*) ' too many levels for input arrays ',nl_in,nl_max
          call wrf_error_fatal ( ' too many levels for input arrays ' )
        end if

!  compute diagnostics,
!  first, convert qv(g/kg) to qv(g/g)

      do k=1,nl
        qv_input(k) = 0.001*qv_input(k)
      enddo

      p_surf = 100.*p_surf  ! convert to pascals
      qvf = 1. + rvovrd*qv_input(1)
      rho_surf = 1./((r/p1000mb)*th_surf*qvf*((p_surf/p1000mb)**cvpm))
      pi_surf = (p_surf/p1000mb)**(r/cp)

      if(debug) then
        write(6,*) ' surface density is ',rho_surf
        write(6,*) ' surface pi is      ',pi_surf
      end if


!  integrate moist sounding hydrostatically, starting from the
!  specified surface pressure
!  -> first, integrate from surface to lowest level

          qvf = 1. + rvovrd*qv_input(1)
          qvf1 = 1. + qv_input(1)
          rho_input(1) = rho_surf
          dz = h_input(1)
          do it=1,10
            pm_input(1) = p_surf &
                    - 0.5*dz*(rho_surf+rho_input(1))*g*qvf1
            rho_input(1) = 1./((r/p1000mb)*th_input(1)*qvf*((pm_input(1)/p1000mb)**cvpm))
          enddo

! integrate up the column

          do k=2,nl
            rho_input(k) = rho_input(k-1)
            dz = h_input(k)-h_input(k-1)
            qvf1 = 0.5*(2.+(qv_input(k-1)+qv_input(k)))
            qvf = 1. + rvovrd*qv_input(k)   ! qv is in g/kg here
 
            do it=1,10
              pm_input(k) = pm_input(k-1) &
                      - 0.5*dz*(rho_input(k)+rho_input(k-1))*g*qvf1
              rho_input(k) = 1./((r/p1000mb)*th_input(k)*qvf*((pm_input(k)/p1000mb)**cvpm))
            enddo
          enddo

!  we have the moist sounding

!  next, compute the dry sounding using p at the highest level from the
!  moist sounding and integrating down.

        p_input(nl) = pm_input(nl)

          do k=nl-1,1,-1
            dz = h_input(k+1)-h_input(k)
            p_input(k) = p_input(k+1) + 0.5*dz*(rho_input(k)+rho_input(k+1))*g
          enddo


        do k=1,nl

          zk(k) = h_input(k)
          p(k) = pm_input(k)
          p_dry(k) = p_input(k)
          theta(k) = th_input(k)
          rho(k) = rho_input(k)
          u(k) = u_input(k)
          v(k) = v_input(k)
          qv(k) = qv_input(k)

        enddo

     if(debug) then
      write(6,*) ' sounding '
      write(6,*) '  k  height(m)  press (Pa) pd(Pa) theta (K) den(kg/m^3)  u(m/s)     v(m/s)    qv(g/g) '
      do k=1,nl
        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)
      enddo

     end if

      end subroutine get_sounding

!-------------------------------------------------------

      subroutine read_sounding( ps,ts,qvs,h,th,qv,u,v,n,nl,debug )
      implicit none
      integer n,nl
      real ps,ts,qvs,h(n),th(n),qv(n),u(n),v(n)
      logical end_of_file
      logical debug

      integer k

      open(unit=10,file='input_sounding',form='formatted',status='old')
      rewind(10)
      read(10,*) ps, ts, qvs
      if(debug) then
        write(6,*) ' input sounding surface parameters '
        write(6,*) ' surface pressure (mb) ',ps
        write(6,*) ' surface pot. temp (K) ',ts
        write(6,*) ' surface mixing ratio (g/kg) ',qvs
      end if

      end_of_file = .false.
      k = 0

      do while (.not. end_of_file)

        read(10,*,end=100) h(k+1), th(k+1), qv(k+1), u(k+1), v(k+1)
        k = k+1
        if(debug) write(6,'(1x,i3,5(1x,e10.3))') k, h(k), th(k), qv(k), u(k), v(k)
        go to 110
 100    end_of_file = .true.
 110    continue
      enddo

      nl = k

      close(unit=10,status = 'keep')

      end subroutine read_sounding

END MODULE module_initialize_ideal



makefile in ./main dir

Code: Select all
#Makefile modified to support the lat-lon and map projections specifications on ideal simulaiton.
#In this case only the quarter_ss init module was modified
#
#Modified 24/06/2012   by A. Pérez Hortal - aperez1@famaf.unc.edu.ar


LN      =       ln -sf
MAKE    =       make -i -r
RM      =       rm -f

MODULES =  module_wrf_top.F


#map_utils modules and dependencies - .o files from WPS compilation
# The WPS_DIR enviromental variable needs to be set as the root dir of the WPS compilation,
#For example WPS_DIR can be set in the $(HOME)/.profile file, and is loaded  when the session starts.
#WPS needs to be compiled with the same compiler used in WRF
MAP_UTILS = $(WPS_DIR)/util/src/module_map_utils.o $(WPS_DIR)/util/src/constants_module.o $(WPS_DIR)/util/src/misc_definitions_module.o $(WPS_DIR)/util/src/module_debug.o $(WPS_DIR)/util/src/cio.o
#*************************************************************************

OBJS    =

LIBPATHS =

include ../configure.wrf

$(SOLVER)_wrf : wrf.o ../main/module_wrf_top.o
   $(RANLIB) $(LIBWRFLIB)
   $(LD) -o wrf.exe $(LDFLAGS) wrf.o ../main/module_wrf_top.o $(LIBWRFLIB) $(LIB)

$(SOLVER)_wrf_SST_ESMF : wrf_ESMFMod.o wrf_SST_ESMF.o ../main/module_wrf_top.o
   $(RANLIB) $(LIBWRFLIB)
   $(LD) -o wrf_SST_ESMF.exe $(LDFLAGS) wrf_SST_ESMF.o wrf_ESMFMod.o ../main/module_wrf_top.o $(LIBWRFLIB) $(LIB)

$(SOLVER)_ideal : module_initialize ideal_$(SOLVER).o
   $(RANLIB) $(LIBWRFLIB)

#*** Conditional compilation if the module quarter_ss is used. map_utils .o files, dependencies and CONFIG_FLAG definition added to the compilation *******
ifeq ($(IDEAL_CASE), quarter_ss)
   $(LD) -o ideal.exe $(LDFLAGS) -DCONFIGFLAG ideal_$(SOLVER).o ../dyn_$(SOLVER)/module_initialize_$(IDEAL_CASE).o $(MAP_UTILS) $(LIBWRFLIB) $(LIB)
else
#Any other case
   $(LD) -o ideal.exe $(LDFLAGS) ideal_$(SOLVER).o ../dyn_$(SOLVER)/module_initialize_$(IDEAL_CASE).o $(LIBWRFLIB) $(LIB)
endif
#**********************************************************************************************************************************************************

$(SOLVER)_real : module_initialize real_$(SOLVER).o ndown_$(SOLVER).o nup_$(SOLVER).o tc_$(SOLVER).o
   $(RANLIB) $(LIBWRFLIB)
   $(LD) -o real.exe $(LDFLAGS) real_$(SOLVER).o ../dyn_$(SOLVER)/module_initialize_$(IDEAL_CASE).o $(LIBWRFLIB) $(LIB)
   $(LD) -o ndown.exe $(LDFLAGS) ndown_$(SOLVER).o  ../dyn_$(SOLVER)/module_initialize_$(IDEAL_CASE).o $(LIBWRFLIB) $(LIB)
   $(LD) -o nup.exe $(LDFLAGS) nup_$(SOLVER).o  ../dyn_$(SOLVER)/module_initialize_$(IDEAL_CASE).o $(LIBWRFLIB) $(LIB)
   $(LD) -o tc.exe $(LDFLAGS) tc_$(SOLVER).o  ../dyn_$(SOLVER)/module_initialize_$(IDEAL_CASE).o $(LIBWRFLIB) $(LIB)

convert_em : convert_em.o
   $(RANLIB) $(LIBWRFLIB)
   $(LD) -o convert_em.exe $(LDFLAGS) convert_em.o $(LIBWRFLIB) $(LIB)

convert_nmm : convert_nmm.o
   $(RANLIB) $(LIBWRFLIB)
   $(FC) -o convert_nmm.exe $(LDFLAGS) convert_nmm.o $(LIBWRFLIB) $(LIB)

real_nmm : real_nmm.o
   ( cd ../dyn_nmm ;  $(MAKE) module_initialize_real.o )
   $(RANLIB) $(LIBWRFLIB)
   $(FC) -o real_nmm.exe $(LDFLAGS) real_nmm.o $(LIBWRFLIB) $(LIB)

module_initialize : ../dyn_$(SOLVER)/module_initialize_$(IDEAL_CASE).o
   ( cd ../dyn_$(SOLVER) ;  $(MAKE) module_initialize_$(IDEAL_CASE).o )

## prevent real being compiled for OMP -- only for regtesting
#$(SOLVER)_real : module_initialize real_$(SOLVER).o
#   $(RANLIB) $(LIBWRFLIB)
#   if [ -z "$(OMP)" ] ; then $(FC) -o real.exe $(LDFLAGS) real_$(SOLVER).o ../dyn_$(SOLVER)/module_initialize_$(IDEAL_CASE).o $(LIBWRFLIB) $(LIB) ; fi
#
## prevent module_initialize being compiled for OMP --remove after IBM debugging
#module_initialize :
#   if [ -z "$(OMP)" ] ; then ( cd ../dyn_$(SOLVER) ;  $(MAKE) module_initialize_$(IDEAL_CASE).o ) ; fi
# end of regtest changes

clean:
   @ echo 'use the clean script'

# DEPENDENCIES : only dependencies after this line (don't remove the word DEPENDENCIES)

convert_nmm.o: \
   ../frame/module_machine.o \
   ../frame/module_domain.o \
   ../frame/module_driver_constants.o \
   ../frame/module_configure.o \
   ../frame/module_timing.o \
   ../frame/module_dm.o \
   ../share/module_bc.o \
   ../share/module_io_domain.o \
   $(ESMF_MOD_DEPENDENCE) 

convert_em.o: \
   ../frame/module_machine.o \
   ../frame/module_domain.o \
   ../frame/module_driver_constants.o \
   ../frame/module_configure.o \
   ../frame/module_timing.o \
   ../frame/module_dm.o \
   ../share/module_bc.o \
   ../share/module_io_domain.o \
   $(ESMF_MOD_DEPENDENCE)

ideal_em.o: \
   ../frame/module_machine.o \
   ../frame/module_domain.o \
   ../frame/module_driver_constants.o \
   ../frame/module_configure.o \
   ../frame/module_timing.o \
   ../frame/module_dm.o \
   ../share/module_io_domain.o \
   ../dyn_$(SOLVER)/$(CASE_MODULE) \
   $(WPS_DIR)/util/src/module_map_utils.o \
   $(ESMF_MOD_DEPENDENCE)


ideal_nmm.o: \
   ../dyn_$(SOLVER)/$(CASE_MODULE) \
   ../share/module_optional_input.o \
   ../share/module_io_domain.o \
   ../share/input_wrf.o

ndown_em.o: \
   ../frame/module_machine.o \
   ../frame/module_domain.o \
   ../frame/module_driver_constants.o \
   ../frame/module_configure.o \
   ../frame/module_timing.o \
   ../frame/module_dm.o \
   ../frame/module_wrf_error.o \
   ../frame/module_integrate.o \
   ../share/module_bc.o \
   ../share/module_io_domain.o \
   ../share/module_get_file_names.o \
   ../share/module_soil_pre.o \
   ../dyn_em/module_initialize_$(IDEAL_CASE).o \
   ../dyn_em/module_big_step_utilities_em.o \
   $(ESMF_MOD_DEPENDENCE)

nup_em.o: \
   ../frame/module_machine.o \
   ../frame/module_domain.o \
   ../frame/module_streams.o \
   ../frame/module_driver_constants.o \
   ../frame/module_configure.o \
   ../frame/module_timing.o \
   ../frame/module_dm.o \
   ../frame/module_wrf_error.o \
   ../frame/module_integrate.o \
   ../share/module_bc.o \
   ../share/module_io_domain.o \
   ../share/module_get_file_names.o \
   ../share/module_soil_pre.o \
   ../dyn_em/module_initialize_real.o \
   ../dyn_em/module_big_step_utilities_em.o \
   $(ESMF_MOD_DEPENDENCE)

# this already built above :../dyn_em/module_initialize.real.o \
real_em.o: \
   ../frame/module_machine.o \
   ../frame/module_domain.o \
   ../frame/module_driver_constants.o \
   ../frame/module_configure.o \
   ../frame/module_timing.o \
   ../frame/module_dm.o \
   ../dyn_em/module_initialize_$(IDEAL_CASE).o \
   ../dyn_em/module_big_step_utilities_em.o \
   ../share/module_io_domain.o \
   ../share/module_date_time.o \
   ../share/module_optional_input.o \
   ../share/module_bc_time_utilities.o \
        ../dyn_em/module_wps_io_arw.o \
   $(ESMF_MOD_DEPENDENCE)
#   ../chem/module_input_chem_data.o \
#   ../chem/module_input_chem_bioemiss.o \


tc_em.o: \
   ../frame/module_machine.o \
   ../frame/module_domain.o \
   ../frame/module_driver_constants.o \
   ../frame/module_configure.o \
   ../frame/module_timing.o \
   ../frame/module_dm.o \
   ../dyn_em/module_initialize_$(IDEAL_CASE).o \
   ../dyn_em/module_big_step_utilities_em.o \
   ../share/module_io_domain.o \
   ../share/module_date_time.o \
   ../share/module_optional_input.o \
   ../share/module_bc_time_utilities.o \
   $(ESMF_MOD_DEPENDENCE)



wrf.o:  ../main/module_wrf_top.o

wrf_ESMFMod.o:  ../main/module_wrf_top.o

wrf_SST_ESMF.o:  wrf_ESMFMod.o

module_wrf_top.o: ../frame/module_machine.o \
                  ../frame/module_domain.o \
                  ../frame/module_integrate.o \
                  ../frame/module_driver_constants.o \
                  ../frame/module_configure.o \
                  ../frame/module_timing.o \
                  ../frame/module_wrf_error.o \
                  ../frame/module_state_description.o \
                  $(ESMF_MOD_DEPENDENCE)

# DO NOT DELETE



the postamble file in arch dir. The lines in this files are included in configure.wrf whe you execute ./configure
Code: Select all
######################
#     Modified 24/06/2012   by A. Pérez Hortal - aperez1@famaf.unc.edu.ar                                ******************************
# POSTAMBLE

FGREP = fgrep -iq

#**  dd CONFIG_FLAG if the quarter_ss case is selected - This case was modified to include lat, long and map projection info    *******
ifeq ($(IDEAL_CASE), quarter_ss)
  HAVE_CONFIG_FLAG= -DCONFIGFLAG
else
  HAVE_CONFIG_FLAG=   
endif
#Then add HAVE_CONFIG_FLAG to ARCHFLAGS
#***************************************************************************************************************************************

ARCHFLAGS       =    $(COREDEFS) -DIWORDSIZE=$(IWORDSIZE) -DDWORDSIZE=$(DWORDSIZE) -DRWORDSIZE=$(RWORDSIZE) -DLWORDSIZE=$(LWORDSIZE) \
                     $(ARCH_LOCAL) \
                     $(DA_ARCHFLAGS) \
                      CONFIGURE_DMPARALLEL \
                      CONFIGURE_STUBMPI \
                      CONFIGURE_NETCDF_FLAG \
                      CONFIGURE_PNETCDF_FLAG \
                      CONFIGURE_ESMF_FLAG \
                      CONFIGURE_GRIB2_FLAG \
                      CONFIGURE_RTTOV_FLAG \
                      CONFIGURE_CRTM_FLAG \
                      CONFIGURE_4DVAR_FLAG \
                      CONFIGURE_WAVELET_FLAG \
                      CONFIGURE_NESTOPT \
                      -DUSE_ALLOCATABLES \
                      -DGRIB1 \
                      -DINTIO \
                      -DLIMIT_ARGS \
                      -DCONFIG_BUF_LEN=$(CONFIG_BUF_LEN) \
                      -DMAX_DOMAINS_F=$(MAX_DOMAINS) \
                      -DMAX_HISTORY=$(MAX_HISTORY) \
                      -DNMM_NEST=$(WRF_NMM_NEST)  \
                      $(HAVE_CONFIG_FLAG)   

 
CFLAGS          =    $(CFLAGS_LOCAL) CONFIGURE_DMPARALLEL CONFIGURE_STUBMPI \
                      -DMAX_HISTORY=$(MAX_HISTORY)
FCFLAGS         =    $(FCOPTIM) $(FCBASEOPTS)
ESMF_LIB_FLAGS  =    ESMFLIBFLAG
# ESMF 5 -- these are defined in esmf.mk, included above
#NOWIN ESMF_IO_LIB     =    ESMFIOLIB
ESMF_IO_LIB_EXT =    ESMFIOEXTLIB
INCLUDE_MODULES =    $(MODULE_SRCH_FLAG) \
                     $(ESMF_MOD_INC) $(ESMF_LIB_FLAGS) \
                      -I$(WRF_SRC_ROOT_DIR)/main \
                      -I$(WRF_SRC_ROOT_DIR)/external/io_netcdf \
                      -I$(WRF_SRC_ROOT_DIR)/external/io_int \
                      -I$(WRF_SRC_ROOT_DIR)/frame \
                      -I$(WRF_SRC_ROOT_DIR)/share \
                      -I$(WRF_SRC_ROOT_DIR)/phys \
                      -I$(WRF_SRC_ROOT_DIR)/chem -I$(WRF_SRC_ROOT_DIR)/inc \
                      -I$(NETCDFPATH)/include \
                      -I$(WPS_DIR)/util/src \
                      CONFIGURE_RTTOV_INC
REGISTRY        =    Registry

#NOWIN LIB_BUNDLED     = \
#NOWIN                      $(WRF_SRC_ROOT_DIR)/external/fftpack/fftpack5/libfftpack.a \
#NOWIN                      $(WRF_SRC_ROOT_DIR)/external/io_grib1/libio_grib1.a \
#NOWIN                      $(WRF_SRC_ROOT_DIR)/external/io_grib_share/libio_grib_share.a \
#NOWIN                      $(WRF_SRC_ROOT_DIR)/external/io_int/libwrfio_int.a \
#NOWIN                      $(ESMF_IO_LIB) \
#NOWIN                      CONFIGURE_COMMS_LIB \
#NOWIN                      $(WRF_SRC_ROOT_DIR)/frame/module_internal_header_util.o \
#NOWIN                      $(WRF_SRC_ROOT_DIR)/frame/pack_utils.o

#NOWIN LIB_EXTERNAL    = \
#NOWIN                      CONFIGURE_NETCDF_LIB_PATH CONFIGURE_PNETCDF_LIB_PATH CONFIGURE_GRIB2_LIB CONFIGURE_ATMOCN_LIB

LIB             =    $(LIB_BUNDLED) $(LIB_EXTERNAL) $(LIB_LOCAL)
LDFLAGS         =    $(OMP) $(FCFLAGS) $(LDFLAGS_LOCAL) CONFIGURE_LDFLAGS
ENVCOMPDEFS     =    CONFIGURE_COMPILEFLAGS
CPPFLAGS        =    $(ARCHFLAGS) $(ENVCOMPDEFS) -I$(LIBINCLUDE) $(TRADFLAG) CONFIGURE_COMMS_INCLUDE
NETCDFPATH      =    CONFIGURE_NETCDF_PATH
PNETCDFPATH     =    CONFIGURE_PNETCDF_PATH

bundled:  wrf_ioapi_includes wrfio_grib_share wrfio_grib1 wrfio_int esmf_time fftpack CONFIGURE_ATMOCN
external:  CONFIGURE_WRFIO_NF CONFIGURE_WRFIO_PNF CONFIGURE_WRFIO_GRIB2 CONFIGURE_COMMS_EXTERNAL $(ESMF_TARGET)

######################
externals: bundled external

gen_comms_serial :
   ( /bin/rm -f $(WRF_SRC_ROOT_DIR)/tools/gen_comms.c )

module_dm_serial :
   ( if [ ! -e module_dm.F ] ; then /bin/cp module_dm_warning module_dm.F ; cat module_dm_stubs.F >> module_dm.F ; fi )

gen_comms_rsllite :
   ( if [ ! -e $(WRF_SRC_ROOT_DIR)/tools/gen_comms.c ] ; then \
          /bin/cp $(WRF_SRC_ROOT_DIR)/tools/gen_comms_warning $(WRF_SRC_ROOT_DIR)/tools/gen_comms.c ; \
          cat $(WRF_SRC_ROOT_DIR)/external/RSL_LITE/gen_comms.c >> $(WRF_SRC_ROOT_DIR)/tools/gen_comms.c ; fi )

module_dm_rsllite :
   ( if [ ! -e module_dm.F ] ; then /bin/cp module_dm_warning module_dm.F ; \
          cat $(WRF_SRC_ROOT_DIR)/external/RSL_LITE/module_dm.F >> module_dm.F ; fi )

wrfio_nf :
   ( cd $(WRF_SRC_ROOT_DIR)/external/io_netcdf ; \
          make NETCDFPATH="$(NETCDFPATH)" RANLIB="$(RANLIB)" CPP="$(CPP)" \
          CC="$(SCC)" CFLAGS="$(CFLAGS)" \
          FC="$(SFC) $(PROMOTION) $(FCFLAGS)" TRADFLAG="$(TRADFLAG)" AR="$(AR)" ARFLAGS="$(ARFLAGS)" )

wrfio_pnf :
   ( cd $(WRF_SRC_ROOT_DIR)/external/io_pnetcdf ; \
          make NETCDFPATH="$(PNETCDFPATH)" RANLIB="$(RANLIB)" CPP="$(CPP) $(ARCHFLAGS)" \
          FC="$(FC) $(PROMOTION) $(FCFLAGS)" TRADFLAG="$(TRADFLAG)" AR="$(AR)" ARFLAGS="$(ARFLAGS)" )

wrfio_grib_share :
   ( cd $(WRF_SRC_ROOT_DIR)/external/io_grib_share ; \
          make CC="$(SCC)" CFLAGS="$(CFLAGS)" RM="$(RM)" RANLIB="$(RANLIB)" CPP="$(CPP)" \
          FC="$(SFC) $(PROMOTION) -I. $(FCDEBUG) $(FCBASEOPTS) $(FCSUFFIX)" TRADFLAG="$(TRADFLAG)" AR="$(AR)" ARFLAGS="$(ARFLAGS)" archive)

wrfio_grib1 :
   ( cd $(WRF_SRC_ROOT_DIR)/external/io_grib1 ; \
          make CC="$(SCC)" CFLAGS="$(CFLAGS)" RM="$(RM)" RANLIB="$(RANLIB)" CPP="$(CPP)" \
          FC="$(SFC) $(PROMOTION) -I. $(FCDEBUG) $(FCBASEOPTS) $(FCSUFFIX)" TRADFLAG="$(TRADFLAG)" AR="$(AR)" ARFLAGS="$(ARFLAGS)" archive)

wrfio_grib2 :
   ( cd $(WRF_SRC_ROOT_DIR)/external/io_grib2 ; \
          make CC="$(SCC)" CFLAGS="$(CFLAGS) CONFIGURE_GRIB2_INC" RM="$(RM)" RANLIB="$(RANLIB)" \
          CPP="$(CPP)" \
          FC="$(SFC) $(PROMOTION) -I. $(FCDEBUG) $(FCBASEOPTS) $(FCSUFFIX)" TRADFLAG="-traditional" AR="$(AR)" ARFLAGS="$(ARFLAGS)" \
          FIXED="$(FORMAT_FIXED)" archive)

wrfio_int :
   ( cd $(WRF_SRC_ROOT_DIR)/external/io_int ; \
          make CC="$(CC)" RM="$(RM)" RANLIB="$(RANLIB)" CPP="$(CPP)" \
          FC="$(SFC) $(PROMOTION) $(FCDEBUG) $(FCBASEOPTS)" FGREP="$(FGREP)" \
          TRADFLAG="$(TRADFLAG)" AR="$(AR)" ARFLAGS="$(ARFLAGS)" ARCHFLAGS="$(ARCHFLAGS)" all )

esmf_time :
   ( cd $(WRF_SRC_ROOT_DIR)/external/esmf_time_f90 ; \
          make FC="$(SFC) $(PROMOTION) $(FCDEBUG) $(FCBASEOPTS)" RANLIB="$(RANLIB)" \
          CPP="$(CPP) -I$(WRF_SRC_ROOT_DIR)/inc -I. $(ARCHFLAGS) $(TRADFLAG)" AR="$(AR)" ARFLAGS="$(ARFLAGS)" )

fftpack :
   ( cd $(WRF_SRC_ROOT_DIR)/external/fftpack/fftpack5 ; \
          make FC="$(SFC)" FFLAGS="$(PROMOTION) $(FCDEBUG) $(FCBASEOPTS)" RANLIB="$(RANLIB)" AR="$(AR)" ARFLAGS="$(ARFLAGS)" )

atm_ocn :
   ( cd $(WRF_SRC_ROOT_DIR)/external/atm_ocn ; \
          make CC="$(SCC)" CFLAGS="$(CFLAGS) " RM="$(RM)" RANLIB="$(RANLIB)" \
          CPP="$(CPP)" \
          FC="$(DM_FC) $(PROMOTION) -I. $(FCDEBUG) $(FCBASEOPTS) $(FCSUFFIX)" TRADFLAG="-traditional" AR="$(AR)" ARFLAGS="$(ARFLAGS)" \
          FIXED="$(FORMAT_FIXED)" )

$(WRF_SRC_ROOT_DIR)/external/RSL_LITE/librsl_lite.a :
   ( cd $(WRF_SRC_ROOT_DIR)/external/RSL_LITE ; make CC="$(CC) $(CFLAGS)" \
          FC="$(FC) $(FCFLAGS) $(PROMOTION) $(BYTESWAPIO)" \
          CPP="$(CPP) -I. $(ARCHFLAGS) $(TRADFLAG)" AR="$(AR)" ARFLAGS="$(ARFLAGS)" ;\
          $(RANLIB) $(WRF_SRC_ROOT_DIR)/external/RSL_LITE/librsl_lite.a )

######################
#   Macros, these should be generic for all machines

LN   =   ln -sf
MAKE   =   make -i -r
RM   =    rm -f


# These sub-directory builds are identical across all architectures

wrf_ioapi_includes :
   ( cd $(WRF_SRC_ROOT_DIR)/external/ioapi_share ; \
          $(MAKE) NATIVE_RWORDSIZE="$(NATIVE_RWORDSIZE)" RWORDSIZE="$(RWORDSIZE)" AR="$(AR)" ARFLAGS="$(ARFLAGS)" )

wrfio_esmf :
   ( cd $(WRF_SRC_ROOT_DIR)/external/io_esmf ; \
          make FC="$(FC) $(PROMOTION) $(FCDEBUG) $(FCBASEOPTS) $(ESMF_MOD_INC)" \
          RANLIB="$(RANLIB)" CPP="$(CPP) $(POUND_DEF) " AR="$(AR)" ARFLAGS="$(ARFLAGS)" )

#   There is probably no reason to modify these rules

.F.i:
   $(RM) $@
   $(CPP) -I$(WRF_SRC_ROOT_DIR)/inc $(CPPFLAGS) $*.F > $@
   mv $*.i $(DEVTOP)/pick/$*.f90
   cp $*.F $(DEVTOP)/pick

.F.o:
   $(RM) $@
   $(CPP) -I$(WRF_SRC_ROOT_DIR)/inc $(CPPFLAGS) $(OMPCPP) $*.F  > $*.bb
   $(SED_FTN) $*.bb | $(CPP) > $*.f90
   $(RM) $*.b $*.bb
   @ if echo $(ARCHFLAGS) | $(FGREP) 'DVAR4D'; then \
          echo COMPILING $*.F for 4DVAR ; \
          $(WRF_SRC_ROOT_DIR)/var/build/da_name_space.pl $*.f90 > $*.f90.tmp ; \
          mv $*.f90.tmp $*.f90 ; \
        fi
   if $(FGREP) '!$$OMP' $*.f90 ; then \
          if [ -n "$(OMP)" ] ; then echo COMPILING $*.F WITH OMP ; fi ; \
     $(FC) -o $@ -c $(FCFLAGS) $(OMP) $(MODULE_DIRS) $(PROMOTION) $(FCSUFFIX) $*.f90 ; \
        else \
          if [ -n "$(OMP)" ] ; then echo COMPILING $*.F WITHOUT OMP ; fi ; \
     $(FC) -o $@ -c $(FCFLAGS) $(MODULE_DIRS) $(PROMOTION) $(FCSUFFIX) $*.f90 ; \
        fi
       

.F.f90:
   $(RM) $@
   $(SED_FTN) $*.F > $*.b
   $(CPP) -I$(WRF_SRC_ROOT_DIR)/inc $(CPPFLAGS) $*.b  > $@
   $(RM) $*.b

.f90.o:
   $(RM) $@
   $(FC) -o $@ -c $(FCFLAGS) $(PROMOTION) $(FCSUFFIX) $*.f90

setfeenv.o : setfeenv.c
   $(RM) $@
   $(CCOMP) -o $@ -c $(CFLAGS) $(OMPCC) $*.c

.c.o:
   $(RM) $@
   $(CC) -o $@ -c $(CFLAGS) $*.c

andresperezcba
 
Posts: 6
Joined: Sun Apr 15, 2012 5:23 pm

Re: Specifying lat/lon on idealized Sea Breeze

Postby andresperezcba » Wed Aug 01, 2012 6:21 pm

I found a bug in the solution that I post If you use a map projection some definitions are wrong. The map scale factors need to be calculated and I think that more variables also. I'm not found a solution to that yet.
Regards
Andres
andresperezcba
 
Posts: 6
Joined: Sun Apr 15, 2012 5:23 pm

Re: Specifying lat/lon on idealized Sea Breeze

Postby varun_murthy » Fri Nov 29, 2013 11:07 pm

A good example to look at for setting the latitudes and longitudes is the held_suarez global idealized case. This idealized case is global, and hence sets the latitude and longitude everywhere. Modify the initialization code accordingly to set latitudes and longitudes in your regional domain.
varun_murthy
 
Posts: 5
Joined: Fri Nov 29, 2013 10:38 pm

Re: Specifying lat/lon on idealized Sea Breeze

Postby Seishou » Mon Jul 07, 2014 2:09 pm

I have used your some parts of code and got the following error:

how to fix that?

module_initialize_tropical_cyclone.f90(96): error #6457: This derived type name has not been declared. [PROJ_INFO]
TYPE(proj_info) :: proj
--------^
module_initialize_tropical_cyclone.f90(185): error #6632: Keyword arguments are invalid without an explicit interface. [PROJ_CODE]
CALL map_set( proj_code=config_flags%map_proj, proj=proj, lat1=config_flags%cen_lat, truelat1=config_flags%truelat1, &
-----------------^
module_initialize_tropical_cyclone.f90(185): error #6632: Keyword arguments are invalid without an explicit interface. [PROJ]
CALL map_set( proj_code=config_flags%map_proj, proj=proj, lat1=config_flags%cen_lat, truelat1=config_flags%truelat1, &
--------------------------------------------------^
module_initialize_tropical_cyclone.f90(185): error #6404: This name does not have a type, and must have an explicit type. [PROJ]
CALL map_set( proj_code=config_flags%map_proj, proj=proj, lat1=config_flags%cen_lat, truelat1=config_flags%truelat1, &
-------------------------------------------------------^
module_initialize_tropical_cyclone.f90(185): error #6632: Keyword arguments are invalid without an explicit interface. [LAT1]
CALL map_set( proj_code=config_flags%map_proj, proj=proj, lat1=config_flags%cen_lat, truelat1=config_flags%truelat1, &
-------------------------------------------------------------^
module_initialize_tropical_cyclone.f90(185): error #6632: Keyword arguments are invalid without an explicit interface. [TRUELAT1]
CALL map_set( proj_code=config_flags%map_proj, proj=proj, lat1=config_flags%cen_lat, truelat1=config_flags%truelat1, &
----------------------------------------------------------------------------------------^
module_initialize_tropical_cyclone.f90(186): error #6632: Keyword arguments are invalid without an explicit interface. [LON1]
lon1=config_flags%cen_lon, knowni=real(nxc), knownj=real(nyc), dx=grid%dx, dy=grid%dy)
-----------------^
module_initialize_tropical_cyclone.f90(186): error #6632: Keyword arguments are invalid without an explicit interface. [KNOWNI]
lon1=config_flags%cen_lon, knowni=real(nxc), knownj=real(nyc), dx=grid%dx, dy=grid%dy)
--------------------------------------------^
module_initialize_tropical_cyclone.f90(186): error #6632: Keyword arguments are invalid without an explicit interface. [KNOWNJ]
lon1=config_flags%cen_lon, knowni=real(nxc), knownj=real(nyc), dx=grid%dx, dy=grid%dy)
--------------------------------------------------------------^
module_initialize_tropical_cyclone.f90(186): error #6632: Keyword arguments are invalid without an explicit interface. [DX]
lon1=config_flags%cen_lon, knowni=real(nxc), knownj=real(nyc), dx=grid%dx, dy=grid%dy)
--------------------------------------------------------------------------------^
module_initialize_tropical_cyclone.f90(186): error #6632: Keyword arguments are invalid without an explicit interface. [DY]
lon1=config_flags%cen_lon, knowni=real(nxc), knownj=real(nyc), dx=grid%dx, dy=grid%dy)
--------------------------------------------------------------------------------------------^
compilation aborted for module_initialize_tropical_cyclone.f90 (code 1)
make[2]: [module_initialize_tropical_cyclone.o] Error 1 (ignored)
make[2]: Leaving directory `/work1/seishou/IdealWRF/WRFV3/dyn_em'
make[1]: *** No rule to make target `/util/src/module_map_utils.o', needed by `ideal_em.o'. Stop.
make[1]: Leaving directory `/work1/seishou/IdealWRF/WRFV3/main'


Code: Select all
!IDEAL:MODEL_LAYER:INITIALIZATION
!

!  This MODULE holds the routines which are used to perform various initializations
!  for the individual domains. 

!  This MODULE CONTAINS the following routines:

!  initialize_field_test - 1. Set different fields to different constant
!                             values.  This is only a test.  If the correct
!                             domain is not found (based upon the "id")
!                             then a fatal error is issued.               

!-----------------------------------------------------------------------

MODULE module_initialize_ideal

   USE module_domain
   USE module_io_domain
   USE module_state_description
   USE module_model_constants
   USE module_bc
   USE module_timing
   USE module_configure
   USE module_init_utilities
   USE module_soil_pre
#ifdef DM_PARALLEL
   USE module_dm
#endif


CONTAINS


!-------------------------------------------------------------------
! this is a wrapper for the solver-specific init_domain routines.
! Also dereferences the grid variables and passes them down as arguments.
! This is crucial, since the lower level routines may do message passing
! and this will get fouled up on machines that insist on passing down
! copies of assumed-shape arrays (by passing down as arguments, the
! data are treated as assumed-size -- ie. f77 -- arrays and the copying
! business is avoided).  Fie on the F90 designers.  Fie and a pox.
! NOTE:  Modified to remove all but arrays of rank 4 or more from the
!        argument list.  Arrays with rank>3 are still problematic due to the
!        above-noted fie- and pox-ities.  TBH 20061129. 

   SUBROUTINE init_domain ( grid , config_flags)

   IMPLICIT NONE

   !  Input data.
   TYPE (domain), POINTER :: grid
   !Output data
   TYPE (grid_config_rec_type)    :: config_flags



   !  Local data.
   INTEGER :: idum1, idum2

   CALL set_scalar_indices_from_config ( head_grid%id , idum1, idum2 )

     CALL init_domain_rk( grid, config_flags &
!
#include <actual_new_args.inc>
!
                        )
   END SUBROUTINE init_domain

!-------------------------------------------------------------------

   SUBROUTINE init_domain_rk ( grid , config_flags &
!
# include <dummy_new_args.inc>
!
)
   IMPLICIT NONE

   !  Input data.
   TYPE (domain), POINTER       :: grid
   TYPE (grid_config_rec_type)  :: config_flags

# include <dummy_new_decl.inc>


   !  Local data
   TYPE(proj_info)              :: proj

   INTEGER                             ::                       &
                                  ids, ide, jds, jde, kds, kde, &
                                  ims, ime, jms, jme, kms, kme, &
                                  its, ite, jts, jte, kts, kte, &
                                  i, j, k
   !  Local data

   INTEGER, PARAMETER :: nl_max = 1000
   REAL, DIMENSION(nl_max) :: zk, p_in, theta, rho, u, v, qv, pd_in
   INTEGER :: nl_in

   INTEGER :: icm,jcm, ii, im1, jj, jm1, loop, error, fid, nxc, nyc, lm
   REAL    :: u_mean,v_mean, f0, p_surf, p_level, qvf, z_at_v, z_at_u
   REAL    :: z_scale, xrad, yrad, zrad, rad, delt, cof1, cof2
!   REAL, EXTERNAL :: interp_0
   REAL    :: pi, rnd

! variables/arrays for analytic vortex:
    integer :: nref,kref,nloop,i1,i2
    real :: r0,zdd,dd1,dd2,xref,vr,fcor,qvs,e1,tx,px,qx,ric,rjc,rr,diff,sst
    real*8 :: rmax,vmax,frac,angle
    real, dimension(:), allocatable :: rref,zref,th0,qv0,thv0,prs0,pi0,rh0
    real, dimension(:,:), allocatable :: vref,piref,pref,thref,thvref,qvref

!  stuff from original initialization that has been dropped from the Registry
   REAL    :: vnu, xnu, xnus, dinit0, cbh, p0_temp, t0_temp, zd, zt
   REAL    :: qvf1, qvf2, pd_surf
   INTEGER :: it
   real :: thtmp, ptmp, temp(3)

   LOGICAL :: moisture_init
   LOGICAL :: stretch_grid, dry_sounding
   character (len=256) :: mminlu2
   
#ifdef DM_PARALLEL
#    include <data_calls.inc>
#endif


   SELECT CASE ( model_data_order )
         CASE ( DATA_ORDER_ZXY )
   kds = grid%sd31 ; kde = grid%ed31 ;
   ids = grid%sd32 ; ide = grid%ed32 ;
   jds = grid%sd33 ; jde = grid%ed33 ;

   kms = grid%sm31 ; kme = grid%em31 ;
   ims = grid%sm32 ; ime = grid%em32 ;
   jms = grid%sm33 ; jme = grid%em33 ;

   kts = grid%sp31 ; kte = grid%ep31 ;   ! note that tile is entire patch
   its = grid%sp32 ; ite = grid%ep32 ;   ! note that tile is entire patch
   jts = grid%sp33 ; jte = grid%ep33 ;   ! note that tile is entire patch
         CASE ( DATA_ORDER_XYZ )
   ids = grid%sd31 ; ide = grid%ed31 ;
   jds = grid%sd32 ; jde = grid%ed32 ;
   kds = grid%sd33 ; kde = grid%ed33 ;

   ims = grid%sm31 ; ime = grid%em31 ;
   jms = grid%sm32 ; jme = grid%em32 ;
   kms = grid%sm33 ; kme = grid%em33 ;

   its = grid%sp31 ; ite = grid%ep31 ;   ! note that tile is entire patch
   jts = grid%sp32 ; jte = grid%ep32 ;   ! note that tile is entire patch
   kts = grid%sp33 ; kte = grid%ep33 ;   ! note that tile is entire patch
         CASE ( DATA_ORDER_XZY )
   ids = grid%sd31 ; ide = grid%ed31 ;
   kds = grid%sd32 ; kde = grid%ed32 ;
   jds = grid%sd33 ; jde = grid%ed33 ;

   ims = grid%sm31 ; ime = grid%em31 ;
   kms = grid%sm32 ; kme = grid%em32 ;
   jms = grid%sm33 ; jme = grid%em33 ;

   its = grid%sp31 ; ite = grid%ep31 ;   ! note that tile is entire patch
   kts = grid%sp32 ; kte = grid%ep32 ;   ! note that tile is entire patch
   jts = grid%sp33 ; jte = grid%ep33 ;   ! note that tile is entire patch

   END SELECT

!-----------------------------------------------------------------------
!               USER SETTINGS

!  Parameters for analytic vortex:
!  Reference:  Rotunno and Emanuel, 1987, JAS, p. 549

    r0     =   412500.0     ! outer radius (m)
    rmax   =    82500.0     ! approximate radius of max winds (m)
    vmax   =       25.0     ! approximate value of max wind speed (m/s)
    zdd    =    20000.0     ! depth of vortex (m)


! other settings:

    fcor   =  5.0e-5        ! Coriolis parameter (1/s)
    sst    =  30.0          ! sea-surface temperature (Celsius)

!-----------------------------------------------------------------------

   stretch_grid = .true.
   delt = 6.
!   z_scale = .50
   z_scale = .40
   pi = 2.*asin(1.0)
   write(6,*) ' pi is ',pi
   nxc = (ide-ids)/2
   nyc = jde/2
   icm = ide/2
! lm is the half width of the land in terms of grid points
   lm = 25
   write(6,*) 'lm,icm-lm,icm+lm = ', lm,icm-lm,icm+lm

!*********Set of lat long and map projection info*******************
    mminlu2 = ' '
    mminlu2(1:4) = 'USGS'
    CALL nl_set_mminlu(grid%id, '    ')
!   CALL nl_set_mminlu(1, 'USGS')
    CALL nl_set_iswater(grid%id,16)
    CALL nl_set_isice(1,3)
    CALL nl_set_cen_lat(grid%id,20.)
    CALL nl_set_cen_lon(grid%id,-105.)
    CALL nl_set_truelat1(grid%id,20.)
    CALL nl_set_truelat2(grid%id,0.)
    CALL nl_set_moad_cen_lat (grid%id,20.)
    CALL nl_set_stand_lon (grid%id,-105.)
    CALL nl_set_pole_lon (grid%id,0.)
    CALL nl_set_pole_lat (grid%id,90.)
    CALL nl_set_map_proj(grid%id,3)
!   CALL model_to_grid_config_rec(1,model_config_rec,config_flags)
    CALL nl_get_iswater(1,grid%iswater)

!****        Save the lat, long and map projections info    in model_config_rec  ***********************
   CALL model_to_grid_config_rec ( grid%id , model_config_rec , config_flags )

!****  Populate the proj data type if undefined values are present   **********************************
   CALL map_set( proj_code=config_flags%map_proj, proj=proj, lat1=config_flags%cen_lat, truelat1=config_flags%truelat1, &
                 lon1=config_flags%cen_lon, knowni=real(nxc), knownj=real(nyc), dx=grid%dx, dy=grid%dy)
!****  Set the xlat and xlong values for all grid points   *********************************************
   DO j = jts, jte
    DO i = its, ite
      CALL ij_to_latlon(proj,real(i) , real(j), grid%xlat(i,j), grid%xlong(i,j))
    END DO
  END DO



! here we check to see if the boundary conditions are set properly

   CALL boundary_condition_check( config_flags, bdyzone, error, grid%id )

   moisture_init = .true.

    grid%itimestep=0

#ifdef DM_PARALLEL
   CALL wrf_dm_bcast_bytes( icm , IWORDSIZE )
   CALL wrf_dm_bcast_bytes( jcm , IWORDSIZE )
#endif


!  here we initialize data that currently is not initialized
!  in the input data


    DO j = jts, jte
      DO i = its, ite
         grid%ht(i,j)       = 0.
         grid%msft(i,j)     = 1.
         grid%msfu(i,j)     = 1.
         grid%msfv(i,j)     = 1.
         grid%msftx(i,j)    = 1.
         grid%msfty(i,j)    = 1.
         grid%msfux(i,j)    = 1.
         grid%msfuy(i,j)    = 1.
         grid%msfvx(i,j)    = 1.
         grid%msfvy(i,j)    = 1.
         grid%msfvx_inv(i,j)= 1.
         grid%sina(i,j)     = 0.
         grid%cosa(i,j)     = 1.
         grid%xlong(i,j)    = 0.0
         grid%e(i,j)        = 0.0
         grid%f(i,j)        = fcor
         grid%xlat(i,j)     = asin(0.5*fcor/EOMEG)/DEGRAD
! Hard-wire the ocean configuration
    grid%xland(i,j)     = 2.
         grid%lu_index(i,j)  = 16
         grid%tsk(i,j) = 273.15 + sst
         ! I think tmn is not used for ocean points, but set a value anyway:
         grid%tmn(i,j) = grid%tsk(i,j) - 10.0
      END DO
   END DO

    print *,'   f = ',grid%f(its,jts)
    print *,'   lat = ',grid%xlat(its,jts)

! for Noah LSM, additional variables need to be initialized

   other_masked_fields : SELECT CASE ( model_config_rec%sf_surface_physics(grid%id) )

      CASE (SLABSCHEME)

      CASE (LSMSCHEME)

        DO j = jts , MIN(jde-1,jte)
           DO i = its , MIN(ide-1,ite)
              IF (grid%xland(i,j) .lt. 1.5) THEN
                 grid%vegfra(i,j) = 0.5
                 grid%canwat(i,j) = 0.
                 grid%ivgtyp(i,j) = 18
                 grid%isltyp(i,j) = 8
                 grid%xice(i,j) = 0.
                 grid%snow(i,j) = 0.
              ELSE
                 grid%vegfra(i,j) = 0.
                 grid%canwat(i,j) = 0.
                 grid%ivgtyp(i,j) = 16
                 grid%isltyp(i,j) = 14
                 grid%xice(i,j) = 0.
                 grid%snow(i,j) = 0.
              ENDIF
           END DO
        END DO

      CASE (RUCLSMSCHEME)

   END SELECT other_masked_fields

   DO j = jts, jte
    DO k = kts, kte
      DO i = its, ite
         grid%ww(i,k,j)     = 0.
      END DO
    END DO
   END DO

   grid%step_number = 0

! Process the soil; note that there are some things hard-wired into share/module_soil_pre.F
      CALL process_soil_ideal(grid%xland,grid%xice,grid%vegfra,grid%snow,grid%canwat, &
                     grid%ivgtyp,grid%isltyp,grid%tslb,grid%smois, &
                     grid%tsk,grid%tmn,grid%zs,grid%dzs,model_config_rec%num_soil_layers, &
                     model_config_rec%sf_surface_physics(grid%id), &
                                   ids,ide, jds,jde, kds,kde,&
                                   ims,ime, jms,jme, kms,kme,&
                                   its,ite, jts,jte, kts,kte )

! set up the grid

   IF (stretch_grid) THEN ! exponential stretch for eta (nearly constant dz)
     DO k=1, kde
      grid%znw(k) = (exp(-(k-1)/float(kde-1)/z_scale) - exp(-1./z_scale))/ &
                                (1.-exp(-1./z_scale))
     ENDDO
   ELSE
     DO k=1, kde
      grid%znw(k) = 1. - float(k-1)/float(kde-1)
     ENDDO
   ENDIF

   DO k=1, kde-1
    grid%dnw(k) = grid%znw(k+1) - grid%znw(k)
    grid%rdnw(k) = 1./grid%dnw(k)
    grid%znu(k) = 0.5*(grid%znw(k+1)+grid%znw(k))
   ENDDO
   DO k=2, kde-1
    grid%dn(k) = 0.5*(grid%dnw(k)+grid%dnw(k-1))
    grid%rdn(k) = 1./grid%dn(k)
    grid%fnp(k) = .5* grid%dnw(k  )/grid%dn(k)
    grid%fnm(k) = .5* grid%dnw(k-1)/grid%dn(k)
   ENDDO

   cof1 = (2.*grid%dn(2)+grid%dn(3))/(grid%dn(2)+grid%dn(3))*grid%dnw(1)/grid%dn(2)
   cof2 =     grid%dn(2)        /(grid%dn(2)+grid%dn(3))*grid%dnw(1)/grid%dn(3)
   grid%cf1  = grid%fnp(2) + cof1
   grid%cf2  = grid%fnm(2) - cof1 - cof2
   grid%cf3  = cof2       

   grid%cfn  = (.5*grid%dnw(kde-1)+grid%dn(kde-1))/grid%dn(kde-1)
   grid%cfn1 = -.5*grid%dnw(kde-1)/grid%dn(kde-1)
   grid%rdx = 1./config_flags%dx
   grid%rdy = 1./config_flags%dy

!  get the sounding from the ascii sounding file, first get dry sounding and
!  calculate base state

  write(6,*) ' getting dry sounding for base state '
  dry_sounding = .true.
  CALL get_sounding( zk, p_in, pd_in, theta, rho, u, v, qv, dry_sounding, nl_max, nl_in )

  write(6,*) ' returned from reading sounding, nl_in is ',nl_in

!  find ptop for the desired ztop (ztop is input from the namelist),
!  and find surface pressure

  grid%p_top = interp_0( p_in, zk, config_flags%ztop, nl_in )

  DO j=jts,jte
  DO i=its,ite  ! flat surface
    grid%phb(i,1,j) = 0.
    grid%php(i,1,j) = 0.
    grid%ph0(i,1,j) = 0.
    grid%ht(i,j) = 0.
  ENDDO
  ENDDO

  DO J = jts, jte
  DO I = its, ite

    p_surf = interp_0( p_in, zk, grid%phb(i,1,j)/g, nl_in )
    grid%mub(i,j) = p_surf-grid%p_top

!  this is dry hydrostatic sounding (base state), so given grid%p (coordinate),
!  interp theta (from interp) and compute 1/rho from eqn. of state

    DO K = 1, kte-1
      p_level = grid%znu(k)*(p_surf - grid%p_top) + grid%p_top
      grid%pb(i,k,j) = p_level
      grid%t_init(i,k,j) = interp_0( theta, p_in, p_level, nl_in ) - t0
      grid%alb(i,k,j) = (r_d/p1000mb)*(grid%t_init(i,k,j)+t0)*(grid%pb(i,k,j)/p1000mb)**cvpm
    ENDDO

!  calc hydrostatic balance (alternatively we could interp the geopotential from the
!  sounding, but this assures that the base state is in exact hydrostatic balance with
!  respect to the model eqns.

    DO k  = 2,kte
      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)
    ENDDO

  ENDDO
  ENDDO

  write(6,*) ' ptop is ',grid%p_top
  write(6,*) ' base state grid%mub(1,1), p_surf is ',grid%mub(1,1),grid%mub(1,1)+grid%p_top

!  calculate full state for each column - this includes moisture.

  write(6,*) ' getting moist sounding for full state '
  dry_sounding = .false.
  CALL get_sounding( zk, p_in, pd_in, theta, rho, u, v, qv, dry_sounding, nl_max, nl_in )

  DO J = jts, min(jde-1,jte)
  DO I = its, min(ide-1,ite)

!  At this point grid%p_top is already set. find the DRY mass in the column
!  by interpolating the DRY pressure. 

   pd_surf = interp_0( pd_in, zk, grid%phb(i,1,j)/g, nl_in )

!  compute the perturbation mass and the full mass

    grid%mu_1(i,j) = pd_surf-grid%p_top - grid%mub(i,j)
    grid%mu_2(i,j) = grid%mu_1(i,j)
    grid%mu0(i,j) = grid%mu_1(i,j) + grid%mub(i,j)

! given the dry pressure and coordinate system, interp the potential
! temperature and qv

    do k=1,kde-1

      p_level = grid%znu(k)*(pd_surf - grid%p_top) + grid%p_top

      moist(i,k,j,P_QV) = interp_0( qv, pd_in, p_level, nl_in )
      grid%t_1(i,k,j)          = interp_0( theta, pd_in, p_level, nl_in ) - t0
      grid%t_2(i,k,j)          = grid%t_1(i,k,j)
     

    enddo

!  integrate the hydrostatic equation (from the RHS of the bigstep
!  vertical momentum equation) down from the top to get grid%p.
!  first from the top of the model to the top pressure

    k = kte-1  ! top level

    qvf1 = 0.5*(moist(i,k,j,P_QV)+moist(i,k,j,P_QV))
    qvf2 = 1./(1.+qvf1)
    qvf1 = qvf1*qvf2

!    grid%p(i,k,j) = - 0.5*grid%mu_1(i,j)/grid%rdnw(k)
    grid%p(i,k,j) = - 0.5*(grid%mu_1(i,j)+qvf1*grid%mub(i,j))/grid%rdnw(k)/qvf2
    qvf = 1. + rvovrd*moist(i,k,j,P_QV)
    grid%alt(i,k,j) = (r_d/p1000mb)*(grid%t_1(i,k,j)+t0)*qvf* &
                (((grid%p(i,k,j)+grid%pb(i,k,j))/p1000mb)**cvpm)
    grid%al(i,k,j) = grid%alt(i,k,j) - grid%alb(i,k,j)

!  down the column

    do k=kte-2,1,-1
      qvf1 = 0.5*(moist(i,k,j,P_QV)+moist(i,k+1,j,P_QV))
      qvf2 = 1./(1.+qvf1)
      qvf1 = qvf1*qvf2
      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)
      qvf = 1. + rvovrd*moist(i,k,j,P_QV)
      grid%alt(i,k,j) = (r_d/p1000mb)*(grid%t_1(i,k,j)+t0)*qvf* &
                  (((grid%p(i,k,j)+grid%pb(i,k,j))/p1000mb)**cvpm)
      grid%al(i,k,j) = grid%alt(i,k,j) - grid%alb(i,k,j)
    enddo

!  this is the hydrostatic equation used in the model after the
!  small timesteps.  In the model, grid%al (inverse density)
!  is computed from the geopotential.


    grid%ph_1(i,1,j) = 0.
    DO k  = 2,kte
      grid%ph_1(i,k,j) = grid%ph_1(i,k-1,j) - (1./grid%rdnw(k-1))*(       &
                   (grid%mub(i,j)+grid%mu_1(i,j))*grid%al(i,k-1,j)+ &
                    grid%mu_1(i,j)*grid%alb(i,k-1,j)  )
                                                   
      grid%ph_2(i,k,j) = grid%ph_1(i,k,j)
      grid%ph0(i,k,j) = grid%ph_1(i,k,j) + grid%phb(i,k,j)
    ENDDO

    if((i==2) .and. (j==2)) then
     write(6,*) ' grid%ph_1 calc ',grid%ph_1(2,1,2),grid%ph_1(2,2,2),&
                              grid%mu_1(2,2)+grid%mub(2,2),grid%mu_1(2,2), &
                              grid%alb(2,1,2),grid%al(1,2,1),grid%rdnw(1)
    endif

  ENDDO
  ENDDO

!-----------------------------------------------------------------------
!  Analytic vortex.
!  Reference:  Rotunno and Emanuel, 1987, JAS, p. 549

    dd2 = 2.0 * rmax / ( r0 + rmax )

    nref = 1 + int( float(ide-ids+1)/2.0 )
    kref = kte-1

    print *,'  ids,ide,kds,kds  = ',ids,ide,kds,kde
    print *,'  its,ite,kts,kts  = ',its,ite,kts,kte
    print *,'  nref,fcor        = ',nref,fcor
    print *,'  r0,rmax,vmax,zdd = ',r0,rmax,vmax,zdd

    allocate(  rref(nref)         )
    allocate(  zref(0:kref+1)     )
    allocate(   th0(0:kref+1)     )
    allocate(   qv0(0:kref+1)     )
    allocate(  thv0(0:kref+1)     )
    allocate(  prs0(0:kref+1)     )
    allocate(   pi0(0:kref+1)     )
    allocate(   rh0(0:kref+1)     )
    allocate(  vref(nref,0:kref+1))
    allocate( piref(nref,0:kref+1))
    allocate(  pref(nref,0:kref+1))
    allocate( thref(nref,0:kref+1))
    allocate(thvref(nref,0:kref+1))
    allocate( qvref(nref,0:kref+1))

    ! get base state:
    print *,'  zref,th0,qv0,thv0:'
    do k=1,kref
      th0(k) = t0+grid%t_1(1,k,1)
      qv0(k) = moist(1,k,1,P_QV)
      thv0(k) = th0(k)*(1.0+(r_v/r_d)*qv0(k))/(1.0+qv0(k))
      zref(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
      print *,k,zref(k),th0(k),qv0(k),thv0(k)
    enddo

    print *,'  prs0,pi0,rh0:'
    do k=1,kref
      prs0(k) = grid%p(1,k,1)+grid%pb(1,k,1)
      pi0(k) = (prs0(k)/p0)**(r_d/cp)
      E1=1000.0*SVP1*EXP(SVP2*(th0(k)*pi0(k)-SVPT0)/(th0(k)*pi0(k)-SVP3))
      qvs = EP_2*E1/(prs0(k)-E1)
      rh0(k) = qv0(k)/qvs
      print *,k,prs0(k),pi0(k),rh0(k)
    enddo

    zref(0) = -zref(1)
    zref(kref+1) = zref(kref)+(zref(kref)-zref(kref-1))

      rref=0.0
      vref=0.0
     piref=0.0
      pref=0.0
     thref=0.0
    thvref=0.0
     qvref=0.0

    do i=1,nref
      rref(i) = config_flags%dx*(float(i-1)+0.5)
    enddo

    print *,'  zref:'
    do k=0,kref+1
      print *,k,zref(k)
    enddo

    print *,'  vref:'
    do k=1,kref
      do i=1,nref
        if(rref(i).lt.r0)then
          dd1 = 2.0 * rmax / ( rref(i) + rmax )
          vr = sqrt( vmax**2 * (rref(i)/rmax)**2     &
          * ( dd1 ** 3 - dd2 ** 3 ) + 0.25*fcor*fcor*rref(i)*rref(i) )   &
                  - 0.5 * fcor * rref(i)
        else
          vr = 0.0
        endif
        if(zref(k).lt.zdd)then
          vref(i,k) = vr * (zdd-zref(k))/(zdd-0.0)
        else
          vref(i,k) = 0.0
        endif
        if(k.eq.1) print *,i,rref(i),vref(i,k)
      enddo
    enddo

    print *,'  Iterate:'
    DO nloop=1,20

      ! get qv and thv from rh and th:
      do k=1,kref
      do i=1,nref
        tx = (pi0(k)+piref(i,k))*(th0(k)+thref(i,k))
        px = p0*((pi0(k)+piref(i,k))**(cp/r_d))
        E1 = 1000.0*SVP1*EXP(SVP2*(tx-SVPT0)/(tx-SVP3))
        qvs = EP_2*E1/(px-E1)
        qvref(i,k) = rh0(k)*qvs
        thvref(i,k)=(th0(k)+thref(i,k))*(1.0+(r_v/r_d)*qvref(i,k))   &
                                           /(1.0+qvref(i,k))
      enddo
      enddo

      ! get nondimensional pressure perturbation (piref):
      do k=1,kref
        piref(nref,k)=0.0
        do i=nref,2,-1
          piref(i-1,k) = piref(i,k)                                       &
       + (rref(i-1)-rref(i))/(cp*0.5*(thvref(i-1,k)+thvref(i,k))) * 0.5 * &
           ( vref(i  ,k)*vref(i  ,k)/rref(i)                              &
            +vref(i-1,k)*vref(i-1,k)/rref(i-1)                            &
             + fcor * ( vref(i,k) + vref(i-1,k) ) )
        enddo
      enddo

      do i=1,nref
        piref(i,   0) = piref(i, 1)
        piref(i,kref+1) = piref(i,kref)
      enddo

      ! get potential temperature perturbation (thref):
      do k=2,kref
      do i=1,nref
        thref(i,k) = 0.5*( cp*0.5*(thvref(i,k)+thvref(i,k+1))*(piref(i,k+1)-piref(i,k))/(zref(k+1)-zref(k))     &
                          +cp*0.5*(thvref(i,k)+thvref(i,k-1))*(piref(i,k)-piref(i,k-1))/(zref(k)-zref(k-1)) )   &
                        *thv0(k)/g
        thref(i,k)=(thv0(k)+thref(i,k))*(1.0+qvref(i,k))/(1.0+(r_v/r_d)*qvref(i,k))-th0(k)
      enddo
      enddo

      k=1
      do i=1,nref
        thref(i,k) = ( cp*0.5*(thvref(i,k)+thvref(i,k+1))*(piref(i,k+1)-piref(i,k))/(zref(k+1)-zref(k)) )   &
                        *thv0(k)/g
        thref(i,k)=(thv0(k)+thref(i,k))*(1.0+qvref(i,k))/(1.0+(r_v/r_d)*qvref(i,k))-th0(k)
      enddo

      print *,'  th,qv,pi = ',nloop,thref(1,1),qvref(1,1),piref(1,1)

    ENDDO   ! enddo for iteration

    ! reference (total) pressure:
    do k=1,kref
    do i=1,nref
      pref(i,k) = p0*( ( pi0(k)+piref(i,k) )**(cp/r_d) )
    enddo
    enddo

    ! analytic axisymmetric vortex is ready ... now interpolate to 3D grid:
    ! (note:  vortex is placed in center of domain)

    ric = float(ide-ids+1)/2.0
    rjc = float(jde-jds+1)/2.0

    print *,'  ids,ide,jds,jde = ',ids,ide,jds,jde
    print *,'  ric,rjc = ',ric,rjc

    print *,'  zk:'
    do k=1,kte
      zk(k) = zref(k)
      print *,k,zk(k)
    enddo

    nl_in = kte-1
    print *,'  nl_in = ',nl_in

    DO J = jts, min(jde-1,jte)
    DO I = its, min(ide-1,ite)
      rr = sqrt( ( (float(i)-ric)*config_flags%dx )**2 + ( (float(j)-rjc)*config_flags%dy )**2 )
      rr = min( rr , rref(nref) )
            diff = -1.0e20
            ii = 0
            do while( diff.lt.0.0 )
              ii = ii + 1
              diff = rref(ii)-rr
            enddo
            i2 = max( ii , 2 )
            i1 = i2-1
            frac = (      rr-rref(i1))   &
                  /(rref(i2)-rref(i1))
            do k=1,kte
              px = p0*( ( pi0(k)+piref(i1,k)+(piref(i2,k)-piref(i1,k))*frac )**(cp/r_d) )
              qx = qvref(i1,k)+(qvref(i2,k)-qvref(i1,k))*frac
              qv(k) = qx
              theta(k) = th0(k)+thref(i1,k)+(thref(i2,k)-thref(i1,k))*frac
              pd_in(k) = px/(1.0+((r_v/r_d)*qx))
            enddo

!  At this point grid%p_top is already set. find the DRY mass in the column
!  by interpolating the DRY pressure. 

   pd_surf = interp_0( pd_in, zk, grid%phb(i,1,j)/g, nl_in )

!  compute the perturbation mass and the full mass

    grid%mu_1(i,j) = pd_surf-grid%p_top - grid%mub(i,j)
    grid%mu_2(i,j) = grid%mu_1(i,j)
    grid%mu0(i,j) = grid%mu_1(i,j) + grid%mub(i,j)

! given the dry pressure and coordinate system, interp the potential
! temperature and qv

    do k=1,kde-1

      p_level = grid%znu(k)*(pd_surf - grid%p_top) + grid%p_top

      moist(i,k,j,P_QV) = interp_0( qv, pd_in, p_level, nl_in )
      grid%t_1(i,k,j)          = interp_0( theta, pd_in, p_level, nl_in ) - t0
      grid%t_2(i,k,j)          = grid%t_1(i,k,j)
     

    enddo



!  integrate the hydrostatic equation (from the RHS of the bigstep
!  vertical momentum equation) down from the top to get grid%p.
!  first from the top of the model to the top pressure

    k = kte-1  ! top level

    qvf1 = 0.5*(moist(i,k,j,P_QV)+moist(i,k,j,P_QV))
    qvf2 = 1./(1.+qvf1)
    qvf1 = qvf1*qvf2

!    grid%p(i,k,j) = - 0.5*grid%mu_1(i,j)/grid%rdnw(k)
    grid%p(i,k,j) = - 0.5*(grid%mu_1(i,j)+qvf1*grid%mub(i,j))/grid%rdnw(k)/qvf2
    qvf = 1. + rvovrd*moist(i,k,j,P_QV)
    grid%alt(i,k,j) = (r_d/p1000mb)*(grid%t_1(i,k,j)+t0)*qvf* &
                (((grid%p(i,k,j)+grid%pb(i,k,j))/p1000mb)**cvpm)
    grid%al(i,k,j) = grid%alt(i,k,j) - grid%alb(i,k,j)

!  down the column

    do k=kte-2,1,-1
      qvf1 = 0.5*(moist(i,k,j,P_QV)+moist(i,k+1,j,P_QV))
      qvf2 = 1./(1.+qvf1)
      qvf1 = qvf1*qvf2
      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)
      qvf = 1. + rvovrd*moist(i,k,j,P_QV)
      grid%alt(i,k,j) = (r_d/p1000mb)*(grid%t_1(i,k,j)+t0)*qvf* &
                  (((grid%p(i,k,j)+grid%pb(i,k,j))/p1000mb)**cvpm)
      grid%al(i,k,j) = grid%alt(i,k,j) - grid%alb(i,k,j)
    enddo

!  this is the hydrostatic equation used in the model after the
!  small timesteps.  In the model, grid%al (inverse density)
!  is computed from the geopotential.


    grid%ph_1(i,1,j) = 0.
    DO k  = 2,kte
      grid%ph_1(i,k,j) = grid%ph_1(i,k-1,j) - (1./grid%rdnw(k-1))*(       &
                   (grid%mub(i,j)+grid%mu_1(i,j))*grid%al(i,k-1,j)+ &
                    grid%mu_1(i,j)*grid%alb(i,k-1,j)  )
                                                   
      grid%ph_2(i,k,j) = grid%ph_1(i,k,j)
      grid%ph0(i,k,j) = grid%ph_1(i,k,j) + grid%phb(i,k,j)
    ENDDO

    ENDDO   ! do loop for i
    ENDDO   ! do loop for j

!-------------------------------------------
! Done with mass fields, now get winds:

! interp v

  DO J = jts, jte
  DO I = its, min(ide-1,ite)

      rr = sqrt( ( (float(i)-ric)*config_flags%dx )**2 + ( (float(j)-0.5-rjc)*config_flags%dy )**2 )
      rr = min( rr , rref(nref) )
            diff = -1.0e20
            ii = 0
            do while( diff.lt.0.0 )
              ii = ii + 1
              diff = rref(ii)-rr
            enddo
            i2 = max( ii , 2 )
            i1 = i2-1
            frac = (      rr-rref(i1))   &
                  /(rref(i2)-rref(i1))
            angle = datan2(dble( (float(j)-0.5-rjc)*config_flags%dy ),   &
                           dble( (float(i)-ric)*config_flags%dx ) )
            do k=1,kte
              v(k) = (vref(i1,k)+( vref(i2,k)- vref(i1,k))*frac )*cos(angle)
              p_in(k) = pref(i1,k)+(pref(i2,k)-pref(i1,k))*frac
            enddo

    IF (j == jds) THEN
      z_at_v = grid%phb(i,1,j)/g
    ELSE IF (j == jde) THEN
      z_at_v = grid%phb(i,1,j-1)/g
    ELSE
      z_at_v = 0.5*(grid%phb(i,1,j)+grid%phb(i,1,j-1))/g
    END IF

    p_surf = interp_0( p_in, zk, z_at_v, nl_in )

    DO K = 1, kte
      p_level = grid%znu(k)*(p_surf - grid%p_top) + grid%p_top
      grid%v_1(i,k,j) = interp_0( v, p_in, p_level, nl_in )
      grid%v_2(i,k,j) = grid%v_1(i,k,j)
    ENDDO

  ENDDO
  ENDDO

! interp u

  DO J = jts, min(jde-1,jte)
  DO I = its, ite

      rr = sqrt( ( (float(i)-ric-0.5)*config_flags%dx )**2 + ( (float(j)-rjc)*config_flags%dy )**2 )
      rr = min( rr , rref(nref) )
            diff = -1.0e20
            ii = 0
            do while( diff.lt.0.0 )
              ii = ii + 1
              diff = rref(ii)-rr
            enddo
            i2 = max( ii , 2 )
            i1 = i2-1
            frac = (      rr-rref(i1))   &
                  /(rref(i2)-rref(i1))
            angle = datan2(dble( (float(j)-rjc)*config_flags%dy ),   &
                           dble( (float(i)-0.5-ric)*config_flags%dx ) )
            do k=1,kte
              u(k) = -(vref(i1,k)+( vref(i2,k)- vref(i1,k))*frac )*sin(angle)
              p_in(k) = pref(i1,k)+(pref(i2,k)-pref(i1,k))*frac
            enddo

    IF (i == ids) THEN
      z_at_u = grid%phb(i,1,j)/g
    ELSE IF (i == ide) THEN
      z_at_u = grid%phb(i-1,1,j)/g
    ELSE
      z_at_u = 0.5*(grid%phb(i,1,j)+grid%phb(i-1,1,j))/g
    END IF

    p_surf = interp_0( p_in, zk, z_at_u, nl_in )

    DO K = 1, kte
      p_level = grid%znu(k)*(p_surf - grid%p_top) + grid%p_top
      grid%u_1(i,k,j) = interp_0( u, p_in, p_level, nl_in )
      grid%u_2(i,k,j) = grid%u_1(i,k,j)
    ENDDO

  ENDDO
  ENDDO

    ! All done ... deallocate arrays:

    deallocate(  rref )
    deallocate(  zref )
    deallocate(   th0 )
    deallocate(   qv0 )
    deallocate(  thv0 )
    deallocate(  prs0 )
    deallocate(   pi0 )
    deallocate(   rh0 )
    deallocate(  vref )
    deallocate( piref )
    deallocate(  pref )
    deallocate( thref )
    deallocate(thvref )
    deallocate( qvref )

    print *,'  completed vortex init successfully '

!-----------------------------------------------------------------------

   if (0.gt.1) then
!#if 0
!  The tropical_cyclone case is adapted from the squall line case
!  so we just turn off the thermal perturbation

!  thermal perturbation to kick off convection
   call random_seed
  write(6,*) ' nxc, nyc for perturbation ',nxc,nyc
  write(6,*) ' delt for perturbation ',delt

  DO J = jts, min(jde-1,jte)
!    yrad = config_flags%dy*float(j-nyc)/4000.
     yrad = 0.
    DO I = its, min(ide-1,ite)
      xrad = config_flags%dx*float(i-nxc)/10000.
!     xrad = 0.
      DO K = 1, 35

!  put in preturbation theta (bubble) and recalc density.  note,
!  the mass in the column is not changing, so when theta changes,
!  we recompute density and geopotential
        zrad = 0.5*(grid%ph_1(i,k,j)+grid%ph_1(i,k+1,j)  &
                   +grid%phb(i,k,j)+grid%phb(i,k+1,j))/g
        zrad = (zrad-1500.)/1500.
        RAD=SQRT(xrad*xrad+yrad*yrad+zrad*zrad)
!        IF(RAD <= 1.) THEN
   call RANDOM_NUMBER(rnd)
     grid%t_1(i,k,j)=grid%t_1(i,k,j)+delt*(rnd-0.5)
         !  grid%t_1(i,k,j)=grid%t_1(i,k,j)+delt*COS(.5*PI*RAD)**2
           grid%t_2(i,k,j)=grid%t_1(i,k,j)
           qvf = 1. + rvovrd*moist(i,k,j,P_QV)
           grid%alt(i,k,j) = (r_d/p1000mb)*(grid%t_1(i,k,j)+t0)*qvf* &
                        (((grid%p(i,k,j)+grid%pb(i,k,j))/p1000mb)**cvpm)
           grid%al(i,k,j) = grid%alt(i,k,j) - grid%alb(i,k,j)
!        ENDIF
      ENDDO

!  rebalance hydrostatically

      DO k  = 2,kte
        grid%ph_1(i,k,j) = grid%ph_1(i,k-1,j) - (1./grid%rdnw(k-1))*(       &
                     (grid%mub(i,j)+grid%mu_1(i,j))*grid%al(i,k-1,j)+ &
                      grid%mu_1(i,j)*grid%alb(i,k-1,j)  )
                                                   
        grid%ph_2(i,k,j) = grid%ph_1(i,k,j)
        grid%ph0(i,k,j) = grid%ph_1(i,k,j) + grid%phb(i,k,j)
      ENDDO

    ENDDO
  ENDDO
   endif
!#endif

   write(6,*) ' grid%mu_1 from comp ', grid%mu_1(1,1)
   write(6,*) ' full state sounding from comp, ph, grid%p, grid%al, grid%t_1, qv '
   do k=1,kde-1
     write(6,'(i3,1x,5(1x,1pe10.3))') k, grid%ph_1(1,k,1)+grid%phb(1,k,1), &
                                      grid%p(1,k,1)+grid%pb(1,k,1), grid%alt(1,k,1), &
                                      grid%t_1(1,k,1)+t0, moist(1,k,1,P_QV)
   enddo

   write(6,*) ' pert state sounding from comp, grid%ph_1, pp, alp, grid%t_1, qv '
   do k=1,kde-1
     write(6,'(i3,1x,5(1x,1pe10.3))') k, grid%ph_1(1,k,1), &
                                      grid%p(1,k,1), grid%al(1,k,1), &
                                      grid%t_1(1,k,1), moist(1,k,1,P_QV)
   enddo

!! interp v
!
!  DO J = jts, jte
!  DO I = its, min(ide-1,ite)
!
!    IF (j == jds) THEN
!      z_at_v = grid%phb(i,1,j)/g
!    ELSE IF (j == jde) THEN
!      z_at_v = grid%phb(i,1,j-1)/g
!    ELSE
!      z_at_v = 0.5*(grid%phb(i,1,j)+grid%phb(i,1,j-1))/g
!    END IF
!
!    p_surf = interp_0( p_in, zk, z_at_v, nl_in )
!
!    DO K = 1, kte
!      p_level = grid%znu(k)*(p_surf - grid%p_top) + grid%p_top
!      grid%v_1(i,k,j) = interp_0( v, p_in, p_level, nl_in )
!      grid%v_2(i,k,j) = grid%v_1(i,k,j)
!    ENDDO
!
!  ENDDO
!  ENDDO
!
!! interp u
!
!  DO J = jts, min(jde-1,jte)
!  DO I = its, ite
!
!    IF (i == ids) THEN
!      z_at_u = grid%phb(i,1,j)/g
!    ELSE IF (i == ide) THEN
!      z_at_u = grid%phb(i-1,1,j)/g
!    ELSE
!      z_at_u = 0.5*(grid%phb(i,1,j)+grid%phb(i-1,1,j))/g
!    END IF
!
!    p_surf = interp_0( p_in, zk, z_at_u, nl_in )
!
!    DO K = 1, kte
!      p_level = grid%znu(k)*(p_surf - grid%p_top) + grid%p_top
!      grid%u_1(i,k,j) = interp_0( u, p_in, p_level, nl_in )
!      grid%u_2(i,k,j) = grid%u_1(i,k,j)
!    ENDDO
!
!  ENDDO
!  ENDDO

!  set w

  DO J = jts, min(jde-1,jte)
  DO K = kts, kte
  DO I = its, min(ide-1,ite)
    grid%w_1(i,k,j) = 0.
    grid%w_2(i,k,j) = 0.
  ENDDO
  ENDDO
  ENDDO

!  set a few more things

  DO J = jts, min(jde-1,jte)
  DO K = kts, kte-1
  DO I = its, min(ide-1,ite)
    grid%h_diabatic(i,k,j) = 0.
  ENDDO
  ENDDO
  ENDDO

  DO k=1,kte-1
    grid%t_base(k) = grid%t_1(1,k,1)
    grid%qv_base(k) = moist(1,k,1,P_QV)
    grid%u_base(k) = grid%u_1(1,k,1)
    grid%v_base(k) = grid%v_1(1,k,1)
    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
  ENDDO

  DO J = jts, min(jde-1,jte)
  DO I = its, min(ide-1,ite)
     thtmp   = grid%t_2(i,1,j)+t0
     ptmp    = grid%p(i,1,j)+grid%pb(i,1,j)
     temp(1) = thtmp * (ptmp/p1000mb)**rcp
     thtmp   = grid%t_2(i,2,j)+t0
     ptmp    = grid%p(i,2,j)+grid%pb(i,2,j)
     temp(2) = thtmp * (ptmp/p1000mb)**rcp
     thtmp   = grid%t_2(i,3,j)+t0
     ptmp    = grid%p(i,3,j)+grid%pb(i,3,j)
     temp(3) = thtmp * (ptmp/p1000mb)**rcp

!     grid%tsk(I,J)=grid%cf1*temp(1)+grid%cf2*temp(2)+grid%cf3*temp(3)
     ! I don't know why this is here, so I have commented it out:
!!!     grid%tmn(I,J)=grid%tsk(I,J)-0.5
  ENDDO
  ENDDO

  RETURN

 END SUBROUTINE init_domain_rk

   SUBROUTINE init_module_initialize
   END SUBROUTINE init_module_initialize

!---------------------------------------------------------------------

!  test driver for get_sounding
!
!      implicit none
!      integer n
!      parameter(n = 1000)
!      real zk(n),p(n),theta(n),rho(n),u(n),v(n),qv(n),pd(n)
!      logical dry
!      integer nl,k
!
!      dry = .false.
!      dry = .true.
!      call get_sounding( zk, p, pd, theta, rho, u, v, qv, dry, n, nl )
!      write(6,*) ' input levels ',nl
!      write(6,*) ' sounding '
!      write(6,*) '  k  height(m)  press (Pa) pd(Pa) theta (K) den(kg/m^3)  u(m/s)     v(m/s)    qv(g/g) '
!      do k=1,nl
!        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)
!      enddo
!      end
!
!---------------------------------------------------------------------------

      subroutine get_sounding( zk, p, p_dry, theta, rho, &
                               u, v, qv, dry, nl_max, nl_in )
      implicit none

      integer nl_max, nl_in
      real zk(nl_max), p(nl_max), theta(nl_max), rho(nl_max), &
           u(nl_max), v(nl_max), qv(nl_max), p_dry(nl_max)
      logical dry

      integer n
      parameter(n=3000)
      logical debug
      parameter( debug = .true.)

! input sounding data

      real p_surf, th_surf, qv_surf
      real pi_surf, pi(n)
      real h_input(n), th_input(n), qv_input(n), u_input(n), v_input(n)

! diagnostics

      real rho_surf, p_input(n), rho_input(n)
      real pm_input(n)  !  this are for full moist sounding

! local data

      real r
      parameter (r = r_d)
      integer k, it, nl
      real qvf, qvf1, dz

!  first, read the sounding

      call read_sounding( p_surf, th_surf, qv_surf, &
                          h_input, th_input, qv_input, u_input, v_input,n, nl, debug )

      if(dry) then
       do k=1,nl
         qv_input(k) = 0.
       enddo
      endif

      if(debug) write(6,*) ' number of input levels = ',nl

        nl_in = nl
        if(nl_in .gt. nl_max ) then
          write(6,*) ' too many levels for input arrays ',nl_in,nl_max
          call wrf_error_fatal ( ' too many levels for input arrays ' )
        end if

!  compute diagnostics,
!  first, convert qv(g/kg) to qv(g/g)

      do k=1,nl
        qv_input(k) = 0.001*qv_input(k)
      enddo

      p_surf = 100.*p_surf  ! convert to pascals
      qvf = 1. + rvovrd*qv_input(1)
      rho_surf = 1./((r/p1000mb)*th_surf*qvf*((p_surf/p1000mb)**cvpm))
      pi_surf = (p_surf/p1000mb)**(r/cp)

      if(debug) then
        write(6,*) ' surface density is ',rho_surf
        write(6,*) ' surface pi is      ',pi_surf
      end if


!  integrate moist sounding hydrostatically, starting from the
!  specified surface pressure
!  -> first, integrate from surface to lowest level

          qvf = 1. + rvovrd*qv_input(1)
          qvf1 = 1. + qv_input(1)
          rho_input(1) = rho_surf
          dz = h_input(1)
          do it=1,10
            pm_input(1) = p_surf &
                    - 0.5*dz*(rho_surf+rho_input(1))*g*qvf1
            rho_input(1) = 1./((r/p1000mb)*th_input(1)*qvf*((pm_input(1)/p1000mb)**cvpm))
          enddo

! integrate up the column

          do k=2,nl
            rho_input(k) = rho_input(k-1)
            dz = h_input(k)-h_input(k-1)
            qvf1 = 0.5*(2.+(qv_input(k-1)+qv_input(k)))
            qvf = 1. + rvovrd*qv_input(k)   ! qv is in g/kg here
 
            do it=1,10
              pm_input(k) = pm_input(k-1) &
                      - 0.5*dz*(rho_input(k)+rho_input(k-1))*g*qvf1
              rho_input(k) = 1./((r/p1000mb)*th_input(k)*qvf*((pm_input(k)/p1000mb)**cvpm))
            enddo
          enddo

!  we have the moist sounding

!  next, compute the dry sounding using p at the highest level from the
!  moist sounding and integrating down.

        p_input(nl) = pm_input(nl)

          do k=nl-1,1,-1
            dz = h_input(k+1)-h_input(k)
            p_input(k) = p_input(k+1) + 0.5*dz*(rho_input(k)+rho_input(k+1))*g
          enddo


        do k=1,nl
          zk(k) = h_input(k)
          p(k) = pm_input(k)
          p_dry(k) = p_input(k)
          theta(k) = th_input(k)
          rho(k) = rho_input(k)
          u(k) = u_input(k)
          v(k) = v_input(k)
          qv(k) = qv_input(k)

        enddo

     if(debug) then
      write(6,*) ' sounding '
      write(6,*) '  k  height(m)  press (Pa) pd(Pa) theta (K) den(kg/m^3)  u(m/s)     v(m/s)    qv(g/g) '
      do k=1,nl
        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)
      enddo

     end if

      end subroutine get_sounding

!-------------------------------------------------------

      subroutine read_sounding( ps,ts,qvs,h,th,qv,u,v,n,nl,debug )
      implicit none
      integer n,nl
      real ps,ts,qvs,h(n),th(n),qv(n),u(n),v(n)
      logical end_of_file
      logical debug

      integer k

      open(unit=10,file='input_sounding',form='formatted',status='old')
      rewind(10)
      read(10,*) ps, ts, qvs
      if(debug) then
        write(6,*) ' input sounding surface parameters '
        write(6,*) ' surface pressure (mb) ',ps
        write(6,*) ' surface pot. temp (K) ',ts
        write(6,*) ' surface mixing ratio (g/kg) ',qvs
      end if

      end_of_file = .false.
      k = 0

      do while (.not. end_of_file)

        read(10,*,end=100) h(k+1), th(k+1), qv(k+1), u(k+1), v(k+1)
        k = k+1
        if(debug) write(6,'(1x,i3,5(1x,e10.3))') k, h(k), th(k), qv(k), u(k), v(k)
        go to 110
 100    end_of_file = .true.
 110    continue
      enddo

      nl = k

      close(unit=10,status = 'keep')

      end subroutine read_sounding

END MODULE module_initialize_ideal




Regards
Peter
Seishou
 
Posts: 1
Joined: Mon Jul 07, 2014 2:04 pm


Return to Idealized Simulations

Who is online

Users browsing this forum: No registered users and 1 guest