external/io_pio/pio_routines.F90

   1 !---------------------------------------------------------------------------
   2 !
   3 ! WRF Parallel I/O
   4 ! Author:  Wei Huang huangwei@ucar.edu
   5 ! Date:    June 01, 2014
   6 !
   7 !---------------------------------------------------------------------------
   8 !$Id$
   9 !---------------------------------------------------------------------------
  10
  11 module pio_routines
  12
  13   use pio_kinds
  14   use pio
  15
  16   use module_domain
  17
  18   use wrf_data_pio
  19
  20   implicit none
  21
  22   include 'mpif.h'
  23
  24   integer(i4) :: nprocs, myrank
  25   integer :: ids, ide, jds, jde, kds, kde, &
  26              ims, ime, jms, jme, kms, kme, &
  27              its, ite, jts, jte, kts, kte
  28
  29 CONTAINS
  30
  31 subroutine allocHandle(DataHandle,DH,Status)
  32   implicit none
  33   include 'wrf_status_codes.h'
  34   integer              ,intent(out) :: DataHandle
  35   type(wrf_data_handle),pointer     :: DH
  36   integer              ,intent(out) :: Status
  37   integer                           :: i, n
  38   integer                           :: stat
  39
  40   do i=1,WrfDataHandleMax
  41     if(WrfDataHandles(i)%Free) then
  42       DH => WrfDataHandles(i)
  43       DataHandle = i
  44       do n = 1, MaxVars
  45          DH%vartype(n) = NOT_LAND_SOIL_VAR
  46       end do
  47       exit
  48     endif
  49     if(i==WrfDataHandleMax) then
  50       Status = WRF_WARN_TOO_MANY_FILES
  51       write(msg,*) 'Warning TOO MANY FILES in ',__FILE__,', line', __LINE__
  52       call wrf_debug ( WARN , TRIM(msg))
  53       write(msg,*) 'Did you call ext_pio_ioinit?'
  54       call wrf_debug ( WARN , TRIM(msg))
  55       return
  56     endif
  57   enddo
  58   DH%Free      =.false.
  59   DH%Write     =.false.
  60   DH%first_operation  = .TRUE.
  61   DH%CurrentVariable = 0
  62   Status = WRF_NO_ERR
  63 end subroutine allocHandle
  64
  65 subroutine deallocHandle(DataHandle, Status)
  66   implicit none
  67   include 'wrf_status_codes.h'
  68   integer              ,intent(in) :: DataHandle
  69   integer              ,intent(out) :: Status
  70   type(wrf_data_handle),pointer     :: DH
  71   integer                           :: i
  72   integer                           :: stat
  73
  74   IF ( DataHandle .GE. 1 .AND. DataHandle .LE. WrfDataHandleMax ) THEN
  75     if(.NOT. WrfDataHandles(DataHandle)%Free) then
  76       DH => WrfDataHandles(DataHandle)
  77       DH%Free      =.TRUE.
  78     endif
  79
  80    !deallocate(DH%iosystem)
  81   ENDIF
  82   Status = WRF_NO_ERR
  83 end subroutine deallocHandle
  84
  85 subroutine GetDH(DataHandle,DH,Status)
  86   implicit none
  87   include 'wrf_status_codes.h'
  88   integer               ,intent(in)     :: DataHandle
  89   type(wrf_data_handle) ,pointer        :: DH
  90   integer               ,intent(out)    :: Status
  91
  92   if(DataHandle < 1 .or. DataHandle > WrfDataHandleMax) then
  93     Status = WRF_WARN_BAD_DATA_HANDLE
  94     return
  95   endif
  96   DH => WrfDataHandles(DataHandle)
  97   if(DH%Free) then
  98     Status = WRF_WARN_BAD_DATA_HANDLE
  99     return
 100   endif
 101   Status = WRF_NO_ERR
 102   return
 103 end subroutine GetDH
 104
 105 subroutine DateCheck(Date,Status)
 106   implicit none
 107   include 'wrf_status_codes.h'
 108   character*(*) ,intent(in)      :: Date
 109   integer       ,intent(out)     :: Status
 110
 111   if(len(Date) /= DateStrLen) then
 112     Status = WRF_WARN_DATESTR_BAD_LENGTH
 113   else
 114     Status = WRF_NO_ERR
 115   endif
 116   return
 117 end subroutine DateCheck
 118
 119 subroutine GetName(Element,Var,Name,Status)
 120   implicit none
 121   include 'wrf_status_codes.h'
 122   character*(*) ,intent(in)     :: Element
 123   character*(*) ,intent(in)     :: Var
 124   character*(*) ,intent(out)    :: Name
 125   integer       ,intent(out)    :: Status
 126   character (VarNameLen)        :: VarName
 127   character (1)                 :: c
 128   integer                       :: i
 129   integer, parameter            ::  upper_to_lower =IACHAR('a')-IACHAR('A')
 130
 131   VarName = Var
 132   Name = 'MD___'//trim(Element)//VarName
 133   do i=1,len(Name)
 134     c=Name(i:i)
 135     if('A'<=c .and. c <='Z') Name(i:i)=achar(iachar(c)+upper_to_lower)
 136     if(c=='-'.or.c==':') Name(i:i)='_'
 137   enddo
 138   Status = WRF_NO_ERR
 139   return
 140 end subroutine GetName
 141
 142 subroutine GetTimeIndex(IO,DataHandle,DateStr,TimeIndex,Status)
 143   implicit none
 144   include 'wrf_status_codes.h'
 145   character (*)         ,intent(in)     :: IO
 146   integer               ,intent(in)     :: DataHandle
 147   character*(*)         ,intent(in)     :: DateStr
 148   integer               ,intent(out)    :: TimeIndex
 149   integer               ,intent(out)    :: Status
 150   type(wrf_data_handle) ,pointer        :: DH
 151   integer                               :: VStart(2)
 152   integer                               :: VCount(2)
 153   integer                               :: stat
 154   integer                               :: i
 155   character(len=DateStrLen)             :: tmpdatestr(1)
 156
 157   if(len(Datestr) == DateStrLen) then
 158     tmpdatestr = DateStr
 159   else
 160     write(unit=0, fmt='(3a,i6)') 'file: ', __FILE__, ', line: ', __LINE__
 161     write(unit=0, fmt='(3a)') 'IO: <', trim(IO), '>'
 162     write(unit=0, fmt='(a,i3)') 'DataHandle = ', DataHandle
 163     write(unit=0, fmt='(3a)') 'DateStr: <', trim(DateStr), '>'
 164     write(unit=0, fmt='(a,i6,a,i6)') 'DateStrLen = ', DateStrLen, &
 165                                      ' did not equal len(DateStr): ', len(DateStr)
 166     Status =  WRF_WARN_DATESTR_ERROR
 167     write(msg,*) 'Warning DATE STRING ERROR in ',__FILE__,', line', __LINE__
 168     call wrf_debug ( WARN , TRIM(msg))
 169     return
 170   end if
 171
 172   DH => WrfDataHandles(DataHandle)
 173   call DateCheck(DateStr,Status)
 174   if(Status /= WRF_NO_ERR) then
 175     Status =  WRF_WARN_DATESTR_ERROR
 176     write(msg,*) 'Warning DATE STRING ERROR in ',__FILE__,', line', __LINE__
 177     call wrf_debug ( WARN , TRIM(msg))
 178     return
 179   endif
 180   if(IO == 'write') then
 181     TimeIndex = DH%TimeIndex
 182     if(TimeIndex <= 0) then
 183       TimeIndex = 1
 184     elseif(DateStr == DH%Times(TimeIndex)) then
 185       Status = WRF_NO_ERR
 186       return
 187     else
 188       TimeIndex = TimeIndex + 1
 189       if(TimeIndex > MaxTimes) then
 190         Status = WRF_WARN_TIME_EOF
 191         write(msg,*) 'Warning TIME EOF in ',__FILE__,', line', __LINE__
 192         call wrf_debug ( WARN , TRIM(msg))
 193         return
 194       endif
 195     endif
 196     DH%TimeIndex        = TimeIndex
 197     DH%Times(TimeIndex) = DateStr
 198     VStart(1) = 1
 199     VStart(2) = TimeIndex
 200     VCount(1) = DateStrLen
 201     VCount(2) = 1
 202    !write(unit=0, fmt='(3a,i6)') 'DateStr: <', trim(DateStr), '>, TimeIndex =', TimeIndex
 203     stat = pio_put_var(DH%file_handle, DH%vtime, VStart, VCount, tmpdatestr)
 204     call netcdf_err(stat,Status)
 205     if(Status /= WRF_NO_ERR) then
 206       write(msg,*) 'NetCDF error in ',__FILE__,', line', __LINE__
 207       call wrf_debug ( WARN , TRIM(msg))
 208       return
 209     endif
 210    !call pio_advanceframe(DH%vtime)
 211   else
 212     do i=1,MaxTimes
 213       if(DH%Times(i)==DateStr) then
 214         Status = WRF_NO_ERR
 215         TimeIndex = i
 216         exit
 217       endif
 218       if(i==MaxTimes) then
 219         Status = WRF_WARN_TIME_NF
 220         write(msg,*) 'Warning TIME ',DateStr,' NOT FOUND in ',__FILE__,', line', __LINE__
 221         call wrf_debug ( WARN , TRIM(msg))
 222         return
 223       endif
 224     enddo
 225   endif
 226   return
 227 end subroutine GetTimeIndex
 228
 229 subroutine GetDim(MemoryOrder,NDim,Status)
 230   implicit none
 231   include 'wrf_status_codes.h'
 232   character*(*) ,intent(in)  :: MemoryOrder
 233   integer       ,intent(out) :: NDim
 234   integer       ,intent(out) :: Status
 235   character*3                :: MemOrd
 236
 237   call LowerCase(MemoryOrder,MemOrd)
 238   select case (MemOrd)
 239     case ('xyz','xzy','yxz','yzx','zxy','zyx','xsz','xez','ysz','yez')
 240       NDim = 3
 241     case ('xy','yx','xs','xe','ys','ye')
 242       NDim = 2
 243     case ('z','c')
 244       NDim = 1
 245     case ('0')  ! NDim=0 for scalars.  TBH:  20060502
 246       NDim = 0
 247     case default
 248       print *, 'memory order = ',MemOrd,'  ',MemoryOrder
 249       Status = WRF_WARN_BAD_MEMORYORDER
 250       return
 251   end select
 252   Status = WRF_NO_ERR
 253   return
 254 end subroutine GetDim
 255
 256 subroutine GetIndices(NDim,Start,End,i1,i2,j1,j2,k1,k2)
 257   implicit none
 258   include 'wrf_status_codes.h'
 259   integer              ,intent(in)  :: NDim
 260   integer ,dimension(*),intent(in)  :: Start,End
 261   integer              ,intent(out) :: i1,i2,j1,j2,k1,k2
 262
 263   i1=1
 264   i2=1
 265   j1=1
 266   j2=1
 267   k1=1
 268   k2=1
 269   if(NDim == 0) return  ! NDim=0 for scalars.  TBH:  20060502
 270   i1 = Start(1)
 271   i2 = End  (1)
 272   if(NDim == 1) return
 273   j1 = Start(2)
 274   j2 = End  (2)
 275   if(NDim == 2) return
 276   k1 = Start(3)
 277   k2 = End  (3)
 278   return
 279 end subroutine GetIndices
 280
 281 logical function ZeroLengthHorzDim(MemoryOrder,Vector,Status)
 282   implicit none
 283   include 'wrf_status_codes.h'
 284   character*(*)              ,intent(in)    :: MemoryOrder
 285   integer,dimension(*)       ,intent(in)    :: Vector
 286   integer                    ,intent(out)   :: Status
 287   integer                                   :: NDim
 288   integer,dimension(NVarDims)               :: temp
 289   character*3                               :: MemOrd
 290   logical zero_length
 291
 292   call GetDim(MemoryOrder,NDim,Status)
 293   temp(1:NDim) = Vector(1:NDim)
 294   call LowerCase(MemoryOrder,MemOrd)
 295   zero_length = .false.
 296   select case (MemOrd)
 297     case ('xsz','xez','ysz','yez','xs','xe','ys','ye','z','c')
 298       continue
 299     case ('0')
 300       continue  ! NDim=0 for scalars.  TBH:  20060502
 301     case ('xzy','yzx')
 302       zero_length = temp(1) .lt. 1 .or. temp(3) .lt. 1
 303     case ('xy','yx','xyz','yxz')
 304       zero_length = temp(1) .lt. 1 .or. temp(2) .lt. 1
 305     case ('zxy','zyx')
 306       zero_length = temp(2) .lt. 1 .or. temp(3) .lt. 1
 307     case default
 308       Status = WRF_WARN_BAD_MEMORYORDER
 309       ZeroLengthHorzDim = .true.
 310       return
 311   end select
 312   Status = WRF_NO_ERR
 313   ZeroLengthHorzDim = zero_length
 314   return
 315 end function ZeroLengthHorzDim
 316
 317 subroutine ExtOrder(MemoryOrder,Vector,Status)
 318   implicit none
 319   include 'wrf_status_codes.h'
 320   character*(*)              ,intent(in)    :: MemoryOrder
 321   integer,dimension(*)       ,intent(inout) :: Vector
 322   integer                    ,intent(out)   :: Status
 323   integer                                   :: NDim
 324   integer,dimension(NVarDims)               :: temp
 325   character*3                               :: MemOrd
 326
 327   call GetDim(MemoryOrder,NDim,Status)
 328   temp(1:NDim) = Vector(1:NDim)
 329   call LowerCase(MemoryOrder,MemOrd)
 330   select case (MemOrd)
 331
 332     case ('xyz','xsz','xez','ysz','yez','xy','xs','xe','ys','ye','z','c')
 333       continue
 334     case ('0')
 335       continue  ! NDim=0 for scalars.  TBH:  20060502
 336     case ('xzy')
 337       Vector(2) = temp(3)
 338       Vector(3) = temp(2)
 339     case ('yxz')
 340       Vector(1) = temp(2)
 341       Vector(2) = temp(1)
 342     case ('yzx')
 343       Vector(1) = temp(3)
 344       Vector(2) = temp(1)
 345       Vector(3) = temp(2)
 346     case ('zxy')
 347       Vector(1) = temp(2)
 348       Vector(2) = temp(3)
 349       Vector(3) = temp(1)
 350     case ('zyx')
 351       Vector(1) = temp(3)
 352       Vector(3) = temp(1)
 353     case ('yx')
 354       Vector(1) = temp(2)
 355       Vector(2) = temp(1)
 356     case default
 357       Status = WRF_WARN_BAD_MEMORYORDER
 358       return
 359   end select
 360   Status = WRF_NO_ERR
 361   return
 362 end subroutine ExtOrder
 363
 364 subroutine ExtOrderStr(MemoryOrder,Vector,ROVector,Status)
 365   implicit none
 366   include 'wrf_status_codes.h'
 367   character*(*)                    ,intent(in)    :: MemoryOrder
 368   character*(*),dimension(*)       ,intent(in)    :: Vector
 369   character(80),dimension(NVarDims),intent(out)   :: ROVector
 370   integer                          ,intent(out)   :: Status
 371   integer                                         :: NDim
 372   character*3                                     :: MemOrd
 373
 374   call GetDim(MemoryOrder,NDim,Status)
 375   ROVector(1:NDim) = Vector(1:NDim)
 376   call LowerCase(MemoryOrder,MemOrd)
 377   select case (MemOrd)
 378
 379     case ('xyz','xsz','xez','ysz','yez','xy','xs','xe','ys','ye','z','c')
 380       continue
 381     case ('0')
 382       continue  ! NDim=0 for scalars.  TBH:  20060502
 383     case ('xzy')
 384       ROVector(2) = Vector(3)
 385       ROVector(3) = Vector(2)
 386     case ('yxz')
 387       ROVector(1) = Vector(2)
 388       ROVector(2) = Vector(1)
 389     case ('yzx')
 390       ROVector(1) = Vector(3)
 391       ROVector(2) = Vector(1)
 392       ROVector(3) = Vector(2)
 393     case ('zxy')
 394       ROVector(1) = Vector(2)
 395       ROVector(2) = Vector(3)
 396       ROVector(3) = Vector(1)
 397     case ('zyx')
 398       ROVector(1) = Vector(3)
 399       ROVector(3) = Vector(1)
 400     case ('yx')
 401       ROVector(1) = Vector(2)
 402       ROVector(2) = Vector(1)
 403     case default
 404       Status = WRF_WARN_BAD_MEMORYORDER
 405       return
 406   end select
 407   Status = WRF_NO_ERR
 408   return
 409 end subroutine ExtOrderStr
 410
 411
 412 subroutine LowerCase(MemoryOrder,MemOrd)
 413   implicit none
 414   include 'wrf_status_codes.h'
 415   character*(*) ,intent(in)  :: MemoryOrder
 416   character*(*) ,intent(out) :: MemOrd
 417   character*1                :: c
 418   integer       ,parameter   :: upper_to_lower =IACHAR('a')-IACHAR('A')
 419   integer                    :: i,N
 420
 421   MemOrd = ' '
 422   N = len(MemoryOrder)
 423   MemOrd(1:N) = MemoryOrder(1:N)
 424   do i=1,N
 425     c = MemoryOrder(i:i)
 426     if('A'<=c .and. c <='Z') MemOrd(i:i)=achar(iachar(c)+upper_to_lower)
 427   enddo
 428   return
 429 end subroutine LowerCase
 430
 431 subroutine UpperCase(MemoryOrder,MemOrd)
 432   implicit none
 433   include 'wrf_status_codes.h'
 434   character*(*) ,intent(in)  :: MemoryOrder
 435   character*(*) ,intent(out) :: MemOrd
 436   character*1                :: c
 437   integer     ,parameter     :: lower_to_upper =IACHAR('A')-IACHAR('a')
 438   integer                    :: i,N
 439
 440   MemOrd = ' '
 441   N = len(MemoryOrder)
 442   MemOrd(1:N) = MemoryOrder(1:N)
 443   do i=1,N
 444     c = MemoryOrder(i:i)
 445     if('a'<=c .and. c <='z') MemOrd(i:i)=achar(iachar(c)+lower_to_upper)
 446   enddo
 447   return
 448 end subroutine UpperCase
 449
 450 subroutine netcdf_err(err,Status)
 451   implicit none
 452   include 'wrf_status_codes.h'
 453   integer  ,intent(in)  :: err
 454   integer  ,intent(out) :: Status
 455   character(len=80)     :: errmsg
 456   integer               :: stat
 457
 458   if(err == PIO_NOERR)then
 459     Status = WRF_NO_ERR
 460   else
 461     write(msg,*) 'NetCDF error: ', 'from PIO'
 462     call wrf_debug ( WARN , TRIM(msg))
 463     Status = WRF_WARN_NETCDF
 464   endif
 465   return
 466 end subroutine netcdf_err
 467
 468 subroutine find_iodesc(DH,MemoryOrder,Stagger,FieldTYpe,whole)
 469   implicit none
 470   type(wrf_data_handle), pointer   :: DH
 471   character*(*),     intent(in)    :: MemoryOrder
 472   character*(*),     intent(in)    :: Stagger
 473   integer,           intent(in)    :: FieldType
 474   logical,           intent(out)   :: whole
 475   character*3                      :: MemOrd
 476   character*1                      :: Stag
 477   integer           ,parameter     :: MaxUpperCase=IACHAR('Z')
 478
 479   whole = .false.
 480
 481   call LowerCase(MemoryOrder,MemOrd)
 482   call LowerCase(Stagger,Stag)
 483
 484   select case (MemOrd)
 485     case ('xzy')
 486       select case (FieldType)
 487         case (WRF_REAL)
 488           select case (Stag)
 489             case ('x')
 490                  DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_u_real
 491             case ('y')
 492                  DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_v_real
 493             case ('z')
 494                  if(LAND_CAT_VAR == DH%vartype(DH%CurrentVariable)) then
 495                     DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_land_real
 496                  else if(SOIL_CAT_VAR == DH%vartype(DH%CurrentVariable)) then
 497                     DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_soil_real
 498                  else if(SOIL_LAYERS_VAR == DH%vartype(DH%CurrentVariable)) then
 499                     DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_soil_layers_real
 500                  else if(MDL_CPL_VAR == DH%vartype(DH%CurrentVariable)) then
 501                     DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_mdl_cpl_real
 502                  else if(ENSEMBLE_VAR == DH%vartype(DH%CurrentVariable)) then
 503                     DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_ensemble_real
 504                  else
 505                     DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_w_real
 506                  endif
 507             case default
 508                  DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_m_real
 509           end select
 510         case (WRF_DOUBLE)
 511           select case (Stag)
 512             case ('x')
 513                  DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_u_double
 514             case ('y')
 515                  DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_v_double
 516             case ('z')
 517                  if(LAND_CAT_VAR == DH%vartype(DH%CurrentVariable)) then
 518                     DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_land_double
 519                  else if(SOIL_CAT_VAR == DH%vartype(DH%CurrentVariable)) then
 520                     DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_soil_double
 521                  else if(SOIL_LAYERS_VAR == DH%vartype(DH%CurrentVariable)) then
 522                     DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_soil_layers_double
 523                  else if(MDL_CPL_VAR == DH%vartype(DH%CurrentVariable)) then
 524                     DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_mdl_cpl_double
 525                  else if(ENSEMBLE_VAR == DH%vartype(DH%CurrentVariable)) then
 526                     DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_ensemble_double
 527                  else
 528                     DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_w_double
 529                  endif
 530             case default
 531                  DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_m_double
 532           end select
 533         case (WRF_INTEGER)
 534           select case (Stag)
 535             case ('x')
 536                  DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_u_int
 537             case ('y')
 538                  DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_v_int
 539             case ('z')
 540                  if(LAND_CAT_VAR == DH%vartype(DH%CurrentVariable)) then
 541                     DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_land_int
 542                  else if(SOIL_CAT_VAR == DH%vartype(DH%CurrentVariable)) then
 543                     DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_soil_int
 544                  else if(SOIL_LAYERS_VAR == DH%vartype(DH%CurrentVariable)) then
 545                     DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_soil_layers_int
 546                  else if(MDL_CPL_VAR == DH%vartype(DH%CurrentVariable)) then
 547                     DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_mdl_cpl_int
 548                  else if(ENSEMBLE_VAR == DH%vartype(DH%CurrentVariable)) then
 549                     DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_ensemble_int
 550                  else
 551                     DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_w_int
 552                  endif
 553             case default
 554                  DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_m_int
 555           end select
 556         case (WRF_LOGICAL)
 557           select case (Stag)
 558             case ('x')
 559                  DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_u_int
 560             case ('y')
 561                  DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_v_int
 562             case ('z')
 563                  if(LAND_CAT_VAR == DH%vartype(DH%CurrentVariable)) then
 564                     DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_land_int
 565                  else if(SOIL_CAT_VAR == DH%vartype(DH%CurrentVariable)) then
 566                     DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_soil_int
 567                  else if(SOIL_LAYERS_VAR == DH%vartype(DH%CurrentVariable)) then
 568                     DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_soil_layers_int
 569                  else if(MDL_CPL_VAR == DH%vartype(DH%CurrentVariable)) then
 570                     DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_mdl_cpl_int
 571                  else if(ENSEMBLE_VAR == DH%vartype(DH%CurrentVariable)) then
 572                     DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_ensemble_int
 573                  else
 574                     DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_w_int
 575                  endif
 576             case default
 577                  DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_m_int
 578           end select
 579         case default
 580           write(msg,*) 'Warning DATA TYPE NOT FOUND in ',__FILE__,', line', __LINE__
 581           call wrf_debug ( WARN , TRIM(msg))
 582           whole = .true.
 583           return
 584       end select
 585     case ('xy')
 586       select case (FieldType)
 587         case (WRF_REAL)
 588           select case (Stag)
 589             case ('x')
 590                  DH%ioVar(DH%CurrentVariable) = DH%iodesc2d_u_real
 591             case ('y')
 592                  DH%ioVar(DH%CurrentVariable) = DH%iodesc2d_v_real
 593             case default
 594                  DH%ioVar(DH%CurrentVariable) = DH%iodesc2d_m_real
 595           end select
 596         case (WRF_DOUBLE)
 597           select case (Stag)
 598             case ('x')
 599                  DH%ioVar(DH%CurrentVariable) = DH%iodesc2d_u_double
 600             case ('y')
 601                  DH%ioVar(DH%CurrentVariable) = DH%iodesc2d_v_double
 602             case default
 603                  DH%ioVar(DH%CurrentVariable) = DH%iodesc2d_m_double
 604           end select
 605         case (WRF_INTEGER)
 606           select case (Stag)
 607             case ('x')
 608                  DH%ioVar(DH%CurrentVariable) = DH%iodesc2d_u_int
 609             case ('y')
 610                  DH%ioVar(DH%CurrentVariable) = DH%iodesc2d_v_int
 611             case default
 612                  DH%ioVar(DH%CurrentVariable) = DH%iodesc2d_m_int
 613           end select
 614         case (WRF_LOGICAL)
 615           select case (Stag)
 616             case ('x')
 617                  DH%ioVar(DH%CurrentVariable) = DH%iodesc2d_u_int
 618             case ('y')
 619                  DH%ioVar(DH%CurrentVariable) = DH%iodesc2d_v_int
 620             case default
 621                  DH%ioVar(DH%CurrentVariable) = DH%iodesc2d_m_int
 622           end select
 623         case default
 624           write(msg,*) 'Warning DATA TYPE NOT FOUND in ',__FILE__,', line', __LINE__
 625           call wrf_debug ( WARN , TRIM(msg))
 626           whole = .true.
 627           return
 628       end select
 629     case ('xsz')
 630       DH%vartype(DH%CurrentVariable) = BDY_VAR
 631       select case (FieldType)
 632         case (WRF_REAL)
 633           select case (Stag)
 634             case ('x')
 635                  DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_xsz_u_real
 636             case ('y')
 637                  DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_xsz_v_real
 638             case ('z')
 639                  DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_xsz_w_real
 640             case default
 641                  DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_xsz_m_real
 642           end select
 643         case (WRF_DOUBLE)
 644           select case (Stag)
 645             case ('x')
 646                  DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_xsz_u_double
 647             case ('y')
 648                  DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_xsz_v_double
 649             case ('z')
 650                  DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_xsz_w_double
 651             case default
 652                  DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_xsz_m_double
 653           end select
 654         case (WRF_INTEGER)
 655           select case (Stag)
 656             case ('x')
 657                  DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_xsz_u_int
 658             case ('y')
 659                  DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_xsz_v_int
 660             case ('z')
 661                  DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_xsz_w_int
 662             case default
 663                  DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_xsz_m_int
 664           end select
 665         case (WRF_LOGICAL)
 666           select case (Stag)
 667             case ('x')
 668                  DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_xsz_u_int
 669             case ('y')
 670                  DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_xsz_v_int
 671             case ('z')
 672                  DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_xsz_w_int
 673             case default
 674                  DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_xsz_m_int
 675           end select
 676         case default
 677           write(msg,*) 'Warning DATA TYPE NOT FOUND in ',__FILE__,', line', __LINE__
 678           call wrf_debug ( WARN , TRIM(msg))
 679           whole = .true.
 680           return
 681       end select
 682     case ('xez')
 683       DH%vartype(DH%CurrentVariable) = BDY_VAR
 684       select case (FieldType)
 685         case (WRF_REAL)
 686           select case (Stag)
 687             case ('x')
 688                  DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_xez_u_real
 689             case ('y')
 690                  DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_xez_v_real
 691             case ('z')
 692                  DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_xez_w_real
 693             case default
 694                  DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_xez_m_real
 695           end select
 696         case (WRF_DOUBLE)
 697           select case (Stag)
 698             case ('x')
 699                  DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_xez_u_double
 700             case ('y')
 701                  DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_xez_v_double
 702             case ('z')
 703                  DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_xez_w_double
 704             case default
 705                  DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_xez_m_double
 706           end select
 707         case (WRF_INTEGER)
 708           select case (Stag)
 709             case ('x')
 710                  DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_xez_u_int
 711             case ('y')
 712                  DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_xez_v_int
 713             case ('z')
 714                  DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_xez_w_int
 715             case default
 716                  DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_xez_m_int
 717           end select
 718         case (WRF_LOGICAL)
 719           select case (Stag)
 720             case ('x')
 721                  DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_xez_u_int
 722             case ('y')
 723                  DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_xez_v_int
 724             case ('z')
 725                  DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_xez_w_int
 726             case default
 727                  DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_xez_m_int
 728           end select
 729         case default
 730           write(msg,*) 'Warning DATA TYPE NOT FOUND in ',__FILE__,', line', __LINE__
 731           call wrf_debug ( WARN , TRIM(msg))
 732           whole = .true.
 733           return
 734       end select
 735     case ('ysz')
 736       DH%vartype(DH%CurrentVariable) = BDY_VAR
 737       select case (FieldType)
 738         case (WRF_REAL)
 739           select case (Stag)
 740             case ('x')
 741                  DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_ysz_u_real
 742             case ('y')
 743                  DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_ysz_v_real
 744             case ('z')
 745                  DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_ysz_w_real
 746             case default
 747                  DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_ysz_m_real
 748           end select
 749         case (WRF_DOUBLE)
 750           select case (Stag)
 751             case ('x')
 752                  DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_ysz_u_double
 753             case ('y')
 754                  DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_ysz_v_double
 755             case ('z')
 756                  DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_ysz_w_double
 757             case default
 758                  DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_ysz_m_double
 759           end select
 760         case (WRF_INTEGER)
 761           select case (Stag)
 762             case ('x')
 763                  DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_ysz_u_int
 764             case ('y')
 765                  DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_ysz_v_int
 766             case ('z')
 767                  DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_ysz_w_int
 768             case default
 769                  DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_ysz_m_int
 770           end select
 771         case (WRF_LOGICAL)
 772           select case (Stag)
 773             case ('x')
 774                  DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_ysz_u_int
 775             case ('y')
 776                  DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_ysz_v_int
 777             case ('z')
 778                  DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_ysz_w_int
 779             case default
 780                  DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_ysz_m_int
 781           end select
 782         case default
 783           write(msg,*) 'Warning DATA TYPE NOT FOUND in ',__FILE__,', line', __LINE__
 784           call wrf_debug ( WARN , TRIM(msg))
 785           whole = .true.
 786           return
 787       end select
 788     case ('yez')
 789       DH%vartype(DH%CurrentVariable) = BDY_VAR
 790       select case (FieldType)
 791         case (WRF_REAL)
 792           select case (Stag)
 793             case ('x')
 794                  DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_yez_u_real
 795             case ('y')
 796                  DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_yez_v_real
 797             case ('z')
 798                  DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_yez_w_real
 799             case default
 800                  DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_yez_m_real
 801           end select
 802         case (WRF_DOUBLE)
 803           select case (Stag)
 804             case ('x')
 805                  DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_yez_u_double
 806             case ('y')
 807                  DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_yez_v_double
 808             case ('z')
 809                  DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_yez_w_double
 810             case default
 811                  DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_yez_m_double
 812           end select
 813         case (WRF_INTEGER)
 814           select case (Stag)
 815             case ('x')
 816                  DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_yez_u_int
 817             case ('y')
 818                  DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_yez_v_int
 819             case ('z')
 820                  DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_yez_w_int
 821             case default
 822                  DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_yez_m_int
 823           end select
 824         case (WRF_LOGICAL)
 825           select case (Stag)
 826             case ('x')
 827                  DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_yez_u_int
 828             case ('y')
 829                  DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_yez_v_int
 830             case ('z')
 831                  DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_yez_w_int
 832             case default
 833                  DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_yez_m_int
 834           end select
 835         case default
 836           write(msg,*) 'Warning DATA TYPE NOT FOUND in ',__FILE__,', line', __LINE__
 837           call wrf_debug ( WARN , TRIM(msg))
 838           whole = .true.
 839           return
 840       end select
 841     case ('xs')
 842       DH%vartype(DH%CurrentVariable) = BDY_VAR
 843       select case (FieldType)
 844         case (WRF_REAL)
 845              DH%ioVar(DH%CurrentVariable) = DH%iodesc2d_xs_m_real
 846         case (WRF_DOUBLE)
 847              DH%ioVar(DH%CurrentVariable) = DH%iodesc2d_xs_m_double
 848         case (WRF_INTEGER)
 849              DH%ioVar(DH%CurrentVariable) = DH%iodesc2d_xs_m_int
 850         case (WRF_LOGICAL)
 851              DH%ioVar(DH%CurrentVariable) = DH%iodesc2d_xs_m_int
 852         case default
 853           write(msg,*) 'Warning DATA TYPE NOT FOUND in ',__FILE__,', line', __LINE__
 854           call wrf_debug ( WARN , TRIM(msg))
 855           whole = .true.
 856           return
 857       end select
 858     case ('xe')
 859       DH%vartype(DH%CurrentVariable) = BDY_VAR
 860       select case (FieldType)
 861         case (WRF_REAL)
 862              DH%ioVar(DH%CurrentVariable) = DH%iodesc2d_xe_m_real
 863         case (WRF_DOUBLE)
 864              DH%ioVar(DH%CurrentVariable) = DH%iodesc2d_xe_m_double
 865         case (WRF_INTEGER)
 866              DH%ioVar(DH%CurrentVariable) = DH%iodesc2d_xe_m_int
 867         case (WRF_LOGICAL)
 868              DH%ioVar(DH%CurrentVariable) = DH%iodesc2d_xe_m_int
 869         case default
 870           write(msg,*) 'Warning DATA TYPE NOT FOUND in ',__FILE__,', line', __LINE__
 871           call wrf_debug ( WARN , TRIM(msg))
 872           whole = .true.
 873           return
 874       end select
 875     case ('ys')
 876       DH%vartype(DH%CurrentVariable) = BDY_VAR
 877       select case (FieldType)
 878         case (WRF_REAL)
 879              DH%ioVar(DH%CurrentVariable) = DH%iodesc2d_ys_m_real
 880         case (WRF_DOUBLE)
 881              DH%ioVar(DH%CurrentVariable) = DH%iodesc2d_ys_m_double
 882         case (WRF_INTEGER)
 883              DH%ioVar(DH%CurrentVariable) = DH%iodesc2d_ys_m_int
 884         case (WRF_LOGICAL)
 885              DH%ioVar(DH%CurrentVariable) = DH%iodesc2d_ys_m_int
 886         case default
 887           write(msg,*) 'Warning DATA TYPE NOT FOUND in ',__FILE__,', line', __LINE__
 888           call wrf_debug ( WARN , TRIM(msg))
 889           whole = .true.
 890           return
 891       end select
 892     case ('ye')
 893       DH%vartype(DH%CurrentVariable) = BDY_VAR
 894       select case (FieldType)
 895         case (WRF_REAL)
 896              DH%ioVar(DH%CurrentVariable) = DH%iodesc2d_ye_m_real
 897         case (WRF_DOUBLE)
 898              DH%ioVar(DH%CurrentVariable) = DH%iodesc2d_ye_m_double
 899         case (WRF_INTEGER)
 900              DH%ioVar(DH%CurrentVariable) = DH%iodesc2d_ye_m_int
 901         case (WRF_LOGICAL)
 902              DH%ioVar(DH%CurrentVariable) = DH%iodesc2d_ye_m_int
 903         case default
 904           write(msg,*) 'Warning DATA TYPE NOT FOUND in ',__FILE__,', line', __LINE__
 905           call wrf_debug ( WARN , TRIM(msg))
 906           whole = .true.
 907           return
 908       end select
 909     case ('xyz')
 910       select case (Stag)
 911         case ('z')
 912              if(ENSEMBLE_VAR == DH%vartype(DH%CurrentVariable)) then
 913                 select case (FieldType)
 914                   case (WRF_REAL)
 915                        DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_ensemble_real
 916                   case (WRF_DOUBLE)
 917                        DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_ensemble_double
 918                   case (WRF_INTEGER)
 919                        DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_ensemble_int
 920                   case (WRF_LOGICAL)
 921                        DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_ensemble_int
 922                   case default
 923                        write(msg,*) 'Warning DO NOT KNOW HOW TO HANDLE this FieldType in ',__FILE__,', line', __LINE__
 924                        call wrf_debug ( WARN , TRIM(msg))
 925                        DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_m_real
 926                 end select
 927              else
 928                 write(msg,*) 'Warning DO NOT KNOW HOW TO HANDLE THIS VAR KIND in ',__FILE__,', line', __LINE__
 929                 call wrf_debug ( WARN , TRIM(msg))
 930                 DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_w_double
 931              end if
 932         case default
 933              write(msg,*) 'Warning DO NOT KNOW HOW TO HANDLE THIS STAG in ',__FILE__,', line', __LINE__
 934              call wrf_debug ( WARN , TRIM(msg))
 935              DH%ioVar(DH%CurrentVariable) = DH%iodesc3d_m_real
 936       end select
 937    !case ('z','c')
 938    !  whole = .true.
 939     case default
 940       whole = .true.
 941       if(ENSEMBLE_VAR == DH%vartype(DH%CurrentVariable)) then
 942          whole = .false.
 943       end if
 944 #if 0
 945       select case (FieldType)
 946         case (WRF_REAL)
 947              DH%ioVar(DH%CurrentVariable) = DH%iodesc1d_real
 948         case (WRF_DOUBLE)
 949              DH%ioVar(DH%CurrentVariable) = DH%iodesc1d_double
 950         case (WRF_INTEGER)
 951              DH%ioVar(DH%CurrentVariable) = DH%iodesc1d_int
 952         case (WRF_LOGICAL)
 953              DH%ioVar(DH%CurrentVariable) = DH%iodesc1d_int
 954         case default
 955              write(msg,*) 'Warning DATA TYPE NOT FOUND in ',__FILE__,', line', __LINE__
 956              call wrf_debug ( WARN , TRIM(msg))
 957              return
 958       end select
 959 #endif
 960   end select
 961 end subroutine find_iodesc
 962
 963 logical function is_boundary(MemoryOrder)
 964
 965   implicit none
 966
 967   character*(*), intent(in) :: MemoryOrder
 968
 969   logical     :: isbdy
 970   character*3 :: MemOrd
 971
 972   isbdy = .false.
 973
 974   call LowerCase(MemoryOrder,MemOrd)
 975
 976   select case (MemOrd)
 977     case ('xsz', 'xez', 'ysz', 'yez')
 978       isbdy = .true.
 979     case ('xs', 'xe', 'ys', 'ye')
 980       isbdy = .true.
 981     case default
 982       isbdy = .false.
 983   end select
 984
 985   is_boundary = isbdy
 986 end function is_boundary
 987
 988 subroutine FieldIO(IO,DataHandle,DateStr,Dimens,Starts,Counts,Length,MemoryOrder, &
 989                    Stagger,FieldType,Field,Status)
 990   implicit none
 991   include 'wrf_status_codes.h'
 992   character (*)              ,intent(in)    :: IO
 993   integer                    ,intent(in)    :: DataHandle
 994   character*(*)              ,intent(in)    :: DateStr
 995   integer,dimension(NVarDims),intent(inout) :: Dimens
 996   integer,dimension(NVarDims),intent(inout) :: Starts
 997   integer,dimension(NVarDims),intent(inout) :: Counts
 998   integer,dimension(NVarDims),intent(in)    :: Length
 999   character*(*)              ,intent(in)    :: MemoryOrder
1000   character*(*)              ,intent(in)    :: Stagger
1001   integer                    ,intent(in)    :: FieldType
1002   integer,dimension(*)       ,intent(inout) :: Field
1003   integer                    ,intent(out)   :: Status
1004   integer                                   :: TimeIndex
1005   logical                                   :: whole, isbdy
1006   integer                                   :: NDim
1007   integer                                   :: fldsize, datasize
1008   integer                                   :: n
1009   type(wrf_data_handle)      ,pointer       :: DH
1010   integer(KIND=PIO_OFFSET)                  :: pioidx
1011
1012   DH => WrfDataHandles(DataHandle)
1013   call GetTimeIndex(IO,DataHandle,DateStr,TimeIndex,Status)
1014   if(Status /= WRF_NO_ERR) then
1015     write(msg,*) 'Warning in ',__FILE__,', line', __LINE__
1016     call wrf_debug ( WARN , TRIM(msg))
1017     write(msg,*) '  Bad time index for DateStr = ',DateStr
1018     call wrf_debug ( WARN , TRIM(msg))
1019     return
1020   endif
1021   call GetDim(MemoryOrder,NDim,Status)
1022
1023   fldsize = 1
1024   datasize = 1
1025   do n = 1, NDim
1026      fldsize = fldsize * Length(n)
1027      datasize = datasize * Counts(n)
1028   end do
1029
1030   Starts(NDim+1) = TimeIndex
1031   Counts(NDim+1) = 1
1032
1033   call find_iodesc(DH,MemoryOrder,Stagger,FieldTYpe,whole)
1034   isbdy = is_boundary(MemoryOrder)
1035  !isbdy = BDY_VAR == DH%vartype(DH%CurrentVariable)
1036
1037   pioidx = TimeIndex
1038   call pio_setframe(DH%descVar(DH%CurrentVariable), pioidx)
1039  !DH%descVar(DH%CurrentVariable)%rec = TimeIndex
1040
1041  !write(unit=0, fmt='(3a,i6)') 'File: ', __FILE__, ', line: ', __LINE__
1042  !write(unit=0, fmt='(3a,l8)') 'IO = ', trim(IO), ', whole = ', whole
1043  !write(unit=0, fmt='(4a)') 'MemoryOrder = ', trim(MemoryOrder), ', Stagger = ', trim(Stagger)
1044  !write(unit=0, fmt='(a,i4,a,i3)') 'DH%vartype(', DH%CurrentVariable, ') = ', DH%vartype(DH%CurrentVariable)
1045
1046  !if(whole .and. (ENSEMBLE_VAR == DH%vartype(DH%CurrentVariable))) then
1047  !   whole = .false.
1048  !end if
1049
1050   select case (FieldType)
1051     case (WRF_REAL)
1052       if(isbdy .and. (IO == 'read')) then
1053         Dimens(NDim+1) = TimeIndex
1054         call read_bdy_RealFieldIO(DH,NDim,Dimens,Starts,Counts,Field,Status)
1055       else
1056         call ext_pio_RealFieldIO(whole,IO,DH,Starts,Counts,fldsize,datasize,Field,Status)
1057       endif
1058     case (WRF_DOUBLE)
1059       if(isbdy .and. (IO == 'read')) then
1060         Dimens(NDim+1) = TimeIndex
1061         call read_bdy_DoubleFieldIO(DH,NDim,Dimens,Starts,Counts,Field,Status)
1062       else
1063         call ext_pio_DoubleFieldIO(whole,IO,DH,Starts,Counts,fldsize,datasize,Field,Status)
1064       endif
1065     case (WRF_INTEGER)
1066       call ext_pio_IntFieldIO(whole,IO,DH,Starts,Counts,fldsize,datasize,Field,Status)
1067     case (WRF_LOGICAL)
1068       call ext_pio_LogicalFieldIO(whole,IO,DH,Starts,Counts,fldsize,datasize,Field,Status)
1069     case default
1070       Status = WRF_WARN_DATA_TYPE_NOT_FOUND
1071       write(msg,*) 'Warning DATA TYPE NOT FOUND in ',__FILE__,', line', __LINE__
1072       call wrf_debug ( WARN , TRIM(msg))
1073       return
1074   end select
1075
1076   return
1077 end subroutine FieldIO
1078
1079 subroutine FieldBDY(IO,DataHandle,DateStr,NDim,Domains, &
1080                     MemoryStart,MemoryEnd,PatchStart,PatchEnd, &
1081                     FieldType,Field,Status)
1082   implicit none
1083   include 'wrf_status_codes.h'
1084   character (*)              ,intent(in)    :: IO
1085   integer                    ,intent(in)    :: DataHandle,NDim
1086   character*(*)              ,intent(in)    :: DateStr
1087   integer,dimension(*)       ,intent(inout) :: Domains
1088   integer,dimension(*)       ,intent(in)    :: MemoryStart, MemoryEnd
1089   integer,dimension(*)       ,intent(in)    :: PatchStart,  PatchEnd
1090   integer                    ,intent(in)    :: FieldType
1091   integer,dimension(*)       ,intent(inout) :: Field
1092   integer                    ,intent(out)   :: Status
1093   integer                                   :: TimeIndex
1094   type(wrf_data_handle)      ,pointer       :: DH
1095   integer(KIND=PIO_OFFSET)                  :: pioidx
1096
1097   DH => WrfDataHandles(DataHandle)
1098   call GetTimeIndex(IO,DataHandle,DateStr,TimeIndex,Status)
1099   if(Status /= WRF_NO_ERR) then
1100     write(msg,*) 'Warning in ',__FILE__,', line', __LINE__
1101     call wrf_debug ( WARN , TRIM(msg))
1102     write(msg,*) '  Bad time index for DateStr = ',DateStr
1103     call wrf_debug ( WARN , TRIM(msg))
1104     return
1105   endif
1106
1107   pioidx = TimeIndex
1108   call pio_setframe(DH%descVar(DH%CurrentVariable), pioidx)
1109  !DH%descVar(DH%CurrentVariable)%rec = TimeIndex
1110   Domains(NDim+1) = TimeIndex
1111
1112   select case (FieldType)
1113     case (WRF_REAL)
1114          call read_bdy_RealFieldIO(DH,NDim,Domains,MemoryStart,MemoryEnd,PatchStart,PatchEnd,Field,Status)
1115     case (WRF_DOUBLE)
1116          call read_bdy_DoubleFieldIO(DH,NDim,Domains,MemoryStart,MemoryEnd,PatchStart,PatchEnd,Field,Status)
1117     case default
1118          Status = WRF_WARN_DATA_TYPE_NOT_FOUND
1119          write(msg,*) 'Warning DATA TYPE NOT FOUND in ',__FILE__,', line', __LINE__
1120          call wrf_debug ( WARN , TRIM(msg))
1121          return
1122   end select
1123
1124   return
1125 end subroutine FieldBDY
1126
1127 ! Returns .TRUE. iff it is OK to write time-independent domain metadata to the
1128 ! file referenced by DataHandle.  If DataHandle is invalid, .FALSE. is
1129 ! returned.
1130 LOGICAL FUNCTION ncd_ok_to_put_dom_ti( DataHandle )
1131     implicit none
1132     include 'wrf_status_codes.h'
1133     INTEGER, INTENT(IN) :: DataHandle
1134     CHARACTER*80 :: fname
1135     INTEGER :: filestate
1136     INTEGER :: Status
1137     LOGICAL :: dryrun, retval
1138     call ext_pio_inquire_filename( DataHandle, fname, filestate, Status )
1139     IF ( Status /= WRF_NO_ERR ) THEN
1140       write(msg,*) 'Warning Status = ',Status,' in ',__FILE__, &
1141                    ', line', __LINE__
1142       call wrf_debug ( WARN , TRIM(msg) )
1143       retval = .FALSE.
1144     ELSE
1145       dryrun = ( filestate .EQ. WRF_FILE_OPENED_NOT_COMMITTED )
1146       retval = dryrun
1147     ENDIF
1148     ncd_ok_to_put_dom_ti = retval
1149     RETURN
1150 END FUNCTION ncd_ok_to_put_dom_ti
1151
1152 ! Returns .TRUE. iff it is OK to read time-independent domain metadata from the
1153 ! file referenced by DataHandle.  If DataHandle is invalid, .FALSE. is
1154 ! returned.
1155 LOGICAL FUNCTION ncd_ok_to_get_dom_ti( DataHandle )
1156     implicit none
1157     include 'wrf_status_codes.h'
1158     INTEGER, INTENT(IN) :: DataHandle
1159     CHARACTER*80 :: fname
1160     INTEGER :: filestate
1161     INTEGER :: Status
1162     LOGICAL :: dryrun, retval
1163     call ext_pio_inquire_filename( DataHandle, fname, filestate, Status )
1164     IF ( Status /= WRF_NO_ERR ) THEN
1165       write(msg,*) 'Warning Status = ',Status,' in ',__FILE__, &
1166                    ', line', __LINE__
1167       call wrf_debug ( WARN , TRIM(msg) )
1168       retval = .FALSE.
1169     ELSE
1170       dryrun = ( filestate .EQ. WRF_FILE_OPENED_NOT_COMMITTED )
1171       retval = .NOT. dryrun
1172     ENDIF
1173     ncd_ok_to_get_dom_ti = retval
1174     RETURN
1175 END FUNCTION ncd_ok_to_get_dom_ti
1176
1177 subroutine initialize_pio(grid, DH)
1178    implicit none
1179
1180    type(domain)                   :: grid
1181    type(wrf_data_handle), pointer :: DH
1182
1183    integer     :: ierr
1184    integer(i4) :: communicator, pioprocs, piostart, piostride, pioshift
1185
1186    communicator = grid%communicator
1187
1188    if(.not. associated(DH%iosystem)) then
1189       allocate(DH%iosystem)
1190    end if
1191
1192    DH%Write = 0
1193
1194   !call pio_setdebuglevel(1)
1195
1196    call mpi_comm_size(communicator, nprocs, ierr)
1197    call mpi_comm_rank(communicator, myrank, ierr)
1198
1199    if(grid%pioprocs > nprocs) then
1200      !Force pioprocs to be nprocs.
1201       pioprocs = nprocs
1202    else if(grid%pioprocs < 1) then
1203      !Force pioprocs to be 1.
1204       pioprocs = 1
1205    else
1206       pioprocs = grid%pioprocs
1207    endif
1208
1209    piostride = nprocs / grid%pioprocs
1210
1211    if((grid%pioprocs * piostride) < nprocs) then
1212      !We expect that: nprocs = piostride * grid%pioprocs
1213       piostride = piostride + 1
1214    endif
1215
1216    if(piostride /= grid%piostride) then
1217      !We expect that user's piostride equals what we calculated here.
1218      !If not, override it.
1219       write(unit=0, fmt='(3a,i6)') 'file: ', __FILE__, ', line: ', __LINE__
1220       write(unit=0, fmt='(a,i6)') 'Calculated piostride = ', piostride, &
1221                                   'User provided piostride = ', grid%piostride
1222    endif
1223
1224    if(grid%pioshift < 0) then
1225      !pioshift can from 0, but can not less than 0, usually, we
1226       if(grid%piostride > 1) then
1227          pioshift = 1
1228       else
1229          pioshift = 0
1230       endif
1231       write(unit=0, fmt='(3a,i6)') 'file: ', __FILE__, ', line: ', __LINE__
1232       write(unit=0, fmt='(a,i6)') 'PIO has forced pioshift to: ', pioshift
1233    else if(grid%pioshift >= grid%piostride) then
1234      !pioshift can not large then piostride
1235       write(unit=0, fmt='(3a,i6)') 'file: ', __FILE__, ', line: ', __LINE__
1236       write(unit=0, fmt='(a,i6)') 'User provided a pioshift of: ', grid%pioshift
1237       if(grid%piostride > 1) then
1238          pioshift = 1
1239       else
1240          pioshift = 0
1241       endif
1242       write(unit=0, fmt='(3a,i6)') 'file: ', __FILE__, ', line: ', __LINE__
1243       write(unit=0, fmt='(a,i6)') 'PIO has forced pioshift to: ', pioshift
1244    else
1245       pioshift = grid%pioshift
1246    endif
1247
1248    if(grid%piostart < 0) then
1249      !Force piostart from 0
1250       write(unit=0, fmt='(3a,i6)') 'file: ', __FILE__, ', line: ', __LINE__
1251       write(unit=0, fmt='(a,i6)') 'User provided a piostart of: ', grid%piostart
1252       write(unit=0, fmt='(a,i6)') 'PIO has forced piosstart to: ', 0
1253       piostart = 0
1254    else
1255       piostart = grid%piostart
1256    endif
1257
1258   !write(unit=0, fmt='(3a,i6)') 'file: ', __FILE__, ', line: ', __LINE__
1259   !write(unit=0, fmt='(2(a,i6))') 'nprocs = ', nprocs, ', myrank = ', myrank
1260   !write(unit=0, fmt='(4(a,i6))') 'pioprocs = ', pioprocs, &
1261   !                             ', piostride = ', piostride, &
1262   !                             ', piostart = ', piostart, &
1263   !                             ', pioshift = ', pioshift
1264
1265   !call PIO_init to initiate iosystem
1266   !call PIO_init(my_rank, MPI_COMM_WORLD, 4, 0, 4, PIO_rearr_box, iosystem, 1)
1267   !call PIO_init(myrank, MPI_COMM_WORLD, pioprocs, &
1268
1269    call PIO_init(myrank, communicator, pioprocs, &
1270                  piostart, piostride, &
1271                  PIO_rearr_box, DH%iosystem, pioshift)
1272
1273    DH%nprocs = nprocs
1274    DH%myrank = myrank
1275
1276    DH%piostart = piostart
1277    DH%pioshift = pioshift
1278    DH%pioprocs = pioprocs
1279    DH%piostride = piostride
1280
1281    call get_ijk_from_grid(grid,                         &
1282                           ids, ide, jds, jde, kds, kde, &
1283                           ims, ime, jms, jme, kms, kme, &
1284                           its, ite, jts, jte, kts, kte)
1285
1286 end subroutine initialize_pio
1287
1288 subroutine define_pio_iodesc(grid, DH)
1289    implicit none
1290
1291    type(domain)                   :: grid
1292    type(wrf_data_handle), pointer :: DH
1293
1294    integer(i4) :: communicator, myrank
1295    integer(i4) :: iostat
1296
1297    integer(kind=PIO_Offset), &
1298            dimension((ime - ims + 1) * (jme - jms + 1) * (kme - kms + 1)) &
1299            :: compdof_3d
1300    integer(kind=PIO_Offset), &
1301            dimension((ime - ims + 1) * (jme - jms + 1) * grid%num_land_cat) &
1302            :: compdof_3d_land
1303    integer(kind=PIO_Offset), &
1304            dimension((ime - ims + 1) * (jme - jms + 1) * grid%num_soil_cat) &
1305            :: compdof_3d_soil
1306    integer(kind=PIO_Offset), &
1307            dimension((ime - ims + 1) * (jme - jms + 1) * grid%num_soil_layers) &
1308            :: compdof_3d_soil_layers
1309    integer(kind=PIO_Offset), &
1310            dimension((ime - ims + 1) * (jme - jms + 1) * grid%num_ext_model_couple_dom) &
1311            :: compdof_3d_mdl_cpl
1312    integer(kind=PIO_Offset), &
1313            dimension((ime - ims + 1) * (jme - jms + 1) * grid%ensdim) &
1314            :: compdof_3d_ensemble
1315    integer(kind=PIO_Offset), &
1316            dimension((jme - jms + 1) * (kme - kms + 1) * grid%spec_bdy_width ) &
1317            :: compdof_3d_xsz, compdof_3d_xez
1318    integer(kind=PIO_Offset), &
1319            dimension((ime - ims + 1) * (kme - kms + 1) * grid%spec_bdy_width ) &
1320            :: compdof_3d_ysz, compdof_3d_yez
1321    integer(kind=PIO_Offset), &
1322            dimension((jme - jms + 1) * grid%spec_bdy_width ) &
1323            :: compdof_2d_xs, compdof_2d_xe
1324    integer(kind=PIO_Offset), &
1325            dimension((ime - ims + 1) * grid%spec_bdy_width ) &
1326            :: compdof_2d_ys, compdof_2d_ye
1327    integer(kind=PIO_Offset), &
1328            dimension((ime - ims + 1) * (jme - jms + 1)) &
1329            :: compdof_2d
1330    integer :: dims3d(3), dims2d(2), dims2di(3)
1331    integer :: dims3d_xb(3), dims2d_xb(2)
1332    integer :: dims3d_yb(3), dims2d_yb(2)
1333    integer :: dims3d_land(3), dims3d_soil(3), dims3d_soil_layers(3)
1334    integer :: dims3d_mdl_cpl(3)
1335    integer :: dims3d_ensemble(3)
1336    integer :: lite, ljte, lkte
1337    integer :: i, j, k, n, npos
1338
1339    DH%first_operation = .false.
1340    communicator = grid%communicator
1341    myrank = DH%myrank
1342
1343 !--For MASS variables
1344    dims3d(1) = ide - 1
1345    dims3d(2) = jde - 1
1346    dims3d(3) = kde - 1
1347
1348    lite = ite
1349    ljte = jte
1350    lkte = kte
1351
1352    if(lite > dims3d(1)) lite = dims3d(1)
1353    if(ljte > dims3d(2)) ljte = dims3d(2)
1354    if(lkte > dims3d(3)) lkte = dims3d(3)
1355
1356    dims3d_land(1) = dims3d(1)
1357    dims3d_land(2) = dims3d(2)
1358    dims3d_land(3) = grid%num_land_cat
1359
1360    dims3d_soil(1) = dims3d(1)
1361    dims3d_soil(2) = dims3d(2)
1362    dims3d_soil(3) = grid%num_soil_cat
1363
1364    dims3d_soil_layers(1) = dims3d(1)
1365    dims3d_soil_layers(2) = dims3d(2)
1366    dims3d_soil_layers(3) = grid%num_soil_layers
1367
1368    dims3d_mdl_cpl(1) = dims3d(1)
1369    dims3d_mdl_cpl(2) = dims3d(2)
1370    dims3d_mdl_cpl(3) = grid%num_ext_model_couple_dom
1371
1372    dims3d_ensemble(1) = dims3d(1)
1373    dims3d_ensemble(2) = dims3d(2)
1374    dims3d_ensemble(3) = grid%ensdim
1375
1376    dims2d(1) = dims3d(1)
1377    dims2d(2) = dims3d(2)
1378
1379    dims2di(1) = dims3d(1)
1380    dims2di(2) = dims3d(2)
1381    dims2di(3) = 1
1382
1383    dims3d_xb(1) = dims3d(2)
1384    dims3d_xb(2) = dims3d(3)
1385    dims3d_xb(3) = grid%spec_bdy_width
1386
1387    dims3d_yb(1) = dims3d(1)
1388    dims3d_yb(2) = dims3d(3)
1389    dims3d_yb(3) = grid%spec_bdy_width
1390
1391    dims2d_xb(1) = dims2d(2)
1392    dims2d_xb(2) = grid%spec_bdy_width
1393
1394    dims2d_yb(1) = dims2d(1)
1395    dims2d_yb(2) = grid%spec_bdy_width
1396
1397   !write(unit=0, fmt='(3a,i6)') 'file: ', __FILE__, ', line: ', __LINE__
1398   !write(unit=0, fmt='(a, 6i6)') 'dims2d = ', dims2d
1399   !write(unit=0, fmt='(a, 6i6)') 'dims3d = ', dims3d
1400   !write(unit=0, fmt='(a, 6i6)') 'dims3d_land = ', dims3d_land
1401   !write(unit=0, fmt='(a, 6i6)') 'dims3d_soil = ', dims3d_soil
1402   !write(unit=0, fmt='(a, 6i6)') 'grid%num_land_cat = ', grid%num_land_cat
1403   !write(unit=0, fmt='(a, 6i6)') 'grid%num_soil_cat = ', grid%num_soil_cat
1404   !write(unit=0, fmt='(a, 6i6)') 'grid%num_soil_layers = ', grid%num_soil_layers
1405   !write(unit=0, fmt='(a, 6i6)') 'grid%num_ext_model_couple_dom = ', grid%num_ext_model_couple_dom
1406   !write(unit=0, fmt='(a, 6i6)') 'grid%spec_bdy_width = ', grid%spec_bdy_width
1407
1408    do j = jms, jme
1409       do i = ims, ime
1410          npos = (i - ims + 1) + (ime - ims + 1) * (j - jms)
1411          compdof_2d(npos) = 0
1412       enddo
1413
1414       do k = kms, kme
1415       do i = ims, ime
1416          npos = (i - ims + 1) + (ime - ims + 1) * (k - kms + (kme - kms + 1) * (j - jms))
1417          compdof_3d(npos) = 0
1418       enddo
1419       enddo
1420
1421       do k = 1, dims3d_land(3)
1422       do i = ims, ime
1423          npos = (i - ims + 1) + (ime - ims + 1) * (k - 1 + dims3d_land(3) * (j - jms))
1424          compdof_3d_land(npos) = 0
1425       enddo
1426       enddo
1427
1428       do k = 1, dims3d_soil(3)
1429       do i = ims, ime
1430          npos = (i - ims + 1) + (ime - ims + 1) * (k - 1 + dims3d_soil(3) * (j - jms))
1431          compdof_3d_soil(npos) = 0
1432       enddo
1433       enddo
1434
1435       do k = 1, dims3d_soil_layers(3)
1436       do i = ims, ime
1437          npos = (i - ims + 1) + (ime - ims + 1) * (k - 1 + dims3d_soil_layers(3) * (j - jms))
1438          compdof_3d_soil_layers(npos) = 0
1439       enddo
1440       enddo
1441
1442       do k = 1, dims3d_mdl_cpl(3)
1443       do i = ims, ime
1444          npos = (i - ims + 1) + (ime - ims + 1) * (k - 1 + dims3d_mdl_cpl(3) * (j - jms))
1445          compdof_3d_mdl_cpl(npos) = 0
1446       enddo
1447       enddo
1448    enddo
1449
1450    do n = 1, grid%spec_bdy_width
1451       do i = ims, ime
1452          npos = i - ims + 1 + (ime - ims + 1) * (n - 1)
1453          compdof_2d_ys(npos) = 0
1454          compdof_2d_ye(npos) = 0
1455       enddo
1456
1457       do j = jms, jme
1458          npos = j - jms + 1 + (jme - jms + 1) * (n - 1)
1459          compdof_2d_xs(npos) = 0
1460          compdof_2d_xe(npos) = 0
1461       enddo
1462
1463       do k = kms, kme
1464       do i = ims, ime
1465          npos = i - ims + 1 + (ime - ims + 1) * (k - kms + (kme - kms + 1) * (n - 1))
1466          compdof_3d_ysz(npos) = 0
1467          compdof_3d_yez(npos) = 0
1468       enddo
1469       enddo
1470
1471       do k = kms, kme
1472       do j = jms, jme
1473          npos = j - jms + 1 + (jme - jms + 1) * (k - kms + (kme - kms + 1) * (n - 1))
1474          compdof_3d_xsz(npos) = 0
1475          compdof_3d_xez(npos) = 0
1476       enddo
1477       enddo
1478    enddo
1479
1480    do j = jts, ljte
1481       do i = its, lite
1482          npos = (i - ims + 1) + (ime - ims + 1) * (j - jms)
1483          compdof_2d(npos) = i + dims2d(1) * (j - 1)
1484       end do
1485
1486       do k = kts, lkte
1487       do i = its, lite
1488          npos = (i - ims + 1) + (ime - ims + 1) * (k - kms + (kme - kms + 1) * (j - jms))
1489          compdof_3d(npos) = i + dims3d(1) * (j - 1 + dims3d(2) * (k - 1))
1490       enddo
1491       enddo
1492
1493       do k = 1, dims3d_land(3)
1494       do i = its, lite
1495          npos = (i - ims + 1) + (ime - ims + 1) * (k - 1 + dims3d_land(3) * (j - jms))
1496          compdof_3d_land(npos) = i + dims3d_land(1) * (j - 1 + dims3d_land(2) * (k - 1))
1497       enddo
1498       enddo
1499
1500       do k = 1, dims3d_soil(3)
1501       do i = its, lite
1502          npos = (i - ims + 1) + (ime - ims + 1) * (k - 1 + dims3d_soil(3) * (j - jms))
1503          compdof_3d_soil(npos) = i + dims3d_soil(1) * (j - 1 + dims3d_soil(2) * (k - 1))
1504       enddo
1505       enddo
1506
1507       do k = 1, dims3d_soil_layers(3)
1508       do i = its, lite
1509          npos = (i - ims + 1) + (ime - ims + 1) * (k - 1 + dims3d_soil_layers(3) * (j - jms))
1510          compdof_3d_soil_layers(npos) = i + dims3d_soil_layers(1) * (j - 1 + dims3d_soil_layers(2) * (k - 1))
1511       enddo
1512       enddo
1513
1514       do k = 1, dims3d_mdl_cpl(3)
1515       do i = its, lite
1516          npos = (i - ims + 1) + (ime - ims + 1) * (k - 1 + dims3d_mdl_cpl(3) * (j - jms))
1517          compdof_3d_mdl_cpl(npos) = i + dims3d_mdl_cpl(1) * (j - 1 + dims3d_mdl_cpl(2) * (k - 1))
1518       enddo
1519       enddo
1520    enddo
1521
1522    do k = 1, dims3d_ensemble(3)
1523       do j = jms, jme
1524       do i = ims, ime
1525          npos = (i - ims + 1) + (ime - ims + 1) * (j - jms + (jme - jms + 1) * (k - 1))
1526          compdof_3d_ensemble(npos) = 0
1527       enddo
1528       enddo
1529
1530       do j = jts, ljte
1531       do i = its, lite
1532          npos = (i - ims + 1) + (ime - ims + 1) * (j - jms + (jme - jms + 1) * (k - 1))
1533          compdof_3d_ensemble(npos) = i + dims3d_ensemble(1) * (j - 1 + dims3d_ensemble(2) * (k - 1))
1534       enddo
1535       enddo
1536    enddo
1537
1538   !write(unit=0, fmt='(3a,i6)') 'File: ', __FILE__, ', line: ', __LINE__
1539   !write(unit=0, fmt='(4x,a,i6)') 'npos = ', npos
1540   !write(unit=0, fmt='(4x,a,i16)') 'compdof_3d_ensemble(npos) = ', compdof_3d_ensemble(npos)
1541
1542    if(1 == its) then
1543       do n = 1, grid%spec_bdy_width
1544       do j = jts, ljte
1545          npos = j - jms + 1 + (jme - jms + 1) * (n - 1)
1546          compdof_2d_xs(npos) = j + dims2d_xb(1) * (n - 1)
1547       enddo
1548       enddo
1549    endif
1550
1551    if(1 == jts) then
1552       do n = 1, grid%spec_bdy_width
1553       do i = its, lite
1554          npos = i - ims + 1 + (ime - ims + 1) * (n - 1)
1555          compdof_2d_ys(npos) = i + dims2d_yb(1) * (n - 1)
1556       enddo
1557       enddo
1558    endif
1559
1560    if(dims2d(1) == lite) then
1561       do n = 1, grid%spec_bdy_width
1562       do j = jts, ljte
1563          npos = j - jms + 1 + (jme - jms + 1) * (n - 1)
1564          compdof_2d_xe(npos) = j + dims2d_xb(1) * (n - 1)
1565       enddo
1566       enddo
1567    endif
1568
1569    if(dims2d(2) == ljte) then
1570       do n = 1, grid%spec_bdy_width
1571       do i = its, lite
1572          npos = i - ims + 1 + (ime - ims + 1) * (n - 1)
1573          compdof_2d_ye(npos) = i + dims2d_yb(1) * (n - 1)
1574       enddo
1575       enddo
1576    endif
1577
1578    if(1 == its) then
1579       do n = 1, grid%spec_bdy_width
1580       do k = kts, lkte
1581       do j = jts, ljte
1582          npos = j - jms + 1 + (jme - jms + 1) * (k - kms + (kme - kms + 1) * (n - 1))
1583          compdof_3d_xsz(npos) = j + dims3d_xb(1) * (k - 1 + dims3d_xb(2) * (n - 1))
1584       enddo
1585       enddo
1586       enddo
1587    endif
1588
1589    if(1 == jts) then
1590       do n = 1, grid%spec_bdy_width
1591       do k = kts, lkte
1592       do i = its, lite
1593          npos = i - ims + 1 + (ime - ims + 1) * (k - kms + (kme - kms + 1) * (n - 1))
1594          compdof_3d_ysz(npos) = i + dims3d_yb(1) * (k - 1 + dims3d_yb(2) * (n - 1))
1595       enddo
1596       enddo
1597       enddo
1598    endif
1599
1600    if(dims2d(1) == lite) then
1601       do n = 1, grid%spec_bdy_width
1602       do k = kts, lkte
1603       do j = jts, ljte
1604          npos = j - jms + 1 + (jme - jms + 1) * (k - kms + (kme - kms + 1) * (n - 1))
1605          compdof_3d_xez(npos) = j + dims3d_xb(1) * (k - 1 + dims3d_xb(2) * (n - 1))
1606       enddo
1607       enddo
1608       enddo
1609    endif
1610
1611    if(dims2d(2) == ljte) then
1612       do n = 1, grid%spec_bdy_width
1613       do k = kts, lkte
1614       do i = its, lite
1615          npos = i - ims + 1 + (ime - ims + 1) * (k - kms + (kme - kms + 1) * (n - 1))
1616          compdof_3d_yez(npos) = i + dims3d_yb(1) * (k - 1 + dims3d_yb(2) * (n - 1))
1617       enddo
1618       enddo
1619       enddo
1620    endif
1621
1622 !--call init_decomp in order to setup the IO decomposition with PIO
1623   !call pio_setdebuglevel(1)
1624
1625    call PIO_initdecomp(DH%iosystem, PIO_int,    dims3d, compdof_3d, DH%iodesc3d_m_int)
1626    call PIO_initdecomp(DH%iosystem, PIO_real,   dims3d, compdof_3d, DH%iodesc3d_m_real)
1627    call PIO_initdecomp(DH%iosystem, PIO_double, dims3d, compdof_3d, DH%iodesc3d_m_double)
1628
1629    call PIO_initdecomp(DH%iosystem, PIO_int,    dims3d_land, compdof_3d_land, DH%iodesc3d_land_int)
1630    call PIO_initdecomp(DH%iosystem, PIO_real,   dims3d_land, compdof_3d_land, DH%iodesc3d_land_real)
1631    call PIO_initdecomp(DH%iosystem, PIO_double, dims3d_land, compdof_3d_land, DH%iodesc3d_land_double)
1632
1633    call PIO_initdecomp(DH%iosystem, PIO_int,    dims3d_soil, compdof_3d_soil, DH%iodesc3d_soil_int)
1634    call PIO_initdecomp(DH%iosystem, PIO_real,   dims3d_soil, compdof_3d_soil, DH%iodesc3d_soil_real)
1635    call PIO_initdecomp(DH%iosystem, PIO_double, dims3d_soil, compdof_3d_soil, DH%iodesc3d_soil_double)
1636
1637    call PIO_initdecomp(DH%iosystem, PIO_int,    dims3d_soil_layers, compdof_3d_soil_layers, DH%iodesc3d_soil_layers_int)
1638    call PIO_initdecomp(DH%iosystem, PIO_real,   dims3d_soil_layers, compdof_3d_soil_layers, DH%iodesc3d_soil_layers_real)
1639    call PIO_initdecomp(DH%iosystem, PIO_double, dims3d_soil_layers, compdof_3d_soil_layers, DH%iodesc3d_soil_layers_double)
1640
1641    call PIO_initdecomp(DH%iosystem, PIO_int,    dims3d_mdl_cpl, compdof_3d_mdl_cpl, DH%iodesc3d_mdl_cpl_int)
1642    call PIO_initdecomp(DH%iosystem, PIO_real,   dims3d_mdl_cpl, compdof_3d_mdl_cpl, DH%iodesc3d_mdl_cpl_real)
1643    call PIO_initdecomp(DH%iosystem, PIO_double, dims3d_mdl_cpl, compdof_3d_mdl_cpl, DH%iodesc3d_mdl_cpl_double)
1644
1645   !call pio_setdebuglevel(1)
1646    call PIO_initdecomp(DH%iosystem, PIO_int,    dims3d_ensemble, compdof_3d_ensemble, DH%iodesc3d_ensemble_int)
1647    call PIO_initdecomp(DH%iosystem, PIO_real,   dims3d_ensemble, compdof_3d_ensemble, DH%iodesc3d_ensemble_real)
1648    call PIO_initdecomp(DH%iosystem, PIO_double, dims3d_ensemble, compdof_3d_ensemble, DH%iodesc3d_ensemble_double)
1649   !call pio_setdebuglevel(0)
1650
1651 #ifndef INTSPECIAL
1652    call PIO_initdecomp(DH%iosystem, PIO_int,    dims2d, compdof_2d, DH%iodesc2d_m_int)
1653 #else
1654    call PIO_initdecomp(DH%iosystem, PIO_int,    dims2di, compdof_2d, DH%iodesc2d_m_int)
1655 #endif
1656    call PIO_initdecomp(DH%iosystem, PIO_real,   dims2d, compdof_2d, DH%iodesc2d_m_real)
1657    call PIO_initdecomp(DH%iosystem, PIO_double, dims2d, compdof_2d, DH%iodesc2d_m_double)
1658
1659    call PIO_initdecomp(DH%iosystem, PIO_int,    dims3d_xb, compdof_3d_xsz, DH%iodesc3d_xsz_m_int)
1660    call PIO_initdecomp(DH%iosystem, PIO_real,   dims3d_xb, compdof_3d_xsz, DH%iodesc3d_xsz_m_real)
1661    call PIO_initdecomp(DH%iosystem, PIO_double, dims3d_xb, compdof_3d_xsz, DH%iodesc3d_xsz_m_double)
1662
1663    call PIO_initdecomp(DH%iosystem, PIO_int,    dims3d_xb, compdof_3d_xez, DH%iodesc3d_xez_m_int)
1664    call PIO_initdecomp(DH%iosystem, PIO_real,   dims3d_xb, compdof_3d_xez, DH%iodesc3d_xez_m_real)
1665    call PIO_initdecomp(DH%iosystem, PIO_double, dims3d_xb, compdof_3d_xez, DH%iodesc3d_xez_m_double)
1666
1667    call PIO_initdecomp(DH%iosystem, PIO_int,    dims3d_yb, compdof_3d_ysz, DH%iodesc3d_ysz_m_int)
1668    call PIO_initdecomp(DH%iosystem, PIO_real,   dims3d_yb, compdof_3d_ysz, DH%iodesc3d_ysz_m_real)
1669    call PIO_initdecomp(DH%iosystem, PIO_double, dims3d_yb, compdof_3d_ysz, DH%iodesc3d_ysz_m_double)
1670
1671    call PIO_initdecomp(DH%iosystem, PIO_int,    dims3d_yb, compdof_3d_yez, DH%iodesc3d_yez_m_int)
1672    call PIO_initdecomp(DH%iosystem, PIO_real,   dims3d_yb, compdof_3d_yez, DH%iodesc3d_yez_m_real)
1673    call PIO_initdecomp(DH%iosystem, PIO_double, dims3d_yb, compdof_3d_yez, DH%iodesc3d_yez_m_double)
1674
1675    call PIO_initdecomp(DH%iosystem, PIO_int,    dims2d_xb, compdof_2d_xs, DH%iodesc2d_xs_m_int)
1676    call PIO_initdecomp(DH%iosystem, PIO_real,   dims2d_xb, compdof_2d_xs, DH%iodesc2d_xs_m_real)
1677    call PIO_initdecomp(DH%iosystem, PIO_double, dims2d_xb, compdof_2d_xs, DH%iodesc2d_xs_m_double)
1678
1679    call PIO_initdecomp(DH%iosystem, PIO_int,    dims2d_xb, compdof_2d_xe, DH%iodesc2d_xe_m_int)
1680    call PIO_initdecomp(DH%iosystem, PIO_real,   dims2d_xb, compdof_2d_xe, DH%iodesc2d_xe_m_real)
1681    call PIO_initdecomp(DH%iosystem, PIO_double, dims2d_xb, compdof_2d_xe, DH%iodesc2d_xe_m_double)
1682
1683    call PIO_initdecomp(DH%iosystem, PIO_int,    dims2d_yb, compdof_2d_ys, DH%iodesc2d_ys_m_int)
1684    call PIO_initdecomp(DH%iosystem, PIO_real,   dims2d_yb, compdof_2d_ys, DH%iodesc2d_ys_m_real)
1685    call PIO_initdecomp(DH%iosystem, PIO_double, dims2d_yb, compdof_2d_ys, DH%iodesc2d_ys_m_double)
1686
1687    call PIO_initdecomp(DH%iosystem, PIO_int,    dims2d_yb, compdof_2d_ye, DH%iodesc2d_ye_m_int)
1688    call PIO_initdecomp(DH%iosystem, PIO_real,   dims2d_yb, compdof_2d_ye, DH%iodesc2d_ye_m_real)
1689    call PIO_initdecomp(DH%iosystem, PIO_double, dims2d_yb, compdof_2d_ye, DH%iodesc2d_ye_m_double)
1690
1691 !--For X-STAG variables
1692    dims3d(1) = ide
1693    dims3d(2) = jde - 1
1694    dims3d(3) = kde - 1
1695
1696    lite = ite
1697    ljte = jte
1698    lkte = kte
1699
1700    if(lite > dims3d(1)) lite = dims3d(1)
1701    if(ljte > dims3d(2)) ljte = dims3d(2)
1702    if(lkte > dims3d(3)) lkte = dims3d(3)
1703
1704    dims2d(1) = dims3d(1)
1705    dims2d(2) = dims3d(2)
1706
1707    dims3d_xb(1) = dims3d(2)
1708    dims3d_xb(2) = dims3d(3)
1709    dims3d_xb(3) = grid%spec_bdy_width
1710
1711    dims3d_yb(1) = dims3d(1)
1712    dims3d_yb(2) = dims3d(3)
1713    dims3d_yb(3) = grid%spec_bdy_width
1714
1715   !compdof_3d =  0
1716   !compdof_2d =  0
1717
1718    do j = jms, jme
1719       do i = ims, ime
1720          npos = (i - ims + 1) + (ime - ims + 1) * (j - jms)
1721          compdof_2d(npos) = 0
1722       enddo
1723
1724       do k = kms, kme
1725       do i = ims, ime
1726          npos = (i - ims + 1) + (ime - ims + 1) * (k - kms + (kme - kms + 1) * (j - jms))
1727          compdof_3d(npos) = 0
1728       enddo
1729       enddo
1730    enddo
1731
1732    do n = 1, grid%spec_bdy_width
1733       do k = kms, kme
1734       do i = ims, ime
1735          npos = i - ims + 1 + (ime - ims + 1) * (k - kms + (kme - kms + 1) * (n - 1))
1736          compdof_3d_ysz(npos) = 0
1737          compdof_3d_yez(npos) = 0
1738       enddo
1739       enddo
1740
1741       do k = kms, kme
1742       do j = jms, jme
1743          npos = j - jms + 1 + (jme - jms + 1) * (k - kms + (kme - kms + 1) * (n - 1))
1744          compdof_3d_xsz(npos) = 0
1745          compdof_3d_xez(npos) = 0
1746       enddo
1747       enddo
1748    enddo
1749
1750    do j = jts, ljte
1751       do i = its, lite
1752          npos = (i - ims + 1) + (ime - ims + 1) * (j - jms)
1753          compdof_2d(npos) = i + dims2d(1) * (j - 1)
1754       end do
1755
1756       do k = kts, lkte
1757       do i = its, lite
1758          npos = (i - ims + 1) + (ime - ims + 1) * (k - kms + (kme - kms + 1) * (j - jms))
1759          compdof_3d(npos) = i + dims3d(1) * (j - 1 + dims3d(2) * (k - 1))
1760       enddo
1761       enddo
1762    enddo
1763
1764    if(1 == its) then
1765       do n = 1, grid%spec_bdy_width
1766       do k = kts, lkte
1767       do j = jts, ljte
1768          npos = j - jms + 1 + (jme - jms + 1) * (k - kms + (kme - kms + 1) * (n - 1))
1769          compdof_3d_xsz(npos) = j + dims3d_xb(1) * (k - 1 + dims3d_xb(2) * (n - 1))
1770       enddo
1771       enddo
1772       enddo
1773    endif
1774
1775    if(1 == jts) then
1776       do n = 1, grid%spec_bdy_width
1777       do k = kts, lkte
1778       do i = its, lite
1779          npos = i - ims + 1 + (ime - ims + 1) * (k - kms + (kme - kms + 1) * (n - 1))
1780          compdof_3d_ysz(npos) = i + dims3d_yb(1) * (k - 1 + dims3d_yb(2) * (n - 1))
1781       enddo
1782       enddo
1783       enddo
1784    endif
1785
1786    if(dims3d(1) == lite) then
1787       do n = 1, grid%spec_bdy_width
1788       do k = kts, lkte
1789       do j = jts, ljte
1790          npos = j - jms + 1 + (jme - jms + 1) * (k - kms + (kme - kms + 1) * (n - 1))
1791          compdof_3d_xez(npos) = j + dims3d_xb(1) * (k - 1 + dims3d_xb(2) * (n - 1))
1792       enddo
1793       enddo
1794       enddo
1795    endif
1796
1797    if(dims3d(2) == ljte) then
1798       do n = 1, grid%spec_bdy_width
1799       do k = kts, lkte
1800       do i = its, lite
1801          npos = i - ims + 1 + (ime - ims + 1) * (k - kms + (kme - kms + 1) * (n - 1))
1802          compdof_3d_yez(npos) = i + dims3d_yb(1) * (k - 1 + dims3d_yb(2) * (n - 1))
1803       enddo
1804       enddo
1805       enddo
1806    endif
1807
1808 !--call init_decomp in order to setup the IO decomposition with PIO
1809    call PIO_initdecomp(DH%iosystem, PIO_double, dims3d, compdof_3d, DH%iodesc3d_u_double)
1810    call PIO_initdecomp(DH%iosystem, PIO_real,   dims3d, compdof_3d, DH%iodesc3d_u_real)
1811    call PIO_initdecomp(DH%iosystem, PIO_int,    dims3d, compdof_3d, DH%iodesc3d_u_int)
1812
1813    call PIO_initdecomp(DH%iosystem, PIO_double, dims2d, compdof_2d, DH%iodesc2d_u_double)
1814    call PIO_initdecomp(DH%iosystem, PIO_real,   dims2d, compdof_2d, DH%iodesc2d_u_real)
1815    call PIO_initdecomp(DH%iosystem, PIO_int,    dims2d, compdof_2d, DH%iodesc2d_u_int)
1816
1817    call PIO_initdecomp(DH%iosystem, PIO_int,    dims3d_xb, compdof_3d_xsz, DH%iodesc3d_xsz_u_int)
1818    call PIO_initdecomp(DH%iosystem, PIO_real,   dims3d_xb, compdof_3d_xsz, DH%iodesc3d_xsz_u_real)
1819    call PIO_initdecomp(DH%iosystem, PIO_double, dims3d_xb, compdof_3d_xsz, DH%iodesc3d_xsz_u_double)
1820
1821    call PIO_initdecomp(DH%iosystem, PIO_int,    dims3d_xb, compdof_3d_xez, DH%iodesc3d_xez_u_int)
1822    call PIO_initdecomp(DH%iosystem, PIO_real,   dims3d_xb, compdof_3d_xez, DH%iodesc3d_xez_u_real)
1823    call PIO_initdecomp(DH%iosystem, PIO_double, dims3d_xb, compdof_3d_xez, DH%iodesc3d_xez_u_double)
1824
1825    call PIO_initdecomp(DH%iosystem, PIO_int,    dims3d_yb, compdof_3d_ysz, DH%iodesc3d_ysz_u_int)
1826    call PIO_initdecomp(DH%iosystem, PIO_real,   dims3d_yb, compdof_3d_ysz, DH%iodesc3d_ysz_u_real)
1827    call PIO_initdecomp(DH%iosystem, PIO_double, dims3d_yb, compdof_3d_ysz, DH%iodesc3d_ysz_u_double)
1828
1829    call PIO_initdecomp(DH%iosystem, PIO_int,    dims3d_yb, compdof_3d_yez, DH%iodesc3d_yez_u_int)
1830    call PIO_initdecomp(DH%iosystem, PIO_real,   dims3d_yb, compdof_3d_yez, DH%iodesc3d_yez_u_real)
1831    call PIO_initdecomp(DH%iosystem, PIO_double, dims3d_yb, compdof_3d_yez, DH%iodesc3d_yez_u_double)
1832
1833 !--For Y-STAG variables
1834    dims3d(1) = ide - 1
1835    dims3d(2) = jde
1836    dims3d(3) = kde - 1
1837
1838    lite = ite
1839    ljte = jte
1840    lkte = kte
1841
1842    if(lite > dims3d(1)) lite = dims3d(1)
1843    if(ljte > dims3d(2)) ljte = dims3d(2)
1844    if(lkte > dims3d(3)) lkte = dims3d(3)
1845
1846    dims2d(1) = dims3d(1)
1847    dims2d(2) = dims3d(2)
1848
1849    dims3d_xb(1) = dims3d(2)
1850    dims3d_xb(2) = dims3d(3)
1851    dims3d_xb(3) = grid%spec_bdy_width
1852
1853    dims3d_yb(1) = dims3d(1)
1854    dims3d_yb(2) = dims3d(3)
1855    dims3d_yb(3) = grid%spec_bdy_width
1856
1857   !compdof_3d =  0
1858   !compdof_2d =  0
1859
1860    do j = jms, jme
1861       do i = ims, ime
1862          npos = (i - ims + 1) + (ime - ims + 1) * (j - jms)
1863          compdof_2d(npos) = 0
1864       enddo
1865
1866       do k = kms, kme
1867       do i = ims, ime
1868          npos = (i - ims + 1) + (ime - ims + 1) * (k - kms + (kme - kms + 1) * (j - jms))
1869          compdof_3d(npos) = 0
1870       enddo
1871       enddo
1872    enddo
1873
1874    do n = 1, grid%spec_bdy_width
1875       do k = kms, kme
1876       do i = ims, ime
1877          npos = i - ims + 1 + (ime - ims + 1) * (k - kms + (kme - kms + 1) * (n - 1))
1878          compdof_3d_ysz(npos) = 0
1879          compdof_3d_yez(npos) = 0
1880       enddo
1881       enddo
1882
1883       do k = kms, kme
1884       do j = jms, jme
1885          npos = j - jms + 1 + (jme - jms + 1) * (k - kms + (kme - kms + 1) * (n - 1))
1886          compdof_3d_xsz(npos) = 0
1887          compdof_3d_xez(npos) = 0
1888       enddo
1889       enddo
1890    enddo
1891
1892    do j = jts, ljte
1893       do i = its, lite
1894          npos = (i - ims + 1) + (ime - ims + 1) * (j - jms)
1895          compdof_2d(npos) = i + dims2d(1) * (j - 1)
1896       end do
1897
1898       do k = kts, lkte
1899       do i = its, lite
1900          npos = (i - ims + 1) + (ime - ims + 1) * (k - kms + (kme - kms + 1) * (j - jms))
1901          compdof_3d(npos) = i + dims3d(1) * (j - 1 + dims3d(2) * (k - 1))
1902       enddo
1903       enddo
1904    enddo
1905
1906    if(1 == its) then
1907       do n = 1, grid%spec_bdy_width
1908       do k = kts, lkte
1909       do j = jts, ljte
1910          npos = j - jms + 1 + (jme - jms + 1) * (k - kms + (kme - kms + 1) * (n - 1))
1911          compdof_3d_xsz(npos) = j + dims3d_xb(1) * (k - 1 + dims3d_xb(2) * (n - 1))
1912       enddo
1913       enddo
1914       enddo
1915    endif
1916
1917    if(1 == jts) then
1918       do n = 1, grid%spec_bdy_width
1919       do k = kts, lkte
1920       do i = its, lite
1921          npos = i - ims + 1 + (ime - ims + 1) * (k - kms + (kme - kms + 1) * (n - 1))
1922          compdof_3d_ysz(npos) = i + dims3d_yb(1) * (k - 1 + dims3d_yb(2) * (n - 1))
1923       enddo
1924       enddo
1925       enddo
1926    endif
1927
1928    if(dims3d(1) == lite) then
1929       do n = 1, grid%spec_bdy_width
1930       do k = kts, lkte
1931       do j = jts, ljte
1932          npos = j - jms + 1 + (jme - jms + 1) * (k - kms + (kme - kms + 1) * (n - 1))
1933          compdof_3d_xez(npos) = j + dims3d_xb(1) * (k - 1 + dims3d_xb(2) * (n - 1))
1934       enddo
1935       enddo
1936       enddo
1937    endif
1938
1939    if(dims3d(2) == ljte) then
1940       do n = 1, grid%spec_bdy_width
1941       do k = kts, lkte
1942       do i = its, lite
1943          npos = i - ims + 1 + (ime - ims + 1) * (k - kms + (kme - kms + 1) * (n - 1))
1944          compdof_3d_yez(npos) = i + dims3d_yb(1) * (k - 1 + dims3d_yb(2) * (n - 1))
1945       enddo
1946       enddo
1947       enddo
1948    endif
1949
1950 !--call init_decomp in order to setup the IO decomposition with PIO
1951    call PIO_initdecomp(DH%iosystem, PIO_double, dims3d, compdof_3d, DH%iodesc3d_v_double)
1952    call PIO_initdecomp(DH%iosystem, PIO_real,   dims3d, compdof_3d, DH%iodesc3d_v_real)
1953    call PIO_initdecomp(DH%iosystem, PIO_int,    dims3d, compdof_3d, DH%iodesc3d_v_int)
1954
1955    call PIO_initdecomp(DH%iosystem, PIO_double, dims2d, compdof_2d, DH%iodesc2d_v_double)
1956    call PIO_initdecomp(DH%iosystem, PIO_real,   dims2d, compdof_2d, DH%iodesc2d_v_real)
1957    call PIO_initdecomp(DH%iosystem, PIO_int,    dims2d, compdof_2d, DH%iodesc2d_v_int)
1958
1959    call PIO_initdecomp(DH%iosystem, PIO_int,    dims3d_xb, compdof_3d_xsz, DH%iodesc3d_xsz_v_int)
1960    call PIO_initdecomp(DH%iosystem, PIO_real,   dims3d_xb, compdof_3d_xsz, DH%iodesc3d_xsz_v_real)
1961    call PIO_initdecomp(DH%iosystem, PIO_double, dims3d_xb, compdof_3d_xsz, DH%iodesc3d_xsz_v_double)
1962
1963    call PIO_initdecomp(DH%iosystem, PIO_int,    dims3d_xb, compdof_3d_xez, DH%iodesc3d_xez_v_int)
1964    call PIO_initdecomp(DH%iosystem, PIO_real,   dims3d_xb, compdof_3d_xez, DH%iodesc3d_xez_v_real)
1965    call PIO_initdecomp(DH%iosystem, PIO_double, dims3d_xb, compdof_3d_xez, DH%iodesc3d_xez_v_double)
1966
1967    call PIO_initdecomp(DH%iosystem, PIO_int,    dims3d_yb, compdof_3d_ysz, DH%iodesc3d_ysz_v_int)
1968    call PIO_initdecomp(DH%iosystem, PIO_real,   dims3d_yb, compdof_3d_ysz, DH%iodesc3d_ysz_v_real)
1969    call PIO_initdecomp(DH%iosystem, PIO_double, dims3d_yb, compdof_3d_ysz, DH%iodesc3d_ysz_v_double)
1970
1971    call PIO_initdecomp(DH%iosystem, PIO_int,    dims3d_yb, compdof_3d_yez, DH%iodesc3d_yez_v_int)
1972    call PIO_initdecomp(DH%iosystem, PIO_real,   dims3d_yb, compdof_3d_yez, DH%iodesc3d_yez_v_real)
1973    call PIO_initdecomp(DH%iosystem, PIO_double, dims3d_yb, compdof_3d_yez, DH%iodesc3d_yez_v_double)
1974
1975 !--For Z-STAG variables
1976    dims3d(1) = ide - 1
1977    dims3d(2) = jde - 1
1978    dims3d(3) = kde
1979
1980    dims2d(1) = dims3d(1)
1981    dims2d(2) = dims3d(2)
1982
1983    dims3d_xb(1) = dims3d(2)
1984    dims3d_xb(2) = dims3d(3)
1985    dims3d_xb(3) = grid%spec_bdy_width
1986
1987    dims3d_yb(1) = dims3d(1)
1988    dims3d_yb(2) = dims3d(3)
1989    dims3d_yb(3) = grid%spec_bdy_width
1990
1991    lite = ite
1992    ljte = jte
1993    lkte = kte
1994
1995    if(lite > dims3d(1)) lite = dims3d(1)
1996    if(ljte > dims3d(2)) ljte = dims3d(2)
1997    if(lkte > dims3d(3)) lkte = dims3d(3)
1998
1999   !compdof_3d =  0
2000
2001    do j = jms, jme
2002       do k = kms, kme
2003       do i = ims, ime
2004          npos = (i - ims + 1) + (ime - ims + 1) * (k - kms + (kme - kms + 1) * (j - jms))
2005          compdof_3d(npos) = 0
2006       enddo
2007       enddo
2008    enddo
2009
2010    do n = 1, grid%spec_bdy_width
2011       do k = kms, kme
2012       do i = ims, ime
2013          npos = i - ims + 1 + (ime - ims + 1) * (k - kms + (kme - kms + 1) * (n - 1))
2014          compdof_3d_ysz(npos) = 0
2015          compdof_3d_yez(npos) = 0
2016       enddo
2017       enddo
2018
2019       do k = kms, kme
2020       do j = jms, jme
2021          npos = j - jms + 1 + (jme - jms + 1) * (k - kms + (kme - kms + 1) * (n - 1))
2022          compdof_3d_xsz(npos) = 0
2023          compdof_3d_xez(npos) = 0
2024       enddo
2025       enddo
2026    enddo
2027
2028    do j = jts, ljte
2029    do k = kts, lkte
2030    do i = its, lite
2031       npos = (i - ims + 1) + (ime - ims + 1) * (k - kms + (kme - kms + 1) * (j - jms))
2032       compdof_3d(npos) = i + dims3d(1) * (j - 1 + dims3d(2) * (k - 1))
2033    enddo
2034    enddo
2035    enddo
2036
2037    if(1 == its) then
2038       do n = 1, grid%spec_bdy_width
2039       do k = kts, lkte
2040       do j = jts, ljte
2041          npos = j - jms + 1 + (jme - jms + 1) * (k - kms + (kme - kms + 1) * (n - 1))
2042          compdof_3d_xsz(npos) = j + dims3d_xb(1) * (k - 1 + dims3d_xb(2) * (n - 1))
2043       enddo
2044       enddo
2045       enddo
2046    endif
2047
2048    if(1 == jts) then
2049       do n = 1, grid%spec_bdy_width
2050       do k = kts, lkte
2051       do i = its, lite
2052          npos = i - ims + 1 + (ime - ims + 1) * (k - kms + (kme - kms + 1) * (n - 1))
2053          compdof_3d_ysz(npos) = i + dims3d_yb(1) * (k - 1 + dims3d_yb(2) * (n - 1))
2054       enddo
2055       enddo
2056       enddo
2057    endif
2058
2059    if(dims3d(1) == lite) then
2060       do n = 1, grid%spec_bdy_width
2061       do k = kts, lkte
2062       do j = jts, ljte
2063          npos = j - jms + 1 + (jme - jms + 1) * (k - kms + (kme - kms + 1) * (n - 1))
2064          compdof_3d_xez(npos) = j + dims3d_xb(1) * (k - 1 + dims3d_xb(2) * (n - 1))
2065       enddo
2066       enddo
2067       enddo
2068    endif
2069
2070    if(dims3d(2) == ljte) then
2071       do n = 1, grid%spec_bdy_width
2072       do k = kts, lkte
2073       do i = its, lite
2074          npos = i - ims + 1 + (ime - ims + 1) * (k - kms + (kme - kms + 1) * (n - 1))
2075          compdof_3d_yez(npos) = i + dims3d_yb(1) * (k - 1 + dims3d_yb(2) * (n - 1))
2076       enddo
2077       enddo
2078       enddo
2079    endif
2080
2081 !--call init_decomp in order to setup the IO decomposition with PIO
2082    call PIO_initdecomp(DH%iosystem, PIO_double, dims3d, compdof_3d, DH%iodesc3d_w_double)
2083    call PIO_initdecomp(DH%iosystem, PIO_real,   dims3d, compdof_3d, DH%iodesc3d_w_real)
2084    call PIO_initdecomp(DH%iosystem, PIO_int,    dims3d, compdof_3d, DH%iodesc3d_w_int)
2085
2086    call PIO_initdecomp(DH%iosystem, PIO_int,    dims3d_xb, compdof_3d_xsz, DH%iodesc3d_xsz_w_int)
2087    call PIO_initdecomp(DH%iosystem, PIO_real,   dims3d_xb, compdof_3d_xsz, DH%iodesc3d_xsz_w_real)
2088    call PIO_initdecomp(DH%iosystem, PIO_double, dims3d_xb, compdof_3d_xsz, DH%iodesc3d_xsz_w_double)
2089
2090    call PIO_initdecomp(DH%iosystem, PIO_int,    dims3d_xb, compdof_3d_xez, DH%iodesc3d_xez_w_int)
2091    call PIO_initdecomp(DH%iosystem, PIO_real,   dims3d_xb, compdof_3d_xez, DH%iodesc3d_xez_w_real)
2092    call PIO_initdecomp(DH%iosystem, PIO_double, dims3d_xb, compdof_3d_xez, DH%iodesc3d_xez_w_double)
2093
2094    call PIO_initdecomp(DH%iosystem, PIO_int,    dims3d_yb, compdof_3d_ysz, DH%iodesc3d_ysz_w_int)
2095    call PIO_initdecomp(DH%iosystem, PIO_real,   dims3d_yb, compdof_3d_ysz, DH%iodesc3d_ysz_w_real)
2096    call PIO_initdecomp(DH%iosystem, PIO_double, dims3d_yb, compdof_3d_ysz, DH%iodesc3d_ysz_w_double)
2097
2098    call PIO_initdecomp(DH%iosystem, PIO_int,    dims3d_yb, compdof_3d_yez, DH%iodesc3d_yez_w_int)
2099    call PIO_initdecomp(DH%iosystem, PIO_real,   dims3d_yb, compdof_3d_yez, DH%iodesc3d_yez_w_real)
2100    call PIO_initdecomp(DH%iosystem, PIO_double, dims3d_yb, compdof_3d_yez, DH%iodesc3d_yez_w_double)
2101
2102 end subroutine define_pio_iodesc
2103
2104 subroutine reorder (MemoryOrder,MemO)
2105   implicit none
2106   include 'wrf_status_codes.h'
2107   character*(*)     ,intent(in)    :: MemoryOrder
2108   character*3       ,intent(out)   :: MemO
2109   character*3                      :: MemOrd
2110   integer                          :: N,i,i1,i2,i3
2111
2112   MemO = MemoryOrder
2113   N = len_trim(MemoryOrder)
2114   if(N == 1) return
2115   call lowercase(MemoryOrder,MemOrd)
2116 ! never invert the boundary codes
2117   select case ( MemOrd )
2118      case ( 'xsz','xez','ysz','yez' )
2119        return
2120      case default
2121        continue
2122   end select
2123   i1 = 1
2124   i3 = 1
2125   do i=2,N
2126     if(ichar(MemOrd(i:i)) < ichar(MemOrd(i1:i1))) I1 = i
2127     if(ichar(MemOrd(i:i)) > ichar(MemOrd(i3:i3))) I3 = i
2128   enddo
2129   if(N == 2) then
2130     i2=i3
2131   else
2132     i2 = 6-i1-i3
2133   endif
2134   MemO(1:1) = MemoryOrder(i1:i1)
2135   MemO(2:2) = MemoryOrder(i2:i2)
2136   if(N == 3) MemO(3:3) = MemoryOrder(i3:i3)
2137   if(MemOrd(i1:i1) == 's' .or. MemOrd(i1:i1) == 'e') then
2138     MemO(1:N-1) = MemO(2:N)
2139     MemO(N:N  ) = MemoryOrder(i1:i1)
2140   endif
2141   return
2142 end subroutine reorder
2143
2144 end module pio_routines
2145