!WRF:DRIVER_LAYER:TOP

MODULE atm_comp_mct
!<DESCRIPTION>
! This module defines top-level atm_init_mct(), atm_run_mct(), and atm_final_mct()
!
! Juanxing He, 12 Apr 2009: WRF 3.0
! Juanxing He, 04 Dec 2009
! Juanxing He, 10 May 2010: WRF 3.2
! Matt Higgins, 13 Jul 2010: fixed incorrect precipitation coupling
! Matt Higgins, 17 Aug 2010: fixed snowh units to be consistent with WRF registry
! Matt Higgins, 23 Feb 2011: added WRF wind rotation on fields sent to coupler
! Matt Higgins, 04 Mar 2011: complete review of coupling code
! Matt Higgins, 06 Jun 2011: fixed double-counting of precip export
! Matt Higgins, 24 Jun 2011: averaged values, rather than instantaneous values, now exported

!</DESCRIPTION>

#ifdef CCSMCOUPLED
#ifdef SEQ_MCT

  use mct_mod
  use esmf_mod, CCSM_Clock=>ESMF_Clock
  use seq_flds_mod
  use seq_flds_indices
  use seq_cdata_mod
  use seq_infodata_mod
  use seq_timemgr_mod
  use seq_comm_mct, only:seq_comm_iamroot,seq_comm_iamin
  use shr_file_mod     , only: shr_file_getunit, shr_file_freeunit, &
                               shr_file_setLogUnit, shr_file_setLogLevel, &
                               shr_file_getLogUnit, shr_file_getLogLevel, &
		               shr_file_setIO
  use shr_sys_mod      , only: shr_sys_flush, shr_sys_abort

#endif
#endif   
! for WRF
	
   USE module_machine
   USE module_domain
   USE module_domain_type
   USE module_integrate
   USE module_driver_constants
   USE module_configure
   USE module_io_domain
   USE module_io
   USE module_io_wrf
   USE module_streams, only : RESTART_NOW
   USE module_timing
   USE module_utility  ! for using ESMF function in WRF
   USE module_wrf_error

   USE module_atm_communicator, only: mpi_communicator_atm
#ifdef DM_PARALLEL
   USE module_dm
   USE module_comm_dm, ONLY : halo_em_phys_a_sub
#endif

implicit none
!
! coupler
!
#ifdef CCSMCOUPLED
#ifdef SEQ_MCT

!--------------------------------------------------------------------------
! Private data
!--------------------------------------------------------------------------

  type(mct_aVect)   :: a2x_a_SNAP
  type(mct_aVect)   :: a2x_a_SUM
  
  integer, parameter :: r8 = 8

! Time averaged counter for flux fields
  integer :: avg_count

! Time averaged flux fields
  character(*), parameter :: a2x_avg_flds = "Faxa_rainc:Faxa_rainl:Faxa_snowc:Faxa_snowl:Sa_u:Sa_v:Sa_pslv:&
                                             Sa_z:Sa_pbot:Sa_ptem:Sa_tbot:Sa_dens:Sa_shum:Faxa_swnet:Faxa_lwdn:&
                                             Faxa_swndr:Faxa_swvdr:Faxa_swndf:Faxa_swvdf"
!
! Are all surface types present   
!
  logical :: lnd_present ! if true => land is present
  logical :: ocn_present ! if true => ocean is present
     
!  integer :: iulog  ! for I/O
  logical :: exists
!  integer  :: shrlogunit,shrloglev
  integer :: ATMID

  type coupling_state_variables
    real,dimension(:,:),allocatable::tsk,rasdir,rasdif,raldir,raldif
    real::atm_cpl_dt
  end type coupling_state_variables  
 real,   dimension(0:1000 ) :: PSIMTB,PSIHTB

  type (coupling_state_variables) :: ragrid 

  character(len=*), private, parameter :: &
     seq_timemgr_noLeap            = "NO_LEAP"   , &
     seq_timemgr_gregorian         = "GREGORIAN" 

#endif
#endif

!
!  wrf
!
	
   REAL    :: time

   INTEGER :: loop , &
              levels_to_process

   TYPE (domain) , POINTER :: keep_grid, grid_ptr, null_domain
   TYPE (grid_config_rec_type), SAVE :: config_flags
   INTEGER                 :: number_at_same_level
   INTEGER                 :: time_step_begin_restart

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

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

   CHARACTER ( 80)      :: rstname
   CHARACTER ( 80)      :: message

   INTERFACE
     SUBROUTINE Setup_Timekeeping( grid )
      USE module_domain
      TYPE(domain), POINTER :: grid
     END SUBROUTINE Setup_Timekeeping

#if (EM_CORE == 1)
     SUBROUTINE wrf_dfi_write_initialized_state( )
     END SUBROUTINE wrf_dfi_write_initialized_state

     SUBROUTINE wrf_dfi_bck_init( )
     END SUBROUTINE wrf_dfi_bck_init

     SUBROUTINE wrf_dfi_fwd_init( )
     END SUBROUTINE wrf_dfi_fwd_init

     SUBROUTINE wrf_dfi_fst_init( )
     END SUBROUTINE wrf_dfi_fst_init

     SUBROUTINE wrf_dfi_array_reset ( )
     END SUBROUTINE wrf_dfi_array_reset
#endif
   END INTERFACE


CONTAINS


   SUBROUTINE atm_init_mct( &
#ifdef CCSMCOUPLED
#ifdef SEQ_MCT  
   EClock, cdata_a, x2a_a, a2x_a,NLFilename,&
#endif
#endif
   no_init1 )
!<DESCRIPTION>
!     WRF initialization routine.
!</DESCRIPTION>
#ifdef CCSMCOUPLED
#ifdef SEQ_MCT  	
    type(CCSM_Clock), intent(in)                 :: EClock
    type(seq_cdata), intent(inout)              :: cdata_a
    type(mct_aVect), intent(inout)              :: x2a_a
    type(mct_aVect), intent(inout)              :: a2x_a   
    character( *), optional,   intent(IN)    :: NLFilename ! Namelist filename
#endif
#endif
    LOGICAL, OPTIONAL, INTENT(IN) :: no_init1
!
! Locals
!
#ifdef CCSMCOUPLED
#ifdef SEQ_MCT
     type(mct_gsMap), pointer   :: gsMap_atm
     type(mct_gGrid), pointer   :: dom_a
     type(seq_infodata_type),pointer :: infodata
     integer :: mpicom_atm
     integer :: lsize 
     real(r8):: nextsw_cday      ! calendar of next atm shortwave		 
     real(r8):: dayofYear        ! day of year
     integer :: wrfsec           ! wrf time in seconds
     integer :: radsec           ! radiation time in seconds
     integer :: nradt
     integer :: stepno           ! time step			 
     integer :: dtime_sync       ! integer timestep size
     integer :: currentymd       ! current year-month-day
     integer :: dtime            ! time step increment (sec)
     integer :: atm_cpl_dt       ! driver atm coupling time step 
     integer :: start_ymd        ! Start date (YYYYMMDD)
     integer :: start_tod        ! Start time of day (sec)
     integer :: restart_ymd        ! Start date (YYYYMMDD)
     integer :: restart_tod        ! Start time of day (sec)
     integer :: ref_ymd          ! Reference date (YYYYMMDD)
     integer :: ref_tod          ! Reference time of day (sec)
     integer :: stop_ymd         ! Stop date (YYYYMMDD)
     integer :: stop_tod         ! Stop time of day (sec)
     logical :: perpetual_run    ! If in perpetual mode or not
     integer :: perpetual_ymd    ! Perpetual date (YYYYMMDD)
     real(r8) :: eccen       ! Eccentricity of orbit
     real(r8) :: mvelp       ! Locatn of vernal equinox
     real(r8) :: obliqr      ! Obliquity in radians
     real(r8) :: lambm0      ! Long of perh. radians
     real(r8) :: mvelpp      ! Locatn of vernal equinox at perh. 
     integer :: plon,plat    ! global grid number in the lon and lat direction
     character( 80) :: start_type    ! Start type
     character( 80) :: caseid     ! Short case identification
     character( 80) :: ctitle       ! Long case description 
     logical  :: adiabatic ! atm adiabatic mode
     logical  :: ideal_phys! atm idealized-physics mode
     logical   :: aqua_planet   ! aqua_planet mode
     logical   :: brnch_retain_casename
     logical    :: single_column
     real (R8) :: scmlat
     real (R8)  :: scmlon
     integer :: nsrest 
   
     logical :: first_time = .true.
     character( 80) :: calendar  ! Calendar type
     character( 80) :: starttype ! infodata start type

     TYPE(WRFU_Time) :: current_time
     type(WRFU_TimeInterval) :: off
     integer :: ymd              ! WRF current date (YYYYMMDD)
     integer :: yr               ! WRF current year
     integer :: mon              ! WRF current month
     integer :: day              ! WRF current day
     integer :: hour             ! WRF current hour
     integer :: minuate            ! WRF current minuate
     integer :: second             ! WRF current second
     integer :: tod              ! WRF current time of day (sec)
     logical :: restart
     real,dimension(:),allocatable::value
     
     integer :: ids,ide,jds,jde,kds,kde,ims,ime,jms,jme,kms,kme,ips,ipe,jps,jpe,kps,kpe,&
                  start_year,start_month,start_day,start_hour,start_minute,start_second,&
		  i,j,ig
     real::x,y,zoln
     integer::n

      LOGICAL , EXTERNAL      :: wrf_dm_on_monitor	
      TYPE(domain) , POINTER :: tmpgrid	  ! for the temporary integration
      TYPE(domain) , target :: tmpgrid1	  ! for the temporary integration
      INTEGER :: mpicomcart
      INTEGER, DIMENSION(2) :: dims, coords
#endif
#endif
	
     INTEGER myproc, nproc, hostid, loccomm, ierr, buddcounter, mydevice
     INTEGER, ALLOCATABLE :: hostids(:), budds(:)
     CHARACTER*512 hostname
#if defined(DM_PARALLEL) && defined(RUN_ON_GPU)
     include "mpif.h"
#endif
#include "version_decl"


!<DESCRIPTION>
! Program_name, a global variable defined in frame/module_domain.F, is
! set, then a routine <a href=init_modules.html>init_modules</a> is
! called. This calls all the init programs that are provided by the
! modules that are linked into WRF.  These include initialization of
! external I/O packages.   Also, some key initializations for
! distributed-memory parallelism occur here if DM_PARALLEL is specified
! in the compile: setting up I/O quilt processes to act as I/O servers
! and dividing up MPI communicators among those as well as initializing
! external communication packages such as RSL or RSL_LITE.
!
!</DESCRIPTION>

   program_name = "WRF " // TRIM(release_version) // " MODEL"

#ifdef CCSMCOUPLED
#ifdef SEQ_MCT

! 
! Get data from driver routine
!
   
       call seq_cdata_setptrs(cdata_a, ID=ATMID, mpicom=mpicom_atm,&
                          mpicomcart=mpicomcart,mpicomx=local_communicator_x,&
                          mpicomy=local_communicator_y,ntasks=ntasks,ntasks_x=ntasks_x,ntasks_y=ntasks_y,& 
                          gsMap=gsMap_atm, dom=dom_a, infodata=infodata)

    RESTART_NOW = .false.
    if(first_time) then
	    
       mpi_communicator_atm=mpicom_atm  ! get ATM group communicator

       if (seq_comm_iamroot(ATMID)) then
          inquire(file='atm_modelio.nml',exist=exists)
          if (exists) then
             iulog = shr_file_getUnit()
             call shr_file_setIO('atm_modelio.nml',iulog)
          endif
          write(iulog,*) "WRF initialization"
       endif
       
       call shr_file_getLogUnit (shrlogunit)
       call shr_file_getLogLevel(shrloglev)
       call shr_file_setLogUnit (iulog)

!       write(iulog,*)'********ATM communicator*********',mpi_communicator_atm,ATMID
! 
! Get data from infodata object
!
       call seq_infodata_GetData( infodata,                                           &
            case_name=caseid, case_desc=ctitle,                                       &
            start_type=starttype,                                                     &
            atm_adiabatic=adiabatic,                                                  &
            atm_ideal_phys=ideal_phys,                                                &
            aqua_planet=aqua_planet,                                                  &
            brnch_retain_casename=brnch_retain_casename,                              &
            single_column=single_column, scmlat=scmlat, scmlon=scmlon,                &
            orb_eccen=eccen, orb_mvelpp=mvelpp, orb_lambm0=lambm0, orb_obliqr=obliqr, &
            lnd_present=lnd_present, ocn_present=ocn_present,                         & 
            perpetual=perpetual_run, perpetual_ymd=perpetual_ymd)
    
! Get nsrest from startup type methods
!
       if (     trim(starttype) == trim(seq_infodata_start_type_start)) then
          nsrest = 0
       else if (trim(starttype) == trim(seq_infodata_start_type_cont) ) then
          nsrest = 1
       else if (trim(starttype) == trim(seq_infodata_start_type_brnch)) then
          nsrest = 3
       else
          write(iulog,*) 'wrf_init: ERROR: unknown starttype'
          call shr_sys_abort()
       end if
#endif
#endif
!-----------------------------------------------
   
   ! Initialize WRF modules:
   ! Phase 1 returns after MPI_INIT() (if it is called)
   CALL init_modules(1)
   IF ( .NOT. PRESENT( no_init1 ) ) THEN
     ! Initialize utilities (time manager, etc.)
     call seq_timemgr_EClockGetData(EClock, &
                                    calendar=calendar )
     if (trim(calendar) == trim(seq_timemgr_noLeap)) then
        CALL WRFU_Initialize( defaultCalendar=ESMF_CAL_NOLEAP )
     elseif (trim(calendar) == trim(seq_timemgr_gregorian)) then
        CALL WRFU_Initialize( defaultCalendar=WRFU_CAL_GREGORIAN )
     else
        write(iulog,*) 'wrf_init: ERROR: unknown calendar ',trim(calendar)
        call shr_sys_abort('wrf_init: ERROR: unknown calendar')
     endif

   ENDIF
   ! Phase 2 resumes after MPI_INIT() (if it is called)
   CALL init_modules(2)

!<DESCRIPTION>
! The wrf namelist.input file is read and stored in the USE associated
! structure model_config_rec, defined in frame/module_configure.F, by the
! call to <a href=initial_config.html>initial_config</a>.  On distributed
! memory parallel runs this is done only on one processor, and then
! broadcast as a buffer.  For distributed-memory, the broadcast of the
! configuration information is accomplished by first putting the
! configuration information into a buffer (<a
! href=get_config_as_buffer.html>get_config_as_buffer</a>), broadcasting
! the buffer, then setting the configuration information (<a
! href=set_config_as_buffer.html>set_config_as_buffer</a>).
!
!</DESCRIPTION>

#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 )

! ---------------------------------------------------------------------
! added by Juanxiong He
! --------------------------------------------------------------------
#ifdef CCSMCOUPLED
#ifdef SEQ_MCT

      CALL wrf_set_dm_communicator (mpicomcart)
!      write(iulog,*)'local_communicator=',local_communicator
!      write(iulog,*)'local_communicator_x=',local_communicator_x
!      write(iulog,*)'local_communicator_y=',local_communicator_y
      write(iulog,*)'ntasks_x=',ntasks_x
      write(iulog,*)'ntasks_y=',ntasks_y

      CALL mpi_comm_rank( local_communicator, mytask, ierr )
      CALL mpi_cart_coords( local_communicator, mytask, 2, coords, ierr )

      mytask_x = coords(2)   ! col task (x)
      mytask_y = coords(1)   ! row task (y)
      CALL nl_set_nproc_x ( 1, ntasks_x )
      CALL nl_set_nproc_y ( 1, ntasks_y )

#else
#endif
#else
	
    CALL wrf_dm_initialize
#endif
! ---------------------------------------------------------------------
! added by Juanxiong He
! ---------------------------------------------------------------------



#else
   CALL initial_config
#endif

   CALL set_derived_rconfigs

#ifdef RUN_ON_GPU
   CALL wrf_get_myproc( myproc )
   CALL wrf_get_nproc( nproc )
   CALL wrf_get_hostid ( hostid )
#ifdef DM_PARALLEL
   CALL wrf_get_dm_communicator ( loccomm )
   ALLOCATE( hostids(nproc) )
   ALLOCATE( budds(nproc) )
   CALL mpi_allgather( hostid, 1, MPI_INTEGER, hostids, 1, MPI_INTEGER, loccomm, ierr )
   if ( ierr .NE. 0 ) write(0,*)__FILE__,__LINE__,'error in mpi_allgather ',ierr
   budds = -1
   buddcounter = 0
   ! mark the ones i am on the same node with
   DO i = 1, nproc
      IF ( hostid .EQ. hostids(i) ) THEN
         budds(i) = buddcounter
         buddcounter = buddcounter + 1
      ENDIF
   ENDDO
   mydevice = budds(myproc+1)
   DEALLOCATE( hostids )
   DEALLOCATE( budds )
#else
   mydevice = 0
#endif
   CALL wsm5_gpu_init( myproc, nproc, mydevice )
#endif

!<DESCRIPTION>
! Among the configuration variables read from the namelist is
! debug_level. This is retrieved using nl_get_debug_level (Registry
! generated and defined in frame/module_configure.F).  The value is then
! used to set the debug-print information level for use by <a
! href=wrf_debug.html>wrf_debug</a> throughout the code. Debug_level
! of zero (the default) causes no information to be printed when the
! model runs. The higher the number (up to 1000) the more information is
! printed.
!
!</DESCRIPTION>

   CALL nl_get_debug_level ( 1, debug_level )
   CALL set_wrf_debug_level ( debug_level )

   ! allocated and configure the mother domain

   NULLIFY( null_domain )

!<DESCRIPTION>
! RSL is required for WRF nesting options.
! The non-MPI build that allows nesting is only supported on machines
! with the -DSTUBMPI option.  Check to see if the WRF model is being asked
! for a for a multi-domain run (max_dom > 1, from the namelist).  If so,
! then we check to make sure that we are under the parallel
! run option or we are on an acceptable machine.
!</DESCRIPTION>

   CALL nl_get_max_dom( 1, max_dom )
   IF ( max_dom > 1 ) THEN
#if ( ! defined(DM_PARALLEL)  &&   ! defined(STUBMPI) )
   CALL wrf_error_fatal( &
     'nesting requires either an MPI build or use of the -DSTUBMPI option' )
#endif
   END IF

!<DESCRIPTION>
! The top-most domain in the simulation is then allocated and configured
! by calling <a href=alloc_and_configure_domain.html>alloc_and_configure_domain</a>.
! Here, in the case of this root domain, the routine is passed the
! globally accessible pointer to TYPE(domain), head_grid, defined in
! frame/module_domain.F.  The parent is null and the child index is given
! as negative, signifying none.  Afterwards, because the call to
! alloc_and_configure_domain may modify the model's configuration data
! stored in model_config_rec, the configuration information is again
! repacked into a buffer, broadcast, and unpacked on each task (for
! DM_PARALLEL compiles). The call to <a
! href=setup_timekeeping.html>setup_timekeeping</a> for head_grid relies
! on this configuration information, and it must occur after the second
! broadcast of the configuration information.
!
!</DESCRIPTION>
   CALL       wrf_message ( program_name )
   CALL       wrf_debug ( 100 , 'wrf: calling alloc_and_configure_domain ' )
   CALL alloc_and_configure_domain ( domain_id  = 1 ,                  &
                                     grid       = head_grid ,          &
                                     parent     = null_domain ,        &
                                     kid        = -1                   )

   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 )
   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

#if (EM_CORE == 1)
! In case we are doing digital filter initialization, set dfi_stage = DFI_SETUP
!   to indicate in Setup_Timekeeping that we want forecast start and
!   end times at this point
   IF ( head_grid%dfi_opt .NE. DFI_NODFI ) head_grid%dfi_stage = DFI_SETUP
#endif		      
!
! Initialize wrf time manager
!

   CALL Setup_Timekeeping (head_grid)

!<DESCRIPTION>
! The head grid is initialized with read-in data through the call to <a
! href=med_initialdata_input.html>med_initialdata_input</a>, which is
! passed the pointer head_grid and a locally declared configuration data
! structure, config_flags, that is set by a call to <a
! href=model_to_grid_config_rec.html>model_to_grid_config_rec</a>.  It is
! also necessary that the indices into the 4d tracer arrays such as
! moisture be set with a call to <a
! href=set_scalar_indices_from_config.html>set_scalar_indices_from_config</a>
! prior to the call to initialize the domain.  Both of these calls are
! told which domain they are setting up for by passing in the integer id
! of the head domain as <tt>head_grid%id</tt>, which is 1 for the
! top-most domain.
!
! In the case that write_restart_at_0h is set to true in the namelist,
! the model simply generates a restart file using the just read-in data
! and then shuts down. This is used for ensemble breeding, and is not
! typically enabled.
!
!</DESCRIPTION>

      call get_ijk_from_grid (  head_grid ,                   &
                              ids, ide, jds, jde, kds, kde,    &
                              ims, ime, jms, jme, kms, kme,    &
                              ips, ipe, jps, jpe, kps, kpe    )

      allocate(ragrid%tsk(ids:ide,jds:jde))
      allocate(ragrid%rasdir(ids:ide,jds:jde))
      allocate(ragrid%rasdif(ids:ide,jds:jde))
      allocate(ragrid%raldir(ids:ide,jds:jde))
      allocate(ragrid%raldif(ids:ide,jds:jde))

        do n=0,1000
        zoln=-float(n)*0.01
        x=(1-16.*zoln)**0.25
        psimtb(n)=2*alog(0.5*(1+x))+alog(0.5*(1+x*x))- 2.*atan(x)+2.*atan(1.)
        y=(1-16*zoln)**0.5
        psihtb(n)=2*alog(0.5*(1+y))
        enddo

   CALL med_initialdata_input( head_grid , config_flags )
   
   IF ( config_flags%write_restart_at_0h ) THEN
      CALL med_restart_out ( head_grid, config_flags )
#ifndef AUTODOC_BUILD
! prevent this from showing up before the call to integrate in the autogenerated call tree
      CALL wrf_debug ( 0 , ' 0 h restart only wrf: SUCCESS COMPLETE WRF' )
! TBH:  $$$ Unscramble this later...
! TBH:  $$$ Need to add state to avoid calling wrf_finalize() twice when ESMF
! TBH:  $$$ library is used.  Maybe just set clock stop_time=start_time and
! TBH:  $$$ do0com not call wrf_finalize here...
      CALL atm_final_mct( )
#endif
   END IF

! set default values for subtimes
   head_grid%start_subtime = domain_get_start_time ( head_grid )
   head_grid%stop_subtime = domain_get_stop_time ( head_grid )
      
!
! Check consistency of outside clock and WRF clock
! Check consistency of restart time information with input clock
!	
#ifdef CCSMCOUPLED
#ifdef SEQ_MCT
!
! Initialize time manager (outside clock).
!
       call seq_timemgr_EClockGetData(EClock, &
                                      start_ymd=start_ymd, start_tod=start_tod, &
                                      ref_ymd=ref_ymd, ref_tod=ref_tod,         &
                                      stop_ymd=stop_ymd, stop_tod=stop_tod,     &
                                      calendar=calendar )
				      
        call nl_get_restart(1,restart)
	if (restart) then
	 
	    CALL nl_get_start_year(head_grid%id,start_year)
            CALL nl_get_start_month(head_grid%id,start_month)
            CALL nl_get_start_day(head_grid%id,start_day)
            CALL nl_get_start_hour(head_grid%id,start_hour)
            CALL nl_get_start_minute(head_grid%id,start_minute)
            CALL nl_get_start_second(head_grid%id,start_second)    

     end if
!
! Initialize MCT gsMap, domain and attribute vectors
!
       call atm_SetgsMap_mct( head_grid, mpicom_atm, ATMID, gsMap_atm )
       lsize = mct_gsMap_lsize(gsMap_atm, mpicom_atm)
!
! Initialize MCT domain 
!
       call atm_domain_mct( head_grid, lsize, gsMap_atm, dom_a )
!
! Initialize MCT attribute vectors
!
       call mct_aVect_init(a2x_a, rList=seq_flds_a2x_fields, lsize=lsize)
       call mct_aVect_zero(a2x_a)
       
       call mct_aVect_init(x2a_a, rList=seq_flds_x2a_fields, lsize=lsize) 
       call mct_aVect_zero(x2a_a)
       
       call mct_aVect_init(a2x_a_SNAP, rList=a2x_avg_flds, lsize=lsize)
       call mct_aVect_zero(a2x_a_SNAP)
       
       call mct_aVect_init(a2x_a_SUM , rList=a2x_avg_flds, lsize=lsize)
       call mct_aVect_zero(a2x_a_SUM )
!
! Initialize averaging counter
!
       avg_count = 0
!
! Create initial atm export state
!
       call atm_export_mct( head_grid, a2x_a )
!
! Set flag to specify that an extra albedo calculation is to be done (i.e. specify active)
!		
        
       plon=ide-ids  !for global domain
       plat=jde-jds  !for global domain
           
       call seq_infodata_PutData(infodata, atm_prognostic=.true.)
       call seq_infodata_PutData(infodata, atm_nx=plon, atm_ny=plat)
!
! Set time step of radiation computation as the current calday
! This will only be used on the first timestep of an initial run
!

       write(iulog,*)'plon=',plon,mytask 
       write(iulog,*)'plat=',plat,mytask
       write(iulog,*)'nextsw_cday = ', nextsw_cday
         
       call seq_infodata_PutData( infodata, nextsw_cday=nextsw_cday )
       
       first_time = .false.
       write(iulog,*) '******local comunicator*******',local_communicator
	  
       call shr_file_setLogUnit (shrlogunit)
       call shr_file_setLogLevel(shrloglev)

    else
       
! For initial run, nothing wnated to read
! For restart run, read restart x2a_a

!  obtain the begin date of CCSM again

       call shr_file_getLogUnit (shrlogunit)
       call shr_file_getLogLevel(shrloglev)
       call shr_file_setLogUnit (iulog)
       
!       write(iulog,*) '******local comunicator*******',local_communicator
       call seq_timemgr_EClockGetData(EClock,curr_ymd=CurrentYMD, StepNo=StepNo, dtime=DTime_Sync ) 
       call domain_clock_get( head_grid, current_time=current_time)
       call WRFU_TimeGet( current_time, yy=yr,mm=mon,dd=day)   
       ymd=yr*10000+mon*100+day
       if( currentymd.ne.ymd) then
         write(iulog,*) 'ymd=',ymd,'currentymd=',currentymd
         write(iulog,*)'WRF begintime must be the same with CCSM'
         call shr_sys_abort()
       endif
   
!  determine whether be restart       
        call nl_get_restart(1,restart)
        if (.NOT.restart) then       
              tmpgrid=>tmpgrid1
              tmpgrid=head_grid
          write(iulog,*)'wrf: calling import' 	      
          call atm_import_mct( x2a_a, tmpgrid )
          write(iulog,*)'wrf: calling integrate' 
	  write(iulog,*) associated(tmpgrid,head_grid)
#ifdef CCSMCOUPLED
#ifdef SEQ_MCT
          call integrate ( tmpgrid,.true. )
#else
#endif
#else
          call integrate ( tmpgrid)
#endif
          write(iulog,*)'wrf: back from integrate'
          call atm_export_mct( tmpgrid, a2x_a )
	  nullify(tmpgrid)
          write(iulog,*)'wrf: export finish'           
        else
	  call atm_export_mct( head_grid, a2x_a )
        end if

! compute day of next radiation computation, assume wrf and rad 
! timestep divide a year evenly

       ! need rad freq, wrf timestep, coupling period
       nradt = config_flags%radt*60
       dtime = head_grid%dt
       call seq_timemgr_EClockGetData(Eclock,dtime=atm_cpl_dt)

       ! get current day of year
       call domain_clock_get( head_grid, current_time=current_time)
       call WRFU_TimeGet( current_time, dayOfYear_r8=dayofYear)

       ! compute wrf time of year in seconds
       wrfsec = nint(dayofYear*86400.0_r8)

       ! compute a wrf radiation time (either past or next)
       radsec = int(wrfsec/nradt)*nradt + dtime

       ! if radsec is in past, compute next radsec
       if (wrfsec >= radsec) radsec = radsec + nradt
       if (wrfsec >= radsec) call shr_sys_abort('wrf_init nextsw1')

       ! set nextsw_cday if the next coupling period is a radiation ts
       if ((wrfsec + atm_cpl_dt) >= radsec .or. .not.restart) then
          nextsw_cday = mod(real(wrfsec,r8)/86400.0_r8-1.0_r8,365.0_r8)+1.0_r8
       else
          nextsw_cday = -1._r8
       endif

!       write(iulog,*) 'tcx nextsw1 ',nradt,dtime,atm_cpl_dt
!       write(iulog,*) 'tcx nextsw2 ',dayofYear,wrfsec,radsec,nextsw_cday

       call seq_infodata_PutData( infodata, nextsw_cday=nextsw_cday )

! check timestep and couple interval consistency
       
       dtime = head_grid%dt
!       if (atm_cpl_dt*1.0/dtime-int(atm_cpl_dt/dtime) .ne. 0.0) then
!	 write(iulog,*)atm_cpl_dt*1.0/dtime,atm_cpl_dt,dtime
       if (abs(atm_cpl_dt*1.0/dtime-int(atm_cpl_dt/dtime)) > 1.0e-5) then
         write(iulog,*) atm_cpl_dt*1.0/dtime,atm_cpl_dt,dtime, &
            (atm_cpl_dt*1.0/dtime-int(atm_cpl_dt/dtime)) 
         write(iulog,*)'atm_cpl_dt must be integer times of time step'
         call shr_sys_abort()
       endif
       if (dtime > atm_cpl_dt) then
         write(iulog,*)'timestep must be less than or equal to atm_cpl_dt'
	 call shr_sys_abort()
       end if 

    end if

       call shr_file_setLogUnit (shrlogunit)
       call shr_file_setLogLevel(shrloglev)
     
#endif
#endif
   END SUBROUTINE atm_init_mct

   SUBROUTINE atm_run_mct(&
#ifdef CCSMCOUPLED
#ifdef SEQ_MCT   
   EClock, cdata_a, x2a_a, a2x_a &
#endif
#endif
   )
!<DESCRIPTION>
!     WRF run routine.
!</DESCRIPTION>

!<DESCRIPTION>
! Once the top-level domain has been allocated, configured, and
! initialized, the model time integration is ready to proceed.  The start
! and stop times for the domain are set to the start and stop time of the
! model run, and then <a href=integrate.html>integrate</a> is called to
! advance the domain forward through that specified time interval.  On
! return, the simulation is completed.
!
!</DESCRIPTION>

!  The forecast integration for the most coarse grid is now started.  The
!  integration is from the first step (1) to the last step of the simulation.
!-----------------------------------------------------------------------
!
! Arguments
!
#ifdef CCSMCOUPLED
#ifdef SEQ_MCT	
    type(CCSM_Clock)            ,intent(in)    :: EClock
    type(seq_cdata)             ,intent(inout) :: cdata_a
    type(mct_aVect)             ,intent(inout) :: x2a_a
    type(mct_aVect)             ,intent(inout) :: a2x_a
!
! Local variables
!
    type(seq_infodata_type),pointer :: infodata
    integer :: lsize           ! size of attribute vector
    integer :: DTime_Sync      ! integer timestep size
    integer :: iradsw          ! shortwave radation frequency (time steps)
    logical :: dosend          ! true => send data back to driver
    integer :: dtime           ! time step increment (sec)
    integer :: atm_cpl_dt      ! driver atm coupling time step
    integer :: ymd_sync        ! Sync date (YYYYMMDD)
    integer :: yr_sync         ! Sync current year
    integer :: mon_sync        ! Sync current month
    integer :: day_sync        ! Sync current day
    integer :: tod_sync        ! Sync current time of day (sec)
    real(r8):: nextsw_cday     ! calendar of next atm shortwave
    logical :: rstwr           ! .true. ==> write restart file before returning
    logical :: nlend           ! Flag signaling last time-step
    logical :: rstwr_sync      ! .true. ==> write restart file before returning
    logical :: nlend_sync      ! Flag signaling last time-step
    logical :: first_time = .true.    
    real(r8):: dayofYear        ! day of year
    integer :: wrfsec           ! wrf time in seconds
    integer :: radsec           ! radiation time in seconds
    integer :: nradt
    
    TYPE(WRFU_Time) :: current_time
    type(WRFU_TimeInterval) :: off
    integer :: ymd             ! WRF current date (YYYYMMDD)
    integer :: yr              ! WRF current year
    integer :: mon             ! WRF current month
    integer :: day             ! WRF current day
    integer :: hour             ! WRF current hour
    integer :: minuate            ! WRF current minuate
    integer :: second             ! WRF current second
    integer :: tod             ! WRF current time of day (sec)
    integer :: nstep           ! WRF nstep

    integer :: ids,ide,jds,jde,kds,kde,ims,ime,jms,jme,kms,kme,ips,ipe,jps,jpe,kps,kpe,&
		  i,j,ig
    TYPE(domain) , POINTER :: tmpgrid	  ! for the temporary integration
    TYPE(domain) , target :: tmpgrid1	  ! for the temporary integration
    real(r8)::tv
    real :: st
    
#endif
#endif    
!-----------------------------------------------------------------------

#ifdef CCSMCOUPLED
#ifdef SEQ_MCT


! Note that sync clock time should match wrf time at end of time step/loop not beginning

    call seq_cdata_setptrs(cdata_a, infodata=infodata)
    call seq_timemgr_EClockGetData(EClock,curr_ymd=ymd_sync,curr_tod=tod_sync, &
            curr_yr=yr_sync,curr_mon=mon_sync,curr_day=day_sync)

    nlend_sync = seq_timemgr_StopAlarmIsOn(EClock)
    rstwr_sync = seq_timemgr_RestartAlarmIsOn(EClock)

! Map input from mct to wrf data structure

    call atm_import_mct(x2a_a, head_grid )
    
! Cycle over all time steps in a wrf radiation coupling interval

    dosend = .false.
    do while (.not. dosend)
!        write(iulog,*)  'wrf: calling integrate' 
!       call atm_sfcclay_mct(head_grid)
       call domain_clock_get(head_grid, current_time=current_time)
       call WRFU_timeget(current_time, yy=yr, mm=mon, dd=day, h=hour, m=minuate, s=second)
       ymd = yr*10000 + mon*100 + day
       tod = hour*3600+minuate*60+second
!       write(iulog,*) 'b4 integrate',ymd,tod

! Determine if time to write wrf restart and stop
! Need to check time at end of next timestep

       call domain_clock_get(head_grid, current_time=current_time)
       dtime = head_grid%dt
       call WRFU_TimeIntervalSet( off, s=dtime)
       current_time = current_time + off
       call WRFU_timeget(current_time, yy=yr, mm=mon, dd=day, h=hour, m=minuate, s=second)
       ymd = yr*10000 + mon*100 + day
       tod = hour*3600+minuate*60+second
       dosend = (seq_timemgr_EClockDateInSync( EClock, ymd, tod))

       RESTART_NOW = .false.
       if (rstwr_sync .and. dosend) RESTART_NOW = .true.
       nlend = .false.
       if (nlend_sync .and. dosend) nlend = .true.

!       if (dosend) write(iulog,*) 'b4 integrate dosend ',ymd,tod,dosend,RESTART_NOW

#endif
#endif

#ifdef CCSMCOUPLED
#ifdef SEQ_MCT
       call integrate ( head_grid,.false. )
#else
#endif
#else
       call integrate ( head_grid)
#endif

#ifdef CCSMCOUPLED
#ifdef SEQ_MCT

!      write(iulog,*)  'wrf: back from integrate' 
! Determine if dosend
! When time is not updated at the end of the loop - then return only if
! are in sync with clock before time is updated

       call domain_clock_get(head_grid, current_time=current_time)
       call WRFU_timeget(current_time, yy=yr, mm=mon, dd=day, h=hour, m=minuate, s=second)
       ymd = yr*10000 + mon*100 + day
       tod = hour*3600+minuate*60+second
       dosend = (seq_timemgr_EClockDateInSync( EClock, ymd, tod))
!       if (dosend) write(iulog,*) 'af integrate dosend ',ymd,tod,dosend
       
! Map output from wrf to mct data structures

       call atm_export_mct( head_grid, a2x_a )
       
! Compute snapshot attribute vector for accumulation

! don't accumulate on first coupling freq ts1 and ts2
! for consistency with ccsm3 when flxave is off       
       nstep = head_grid%itimestep        
       if (nstep <= 2) then
          call mct_aVect_copy( a2x_a, a2x_a_SUM )
          avg_count = 1
       else
          call mct_aVect_copy( a2x_a, a2x_a_SNAP )
          call mct_aVect_accum( aVin=a2x_a_SNAP, aVout=a2x_a_SUM )
          avg_count = avg_count + 1
       endif

    end do

! Finish accumulation of attribute vector and average and copy accumulation
! field into output attribute vector

    call mct_aVect_avg ( a2x_a_SUM, avg_count)
    call mct_aVect_copy( a2x_a_SUM, a2x_a )
    call mct_aVect_zero( a2x_a_SUM)
    avg_count = 0

! Get time of next radiation calculation - albedos will need to be
! calculated by each surface model at this time

! compute day of next radiation computation, assume wrf and rad 
! timestep divide a year evenly

       ! need rad freq, wrf timestep, coupling period
       nradt = config_flags%radt*60
       dtime = head_grid%dt
       call seq_timemgr_EClockGetData(Eclock,dtime=atm_cpl_dt)

       ! get current day of year
       call domain_clock_get( head_grid, current_time=current_time)
       call WRFU_TimeGet( current_time, dayOfYear_r8=dayofYear)

       ! compute wrf time of year in seconds
       wrfsec = nint(dayofYear*86400.0_r8)

       ! compute a wrf radiation time (either past or next)
       radsec = int(wrfsec/nradt)*nradt + dtime

       ! if radsec is in past, compute next radsec
       if (wrfsec >= radsec) radsec = radsec + nradt
       if (wrfsec >= radsec) call shr_sys_abort('wrf_init nextsw1')

       ! set nextsw_cday if the next coupling period is a radiation ts
       if ((wrfsec + atm_cpl_dt) >= radsec) then
          nextsw_cday = mod(real(wrfsec,r8)/86400.0_r8-1.0_r8,365.0_r8)+1.0_r8
       else
          nextsw_cday = -1._r8
       endif

!       write(iulog,*) 'tcx nextsw1 ',nradt,dtime,atm_cpl_dt
!       write(iulog,*) 'tcx nextsw2 ',dayofYear,wrfsec,radsec,nextsw_cday

       call seq_infodata_PutData( infodata, nextsw_cday=nextsw_cday )

! Write merged surface data restart file if appropriate

!    if (rstwr_sync) then
!       call atm_write_srfrest_mct( head_grid,  config_flags)
!    end if

! Check for consistency of internal wrf clock with master sync clock

    call domain_clock_get(head_grid, current_time=current_time)    
    call WRFU_timeget(current_time, yy=yr, mm=mon, dd=day, h=hour, m=minuate, s=second)
    ymd = yr*10000 + mon*100 + day
    tod = hour*3600+minuate*60+second
       
    if ( .not. seq_timemgr_EClockDateInSync( EClock, ymd, tod ) )then
       call seq_timemgr_EClockGetData(EClock, curr_ymd=ymd_sync, curr_tod=tod_sync )
       write(iulog,*)  ' wrf ymd=',ymd     ,'  wrf tod= ',tod 
       write(iulog,*)  ' sync ymd=',ymd_sync,' sync tod= ',tod_sync
       call shr_sys_abort()
    end if
#endif
#endif
   END SUBROUTINE atm_run_mct

   SUBROUTINE atm_final_mct( no_shutdown )
!<DESCRIPTION>
!     WRF finalize routine.
!</DESCRIPTION>

!<DESCRIPTION>
! A Mediation Layer-provided
! subroutine, <a href=med_shutdown_io.html>med_shutdown_io</a> is called
! to allow the the model to do any I/O specific cleanup and shutdown, and
! then the WRF Driver Layer routine <a
! href=wrf_shutdown.html>wrf_shutdown</a> (quilt servers would be
! directed to shut down here) is called to properly end the run,
! including shutting down the communications (for example, most comm
! layers would call MPI_FINALIZE at this point if they're using MPI).
!
!</DESCRIPTION>
     LOGICAL, OPTIONAL, INTENT(IN) :: no_shutdown


#if CCSMCOUPLED
#if SEQ_MCT
    deallocate(ragrid%tsk)
    deallocate(ragrid%rasdir)
    deallocate(ragrid%raldir)
    deallocate(ragrid%rasdif)
    deallocate(ragrid%raldif)
#endif
#endif

   ! shut down I/O
   CALL med_shutdown_io ( head_grid , config_flags )
   CALL       wrf_debug ( 100 , 'wrf: back from med_shutdown_io' )

   CALL       wrf_debug (   0 , 'wrf: SUCCESS COMPLETE WRF' )

   ! Call wrf_shutdown() (which calls MPI_FINALIZE()
   ! for DM parallel runs).
   IF ( PRESENT( no_shutdown ) ) THEN
     ! Finalize time manager
#if CCSMCOUPLED
     ! do not shutdown mpi here if coupled in ccsm
#else
     CALL WRFU_Finalize
     CALL wrf_shutdown
#endif
   ENDIF

   END SUBROUTINE atm_final_mct

   SUBROUTINE wrf_dfi()

#if CCSMCOUPLED
#if SEQ_MCT
#else
#endif
#else
   
!<DESCRIPTION>
! Runs a digital filter initialization procedure.
!</DESCRIPTION>
      IMPLICIT NONE

#if (EM_CORE == 1)
      ! Initialization procedure
      IF ( config_flags%dfi_opt .NE. DFI_NODFI ) THEN

         SELECT CASE ( config_flags%dfi_opt )

            CASE (DFI_DFL)
               wrf_err_message = 'Initializing with DFL'
               CALL wrf_message(TRIM(wrf_err_message))

               wrf_err_message = '   Filtering forward in time'
               CALL wrf_message(TRIM(wrf_err_message))

               CALL wrf_dfi_fwd_init()
               CALL atm_run_mct()

               CALL wrf_dfi_array_reset()

               CALL wrf_dfi_fst_init()

               IF ( config_flags%dfi_write_filtered_input ) THEN
                  CALL wrf_dfi_write_initialized_state()
               END IF

            CASE (DFI_DDFI)
               wrf_err_message = 'Initializing with DDFI'
               CALL wrf_message(TRIM(wrf_err_message))

               wrf_err_message = '   Integrating backward in time'
               CALL wrf_message(TRIM(wrf_err_message))

               CALL wrf_dfi_bck_init()
               CALL atm_run_mct()

               wrf_err_message = '   Filtering forward in time'
               CALL wrf_message(TRIM(wrf_err_message))

               CALL wrf_dfi_fwd_init()
               CALL atm_run_mct()

               CALL wrf_dfi_array_reset()

               CALL wrf_dfi_fst_init()

               IF ( config_flags%dfi_write_filtered_input ) THEN
                  CALL wrf_dfi_write_initialized_state()
               END IF

            CASE (DFI_TDFI)
               wrf_err_message = 'Initializing with TDFI'
               CALL wrf_message(TRIM(wrf_err_message))

               wrf_err_message = '   Integrating backward in time'
               CALL wrf_message(TRIM(wrf_err_message))

               CALL wrf_dfi_bck_init()
               CALL atm_run_mct()

               CALL wrf_dfi_array_reset()

               wrf_err_message = '   Filtering forward in time'
               CALL wrf_message(TRIM(wrf_err_message))

               CALL wrf_dfi_fwd_init()
               CALL atm_run_mct()

               CALL wrf_dfi_array_reset()

               CALL wrf_dfi_fst_init()

               IF ( config_flags%dfi_write_filtered_input ) THEN
                  CALL wrf_dfi_write_initialized_state()
               END IF

            CASE DEFAULT
               wrf_err_message = 'Unrecognized DFI_OPT in namelist'
               CALL wrf_error_fatal(TRIM(wrf_err_message))

         END SELECT

      END IF
#endif
#endif
   END SUBROUTINE wrf_dfi

   SUBROUTINE set_derived_rconfigs
!<DESCRIPTION>
! Some derived rconfig entries need to be set based on the value of other,
! non-derived entries before package-dependent memory allocation takes place.
! This might be employed when, for example, we want to allocate arrays in
! a package that depends on the setting of two or more namelist variables.
! In this subroutine, we do just that.
!</DESCRIPTION>

      IMPLICIT NONE

      INTEGER :: i


#if (EM_CORE == 1)
      IF ( model_config_rec % dfi_opt .EQ. DFI_NODFI ) THEN
        DO i = 1, model_config_rec % max_dom
           model_config_rec % mp_physics_dfi(i) = -1
        ENDDO
      ELSE
        DO i = 1, model_config_rec % max_dom
           model_config_rec % mp_physics_dfi(i) = model_config_rec % mp_physics(i)
        ENDDO
      END IF
#endif

#if (DA_CORE == 1)
      IF ( model_config_rec % dyn_opt .EQ. 2 ) THEN
        DO i = 1, model_config_rec % max_dom
           model_config_rec % mp_physics_4dvar(i) = -1
        ENDDO
      ELSE
        DO i = 1, model_config_rec % max_dom
           model_config_rec % mp_physics_4dvar(i) = model_config_rec % mp_physics(i)
        ENDDO
      END IF
#endif

   END SUBROUTINE set_derived_rconfigs
   
!================================================================================
#ifdef CCSMCOUPLED
#ifdef SEQ_MCT

    subroutine atm_SetgsMap_mct( grid, mpicom_atm, ATMID, GSMap_atm )

!-------------------------------------------------------------------
!
! Arguments
!
    TYPE(domain) , POINTER :: grid
    integer        , intent(in)  :: mpicom_atm
    integer        , intent(in)  :: ATMID
    type(mct_gsMap), intent(out) :: GSMap_atm
!
! Local variables
!
    integer, allocatable :: gindex(:)
    integer :: i, j, k, n,lsize,gsize
    integer :: ids,ide,jds,jde,kds,kde,ims,ime,jms,jme,kms,kme,ips,ipe,jps,jpe,kps,kpe
    integer :: ier            ! error status
!-------------------------------------------------------------------
! Build the atmosphere grid numbering for MCT
! NOTE:  Numbering scheme is: West to East and South to North
! starting at south pole.  Should be the same as what's used in SCRIP
! Determine global seg map
    call get_ijk_from_grid (  grid ,                   &
                              ids, ide, jds, jde, kds, kde,    &
                              ims, ime, jms, jme, kms, kme,    &
                              ips, ipe, jps, jpe, kps, kpe    )
    if(jpe.eq.jde) jpe=jpe-1
    if(ipe.eq.ide) ipe=ipe-1
			      
    lsize=0
    do j=jps, jpe             
       do i=ips,ipe
             lsize = lsize+1  !local index
       end do
    end do

!    gsize=(ide-ids+1)*(jde-jds+1)
    gsize=(ide-ids)*(jde-jds)
    allocate(gindex(lsize))

    n=0
    do j=jps, jpe             
       do i=ips,ipe
          n=n+1
!          gindex(n) = (j-1)*(ide-ids+1)+i  ! global index
          gindex(n) = (j-1)*(ide-ids)+i  ! global index
       end do
    end do

    call mct_gsMap_init( gsMap_atm, gindex, mpicom_atm, ATMID, lsize, gsize)

    deallocate(gindex)

  end subroutine atm_SetgsMap_mct

!===============================================================================

  subroutine atm_import_mct( x2a_a, grid )

!-----------------------------------------------------------------------
!
! Arguments
!
    TYPE(domain) , POINTER:: grid
    type(mct_aVect)   , intent(inout) :: x2a_a

!
! Local variables
!		
    integer  :: mype            ! processor ID
    integer  :: m             ! indices
    integer  :: i,j,k,ig  ! indices    
    integer :: ids,ide,jds,jde,kds,kde,ims,ime,jms,jme,kms,kme,ips,ipe,jps,jpe,kps,kpe    
    logical :: intp
    logical,save :: first_time=.true.
    real :: tau, taux, tauy
!-----------------------------------------------------------------------
!
       call shr_file_getLogUnit (shrlogunit)
       call shr_file_getLogLevel(shrloglev)
       call shr_file_setLogUnit (iulog)
       
    call get_ijk_from_grid (  grid ,                   &
                              ids, ide, jds, jde, kds, kde,    &
                              ims, ime, jms, jme, kms, kme,    &
                              ips, ipe, jps, jpe, kps, kpe    )

!    write(iulog,*)local_communicator
      
    if(jpe.eq.jde) jpe=jpe-1
    if(ipe.eq.ide) ipe=ipe-1

!     write(iulog,*)ips,ipe,jps,jpe,kps,kpe
    ig=1
    do j=jps,jpe                                                        
       do i =ips,ipe                                                                   

! WRF expects latent heat, sensible heat, upward longwave, and upward moisture flux
! to be positive upward (reversed from what coupler expects). CPL7 expects fluxes to
! be positive downards. Therefore, must multiply by -1 for these fluxes:
!
! surface latent heat flux in W m-2 mapped to -1 * WRF's latent heat flux at surface in W m-2
          grid%lh(i,j)       = -x2a_a%rAttr(index_x2a_Faxx_lat, ig)

! surface sensible heat flux in W m-2 mapped to -1 * WRF's upwards heat flux at the surface in W m-2
          grid%hfx(i,j)       = -x2a_a%rAttr(index_x2a_Faxx_sen, ig)

! surface upwards longwave heat flux in W m-2 mapped to -1 * WRF's surface upward longwave 
! radiation in W m-2. lwups is not a field normally in WRF's Registry. Instead,
! this has been added to the WRF Registry and the WRF radiation module has been 
! modified to correctly utilize lwups.
          grid%lwups(i,j)      = -x2a_a%rAttr(index_x2a_Faxx_lwup,ig)

! evaporation water flux in kg m-2 s-1 mapped to -1 * WRF's upward moisture flux at the surface
! in kg m-2 s-1
          grid%qfx(i,j)    = -x2a_a%rAttr(index_x2a_Faxx_evap,ig)

! The following four albedo terms are not part of stand-alone WRF and have been
! added to the WRF Regsitry. The WRF radiation module has been modified to correctly
! utilize the terms.
! direct albedo (visible radiation); unitless
          grid%rasdir(i,j)     =  x2a_a%rAttr(index_x2a_Sx_avsdr, ig)

! direct albedo (near-infrared radiation); unitless
          grid%raldir(i,j)     =  x2a_a%rAttr(index_x2a_Sx_anidr, ig)

! diffuse albedo (visible radiation); unitless 
          grid%rasdif(i,j)     =  x2a_a%rAttr(index_x2a_Sx_avsdf, ig)

! diffuse albedo (near-infrared radiation); unitless
          grid%raldif(i,j)     =  x2a_a%rAttr(index_x2a_Sx_anidf, ig)

! snow water equivalent over land imported from coupler in m and immediately 
! converted to kg/m2 as WRF expects snow water equiv to be in kg/m2. Note that 
! the coupler field description implies this field is physical snow depth. However,
! according to Tony Craig, it is snow water equivalent.
          grid%snow(i,j) =  1000.*x2a_a%rAttr(index_x2a_Sl_snowh, ig)

! Note that, despite being a coupler field, x2a_si_snowh (snow on ice) is NOT 
! used by WRF (or CAM). Therefore, that field is not imported.
     
! reference temperature at 2 m in K mapped to WRF's temperature at 2 m in K
          grid%t2(i,j)      =  x2a_a%rAttr(index_x2a_Sx_tref,  ig)

! reference specific humidity at 2 m in kg/kg mapped to WRF's specific humidity at 2 m in kg/kg
          grid%q2(i,j)      =  x2a_a%rAttr(index_x2a_Sx_qref,  ig)

! surface ice fraction mapped to WRF's ice fraction
          grid%xice(i,j)   =  x2a_a%rAttr(index_x2a_Sx_ifrac, ig)

! surface land fraction mapped to WRF's land mask. In WRF, 1 is land and 2 is water
          grid%xland(i,j)  =  2-x2a_a%rAttr(index_x2a_Sx_lfrac, ig)


! zonal and meridional surface stress in N/m2
! WRF does not need these fields (x2a_Faxx_taux and x2a_Faxx_tauy) as long as it has u*
! Therefore, these fields are not imported. Previous versions of RACM imported these 
! stresses for the sole purpose of calculating u*.

! friction velocity (u*) in m/s
! u* is calculated from input values of wind stress, rather than the input values of u*. This
! is because while both the land model and the ocean model output u*, the ice model does
! not. Wind stresses, however, are combined by the coupler from all three components. Perhaps
! a future version of CPL7 will output a combined u* or at least a u* from the ice model. 
! Until then, calculate u* here for use by WRF:
!
! a) get incoming wind stresses:
          taux = -x2a_a%rAttr(index_x2a_Faxx_taux,ig)     
          tauy = -x2a_a%rAttr(index_x2a_Faxx_tauy,ig) 

! b) calculate stress from component vectors. Note: we do not rotate incoming vectors to
! WRF grid because all we are interested in is the magnitude of the stress
          tau = sqrt(taux*taux + tauy*tauy)

! c) calculate u* as sqrt (tau / rho)
! rho is taken as density at lowest model level
          grid%ust(i,j) = sqrt(tau / (1.0/grid%alt(i,1,j)*(1+grid%moist(i,1,j,P_QV))))

! surface temperature mapped to WRF's surface skin temperature, both in K
          grid%tsk(i,j)       =  x2a_a%rAttr(index_x2a_Sx_t,ig)

! sea surface temperature mapped to WRF's sea surface skin temperature, both in K
          grid%sst(i,j)       =  x2a_a%rAttr(index_x2a_So_t,ig)

! put WRF's imported sea surface temperature into skin temp variable if we are over
! an ocean and do not have a value for skin temp
          if(grid%tsk(i,j).eq.0.and.grid%xland(i,j).eq.2) grid%tsk(i,j)=grid%sst(i,j) 

! Fix several defaults for snow cover and energy fluxes; CPL7 default is 1.0e30 and this
! will cause problems for WRF.
          if(grid%snow(i,j).eq.1.0e30) grid%snow(i,j)=0.0 
          if(grid%hfx(i,j).eq.1.0e30)  grid%hfx(i,j)=0.0
          if(grid%lh(i,j).eq.1.0e30)   grid%lh(i,j)=0.0
          if(grid%qfx(i,j).eq.1.0e30)  grid%qfx(i,j)=0.0
          
          ig=ig+1
       end do
      end do
      
      if (grid%itimestep.eq.0.and.first_time) then     
         do j=jps,jpe	
            do i=ips,ipe	
               grid%lwups(i,j) = 5.67e-8*(grid%tsk(i,j)**4)
               grid%pblh(i,j)= 2000.0 ! need an initial value
               grid%ust(i,j) = 0.3    ! need an initial value
               first_time=.false.
            enddo
         enddo     
      end if  
      
      call shr_file_setLogUnit (shrlogunit)
      call shr_file_setLogLevel(shrloglev)
     
  end subroutine atm_import_mct

!===============================================================================

  subroutine atm_export_mct( grid, a2x_a )
  
  USE module_state_description
!-------------------------------------------------------------------
!
! Arguments
!
    TYPE(domain) , POINTER:: grid
    type(mct_aVect)    , intent(out) :: a2x_a
!
! Local variables
!
    real :: u_phy,v_phy
    real :: earth_u, earth_v
    real(r8)::z,tv,rd,g,xlapse,alpha,tstar,tt0,alph,beta,psfc,pslv,p
    integer :: i,j,k,ig       ! indices    
    integer :: ids,ide,jds,jde,kds,kde,ims,ime,jms,jme,kms,kme,ips,ipe,jps,jpe,kps,kpe
    integer :: xtime_week
    real :: uv_factor
    real :: xtime_week_rl
    character (256) :: msg
!-----------------------------------------------------------------------
! Copy from component arrays into chunk array data structure
! Rearrange data from chunk structure into lat-lon buffer and subsequently
! create attribute vector
    rd=287.0_r8
    g=9.81_r8
    xlapse= 6.5e-3_r8
    
    call shr_file_getLogUnit (shrlogunit)
    call shr_file_getLogLevel(shrloglev)
    call shr_file_setLogUnit (iulog)

    call get_ijk_from_grid (  grid ,                   &
                           ids, ide, jds, jde, kds, kde,    &
                           ims, ime, jms, jme, kms, kme,    &
                           ips, ipe, jps, jpe, kps, kpe    )

#ifdef DM_PARALLEL
#     include "HALO_EM_PHYS_A.inc"
#endif	
!    write(iulog,*) local_communicator 
    if(jpe.eq.jde) jpe=jpe-1
    if(ipe.eq.ide) ipe=ipe-1

    ig=1
    do j=jps, jpe             
       do i=ips,ipe 

! Sea level pressure, Pa
! First, calculate virtual temperature at lowest model level in K.
	 tv=(grid%t_2(i,1,j)+300)*&  ! This is potential tempearture at lowest model level.
	       (((grid%p(i,1,j) + grid%pb(i,1,j))/100000.0)**(287.0/1004.6))*& ! mult by inverted Poisson's equation
	       (1+grid%moist(i,1,j,P_QV)*0.61)*1.0_r8	     ! multiply by (1+0.61*mixing ratio) to get virtual temp

! Calculate geopotential at surface in m2 s-2
         z=grid%ht(i,j) *9.81_r8 ! terrain height in m
	 	 
	 alpha = rd*xlapse/g ! 0.1903, no units

! Get surface pressure in Pa
	 psfc= grid%psfc(i,j)*1.0_r8

! Calculate pressure in WRF at lowest model level in Pa
	 p=(grid%p(i,1,j) + grid%pb(i,1,j))*1.0_r8

! Use surface pressure for SLP if we are at sea level.
         if ( abs(z/g) < 1.e-4_r8 )then
            pslv=psfc
! Othewise, calculate based on ECMWF method (modified from vintp2p_ecmwf.F)
         else
            tstar=tv*(1._r8+alpha*(psfc/p-1._r8)) 
            tt0=tstar + xlapse*z/g

            if ( tstar<=290.5_r8 .and. tt0>290.5_r8 ) then     
               alph=rd/z*(290.5_r8-tstar)  
            else if (tstar>290.5_r8  .and. tt0>290.5_r8) then  
               alph=0._r8
               tstar= 0.5_r8 * (290.5_r8 + tstar)  
            else  
               alph=alpha  
               if (tstar<255._r8) then  
                  tstar= 0.5_r8 * (255._r8 + tstar)             
               endif
            endif
            beta = z/(rd*tstar)
            pslv=psfc*exp( beta*(1._r8-alph*beta/2._r8+((alph*beta)**2)/3._r8))
         end if
      
! Finally, output sea level pressure to coupler.
         tv=pslv  
         a2x_a%rAttr(index_a2x_Sa_pslv   ,ig) = tv 
	  
! bottom atm level height (depth), m
         tv=(grid%phb(i,1,j)+grid%ph_2(i,1,j)+grid%phb(i,2,j)+grid%ph_2(i,2,j))*0.5_r8/9.81_r8-z/9.81_r8
         a2x_a%rAttr(index_a2x_Sa_z      ,ig) = tv 

! zonal wind at lowest model level, m/s
! meridional wind at lowest model level, m/s
! First, get winds from WRF at lowest model level
         u_phy=(grid%u_2(i,1,j)+grid%u_2(i+1,1,j))/2.0
         v_phy=(grid%v_2(i,1,j)+grid%v_2(i,1,j+1))/2.0
! Now rotate winds from WRF grid to true grid
         earth_u = u_phy*grid%cosa(i,j)-v_phy*grid%sina(i,j)
         earth_v = v_phy*grid%cosa(i,j)+u_phy*grid%sina(i,j)
! Output rotated winds to coupler.
! U
         tv=earth_u*1.0_r8
         a2x_a%rAttr(index_a2x_Sa_u      ,ig) = tv
! V
         tv=earth_v*1.0_r8
         a2x_a%rAttr(index_a2x_Sa_v      ,ig) = tv

! pressure at the lowest model level, Pa  
         u_phy=grid%pb(i,1,j) + grid%p(i,1,j)
         tv=u_phy*1.0_r8
         a2x_a%rAttr(index_a2x_Sa_pbot   ,ig) = tv 

! potential temperature at the lowest model level, K
         u_phy=grid%t_2(i,1,j)+300.0
         tv=u_phy*1.0_r8
         a2x_a%rAttr(index_a2x_Sa_ptem   ,ig) = tv 

! temperature at the lowest model level, K
         u_phy=(grid%t_2(i,1,j)+300.0)*((grid%pb(i,1,j) + grid%p(i,1,j))/100000.0)**0.2854
         tv=u_phy*1.0_r8
         a2x_a%rAttr(index_a2x_Sa_tbot   ,ig) = tv

! air density at the lowest model level, kg/m3
! this includes moisture; density in WRF is dry air density
          if(grid%moist(i,1,j,P_QV).gt.0) then 
             u_phy=1.0/grid%alt(i,1,j)*(1+grid%moist(i,1,j,P_QV))
          else
             u_phy=1.0/grid%alt(i,1,j)
          endif
	  tv=u_phy*1.0_r8      
	  a2x_a%rAttr(index_a2x_Sa_dens   ,ig) =tv

! specific humidity at the lowest model level, kg/kg
! calculated from mixing ratio in WRF
          if(grid%moist(i,1,j,P_QV).gt.0) then
	  u_phy=grid%moist(i,1,j,P_QV)/(1+grid%moist(i,1,j,P_QV)) 
	  else
          u_phy=0.0
          endif 
	  tv=u_phy*1.0_r8      
          a2x_a%rAttr(index_a2x_Sa_shum ,ig) = tv
	  
! net shortwave radiation, W/m2
	  tv=grid%gsw(i,j)*1.0_r8
	  a2x_a%rAttr(index_a2x_Faxa_swnet,ig) = tv

! downward longwave heat flux, W/m2
	  tv=grid%glw(i,j)*1.0_r8
          a2x_a%rAttr(index_a2x_Faxa_lwdn ,ig) = tv

! convective precipitation rate in kg m2-1 s-1
! Note that the above coupler description is very misleading. The coupler                                             
! wants non-frozen convective precip here, not all precip. However, WRF's                                             
! cumulus scheme (Grell, module_cu_gd.F specifically) does not output frozen                                          
! and non-frozen convective precip seperately, and so as we are exporting                                             
! to the coupler, this variable contains both frozen and unfrozen. We could                                           
! To calculate, divide WRF's time-step cumulus precipitation in mm by the time step,                                  
! multiply by 1000 kg/m3 and divide by 1000 mm/m 
	  tv=grid%raincv(i,j)/grid%dt*1.0_r8
          a2x_a%rAttr(index_a2x_Faxa_rainc,ig) = tv

! convective frozen precip rate (water equivalent) in kg m-2 s-1                                                      
! WRF's cumulus schemes do not output convective frozen precip variables at all, and                                  
! there is no way to get at this information. So, we set this equal to zero here. Therefore,                          
! important to note that the convective non-frozen precip is actually frozen+non-frozen.                              
          tv = 0.*grid%dt*1.0_r8
          a2x_a%rAttr(index_a2x_Faxa_snowc,ig) = tv

! large scale (stable) precipitation rate in kg m-2 s-1                                                               
! Note that the above coupler description is very misleading. The coupler wants                                       
! non-frozen nonconvective precip, not total nonconvective precip. To calculate this,                                 
! multiply WRF's total nonconvective precip by the fraction of liquid precip (1-SR).                                  
! SR is the fraction of frozen, non-convective precipitation and is calculated in                                     
! the microphysics scheme (module_mp_gsfcgce.F). Then, divide by the WRF time step,                                   
! multiply by 1000 kg/m3 and divide by 1000 mm/m.                                                                     
          tv=(grid%rainncv(i,j)*(1-grid%sr(i,j)))/grid%dt*1.0_r8
          a2x_a%rAttr(index_a2x_Faxa_rainl,ig) =  tv

! large-scale (stable) frozen precip rate (water equivalent) in kg m-2 s-1                                            
! Add WRF's snowncv, graupelncv, and hailncv (time-step nonconvective snow and ice, graupel,                          
! and hail all in mm) and then divide by the time step, multiply by 1000 kg/m3 and                                    
! divide by 1000 mm/m.                                                                                                
          tv=(grid%snowncv(i,j)+grid%graupelncv(i,j)+grid%hailncv(i,j))/grid%dt*1.0_r8                                
          a2x_a%rAttr(index_a2x_Faxa_snowl,ig) = tv

! The following four radiation fields are not normally part of WRF's Registry.
! Instead, the fields have been added to the WRF Registry and the WRF 
! radiation module has been modified to export these fields to the Registry
! so they may then be exported here to the coupler.

! direct near-infrared incident solar radiation in W/m2
	  tv=grid%swndr(i,j)*1.0_r8
          a2x_a%rAttr(index_a2x_Faxa_swndr,ig) = tv

! direct visible incident solar radiation in W/m2
	  tv=grid%swvdr(i,j)*1.0_r8
          a2x_a%rAttr(index_a2x_Faxa_swvdr,ig) = tv

! diffuse near-infrared incident solar radiation in W/m2
	  tv=grid%swndf(i,j)*1.0_r8
          a2x_a%rAttr(index_a2x_Faxa_swndf,ig) = tv

! diffuse visible incident solar radiation in W/m2
	  tv=grid%swvdf(i,j)*1.0_r8
          a2x_a%rAttr(index_a2x_Faxa_swvdf,ig) = tv

          if (index_a2x_Sa_co2prog /= 0) then
             a2x_a%rAttr(index_a2x_Sa_co2prog,ig) = 3.55e-4_r8 ! atm prognostic co2
          end if
          if (index_a2x_Sa_co2diag /= 0) then
             a2x_a%rAttr(index_a2x_Sa_co2diag,ig) = 3.55e-4_r8 ! atm diagnostic co2
          end if
          ig=ig+1
       end do
    end do
       call shr_file_setLogUnit (shrlogunit)
       call shr_file_setLogLevel(shrloglev)    

  end subroutine atm_export_mct

!===============================================================================

  subroutine atm_domain_mct( grid, lsize, gsMap_a, dom_a )

!-------------------------------------------------------------------
! Arguments
!
    TYPE(domain) , POINTER:: grid
    integer        , intent(in)   :: lsize
    type(mct_gsMap), intent(in)   :: gsMap_a
    type(mct_ggrid), intent(inout):: dom_a    
!
! Local Variables
!
    integer  :: i,j,k,m,n           ! indices	
    integer :: ids,ide,jds,jde,kds,kde,ims,ime,jms,jme,kms,kme,ips,ipe,jps,jpe,kps,kpe
    real(r8)::del_phi,del_theta
    real(r8),dimension(:,:,:),allocatable :: ew_vert,ns_vert,xcenterlat,xcenterlon ! vertex of cell    
    real(r8),dimension(:,:),allocatable :: area          ! area in radians squared for each grid point
    real(r8), pointer  :: data(:)     ! temporary
    integer , pointer  :: idata(:)    ! temporary
    
  TYPE proj_info
 
      INTEGER          :: code     ! Integer code for projection TYPE
      INTEGER          :: nlat     ! For Gaussian -- number of latitude points 
                                   !  north of the equator 
      INTEGER          :: nlon     !
                                   !
      INTEGER          :: ixdim    ! For Rotated Lat/Lon -- number of mass points
                                   !  in an odd row
      INTEGER          :: jydim    ! For Rotated Lat/Lon -- number of rows
      INTEGER          :: stagger  ! For Rotated Lat/Lon -- mass or velocity grid 
      REAL             :: phi      ! For Rotated Lat/Lon -- domain half-extent in 
                                   !  degrees latitude
      REAL             :: lambda   ! For Rotated Lat/Lon -- domain half-extend in
                                   !  degrees longitude
      REAL             :: lat1     ! SW latitude (1,1) in degrees (-90->90N)
      REAL             :: lon1     ! SW longitude (1,1) in degrees (-180->180E)
      REAL             :: lat0     ! For Cassini, latitude of projection pole
      REAL             :: lon0     ! For Cassini, longitude of projection pole
      REAL             :: dx       ! Grid spacing in meters at truelats, used
                                   !  only for ps, lc, and merc projections
      REAL             :: dy       ! Grid spacing in meters at truelats, used
                                   !  only for ps, lc, and merc projections
      REAL             :: latinc   ! Latitude increment for cylindrical lat/lon
      REAL             :: loninc   ! Longitude increment for cylindrical lat/lon
                                   !  also the lon increment for Gaussian grid
      REAL             :: dlat     ! Lat increment for lat/lon grids
      REAL             :: dlon     ! Lon increment for lat/lon grids
      REAL             :: stdlon   ! Longitude parallel to y-axis (-180->180E)
      REAL             :: truelat1 ! First true latitude (all projections)
      REAL             :: truelat2 ! Second true lat (LC only)
      REAL             :: hemi     ! 1 for NH, -1 for SH
      REAL             :: cone     ! Cone factor for LC projections
      REAL             :: polei    ! Computed i-location of pole point
      REAL             :: polej    ! Computed j-location of pole point
      REAL             :: rsw      ! Computed radius to SW corner
      REAL             :: rebydx   ! Earth radius divided by dx
      REAL             :: knowni   ! X-location of known lat/lon
      REAL             :: knownj   ! Y-location of known lat/lon
      REAL             :: re_m     ! Radius of spherical earth, meters
      REAL             :: rho0     ! For Albers equal area
      REAL             :: nc       ! For Albers equal area
      REAL             :: bigc     ! For Albers equal area
      LOGICAL          :: init     ! Flag to indicate if this struct is 
                                   !  ready for use
      LOGICAL          :: wrap     ! For Gaussian -- flag to indicate wrapping 
                                   !  around globe?
      LOGICAL          :: comp_ll  ! Work in computational lat/lon space for Cassini
      REAL, POINTER, DIMENSION(:) :: gauss_lat  ! Latitude array for Gaussian grid
 
   END TYPE proj_info    

      TYPE(proj_info) :: proj

      INTEGER          :: code     ! Integer code for projection TYPE
      INTEGER          :: nlat     ! For Gaussian -- number of latitude points 
                                   !  north of the equator 
      INTEGER          :: nlon     !
                                   !
      INTEGER          :: ixdim    ! For Rotated Lat/Lon -- number of mass points
                                   !  in an odd row
      INTEGER          :: jydim    ! For Rotated Lat/Lon -- number of rows
      INTEGER          :: stagger  ! For Rotated Lat/Lon -- mass or velocity grid 
      REAL             :: phi      ! For Rotated Lat/Lon -- domain half-extent in 
                                   !  degrees latitude
      REAL             :: lambda   ! For Rotated Lat/Lon -- domain half-extend in
                                   !  degrees longitude
      REAL             :: lat1     ! SW latitude (1,1) in degrees (-90->90N)
      REAL             :: lon1     ! SW longitude (1,1) in degrees (-180->180E)
      REAL             :: lat0     ! For Cassini, latitude of projection pole
      REAL             :: lon0     ! For Cassini, longitude of projection pole
      REAL             :: dx       ! Grid spacing in meters at truelats, used
                                   !  only for ps, lc, and merc projections
      REAL             :: dy       ! Grid spacing in meters at truelats, used
                                   !  only for ps, lc, and merc projections
      REAL             :: latinc   ! Latitude increment for cylindrical lat/lon
      REAL             :: loninc   ! Longitude increment for cylindrical lat/lon
                                   !  also the lon increment for Gaussian grid
      REAL             :: dlat     ! Lat increment for lat/lon grids
      REAL             :: dlon     ! Lon increment for lat/lon grids
      REAL             :: stdlon   ! Longitude parallel to y-axis (-180->180E)
      REAL             :: truelat1 ! First true latitude (all projections)
      REAL             :: truelat2 ! Second true lat (LC only)
      REAL             :: hemi     ! 1 for NH, -1 for SH
      REAL             :: cone     ! Cone factor for LC projections
      REAL             :: polei    ! Computed i-location of pole point
      REAL             :: polej    ! Computed j-location of pole point
      REAL             :: rsw      ! Computed radius to SW corner
      REAL             :: rebydx   ! Earth radius divided by dx
      REAL             :: knowni   ! X-location of known lat/lon
      REAL             :: knownj   ! Y-location of known lat/lon
      REAL             :: re_m     ! Radius of spherical earth, meters
      REAL             :: rho0     ! For Albers equal area
      REAL             :: nc       ! For Albers equal area
      REAL             :: bigc     ! For Albers equal area
      LOGICAL          :: init     ! Flag to indicate if this struct is 
                                   !  ready for use
      LOGICAL          :: wrap     ! For Gaussian -- flag to indicate wrapping 
                                   !  around globe?
      LOGICAL          :: comp_ll  ! Work in computational lat/lon space for Cassini
      REAL, POINTER, DIMENSION(:) :: gauss_lat  ! Latitude array for Gaussian grid
 
      REAL :: dummy_lon1, dummy_lon0, dummy_stdlon
      REAL :: ala, ala1, alo1, reflon, scale_top, alo, rm, xx, yy, gi2, r2, &
                  arccos, lat, lon, ref_lat, ref_lon, stand_lon
		  
   real, parameter :: PI = 3.141592653589793
   real, parameter :: OMEGA_E = 7.292e-5 ! Angular rotation rate of the earth

   real, parameter :: DEG_PER_RAD = 180./PI
   real, parameter :: RAD_PER_DEG = PI/180.
 
   real, parameter :: EARTH_RADIUS_M = 6370000.   ! same as MM5 system
    
!-------------------------------------------------------------------
! Initialize map projection
!

      proj%lat1     = -999.9
      proj%lon1     = -999.9
      proj%lat0     = -999.9
      proj%lon0     = -999.9
      proj%dx       = -999.9
      proj%dy       = -999.9
      proj%latinc   = -999.9
      proj%loninc   = -999.9
      proj%stdlon   = -999.9
      proj%truelat1 = -999.9
      proj%truelat2 = -999.9
      proj%phi      = -999.9
      proj%lambda   = -999.9
      proj%ixdim    = -999
      proj%jydim    = -999
      proj%stagger  = 0
      proj%nlat     = 0
      proj%nlon     = 0
      proj%hemi     = 0.0
      proj%cone     = -999.9
      proj%polei    = -999.9
      proj%polej    = -999.9
      proj%rsw      = -999.9
      proj%knowni   = -999.9
      proj%knownj   = -999.9
      proj%re_m     = EARTH_RADIUS_M
      proj%init     = .FALSE.
      proj%wrap     = .FALSE.
      proj%rho0     = 0.
      proj%nc       = 0.
      proj%bigc     = 0.
      proj%comp_ll  = .FALSE.
      nullify(proj%gauss_lat)
      
      dummy_lon1=-999.0
      dummy_lon0=-999.0
      dummy_stdlon=-999.0      
		  
      dx=25000.0
      dy=25000.0
      ref_lat   =   80.5
      ref_lon   = 160.0
      truelat1  =   80.5
      stand_lon = -114.0
      knowni = 277
      knownj  = 207
      lat1=80.5
      lon1=160
	
      ! Given a partially filled proj_info structure, this routine computes
      ! polei, polej, rsw, and cone (if LC projection) to complete the 
      ! structure.  This allows us to eliminate redundant calculations when
      ! calling the coordinate conversion routines multiple times for the
      ! same map.
      ! This will generally be the first routine called when a user wants
      ! to be able to use the coordinate conversion routines, and it
      ! will call the appropriate subroutines based on the 
      ! proj%code which indicates which projection type this is.  
      
      proj%dx       = dx
      proj%dy       = dy
      proj%stdlon   = stand_lon
      proj%truelat1 = truelat1
      proj%truelat2 = truelat2
      proj%knowni   = knowni
      proj%knownj   = knownj
      proj%lat1   = lat1
      proj%lon1   = lon1
      
            IF (truelat1 .LT. 0.) THEN
               proj%hemi = -1.0 
            ELSE
               proj%hemi = 1.0
            ENDIF
            proj%rebydx = proj%re_m / dx

      ! Initializes a polar-stereographic map projection from the partially
      ! filled proj structure. This routine computes the radius to the
      ! southwest corner and computes the i/j location of the pole for use
      ! in llij_ps and ijll_ps.
  
      ! Executable code
      reflon = proj%stdlon + 90.
  
      ! Compute numerator term of map scale factor
      scale_top = 1. + proj%hemi * SIN(proj%truelat1 * rad_per_deg)
  
      ! Compute radius to lower-left (SW) corner
      ala1 = proj%lat1 * rad_per_deg
      proj%rsw = proj%rebydx*COS(ala1)*scale_top/(1.+proj%hemi*SIN(ala1))
  
      ! Find the pole point
      alo1 = (proj%lon1 - reflon) * rad_per_deg
      proj%polei = proj%knowni - proj%rsw * COS(alo1)
      proj%polej = proj%knownj - proj%hemi * proj%rsw * SIN(alo1)

!
! Initialize mct atm domain
!
    call mct_gGrid_init( GGrid=dom_a, CoordChars=trim(seq_flds_dom_coord), OtherChars=trim(seq_flds_dom_other), lsize=lsize )
!
! Allocate memory
!
    allocate(data(lsize))
!
! Initialize attribute vector with special value
!
    call mct_gsMap_orderedPoints(gsMap_a, mytask, idata)
    call mct_gGrid_importIAttr(dom_a,'GlobGridNum',idata,lsize)
!
! Determine domain (numbering scheme is: West to East and South to North to South pole)
! Initialize attribute vector with special value
!
    data(:) = -9999.0_R8 
    call mct_gGrid_importRAttr(dom_a,"lat"  ,data,lsize) 
    call mct_gGrid_importRAttr(dom_a,"lon"  ,data,lsize) 
    call mct_gGrid_importRAttr(dom_a,"area" ,data,lsize) 
    call mct_gGrid_importRAttr(dom_a,"aream",data,lsize) 
    data(:) = 0.0_R8     
    call mct_gGrid_importRAttr(dom_a,"mask" ,data,lsize) 
    data(:) = 1.0_R8
    call mct_gGrid_importRAttr(dom_a,"frac" ,data,lsize)

!
! Fill in correct values for domain components
!

    call get_ijk_from_grid (grid ,                   &
                            ids, ide, jds, jde, kds, kde,    &
                            ims, ime, jms, jme, kms, kme,    &    
			    ips, ipe, jps, jpe, kps, kpe    )
    
    if(jpe.eq.jde) jpe=jpe-1
    if(ipe.eq.ide) ipe=ipe-1
    m=2*(ide-ids)+3
    n=2*(jde-jds)+3
	
       allocate(area(ips:ipe,jps:jpe))
        allocate(xcenterlat(1:m,1:n,1:1))
        allocate(xcenterlon(1:m,1:n,1:1))
   	allocate(ew_vert(4,ids:ide,jds:jde))
	allocate(ns_vert(4,ids:ide,jds:jde))

! lat	
    n=0
    do j=jps, jpe             
       do i=ips,ipe
          n = n+1
          data(n) = grid%xlat(i,j)
       end do
    end do
    call mct_gGrid_importRAttr(dom_a,"lat",data,lsize) 

! lon
    n=0
    do j=jps, jpe             
       do i=ips,ipe
          n = n+1
          data(n) = grid%xlong(i,j)
       end do
    end do
    call mct_gGrid_importRAttr(dom_a,"lon",data,lsize) 

! the corner of the gridcells
    
       do j=1,2*(jde-jds)+3
       do i=1,2*(ide-ids)+3
       
      ! Given latitude (-90 to 90), longitude (-180 to 180), and the
      ! standard polar-stereographic projection information via the 
      ! public proj structure, this routine returns the i/j indices which
      ! if within the domain range from 1->nx and 1->ny, respectively.

      ! BEGIN CODE
    
!      reflon = proj%stdlon + 90.
     
      ! Compute numerator term of map scale factor
  
!      scale_top = 1. + proj%hemi * SIN(proj%truelat1 * rad_per_deg)
  
      ! Find radius to desired point
!      ala = lat * rad_per_deg
!      rm = proj%rebydx * COS(ala) * scale_top/(1. + proj%hemi *SIN(ala))
!      alo = (lon - reflon) * rad_per_deg
!      i = proj%polei + rm * COS(alo)
!      j = proj%polej + proj%hemi * rm * SIN(alo)
 
      ! This is the inverse subroutine of llij_ps.  It returns the 
      ! latitude and longitude of an i/j point given the projection info 
      ! structure.  
    
      ! Begin Code
  
      ! Compute the reference longitude by rotating 90 degrees to the east
      ! to find the longitude line parallel to the positive x-axis.
      reflon = proj%stdlon + 90.
     
      ! Compute numerator term of map scale factor
      scale_top = 1. + proj%hemi * SIN(proj%truelat1 * rad_per_deg)
  
      ! Compute radius to point of interest
      xx = i - proj%polei
      yy = (j - proj%polej) * proj%hemi
      r2 = xx**2 + yy**2
  
      ! Now the magic code
      IF (r2 .EQ. 0.) THEN 
         lat = proj%hemi * 90.
         lon = reflon
      ELSE
         gi2 = (proj%rebydx * scale_top)**2.
         lat = deg_per_rad * proj%hemi * ASIN((gi2-r2)/(gi2+r2))
!         arccos = ACOS(xx/SQRT(r2))
         if( xx/SQRT(r2) <= -1. ) then
            arccos = acos(-1.)
         else if( xx/SQRT(r2) >= 1. ) then
            arccos = acos(1.)
         else
            arccos = ACOS(xx/SQRT(r2))
         endif
         IF (yy .GT. 0) THEN
            lon = reflon + deg_per_rad * arccos
         ELSE
            lon = reflon - deg_per_rad * arccos
         ENDIF
      ENDIF
    
      ! Convert to a -180 -> 180 East convention
      IF (lon .GT. 180.) lon = lon - 360.
      IF (lon .LT. -180.) lon = lon + 360.
      
      if (lon.lt.0) lon=lon+360
      
      xcenterlat(i,j,1)=lat
      xcenterlon(i,j,1)=lon
      
      enddo
      enddo
      
       do j=2,jde-jds+1
         do i=2,ide-ids+1

          ns_vert(1,i-1,j-1) = xcenterlat(2*i-2,2*j-2,1)
          ns_vert(2,i-1,j-1) = xcenterlat(2*i,2*j-2,1)
          ns_vert(3,i-1,j-1) = xcenterlat(2*i,2*j,1)
	  ns_vert(4,i-1,j-1) = xcenterlat(2*i-2,2*j,1)
	  
          if(xcenterlon(2*i-2,2*j-2,1).lt.0) xcenterlon(2*i-2,2*j-2,1)=xcenterlon(2*i-2,2*j-2,1)+360.0
          if(xcenterlon(2*i-2,2*j,1).lt.0) xcenterlon(2*i-2,2*j,1)=xcenterlon(2*i-2,2*j,1)+360.0
          if(xcenterlon(2*i,2*j,1).lt.0) xcenterlon(2*i,2*j,1)=xcenterlon(2*i,2*j,1)+360.0
          if(xcenterlon(2*i,2*j-2,1).lt.0) xcenterlon(2*i,2*j-2,1)=xcenterlon(2*i,2*j-2,1)+360.0
          if(xcenterlon(2*i-1,2*j-1,1).lt.0) xcenterlon(2*i-1,2*j-1,1)=xcenterlon(2*i-1,2*j-1,1)+360.0
	  
	    ew_vert(1,i-1,j-1) = xcenterlon(2*i-2,2*j-2,1)    
            ew_vert(2,i-1,j-1) = xcenterlon(2*i,2*j-2,1)	    
	    ew_vert(3,i-1,j-1) = xcenterlon(2*i,2*j,1)
            ew_vert(4,i-1,j-1) = xcenterlon(2*i-2,2*j,1)
         
         enddo
	 enddo

! area

	 do j=jps,jpe
	 do i=ips,ipe
! segement interation
	   area(i,j)=0.0
	   do k=1,4
	   m=k+1
	   if(m.eq.5) m=1
	   del_phi = (dsin( ns_vert(k,i,j)*RAD_PER_DEG )+ dsin( ns_vert(m,i,j)*RAD_PER_DEG))/2.0
	   if(abs(ew_vert(k,i,j)-ew_vert(m,i,j)).gt.190) then  !gridcell cross 0 longitude
	    if(ew_vert(k,i,j).gt.ew_vert(m,i,j)) then  
	     del_theta = ( ew_vert(k,i,j)-ew_vert(m,i,j)-360.0)*RAD_PER_DEG
	    else 
             del_theta = (ew_vert(k,i,j)+360.0-ew_vert(m,i,j) )*RAD_PER_DEG
	    endif
	   else
	    del_theta = (ew_vert(k,i,j)-ew_vert(m,i,j) )*RAD_PER_DEG
           endif
           area(i,j) = area(i,j)+del_theta*del_phi
	   enddo
	   area(i,j)=dabs(area(i,j))

          if(area(i,j).gt.1) then  ! cover the pole
           area(i,j)=0.0
	   do k=1,4
	   m=k+1
	   if(m.eq.5) m=1
	   del_phi = (dsin( ns_vert(k,i,j)*RAD_PER_DEG)+ dsin( ns_vert(m,i,j)*RAD_PER_DEG))/2.0
	   if(abs(ew_vert(k,i,j)-ew_vert(m,i,j)).gt.190) then  !gridcell cross 0 longitude
	    if(ew_vert(k,i,j).gt.ew_vert(m,i,j)) then  
	     del_theta = ( ew_vert(m,i,j)+360-ew_vert(k,i,j) )*RAD_PER_DEG
	    else 
             del_theta = (ew_vert(k,i,j)-360.0-ew_vert(m,i,j)  )*RAD_PER_DEG
	    endif
	   else
	    del_theta = (ew_vert(k,i,j)-ew_vert(m,i,j)  )*RAD_PER_DEG
           endif
           area(i,j) = area(i,j)+del_theta*(1.0-abs(del_phi))
	   enddo
	   area(i,j)=dabs(area(i,j))
	  endif
	  
	 enddo
      enddo

    n=0
    do j=jps, jpe             
       do i=ips,ipe
          n = n+1
          data(n) = area(i,j) 
       end do
    end do
    call mct_gGrid_importRAttr(dom_a,"area",data,lsize) 

    n=0
    do j=jps, jpe             
       do i=ips,ipe
          n = n+1
          data(n) = 1._r8 ! mask
       end do
    end do
    call mct_gGrid_importRAttr(dom_a,"mask"   ,data,lsize) 
    
	deallocate(data)
	deallocate(area)
	deallocate(xcenterlat)
	deallocate(xcenterlon)
	deallocate(ew_vert)
	deallocate(ns_vert)

  end subroutine atm_domain_mct
 
#endif
#endif
  
END MODULE atm_comp_mct
