utilities_mod sourceΒΆ

   1module utilities_mod
   2  ! MODULE utilities_mod (prefix='utl' category='8. Low-level utilities and constants')
   3  !
   4  !:Purpose: A place to collect numerous simple utility routines.
   5  !
   6  use clibInterfaces_mod
   7  use randomNumber_mod
   8
   9  implicit none
  10  save
  11  private
  12
  13  ! public procedures
  14  public :: utl_fstlir,  utl_fstlir_r4, utl_fstecr
  15  public :: utl_matSqrt, utl_matInverse, utl_eigenDecomp
  16  public :: utl_pseudo_inverse
  17  public :: utl_writeStatus, utl_getfldprm, utl_abort, utl_checkAllocationStatus
  18  public :: utl_stopAndWait4Debug
  19  public :: utl_open_asciifile, utl_stnid_equal, utl_resize, utl_str
  20  public :: utl_get_stringId, utl_get_Id, utl_isNamelistPresent
  21  public :: utl_readFstField
  22  public :: utl_varNamePresentInFile
  23  public :: utl_reAllocate
  24  public :: utl_heapsort2d
  25  public :: utl_combineString, utl_splitString, utl_removeEmptyStrings
  26  public :: utl_stringArrayToIntegerArray, utl_parseColumns
  27  public :: utl_copyFile, utl_allReduce, utl_findloc, utl_findlocs
  28  public :: utl_randomOrderInt
  29  public :: utl_tmg_start, utl_tmg_stop, utl_medianIndex
  30
  31  ! module interfaces
  32  ! -----------------
  33
  34  ! interface for resizing arrays
  35  interface utl_resize
  36    module procedure utl_resize_1d_real
  37    module procedure utl_resize_1d_int
  38    module procedure utl_resize_1d_str
  39    module procedure utl_resize_2d_real
  40    module procedure utl_resize_3d_real
  41  end interface utl_resize
  42
  43  ! interface for conversion to a left-justified string (useful for calls to utl_abort)
  44  interface utl_str
  45    module procedure utl_int2str
  46    module procedure utl_float2str
  47  end interface utl_str
  48
  49  interface utl_reAllocate
  50    module procedure utl_reAllocate_char_1d
  51    module procedure utl_reAllocate_char_2d
  52    module procedure utl_reAllocate_char_3d
  53    module procedure utl_reAllocate_log_1d
  54    module procedure utl_reAllocate_log_2d
  55    module procedure utl_reAllocate_log_3d
  56    module procedure utl_reAllocate_int_1d
  57    module procedure utl_reAllocate_int_2d
  58    module procedure utl_reAllocate_int_3d
  59    module procedure utl_reAllocate_r4_1d
  60    module procedure utl_reAllocate_r8_1d
  61    module procedure utl_reAllocate_r4_2d
  62    module procedure utl_reAllocate_r8_2d
  63    module procedure utl_reAllocate_r4_3d
  64    module procedure utl_reAllocate_r8_3d
  65    module procedure utl_reAllocate_r4_4d
  66    module procedure utl_reAllocate_r8_4d
  67    module procedure utl_reAllocate_r4_5d
  68    module procedure utl_reAllocate_r8_5d
  69  end interface utl_reAllocate
  70
  71  interface utl_findloc
  72    module procedure utl_findloc_char
  73    module procedure utl_findloc_int
  74  end interface utl_findloc
  75
  76  interface utl_findlocs
  77    module procedure utl_findlocs_char
  78  end interface utl_findlocs
  79
  80contains
  81
  82  function utl_fstlir(fld8, iun, ni, nj, nk, datev, etiket, &
  83                      ip1, ip2, ip3, typvar, nomvar) result(vfstlir)
  84    implicit none
  85
  86    ! Arguments:
  87    real(8),          intent(inout) :: fld8(*)
  88    integer,          intent(in)    :: iun
  89    integer,          intent(in)    :: ni
  90    integer,          intent(in)    :: nj
  91    integer,          intent(in)    :: nk
  92    integer,          intent(in)    :: datev
  93    integer,          intent(in)    :: ip1
  94    integer,          intent(in)    :: ip2
  95    integer,          intent(in)    :: ip3
  96    character(len=*), intent(in)    :: etiket
  97    character(len=*), intent(in)    :: nomvar
  98    character(len=*), intent(in)    :: typvar
  99    ! Result:
 100    integer :: vfstlir
 101
 102    ! Locals:
 103    integer :: key1,key2, ilen, jk1, jk2, jk3, la
 104    real(4), allocatable :: buffer4(:)
 105    integer :: fstluk, fstinf
 106
 107    !     Get field dimensions and allow memory for REAL copy of fld8.
 108    key1 = fstinf(iun, ni, nj, nk, datev, etiket, &
 109         ip1, ip2, ip3, typvar, nomvar)
 110
 111    if(key1 >= 0) then
 112       ilen = ni*nj*nk
 113       allocate(buffer4(ilen))
 114       !     Read field
 115       key2 = fstluk(buffer4, key1, ni, nj, nk)
 116       if(key2 >= 0) then
 117          do jk3 = 1,nk
 118             do jk2 = 1,nj
 119                do jk1 = 1,ni
 120                   la=jk1+(jk2-1)*ni+(jk3-1)*ni*nj
 121                   fld8(la) = buffer4(la)
 122                end do
 123             end do
 124          end do
 125       end if
 126
 127       deallocate(buffer4)
 128    end if
 129
 130    vfstlir=key1
 131
 132  end function utl_fstlir
 133
 134  function utl_fstlir_r4(fld_r4, iun, ni, nj, nk, datev, etiket, &
 135                         ip1, ip2, ip3, typvar, nomvar) result(vfstlir)
 136    implicit none
 137
 138    ! Arguments:
 139    real(4),          intent(inout) :: fld_r4(*)
 140    integer,          intent(in)    :: iun
 141    integer,          intent(in)    :: ni
 142    integer,          intent(in)    :: nj
 143    integer,          intent(in)    :: nk
 144    integer,          intent(in)    :: datev
 145    integer,          intent(in)    :: ip1
 146    integer,          intent(in)    :: ip2
 147    integer,          intent(in)    :: ip3
 148    character(len=*), intent(in)    :: etiket
 149    character(len=*), intent(in)    :: nomvar
 150    character(len=*), intent(in)    :: typvar
 151    ! Result:
 152    integer :: vfstlir
 153
 154    ! Locals:
 155    integer :: key1,key2, ilen, jk1, jk2, jk3, la
 156    real(4), allocatable :: buffer_r4(:)
 157    integer :: fstluk, fstinf
 158
 159    !     Get field dimensions.
 160    key1 = fstinf(iun, ni, nj, nk, datev, etiket, &
 161         ip1, ip2, ip3, typvar, nomvar)
 162
 163    if(key1 >= 0) then
 164       ilen = ni*nj*nk
 165       allocate(buffer_r4(ilen))
 166       !     Read field
 167       key2 = fstluk(buffer_r4, key1, ni, nj, nk)
 168       if(key2 >= 0) then
 169          do jk3 = 1,nk
 170             do jk2 = 1,nj
 171                do jk1 = 1,ni
 172                   la=jk1+(jk2-1)*ni+(jk3-1)*ni*nj
 173                   fld_r4(la) = buffer_r4(la)
 174                end do
 175             end do
 176          end do
 177       end if
 178
 179       deallocate(buffer_r4)
 180    end if
 181
 182    vfstlir=key1
 183
 184  end function utl_fstlir_r4
 185
 186  function utl_fstecr(fld8, npak, iun, dateo, deet, &
 187                      npas, ni, nj, nk, ip1, ip2, ip3, typvar, &
 188                      nomvar, etiket, grtyp, ig1, ig2, ig3, ig4, & 
 189                      datyp, rewrit) result(vfstecr)
 190    implicit none
 191
 192    ! Arguments:
 193    integer,          intent(in) :: ni
 194    integer,          intent(in) :: nj
 195    integer,          intent(in) :: nk
 196    real(8),          intent(in) :: fld8(ni,nj,nk)
 197    integer,          intent(in) :: iun
 198    integer,          intent(in) :: ip1
 199    integer,          intent(in) :: ip2
 200    integer,          intent(in) :: ip3
 201    integer,          intent(in) :: ig1
 202    integer,          intent(in) :: ig2
 203    integer,          intent(in) :: ig3
 204    integer,          intent(in) :: ig4
 205    integer,          intent(in) :: npak
 206    integer,          intent(in) :: dateo
 207    integer,          intent(in) :: deet
 208    integer,          intent(in) :: npas
 209    integer,          intent(in) :: datyp
 210    logical,          intent(in) :: rewrit  
 211    character(len=*), intent(in) :: etiket 
 212    character(len=*), intent(in) :: typvar
 213    character(len=*), intent(in) :: grtyp 
 214    character(len=*), intent(in) :: nomvar            
 215    ! Result:
 216    integer :: vfstecr
 217
 218    ! Locals:
 219    real(4) :: work
 220    integer :: ikey, jk1, jk2, jk3
 221    real(4), allocatable :: buffer4(:,:,:)
 222    integer :: fstecr
 223
 224    allocate(buffer4(ni,nj,nk))
 225
 226    do jk3 = 1,nk
 227      do jk2 = 1,nj
 228        do jk1 = 1,ni
 229          buffer4(jk1,jk2,jk3) = fld8(jk1,jk2,jk3)
 230        end do
 231      end do
 232    end do
 233
 234    ikey = fstecr(buffer4, work, npak, iun, dateo, deet, &
 235         npas, ni, nj, nk, ip1, ip2, ip3, typvar, nomvar, & 
 236         etiket, grtyp, ig1, ig2, ig3, ig4, datyp, rewrit)
 237
 238    deallocate(buffer4)
 239
 240    vfstecr=ikey
 241
 242  end function utl_fstecr
 243
 244  subroutine utl_matsqrt(matrix, rank, exponentSign, printInformation_opt )
 245    ! 
 246    !:Purpose: Calculate square root of an error covariance matrix
 247    !
 248    implicit none
 249
 250    ! Arguments:
 251    integer,           intent(in)    :: rank
 252    real(8),           intent(inout) :: matrix(rank,rank)
 253    real(8),           intent(in)    :: exponentSign
 254    logical, optional, intent(in)    :: printInformation_opt ! switch to print be more verbose
 255
 256    ! Locals:
 257    real(8), allocatable :: eigenValues(:)
 258    real(8), allocatable :: work(:)
 259    real(8), allocatable :: eigenVectors(:,:)
 260    integer :: sizework, info, index, index1, index2 
 261    logical :: printInformation
 262
 263    if (present(printInformation_opt)) then
 264       printInformation = printInformation_opt
 265    else
 266       printInformation = .false.
 267    end if
 268
 269    if (printInformation) then
 270      write(*,*)
 271      write(*,*) 'utl_matsqrt: Starting...'
 272    end if
 273
 274    sizework = 64 * rank
 275    allocate(work(sizework))
 276
 277    allocate(eigenValues (rank))
 278    allocate(eigenVectors(rank,rank))
 279
 280    !- Calculate EigenVectors (V) and EigenValues (D) of B matrix
 281    eigenVectors(:,:) = matrix(:,:)
 282
 283    call dsyev('V','U',rank,   & ! IN
 284               eigenVectors,   & ! INOUT
 285               rank,           & ! IN
 286               eigenValues,    & ! OUT
 287               work, sizework, & ! IN
 288               info )            ! OUT
 289
 290    if ( info /= 0 ) then
 291      write(*,*)
 292      write(*,*) 'dsyev: ',info
 293      call utl_abort('utl_matsqrt: DSYEV failed!')
 294    end if
 295
 296    if (printInformation) then
 297      write(*,*)
 298      write(*,'(1x,"Original EIGEN VALUES: ")')
 299      write(*,'(1x,10f7.3)') (eigenValues(index),index=1,rank)
 300      if (exponentSign < 0.d0) then
 301        write(*,*)
 302        write(*,'(A,1x,e14.6)') "Condition number:", &
 303             maxval(eigenValues(:))/minval(eigenValues(:))
 304      end if
 305    end if
 306
 307    !- Calculate Matrix^0.5 = V D^0.5 V^t
 308    where(eigenValues(:) < 0.d0)
 309      eigenValues = 0.d0
 310    end where
 311
 312    eigenValues(:) = eigenValues(:)**(0.5d0*exponentSign)
 313
 314    do index1 = 1, rank
 315      do index2 = 1, rank
 316        matrix(index1,index2) = sum ( eigenVectors (index1,1:rank)   &
 317                                    * eigenVectors (index2,1:rank)   &
 318                                    * eigenValues(1:rank) )
 319      end do
 320    end do
 321
 322    deallocate(eigenVectors)
 323    deallocate(eigenValues)
 324    deallocate(work)
 325    
 326    if (printInformation) then
 327      write(*,*)
 328      write(*,*) 'utl_matsqrt: Ending...'
 329    end if
 330
 331  end subroutine utl_matsqrt
 332
 333  !--------------------------------------------------------------------------
 334  ! utl_matInverse
 335  !--------------------------------------------------------------------------
 336  subroutine utl_matInverse(matrix, rank, inverseSqrt_opt, printInformation_opt)
 337    !
 338    !:Purpose: Calculate the inverse of a covariance matrix 
 339    !          and, optionally, also the inverse square-root.
 340    !
 341    implicit none
 342
 343    ! Arguments:
 344    integer,           intent(in)    :: rank                 ! order of the matrix
 345    real(8),           intent(inout) :: matrix(:,:)          ! on entry, the original matrix; on exit, the inverse
 346    real(8), optional, intent(inout) :: inverseSqrt_opt(:,:) ! if present, the inverse sqrt matrix on exit 
 347    logical, optional, intent(in)    :: printInformation_opt ! switch to print be more verbose
 348
 349    ! Locals:
 350    integer :: index1, index2, info, sizework
 351    real(8) :: sizework_r8
 352    real(8), allocatable :: work(:), eigenVectors(:,:), eigenValues(:)
 353    logical :: printInformation
 354
 355    call utl_tmg_start(180,'low-level--utl_matInverse')
 356
 357    if (present(printInformation_opt)) then
 358      printInformation = printInformation_opt
 359    else
 360      printInformation = .false.
 361    end if
 362
 363    if (printInformation) then
 364      write(*,*)' utl_matInvers: Inverse matrix of a symmetric matrix'
 365    end if
 366
 367    !     1. Computation of eigenvalues and eigenvectors
 368
 369    allocate(eigenVectors(rank,rank))
 370    allocate(eigenValues(rank))
 371
 372    do index2=1,rank
 373      do index1=1,rank
 374        eigenVectors(index1,index2)=matrix(index1,index2)
 375      end do
 376    end do
 377
 378    ! query the size of the 'work' vector by calling 'DSYEV' with 'sizework=-1'
 379    sizework = -1
 380    info = -1
 381    call dsyev('V','U',rank, eigenVectors, rank, eigenValues, sizework_r8, sizework, info)
 382
 383    ! compute the eigenvalues
 384    sizework=int(sizework_r8)
 385    allocate(work(sizework))
 386    call dsyev('V','U',rank, eigenVectors,rank, eigenValues,work, sizework, info)
 387    deallocate(work)
 388
 389    if (printInformation) then
 390      write(*,'(1x,"Original eigen values: ")')
 391      write(*,'(1x,10f7.3)') (eigenValues(index1),index1=1,rank)
 392
 393      if(minval(eigenValues) > 1.0d-10) then
 394        write(*,'(A,1x,e14.6)') "Condition number:", &
 395             maxval(eigenValues)/minval(eigenValues)
 396      end if
 397    end if
 398
 399    !     2.  Take inverse of eigenvalues
 400
 401    do index1=1,rank
 402      if(eigenValues(index1) > 1.0d-10) then
 403        eigenValues(index1)= 1.0d0/eigenValues(index1)
 404      else
 405        write(*,*) 'utl_matInverse: WARNING eigenvalue is too small = ', index1, eigenValues(index1)
 406        eigenValues(index1) = 0.0d0
 407      end if
 408    end do
 409
 410    if (printInformation) then
 411      write(*,'(1x,"Inverse of original eigen values: ")')
 412      write(*,'(1x,10f7.3)') (eigenValues(index1),index1=1,rank)
 413    end if
 414
 415    !     3.  Compute the inverse matrix
 416
 417    do index2 = 1, rank
 418      do index1 = 1, rank
 419        matrix(index1,index2) = sum ( eigenVectors (index1,1:rank)   &
 420                                    * eigenVectors (index2,1:rank)   &
 421                                    * eigenValues(1:rank) )
 422      end do
 423    end do
 424
 425    !     4.  If requested, computed the inverse square-root also
 426
 427    if (present(inverseSqrt_opt)) then
 428      do index1=1,rank
 429        if(eigenValues(index1) > 1.0d-10) then
 430          eigenValues(index1)= sqrt(eigenValues(index1))
 431        else
 432          eigenValues(index1) = 0.0d0
 433        end if
 434      end do
 435      do index2 = 1, rank
 436        do index1 = 1, rank
 437          inverseSqrt_opt(index1,index2) = sum ( eigenVectors (index1,1:rank)   &
 438                                               * eigenVectors (index2,1:rank)   &
 439                                               * eigenValues(1:rank) )
 440        end do
 441      end do
 442    end if
 443
 444    !     5. Deallocate local arrays
 445    deallocate(eigenVectors,eigenValues)
 446
 447    if (printInformation) then
 448      write(*,*) 'utl_matInverse: done'
 449      write(*,*) ' '
 450    end if
 451
 452    call utl_tmg_stop(180)
 453
 454  end subroutine utl_matInverse
 455
 456  !--------------------------------------------------------------------------
 457  ! utl_eigenDecomp
 458  !--------------------------------------------------------------------------
 459  subroutine utl_eigenDecomp(matrix, eigenValues, eigenVectors, tolerance, numReturned, printInformation_opt)
 460    !
 461    !:Purpose: Calculate eigenValues/Vectors and return only those with eigenValues
 462    !          whose magnitude is greater than the specified tolerance.
 463    !
 464    implicit none
 465
 466    ! Arguments:
 467    real(8),           intent(inout) :: matrix(:,:)          ! on entry, the original matrix; on exit, the inverse
 468    real(8),           intent(out)   :: eigenValues(:)       ! computed eigenValues
 469    real(8),           intent(out)   :: eigenVectors(:,:)    ! computed eigenVectors
 470    real(8),           intent(in)    :: tolerance            ! threshold for eigenValue magnitude to be returned
 471    integer,           intent(out)   :: numReturned          ! number of eigenValues/Vectors returned
 472    logical, optional, intent(in)    :: printInformation_opt ! switch to print be more verbose
 473
 474    ! Locals:
 475    integer :: rank, index1, index2, info, sizework
 476    real(8) :: sizework_r8
 477    real(8), allocatable :: work(:), eigenVectorsOrig(:,:), eigenValuesOrig(:)
 478    logical :: printInformation
 479
 480    if (present(printInformation_opt)) then
 481      printInformation = printInformation_opt
 482    else
 483      printInformation = .false.
 484    end if
 485
 486    if (printInformation) then
 487      write(*,*)' utl_eigenDecomp: Eigen decomposition of a symmetric matrix'
 488    end if
 489
 490    !     1. Computation of eigenvalues and eigenvectors
 491
 492    rank = size(matrix,1)
 493    allocate(eigenVectorsOrig(rank,rank))
 494    allocate(eigenValuesOrig(rank))
 495
 496    do index2 = 1, rank
 497      do index1 = 1, rank
 498        eigenVectorsOrig(index1,index2)=matrix(index1,index2)
 499      end do
 500    end do
 501
 502    ! Query the size of the 'work' vector by calling 'DSYEV' with 'sizework=-1'
 503    sizework = -1
 504    info = -1
 505    call dsyev('V', 'U', rank, eigenVectorsOrig, rank, eigenValuesOrig,  &
 506               sizework_r8, sizework, info)
 507
 508    ! Compute the eigenvalues/vectors
 509    sizework = int(sizework_r8)
 510    allocate(work(sizework))
 511    call dsyev('V', 'U', rank, eigenVectorsOrig, rank, eigenValuesOrig,  &
 512               work, sizework, info)
 513    deallocate(work)
 514
 515    if (printInformation) then
 516      write(*,'(1x,"Original eigen values: ")')
 517      write(*,'(1x,10f7.3)') (eigenValuesOrig(index1),index1=1,rank)
 518
 519      if(minval(eigenValuesOrig) > tolerance) then
 520        write(*,'(A,1x,e14.6)') "Condition number:", &
 521             maxval(eigenValuesOrig)/minval(eigenValuesOrig)
 522      end if
 523    end if
 524
 525    !     2.  Determine which eigen values/vectors to return
 526
 527    numReturned = 0
 528    do index1 = rank, 1, -1
 529      if (eigenValuesOrig(index1) > tolerance) then
 530        numReturned = numReturned + 1
 531      else
 532        exit
 533      end if
 534    end do
 535
 536    if (printInformation) then
 537      write(*,*) 'Number of eigen values returned =', numReturned, ' out of', rank
 538    end if
 539
 540    !     3.  Copy eigenValues/Vectors into output arrays with reversed order
 541
 542    do index1 = 1, numReturned
 543      ! And set negative values to zero
 544      eigenValues(index1) = max(0.0D0,eigenValuesOrig(rank-index1+1))
 545    end do
 546    do index1 = numReturned+1, rank
 547      eigenValues(index1) = 0.0D0
 548    end do
 549
 550    do index2 = 1, numReturned
 551      do index1 = 1, rank
 552        eigenVectors(index1,index2) = eigenVectorsOrig(index1,rank-index2+1)
 553      end do
 554    end do
 555    do index2 = numReturned+1, rank
 556      do index1 = 1, rank
 557        eigenVectors(index1,index2) = 0.0D0
 558      end do
 559    end do
 560
 561    !     4. Deallocate local arrays
 562    deallocate(eigenVectorsOrig,eigenValuesOrig)
 563
 564    if (printInformation) then
 565      write(*,*) 'utl_eigenDecomp: done'
 566      write(*,*) ' '
 567    end if
 568
 569  end subroutine utl_eigenDecomp
 570
 571  !-----------------------------------------
 572  ! utl_pseudo_inverse
 573  !-----------------------------------------
 574  subroutine utl_pseudo_inverse(inputMatrix, pseudoInverse, threshold_opt)
 575    !
 576    !:Purpose: to calculate the More-Penrose pseudo inverse of the matrix inputMatrix
 577    !
 578    implicit none
 579
 580    ! Arguments:
 581    real(8),           intent(in)  :: inputMatrix(:,:)   ! Input Matrix
 582    real(8),           intent(out) :: pseudoInverse(:,:) ! its Moore Penrose Pseudo-Inverse
 583    real(8), optional, intent(in)  :: threshold_opt 
 584
 585    ! Locals:
 586    real(8), allocatable :: copyMatrix(:,:), leftSingularVector(:,:), rightSingularVectorT(:,:)
 587    real(8), allocatable :: singularValues(:)
 588    integer :: info, lwork, lineIndex
 589    integer :: lineDim, columnDim, minDim
 590    real(8) :: thresh
 591    real(8), allocatable :: work(:)
 592    character(len=80) :: errorMessage
 593
 594    lineDim = size(inputMatrix, dim=1)
 595    columnDim = size(inputMatrix, dim=2)
 596    minDim = min(lineDim, columnDim)
 597
 598    allocate( copyMatrix(lineDim,columnDim), leftSingularVector(lineDim,lineDim) )
 599    allocate( rightSingularVectorT(columnDim,columnDim), singularValues(minDim) )
 600
 601    copyMatrix(:,:) = inputMatrix(:,:) ! Work with a copy because copyMatrix will be overwriten
 602    lwork = max(10000, max(1, 3 * min(lineDim,columnDim) + max(lineDim,columnDim), 5 * minDim ))
 603    allocate(work(lwork))
 604    call dgesvd("A", "A", lineDim, columnDim, copyMatrix, lineDim, singularValues, &
 605         leftSingularVector, lineDim, rightSingularVectorT, columnDim, work, lwork, info ) 
 606   
 607    if (info /= 0) then
 608      write(errorMessage,*) "utl_pseudo_inverse: Problem in DGESVD ! ",info
 609      call utl_abort(errorMessage)
 610    end if
 611
 612    deallocate(work)
 613
 614    if (present(threshold_opt)) then
 615      thresh = threshold_opt
 616    else
 617      !according to wikipedia... as in matlab or numpy
 618      thresh = epsilon(thresh) * max(lineDim, columnDim) * maxval(singularValues)
 619    end if
 620    print *,"utl_pseudo_inverse: threshold= ",thresh
 621
 622    pseudoInverse(:,:)=0.d0
 623    do lineIndex = 1, minDim
 624      If (singularValues(lineIndex) > thresh) then
 625        pseudoInverse(lineIndex,:) = ( 1.d0 / singularValues(lineIndex) ) * leftSingularVector(:,lineIndex)
 626      end if
 627    end do
 628
 629    pseudoInverse = matmul( transpose(rightSingularVectorT), pseudoInverse)
 630
 631    deallocate( singularValues, rightSingularVectorT )
 632    deallocate( leftSingularVector, copyMatrix )
 633
 634  end subroutine utl_pseudo_inverse
 635
 636
 637  !--------------------------------------------------------------------------
 638  ! utl_writeStatus
 639  !--------------------------------------------------------------------------
 640  subroutine utl_writeStatus(cmsg)
 641    implicit none
 642
 643    ! Arguments:
 644    character(len=*), intent(in) :: cmsg
 645
 646    ! Locals:
 647    INTEGER :: iulstatus,fnom,fclos, ierr
 648    character(len=22):: clmsg
 649
 650    clmsg='VAR3D_STATUS='//cmsg
 651    iulstatus = 0
 652    IERR =  FNOM(iulstatus,'VAR3D_STATUS.dot','SEQ+FMT',0)
 653    rewind (iulstatus)
 654    WRITE(iulstatus,'(a22)') clmsg
 655    ierr = fclos(iulstatus)
 656
 657  end subroutine utl_writeStatus
 658
 659
 660  subroutine utl_getfldprm(kip1s,kip2,kip3,knlev,cdetiket,cdtypvar,kgid, &
 661                           cdvar,kstampv,knmaxlev,kinmpg,kip1style,kip1kind, &
 662                           ktrials,koutmpg)
 663    !
 664    !:Purpose:  Get 3D grid parameters for a specific trial field
 665    !           and check for consitancies between grid parameters
 666    !           of the levels.
 667    !
 668    implicit none
 669
 670    ! Arguments:
 671    integer,          intent(out) :: kstampv
 672    integer,          intent(in)  :: knmaxlev
 673    integer,          intent(out) :: knlev
 674    integer,          intent(out) :: kgid
 675    integer,          intent(out) :: kip1s(knmaxlev)
 676    integer,          intent(out) :: kip1style
 677    integer,          intent(out) :: kip1kind
 678    integer,          intent(out) :: kip2
 679    integer,          intent(out) :: kip3
 680    integer,          intent(in)  :: ktrials
 681    integer,          intent(out) :: koutmpg
 682    integer,          intent(in)  :: kinmpg(ktrials)
 683    character(len=*), intent(out) :: cdtypvar
 684    character(len=*), intent(out) :: cdvar
 685    character(len=*), intent(out) :: cdetiket
 686
 687    ! Locals:
 688    integer :: fstinl,fstprm,ezqkdef,newdate
 689    integer :: ini,inj,ink,jlev,ier
 690    integer :: idateo, idateo2, idatyp, idatyp2, ideet, ideet2, idltf
 691    integer :: iextra1, iextra2, iextra3, iig12, iig22
 692    integer :: iig32, iig42, ilng, inbits,iig1,iig2,iig3,iig4
 693    integer :: inpas,inpas2, iswa, iubc, iip2, iip3
 694    integer :: ipmode,idate2,idate3,idatefull
 695    integer :: k,ier1 
 696    real(4) :: zlev_r4
 697    character(len=12) :: cletiket
 698    character(len=4) :: clnomvar
 699    character(len=3) :: clnomvar_3
 700    character(len=2) :: cltypvar
 701    character(len=1) :: clgrtyp2,clgrtyp,clstring
 702    logical :: llflag
 703    integer :: ikeys(knmaxlev)
 704
 705    knlev = 0
 706
 707    do k=1,ktrials
 708       if(cdvar.eq.'U1') then
 709          clnomvar_3='UT1'
 710          ier = fstinl(kinmpg(k),INI,INJ, INK, kstampv, ' ', -1, -1, -1, &
 711               ' ',clnomvar_3,ikeys, knlev, knmaxlev)
 712       else if(cdvar.eq.'V1') then
 713          clnomvar_3='VT1'
 714          ier = fstinl(kinmpg(k),ini,inj, ink, kstampv, ' ', -1, -1, -1, &
 715               ' ',clnomvar_3,ikeys, knlev, knmaxlev)
 716       else
 717          ier = fstinl(kinmpg(k),INI,INJ, INK, kstampv, ' ', -1, -1, -1, &
 718               ' ',cdvar,IKEYS, KNLEV, knmaxlev)
 719       end if
 720       !
 721       if(knlev > 0 ) then
 722          ier1   = newdate(kstampv,idate2,idate3,-3)
 723
 724          idatefull = idate2*100 + idate3/1000000
 725          idateo = -9999
 726          ideet = -9999
 727          inpas = -9999
 728          cdetiket = '-9999999'
 729          clgrtyp = '-'
 730          kip2 = -9999
 731          kip3 = -9999
 732          cdtypvar = '-'
 733          idatyp = -9999
 734          iig1 = -9999
 735          iig2 = -9999
 736          iig3 = -9999
 737          iig4 = -9999
 738          llflag = .true.
 739          koutmpg = kinmpg(k) 
 740          exit 
 741       end if
 742    end do ! End of loop k   
 743    !
 744    if (knlev.gt.0) then
 745       do jlev = 1, knlev
 746          ier = fstprm(ikeys(jlev), idateo2, ideet2, inpas2, ini, inj, &
 747               ink,inbits,idatyp2, kip1s(jlev),iip2, iip3, &
 748               cltypvar,clnomvar,cletiket,clgrtyp2, iig12, iig22,iig32 &
 749               ,iig42,iswa,ilng,idltf,iubc,iextra1, iextra2, iextra3)
 750          llflag = (llflag.and.(idateo.eq.idateo2.or.idateo.eq.-9999))
 751          llflag = (llflag.and.(ideet.eq.ideet2.or.ideet.eq.-9999))
 752          llflag = (llflag.and.(inpas.eq.inpas2.or.inpas.eq.-9999))
 753          !          llflag = (llflag.and.(cdetiket.eq.cletiket.or.cdetiket.eq.
 754          !     &         '-9999999'))
 755          llflag = (llflag.and.(clgrtyp.eq.clgrtyp2.or.clgrtyp.eq.'-'))
 756          llflag = (llflag.and.(kip2.eq.iip2.or.kip2.eq.-9999))
 757          llflag = (llflag.and.(kip3.eq.iip3.or.kip3.eq.-9999))
 758          llflag = (llflag.and.(cdtypvar.eq.cltypvar.or.cdtypvar.eq.'-'))
 759          llflag = (llflag.and.(idatyp.eq.idatyp2.or.idatyp.eq.-9999))
 760          llflag = (llflag.and.(iig1.eq.iig12.or.iig1.eq.-9999))
 761          llflag = (llflag.and.(iig2.eq.iig22.or.iig2.eq.-9999))
 762          llflag = (llflag.and.(iig3.eq.iig32.or.iig3.eq.-9999))
 763          llflag = (llflag.and.(iig4.eq.iig42.or.iig4.eq.-9999))
 764          if (llflag) then
 765             idateo = idateo2
 766             ideet = ideet2
 767             inpas = inpas2
 768             cdetiket = cletiket
 769             clgrtyp = clgrtyp2
 770             kip2 = iip2
 771             kip3 = iip3
 772             cdtypvar = cltypvar
 773             idatyp = idatyp2
 774             iig1 = iig12
 775             iig2 = iig22
 776             iig3 = iig32
 777             iig4 = iig42
 778          else
 779             write(*,*) &
 780                  '****** Unit ', kinmpg &
 781                  ,' contains mixed dateo,deet,npas,etiket,grtyp,ip2,ip3' &
 782                  ,',typvar,datyp,ig1,ig2,ig3 and/or ig4 ' &
 783                  ,'for variable ',cdvar,' and datev, ',kstampv
 784             call utl_abort('GETFLDPRM2')
 785          end if
 786       end do
 787       !
 788       kgid = ezqkdef(ini,inj,clgrtyp,iig1,iig2,iig3,iig4,koutmpg)
 789       !
 790       !-------Determine the style in which ip1 is encoded (15bits or 31 bits)
 791       !       A value <= 32767 (2**16 -1)  means that ip1 is compacted in 15 bits
 792       !       Determine the type of P which was encoded in IP1
 793       !
 794       if(kip1s(1) .le. 32767) then
 795          kip1style = 3
 796       else
 797          kip1style = 2
 798       end if
 799       !
 800       !-------Determine the type of P  (see doc. of convip)
 801       !
 802       ipmode = -1
 803       call CONVIP(kip1s(1),zlev_r4,KIP1KIND, &
 804            ipmode,clstring, .false. )
 805    else
 806       do k=1,ktrials
 807          ier = fstinl(kinmpg(k),ini,inj, ink, -1, ' ', -1, -1, -1, &
 808               ' ',cdvar,ikeys, knlev, knmaxlev)
 809       end do
 810       write(*,*) 'Error - getfldprm2: no record found at time ' &
 811            ,idatefull,' for field ',cdvar,' but',knlev, &
 812            ' records found in unit ',kinmpg(k)
 813       call utl_abort('GETFLDPRM2')
 814    end if
 815    !
 816  end subroutine utl_getfldprm
 817
 818
 819  subroutine utl_abort(message)
 820    implicit none
 821
 822    ! Arguments:
 823    character(len=*), intent(in) :: message
 824
 825    ! Locals:
 826    integer :: comm, ierr, rpn_comm_comm
 827    
 828    write(6,9000) message
 8299000 format(//,4X,"!!!---ABORT---!!!",/,8X,"MIDAS stopped in ",A)
 830    call flush(6)
 831
 832    comm = rpn_comm_comm("WORLD")
 833    call mpi_abort( comm, 1, ierr )
 834
 835  end subroutine utl_abort
 836
 837  !--------------------------------------------------------------------------
 838  ! utl_stopAndWait4Debug
 839  !--------------------------------------------------------------------------
 840  subroutine utl_stopAndWait4Debug(message)
 841    !
 842    !:Purpose: Stop the execution for the process reaching a call to the
 843    !          subroutine, then wait until all MPI processes reached such a
 844    !          call to utl_stopAndWait4Debug.
 845    !          Intended **for debugging puposes only** since it can cause
 846    !          unwanted MPI deadlocks - processes waiting infinitely because
 847    !          not all MPI processes will ever reach a call to
 848    !          utl_stopAndWait4Debug.
 849    !
 850    implicit none
 851
 852    ! Arguments:
 853    character(len=*), intent(in) :: message
 854
 855    ! Locals:
 856    integer :: comm, ierr, rpn_comm_comm
 857
 858    write(6,9000) message
 8599000 format(//,4X,"!!!---ALL STOP---!!!",/,8X,"Debugging message: ",A)
 860    call flush(6)
 861
 862    call rpn_comm_barrier('WORLD', ierr)
 863    comm = rpn_comm_comm("WORLD")
 864    call mpi_abort( comm, 1, ierr )
 865  end subroutine utl_stopAndWait4Debug
 866
 867  subroutine utl_open_asciifile(filename,unit)
 868    ! 
 869    !:Purpose: Opens an ascii file for output 
 870    !
 871    implicit none
 872
 873    ! Arguments:
 874    character(len=*), intent(in)  :: filename
 875    integer,          intent(out) :: unit
 876
 877    ! Locals:
 878    logical :: file_exists
 879    integer :: ier
 880    character(len=20) :: mode
 881    
 882    inquire(file=trim(filename), exist=file_exists)
 883    
 884    if (file_exists) then
 885       mode = 'FTN+APPEND+R/W'
 886    else
 887       mode = 'FTN+R/W'
 888    end if
 889
 890    unit=0
 891    
 892    ier = utl_open_file(unit,trim(filename),trim(mode))
 893
 894    if (ier.ne.0) call utl_abort('utl_open_messagefile: Error associating unit number')
 895
 896  end subroutine utl_open_asciifile
 897
 898
 899  function utl_open_file(unit,filename,mode) result(ier)
 900    ! 
 901    !:Purpose: This is a temporary subroutine to open a file with fnom that is needed due to
 902    !          a bug in fnom that does not allow an ascii file to be opened in 'APPEND' mode.  
 903    !
 904    implicit none
 905
 906    ! Arguments:
 907    integer,          intent(inout) :: unit
 908    character(len=*), intent(in)    :: filename
 909    character(len=*), intent(in)    :: mode
 910    ! Result:
 911    integer :: ier
 912
 913    ! Locals:
 914    character(len=10) :: position,action
 915    integer :: fnom
 916
 917    if (index(mode,'APPEND').gt.0) then
 918       position = 'APPEND'
 919    else
 920       position = 'ASIS'
 921    end if
 922    
 923    if (index(mode,'R/W').gt.0) then
 924       action = 'READWRITE'
 925    else
 926       action = 'READ'
 927    end if
 928
 929    ier = fnom(unit,filename,mode,0)
 930    
 931    close(unit=unit)
 932    open(unit=unit, file=filename, position=position, action=action)
 933
 934  end function utl_open_file
 935    
 936
 937  function utl_stnid_equal(id1,id2) result(same)
 938    !
 939    !:Purpose: Compares STNID values allowing for * as wildcards and trailing blanks 
 940    !
 941    !:Arguments:
 942    !           :id1: reference stnid
 943    !           :id2: stnid being verified
 944    !           :same: logical indicating if id1 and id2 match
 945    !     
 946    implicit none
 947
 948    ! Arguments:
 949    character(len=*), intent(in) :: id1
 950    character(len=*), intent(in) :: id2
 951    ! Result:
 952    logical :: same
 953
 954    ! Locals:
 955    integer :: ilen1,ilen2,ji
 956
 957    same=.true.
 958    ilen1=len_trim(id1)
 959    ilen2=len_trim(id2)  
 960              
 961    do ji=1,min(ilen1,ilen2)
 962       if ( id1(ji:ji).ne.'*' .and. id2(ji:ji).ne.'*' .and. id2(ji:ji).ne.id1(ji:ji) ) then
 963          same = .false.
 964          exit
 965       end if
 966    end do
 967    
 968    if (same.and.ilen1.gt.ilen2) then
 969       do ji=ilen2+1,ilen1
 970          if (id1(ji:ji).ne.'*') then
 971              same=.false.
 972              exit
 973          end if
 974       end do
 975    else if (same.and.ilen2.gt.ilen1) then
 976       do ji=ilen1+1,ilen2
 977          if (id2(ji:ji).ne.'*') then
 978              same=.false.
 979              exit
 980          end if
 981       end do
 982    end if
 983        
 984  end function utl_stnid_equal
 985 
 986
 987  character(len=20) function utl_int2str(i)
 988    !
 989    !:Purpose: Function for integer to string conversion. Helpful when calling subroutine utl_abort. 
 990    !
 991    implicit none
 992
 993    ! Arguments:
 994    integer, intent(in) :: i
 995    
 996    write(utl_int2str,*) i
 997    utl_int2str = adjustl(utl_int2str)
 998    
 999  end function utl_int2str
1000            
1001
1002  character(len=20) function utl_float2str(x)
1003    !
1004    !:Purpose: Function for integer to string conversion. Helpful when calling subroutine utl_abort.
1005    !
1006    implicit none
1007
1008    ! Arguments:
1009    real(8), intent(in) :: x
1010
1011    write(utl_float2str,*) x
1012    utl_float2str = adjustl(utl_float2str)
1013
1014  end function utl_float2str
1015
1016
1017  subroutine utl_resize_1d_real(arr,dim1)
1018    !
1019    !:Purpose: Resize 1D array
1020    !
1021    implicit none
1022
1023    ! Arguments:
1024    real(8), pointer, intent(inout) :: arr(:)
1025    integer,          intent(in)    :: dim1
1026
1027    ! Locals:
1028    real(8), pointer :: tmp(:)
1029    integer :: dim1_in,d1
1030
1031    dim1_in = size(arr)
1032    d1 = min(dim1_in, dim1)
1033
1034    allocate(tmp(dim1))
1035    tmp(1:d1) = arr(1:d1)
1036
1037    if (dim1.gt.dim1_in) tmp(d1+1:dim1) = 0.0D0
1038    
1039    deallocate(arr)
1040
1041    arr => tmp
1042    
1043    nullify(tmp)
1044
1045  end subroutine utl_resize_1d_real
1046
1047
1048  subroutine utl_resize_1d_int(arr,dim1)
1049    !
1050    !:Purpose: Resize 1D array 
1051    !
1052    implicit none
1053
1054    ! Arguments:
1055    integer, pointer, intent(inout) :: arr(:)
1056    integer,          intent(in)    :: dim1
1057
1058    ! Locals:
1059    integer, pointer :: tmp(:)
1060    integer :: dim1_in,d1
1061
1062    dim1_in = size(arr)
1063    d1 = min(dim1_in, dim1)
1064
1065    allocate(tmp(dim1))
1066    tmp(1:d1) = arr(1:d1)
1067
1068    if (dim1.gt.dim1_in) tmp(d1+1:dim1) = 0
1069    
1070    deallocate(arr)
1071
1072    arr => tmp
1073    
1074    nullify(tmp)
1075
1076  end subroutine utl_resize_1d_int
1077
1078 
1079  subroutine utl_resize_1d_str(arr,dim1)
1080    !
1081    !:Purpose: Resize 1D array
1082    !
1083    implicit none
1084
1085    ! Arguments:
1086    character(len=*), pointer, intent(inout) :: arr(:)
1087    integer,                   intent(in)    :: dim1
1088
1089    ! Locals:
1090    character(len=len(arr(1))), pointer :: tmp(:)
1091    integer :: dim1_in,d1
1092
1093    dim1_in = size(arr)
1094    d1 = min(dim1_in, dim1)
1095
1096    allocate(tmp(dim1))
1097    tmp(1:d1) = arr(1:d1)
1098
1099    if (dim1.gt.dim1_in) tmp(d1+1:dim1) = ""
1100    
1101    deallocate(arr)
1102    arr => tmp
1103    nullify(tmp)
1104
1105  end subroutine utl_resize_1d_str
1106
1107
1108  subroutine utl_resize_2d_real(arr,dim1,dim2)
1109    !
1110    !:Purpose: Resize 2D array
1111    !
1112    implicit none
1113
1114    ! Arguments:
1115    real(8), pointer, intent(inout) :: arr(:,:)
1116    integer,          intent(in)    :: dim1
1117    integer,          intent(in)    :: dim2
1118
1119    ! Locals:
1120    real(8), pointer :: tmp(:,:)
1121    integer :: dim1_in,dim2_in,d1,d2
1122
1123    dim1_in = size(arr,dim=1)
1124    dim2_in = size(arr,dim=2)
1125    d1 = min(dim1_in, dim1)
1126    d2 = min(dim2_in, dim2)
1127
1128    allocate(tmp(dim1,dim2))
1129    tmp(1:d1,1:d2) = arr(1:d1,1:d2)
1130
1131    if (dim1.gt.dim1_in) tmp(d1+1:dim1,:) = 0.0D0
1132    if (dim2.gt.dim2_in) tmp(:,d2+1:dim2) = 0.0D0
1133      
1134    deallocate(arr)
1135
1136    arr => tmp
1137    
1138    nullify(tmp)
1139
1140  end subroutine utl_resize_2d_real
1141
1142
1143  subroutine utl_resize_3d_real(arr,dim1,dim2,dim3)
1144    !
1145    !:Purpose: Resize 3D array
1146    !
1147    implicit none
1148
1149    ! Arguments:
1150    real(8), pointer, intent(inout) :: arr(:,:,:)
1151    integer,          intent(in)    :: dim1
1152    integer,          intent(in)    :: dim2
1153    integer,          intent(in)    :: dim3
1154
1155    ! Locals:
1156    real(8), pointer :: tmp(:,:,:)
1157    integer :: dim1_in,dim2_in,dim3_in,d1,d2,d3
1158
1159    dim1_in = size(arr,dim=1)
1160    dim2_in = size(arr,dim=2)
1161    dim3_in = size(arr,dim=3)
1162    d1 = min(dim1_in, dim1)
1163    d2 = min(dim2_in, dim2)
1164    d3 = min(dim3_in, dim3)
1165
1166    allocate(tmp(dim1,dim2,dim3))
1167    tmp(1:d1,1:d2,1:d3) = arr(1:d1,1:d2,1:d3)
1168
1169    if (dim1.gt.dim1_in) tmp(d1+1:dim1,:,:) = 0.0D0
1170    if (dim2.gt.dim2_in) tmp(:,d2+1:dim2,:) = 0.0D0
1171    if (dim3.gt.dim3_in) tmp(:,:,d3+1:dim3) = 0.0D0
1172    
1173    deallocate(arr)
1174
1175    arr => tmp
1176    
1177    nullify(tmp)
1178
1179  end subroutine utl_resize_3d_real
1180
1181
1182  subroutine utl_get_stringId(cstringin,nobslev,CList,NListSize,Nmax,elemId)
1183    !
1184    !:Purpose: Get element ID from a list of accumulating character strings (e.g. stnids). 
1185    !          Called by filt_topoChm in filterobs_mod.ftn90
1186    !
1187    implicit none
1188
1189    ! Arguments:
1190    integer,           intent(in)    :: Nmax
1191    integer,           intent(in)    :: nobslev
1192    integer,           intent(inout) :: NListSize
1193    integer,           intent(out)   :: elemId
1194    character(len=*),  intent(in)    :: cstringin
1195    character(len=*),  intent(inout) :: CList(Nmax)
1196
1197    ! Locals:
1198    integer :: i
1199    character(len=120) :: cstring
1200    
1201    elemId=0
1202    if (NListSize.gt.Nmax-1) then
1203       call utl_abort('utl_get_stringId: Dimension error, NListSize > Nmax-1.')     
1204    else if (NListSize.gt.0) then
1205       if (nobslev.eq.1) then 
1206          cstring=trim(cstringin)//'U'
1207          do i=1,NListSize
1208             if (trim(cstring).eq.trim(CList(i))) then
1209                 elemId=i
1210                 exit
1211             end if
1212          end do
1213       else 
1214          cstring=trim(cstringin)       
1215          do i=1,NListSize
1216             if (trim(cstring).eq.trim(CList(i))) then
1217                 elemId=i
1218                 exit
1219             end if
1220          end do
1221       end if
1222       
1223       if (elemId.eq.0) then
1224          do i=1,NListSize
1225             if (utl_stnid_equal(trim(CList(i)),trim(cstring))) then
1226                elemId=i
1227                exit
1228             end if
1229          end do
1230       end if
1231    end if
1232
1233    if (elemID.eq.0) then
1234        NListSize=NListSize+1
1235        elemId=NListSize
1236        if (nobslev.eq.1) then
1237           CList(NListSize)=trim(cstringin)//'U'
1238        else
1239           CList(NListSize)=trim(cstringin)
1240        end if
1241    end if
1242    
1243  end subroutine utl_get_stringId
1244
1245
1246  subroutine utl_get_Id(id,IdList,NListSize,Nmax,elemId)
1247    !
1248    !:Purpose: Get element ID from list of accumulating integer IDs.
1249    !
1250    implicit none
1251
1252    ! Arguments:
1253    integer, intent(in)    :: Nmax
1254    integer, intent(in)    :: id
1255    integer, intent(inout) :: NListSize
1256    integer, intent(inout) :: IdList(Nmax)
1257    integer, intent(out)   :: elemId
1258
1259    ! Locals:
1260    integer :: i
1261    
1262    elemId=0
1263    if (NListSize.gt.Nmax-1) then
1264       call utl_abort('utl_get_Id: Dimension error, NListSize > Nmax-1.')     
1265    else if (NListSize.gt.0) then
1266       do i=1,NListSize
1267          if (id.eq.IdList(i)) then
1268              elemId=i
1269              exit
1270          end if
1271       end do
1272    end if
1273
1274    if (elemID.eq.0) then
1275        NListSize=NListSize+1
1276        elemId=NListSize
1277        IdList(NListSize)=id
1278    end if
1279    
1280    
1281  end subroutine utl_get_Id
1282  
1283
1284  subroutine utl_readFstField( fname, varName, iip1, iip2, iip3, etiketi, &
1285                               ni, nj, nkeys, array, xlat_opt, xlong_opt, lvls_opt, kind_opt )
1286    !
1287    !:Purpose:  Read specified field from standard RPN/fst file. Could be one
1288    !           to all levels depending on the input iip1,iip2,iip3 values.
1289    !
1290    !           Currently assumes lat/long (or Gaussian) type grids.
1291    !           See hco_SetupFromFile for example toward future generalizations.
1292    !           Generalization would require having xlat and xlong being 2D.
1293    !
1294    !:Arguments:
1295    !                 :fname: input filename
1296    !                 :varName:  search nomvar
1297    !                 :iip1: search ip1
1298    !                 :iip2: search ip2
1299    !                 :iip3: search ip3
1300    !                 :etiketi: search etiket
1301    !                 :ni: ni values
1302    !                 :nj: nj values
1303    !                 :nkeys: number of records satisfying search criteria
1304    !                 :array: data arrray
1305    !                 :xlat_opt: 1D latitude array (optional)
1306    !                 :xlong_opt: 1D longitude array (optional)
1307    !                 :lvls_opt: 1D vertical coordinate array (optional)
1308    !                 :kind_opt: vertical coordinate type according to convip (optional)
1309    !
1310    implicit none
1311
1312    ! Arguments:
1313    integer,                        intent(in)  :: iip1
1314    integer,                        intent(in)  :: iip2
1315    integer,                        intent(in)  :: iip3
1316    character(len=*),               intent(in)  :: varName
1317    character(len=*),               intent(in)  :: fname
1318    character(len=*),               intent(in)  :: etiketi
1319    integer,                        intent(out) :: ni
1320    integer,                        intent(out) :: nj
1321    integer,                        intent(out) :: nkeys
1322    integer, optional,              intent(out) :: kind_opt
1323    real(8), allocatable,           intent(out) :: array(:,:,:)
1324    real(8), allocatable, optional, intent(out) :: lvls_opt(:)
1325    real(8), allocatable, optional, intent(out) :: xlat_opt(:)
1326    real(8), allocatable, optional, intent(out) :: xlong_opt(:)
1327
1328    ! Locals:
1329    integer, external :: fnom,fclos,fstouv,fstfrm,fstinl,fstlir,fstluk,fstprm
1330    real(4) :: lvl_r4
1331    logical :: Exists
1332    character(len=1) :: string
1333    integer, parameter :: iun=0
1334    integer :: i,ier, kindi
1335    integer, parameter :: maxkeys=1000
1336    integer :: keys(maxkeys),ini,inj,nk
1337    integer :: dateo, deet, npas, nbits, datyp
1338    integer :: ip1, ip2, ip3, swa, lng, dltf, ubc
1339    integer :: extra1, extra2, extra3
1340    integer :: ig1, ig2, ig3, ig4  
1341    character*1 clgrtyp
1342    character*2 cltypvar
1343    character*4 nomvar
1344    character*12 cletiket
1345    real(4), allocatable :: buffer(:,:)
1346    real :: xlat1_4, xlon1_4, xlat2_4, xlon2_4
1347    
1348    ! Open file
1349    inquire(file=trim(fname),exist=Exists)
1350    if(.not.Exists) then
1351      write(*,*) 'File missing=',fname
1352      call utl_abort('utl_read_fst_field: did not find file.')
1353    else
1354      ier=fnom(iun,trim(fname),'RND+OLD+R/O',0)
1355      ier=fstouv(iun,'RND+OLD')
1356    end if
1357
1358    ! Find reports in file for specified varName and iip*.
1359    ier = fstinl(iun,ni,nj,nk,-1,etiketi,iip1,iip2,iip3,'',varName,keys,nkeys,maxkeys) 
1360
1361    if(ier.lt.0.or.nkeys.eq.0) then
1362      write(*,*) 'Search field missing ',varName, ' from file ',fname
1363      call utl_abort('utl_read_fst_field: did not find field.')
1364    else if (nk.gt.1) then
1365      write(*,*) 'Unexpected size nk ',nk,' for ',varName,' of file ',fname 
1366      call utl_abort('utl_read_fst_field')      
1367    end if
1368
1369    if (present(xlat_opt).and.present(xlong_opt)) then
1370
1371       !  Get lat and long if available.
1372
1373       if (allocated(xlat_opt)) deallocate(xlat_opt,xlong_opt)
1374       allocate(xlat_opt(nj),xlong_opt(ni),buffer(ni*nj,1))
1375       xlat_opt(:)=-999.
1376       xlong_opt(:)=-999.
1377
1378       ier = fstprm(keys(1),dateo, deet, npas, ni, nj, nk, nbits,    &         
1379                    datyp, ip1, ip2, ip3, cltypvar, nomvar, cletiket, &
1380                    clgrtyp, ig1, ig2, ig3,                           &
1381                    ig4, swa, lng, dltf, ubc, extra1, extra2, extra3)  
1382    
1383       if (ni.gt.1) then
1384          ier=fstlir(buffer,iun,ni,inj,nk,-1,'',ig1,ig2,ig3,'','>>')
1385          if (ier.ge.0) xlong_opt(:)=buffer(1:ni,1) 
1386       end if
1387       if (nj.gt.1) then
1388          ier=fstlir(buffer,iun,ini,nj,nk,-1,'',ig1,ig2,ig3,'','^^')
1389          if (ier.ge.0) xlat_opt(:)=buffer(1:nj,1)     
1390       end if
1391       deallocate(buffer)
1392
1393       if ( trim(clgrtyp) == 'Z' ) then
1394
1395          ! Check for rotated grid
1396
1397          ier = fstprm(ier,                                         & ! IN
1398                  dateo, deet, npas, ni, nj, nk, nbits,             & ! OUT
1399                  datyp, ip1, ip2, ip3, cltypvar, nomvar, cletiket, & ! OUT
1400                  clgrtyp, ig1, ig2, ig3,                           & ! OUT
1401                  ig4, swa, lng, dltf, ubc, extra1, extra2, extra3 ) ! OUT
1402
1403          call cigaxg (clgrtyp,                                     & ! IN
1404                     xlat1_4, xlon1_4, xlat2_4, xlon2_4,            & ! OUT
1405                     ig1, ig2, ig3, ig4 ) ! IN
1406
1407          if ( xlat1_4 /= xlat2_4 .or. xlon1_4 /= xlon2_4 ) &
1408             call utl_abort('utl_readFstField: Cannot currently handle rotated grid')
1409          
1410       else if (trim(clgrtyp) /= 'G') then
1411
1412         call utl_abort('utl_readFstField: Cannot currently handle grid type ' // trim(clgrtyp) )
1413
1414       end if
1415       
1416    end if 
1417    
1418    ! Get vertical coordinate
1419    
1420    if (present(lvls_opt)) then
1421       if (allocated(lvls_opt)) deallocate(lvls_opt)
1422       allocate(lvls_opt(nkeys))
1423
1424       do i=1,nkeys
1425          ier = fstprm(keys(i),dateo, deet, npas, ni, nj, nk, nbits,    &         
1426                    datyp, ip1, ip2, ip3, cltypvar, nomvar, cletiket, &
1427                    clgrtyp, ig1, ig2, ig3,                           &
1428                    ig4, swa, lng, dltf, ubc, extra1, extra2, extra3)  
1429          call convip(ip1,lvl_r4,kindi,-1,string,.false.)
1430          lvls_opt(i)=lvl_r4
1431       end do
1432    end if
1433     
1434    if (present(kind_opt)) then
1435        if (present(lvls_opt)) then         
1436            kind_opt=kindi
1437        else
1438            kind_opt=-1
1439        end if
1440    end if
1441
1442    ! Get field
1443    
1444    allocate(array(ni,nj,nkeys),buffer(ni,nj))
1445
1446    do i=1,nkeys
1447       ier=fstluk(buffer,keys(i),ni,nj,nk)
1448       array(:,:,i)=buffer(:,:)
1449    end do
1450    
1451    deallocate(buffer)
1452    
1453    ier=fstfrm(iun)  
1454    ier=fclos(iun)  
1455
1456  end subroutine utl_readFstField
1457
1458
1459  subroutine utl_checkAllocationStatus(status, message, alloc_opt)
1460    implicit none
1461
1462    ! Arguments:
1463    character(len=*),  intent(in) :: message
1464    integer,           intent(in) :: status(:)
1465    logical, optional, intent(in) :: alloc_opt 
1466
1467    ! Locals:
1468    logical :: flag
1469
1470    if ( present(alloc_opt) ) then
1471       flag = alloc_opt
1472    else
1473       flag = .true.
1474    end if
1475
1476    if (any(status /= 0)) then
1477       if (flag) then
1478          Write(6,'(A)') "Memory allocation failure !"
1479       else
1480          Write(6,'(A)') "Memory deallocation failure !"
1481       end if
1482       Write(6,'(64(i8,1x))') status
1483       call utl_abort(message)
1484    end if
1485
1486  end subroutine utl_checkAllocationStatus
1487
1488
1489  function utl_varNamePresentInFile(varName, fileName_opt, fileUnit_opt, typvar_opt) result(found)
1490    implicit none
1491
1492    ! Arguments:
1493    character(len=*),           intent(in) :: varName
1494    character(len=*), optional, intent(in) :: fileName_opt
1495    integer,          optional, intent(in) :: fileUnit_opt
1496    character(len=*), optional, intent(in) :: typvar_opt
1497    ! Result:
1498    logical :: found
1499
1500    ! Locals:
1501    integer :: fnom, fstouv, fstfrm, fclos, fstinf
1502    integer :: ni, nj, nk, key, ierr
1503    integer :: unit
1504    character(len=128) :: fileName
1505    character(len=2)   :: typvar
1506    logical :: openFile
1507
1508    if ( present(fileUnit_opt) ) then
1509      unit = fileUnit_opt
1510      openFile = .false.
1511    else
1512      unit = 0
1513      openFile = .true.
1514      if ( present(fileName_opt) ) then
1515        fileName = fileName_opt
1516      else
1517        call utl_abort('utl_varNamePresentInFile: please provide and file name or unit')
1518      end if
1519    end if
1520
1521    if ( present(typvar_opt) ) then
1522      typvar = trim(typvar_opt)
1523    else
1524      typvar = ' '
1525    end if
1526
1527    if (openFile) then
1528      ierr = fnom(unit,fileName,'RND+OLD+R/O',0)
1529      ierr = fstouv(unit,'RND+OLD')
1530    end if
1531
1532    key = fstinf(unit, ni, nj, nk, -1 ,' ', -1, -1, -1, typvar, trim(varName))
1533    
1534    if ( key > 0 )  then
1535      found = .true.
1536    else
1537      found = .false.
1538    end if
1539
1540    if (openFile) then
1541      ierr =  fstfrm(unit)
1542      ierr =  fclos (unit)
1543    end if
1544
1545  end function utl_varNamePresentInFile
1546
1547
1548  subroutine utl_reAllocate_char_1d(array,dim1)
1549    implicit none
1550
1551    ! Arguments:
1552    character(len=128), allocatable, intent(inout) :: array(:)
1553    integer,                         intent(in) :: dim1
1554
1555    if( allocated(array) ) then
1556      if ( size(array) == dim1 ) then
1557        return
1558      else
1559        deallocate(array)
1560      end if
1561    end if
1562
1563    allocate(array(dim1))
1564
1565  end subroutine utl_reAllocate_char_1d
1566
1567  subroutine utl_reAllocate_char_2d(array,dim1,dim2)
1568    implicit none
1569
1570    ! Arguments:
1571    character(len=128), allocatable, intent(inout) :: array(:,:)
1572    integer,                         intent(in)    :: dim1
1573    integer,                         intent(in)    :: dim2
1574
1575    if( allocated(array) ) then
1576      if ( size(array) == dim1*dim2 ) then
1577        return
1578      else
1579        deallocate(array)
1580      end if
1581    end if
1582
1583    allocate(array(dim1,dim2))
1584
1585  end subroutine utl_reAllocate_char_2d
1586
1587  subroutine utl_reAllocate_char_3d(array,dim1,dim2,dim3)
1588    implicit none
1589
1590    ! Arguments:
1591    character(len=128), allocatable, intent(inout) :: array(:,:,:)
1592    integer,                         intent(in)    :: dim1
1593    integer,                         intent(in)    :: dim2
1594    integer,                         intent(in)    :: dim3
1595
1596    if( allocated(array) ) then
1597      if ( size(array) == dim1*dim2*dim3 ) then
1598        return
1599      else
1600        deallocate(array)
1601      end if
1602    end if
1603
1604    allocate(array(dim1,dim2,dim3))
1605
1606  end subroutine utl_reAllocate_char_3d
1607
1608  subroutine utl_reAllocate_log_1d(array,dim1)
1609    implicit none
1610
1611    ! Arguments:
1612    logical, allocatable, intent(inout) :: array(:)
1613    integer,              intent(in)    :: dim1
1614
1615    if( allocated(array) ) then
1616      if ( size(array) == dim1 ) then
1617        return
1618      else
1619        deallocate(array)
1620      end if
1621    end if
1622
1623    allocate(array(dim1))
1624    array(:) = .true.
1625
1626  end subroutine utl_reAllocate_log_1d
1627
1628  subroutine utl_reAllocate_log_2d(array,dim1,dim2)
1629    implicit none
1630
1631    ! Arguments:
1632    logical, allocatable, intent(inout) :: array(:,:)
1633    integer,              intent(in)    :: dim1
1634    integer,              intent(in)    :: dim2
1635
1636    if( allocated(array) ) then
1637      if ( size(array) == dim1*dim2 ) then
1638        return
1639      else
1640        deallocate(array)
1641      end if
1642    end if
1643
1644    allocate(array(dim1,dim2))
1645    array(:,:) = .true.
1646
1647  end subroutine utl_reAllocate_log_2d
1648
1649  subroutine utl_reAllocate_log_3d(array,dim1,dim2,dim3)
1650    implicit none
1651
1652    ! Arguments:
1653    logical, allocatable, intent(inout) :: array(:,:,:)
1654    integer,              intent(in)    :: dim1
1655    integer,              intent(in)    :: dim2
1656    integer,              intent(in)    :: dim3
1657
1658    if( allocated(array) ) then
1659      if ( size(array) == dim1*dim2*dim3 ) then
1660        return
1661      else
1662        deallocate(array)
1663      end if
1664    end if
1665
1666    allocate(array(dim1,dim2,dim3))
1667    array(:,:,:) = .true.
1668
1669  end subroutine utl_reAllocate_log_3d
1670
1671  subroutine utl_reAllocate_int_1d(array,dim1)
1672    implicit none
1673
1674    ! Arguments:
1675    integer, allocatable, intent(inout) :: array(:)
1676    integer,              intent(in)    :: dim1
1677
1678    if( allocated(array) ) then
1679      if ( size(array) == dim1 ) then
1680        return
1681      else
1682        deallocate(array)
1683      end if
1684    end if
1685
1686    allocate(array(dim1))
1687    array(:) = 0d0
1688
1689  end subroutine utl_reAllocate_int_1d
1690
1691  subroutine utl_reAllocate_int_2d(array,dim1,dim2)
1692    implicit none
1693
1694    ! Arguments:
1695    integer, allocatable, intent(inout) :: array(:,:)
1696    integer,              intent(in)    :: dim1
1697    integer,              intent(in)    :: dim2
1698
1699    if( allocated(array) ) then
1700      if ( size(array) == dim1*dim2 ) then
1701        return
1702      else
1703        deallocate(array)
1704      end if
1705    end if
1706
1707    allocate(array(dim1,dim2))
1708    array(:,:) = 0d0
1709
1710  end subroutine utl_reAllocate_int_2d
1711
1712  subroutine utl_reAllocate_int_3d(array,dim1,dim2,dim3)
1713    implicit none
1714
1715    ! Arguments:
1716    integer, allocatable, intent(inout) :: array(:,:,:)
1717    integer,              intent(in)    :: dim1
1718    integer,              intent(in)    :: dim2
1719    integer,              intent(in)    :: dim3
1720
1721    if( allocated(array) ) then
1722      if ( size(array) == dim1*dim2*dim3 ) then
1723        return
1724      else
1725        deallocate(array)
1726      end if
1727    end if
1728
1729    allocate(array(dim1,dim2,dim3))
1730    array(:,:,:) = 0d0
1731
1732  end subroutine utl_reAllocate_int_3d
1733
1734  subroutine utl_reAllocate_r4_1d(array,dim1)
1735    implicit none
1736
1737    ! Arguments:
1738    real(4), allocatable, intent(inout) :: array(:)
1739    integer,              intent(in)    :: dim1
1740
1741    if( allocated(array) ) then
1742      if ( size(array) == dim1 ) then
1743        return
1744      else
1745        deallocate(array)
1746      end if
1747    end if
1748
1749    allocate(array(dim1))
1750    array(:) = 0.0d0
1751
1752  end subroutine utl_reAllocate_r4_1d
1753
1754  subroutine utl_reAllocate_r8_1d(array,dim1)
1755    implicit none
1756
1757    ! Arguments:
1758    real(8), allocatable, intent(inout) :: array(:)
1759    integer,              intent(in)    :: dim1
1760
1761    if( allocated(array) ) then
1762      if ( size(array) == dim1 ) then
1763        return
1764      else
1765        deallocate(array)
1766      end if
1767    end if
1768
1769    allocate(array(dim1))
1770    array(:) = 0.0d0
1771
1772  end subroutine utl_reAllocate_r8_1d
1773
1774  subroutine utl_reAllocate_r4_2d(array,dim1,dim2)
1775    implicit none
1776
1777    ! Arguments:
1778    real(4), allocatable, intent(inout) :: array(:,:)
1779    integer,              intent(in)    :: dim1
1780    integer,              intent(in)    :: dim2
1781
1782    if( allocated(array) ) then
1783      if ( size(array) == dim1*dim2 ) then
1784        return
1785      else
1786        deallocate(array)
1787      end if
1788    end if
1789
1790    allocate(array(dim1,dim2))
1791    array(:,:) = 0.0d0
1792
1793  end subroutine utl_reAllocate_r4_2d
1794
1795  subroutine utl_reAllocate_r8_2d(array,dim1,dim2)
1796    implicit none
1797
1798    ! Arguments:
1799    real(8), allocatable, intent(inout) :: array(:,:)
1800    integer,              intent(in)    :: dim1
1801    integer,              intent(in)    :: dim2
1802
1803    if( allocated(array) ) then
1804      if ( size(array) == dim1*dim2 ) then
1805        return
1806      else
1807        deallocate(array)
1808      end if
1809    end if
1810
1811    allocate(array(dim1,dim2))
1812    array(:,:) = 0.0d0
1813
1814  end subroutine utl_reAllocate_r8_2d
1815
1816  subroutine utl_reAllocate_r4_3d(array,dim1,dim2,dim3)
1817    implicit none
1818
1819    ! Arguments:
1820    real(4), allocatable, intent(inout) :: array(:,:,:)
1821    integer,              intent(in)    :: dim1
1822    integer,              intent(in)    :: dim2
1823    integer,              intent(in)    :: dim3
1824
1825    if( allocated(array) ) then
1826      if ( size(array) == dim1*dim2*dim3 ) then
1827        return
1828      else
1829        deallocate(array)
1830      end if
1831    end if
1832
1833    allocate(array(dim1,dim2,dim3))
1834    array(:,:,:) = 0.0d0
1835
1836  end subroutine utl_reAllocate_r4_3d
1837
1838
1839  subroutine utl_reAllocate_r8_3d(array,dim1,dim2,dim3)
1840    implicit none
1841
1842    ! Arguments:
1843    real(8), allocatable, intent(inout) :: array(:,:,:)
1844    integer,              intent(in)    :: dim1
1845    integer,              intent(in)    :: dim2
1846    integer,              intent(in)    :: dim3
1847
1848    if( allocated(array) ) then
1849      if ( size(array) == dim1*dim2*dim3 ) then
1850        return
1851      else
1852        deallocate(array)
1853      end if
1854    end if
1855
1856    allocate(array(dim1,dim2,dim3))
1857    array(:,:,:) = 0.0d0
1858
1859  end subroutine utl_reAllocate_r8_3d
1860
1861
1862  subroutine utl_reAllocate_r4_4d(array,dim1,dim2,dim3,dim4)
1863    implicit none
1864
1865    ! Arguments:
1866    real(4), allocatable, intent(inout) :: array(:,:,:,:)
1867    integer,              intent(in)    :: dim1
1868    integer,              intent(in)    :: dim2
1869    integer,              intent(in)    :: dim3
1870    integer,              intent(in)    :: dim4
1871
1872    if( allocated(array) ) then
1873      if ( size(array) == dim1*dim2*dim3*dim4 ) then
1874        return
1875      else
1876        deallocate(array)
1877      end if
1878    end if
1879
1880    allocate(array(dim1,dim2,dim3,dim4))
1881    array(:,:,:,:) = 0.0d0
1882
1883  end subroutine utl_reAllocate_r4_4d
1884
1885
1886  subroutine utl_reAllocate_r8_4d(array,dim1,dim2,dim3,dim4)
1887    implicit none
1888
1889    ! Arguments:
1890    real(8), allocatable, intent(inout) :: array(:,:,:,:)
1891    integer,              intent(in)    :: dim1
1892    integer,              intent(in)    :: dim2
1893    integer,              intent(in)    :: dim3
1894    integer,              intent(in)    :: dim4
1895
1896    if( allocated(array) ) then
1897      if ( size(array) == dim1*dim2*dim3*dim4 ) then
1898        return
1899      else
1900        deallocate(array)
1901      end if
1902    end if
1903
1904    allocate(array(dim1,dim2,dim3,dim4))
1905    array(:,:,:,:) = 0.0d0
1906
1907  end subroutine utl_reAllocate_r8_4d
1908
1909
1910  subroutine utl_reAllocate_r4_5d(array,dim1,dim2,dim3,dim4,dim5)
1911    implicit none
1912
1913    ! Arguments:
1914    real(4), allocatable, intent(inout) :: array(:,:,:,:,:)
1915    integer,              intent(in)    :: dim1
1916    integer,              intent(in)    :: dim2
1917    integer,              intent(in)    :: dim3
1918    integer,              intent(in)    :: dim4
1919    integer,              intent(in)    :: dim5
1920
1921    if( allocated(array) ) then
1922      if ( size(array) == dim1*dim2*dim3*dim4*dim5 ) then
1923        return
1924      else
1925        deallocate(array)
1926      end if
1927    end if
1928
1929    allocate(array(dim1,dim2,dim3,dim4,dim5))
1930    array(:,:,:,:,:) = 0.0d0
1931
1932  end subroutine utl_reAllocate_r4_5d
1933
1934
1935  subroutine utl_reAllocate_r8_5d(array,dim1,dim2,dim3,dim4,dim5)
1936    implicit none
1937
1938    ! Arguments:
1939    real(8), allocatable, intent(inout) :: array(:,:,:,:,:)
1940    integer,              intent(in)    :: dim1
1941    integer,              intent(in)    :: dim2
1942    integer,              intent(in)    :: dim3
1943    integer,              intent(in)    :: dim4
1944    integer,              intent(in)    :: dim5
1945
1946    if( allocated(array) ) then
1947      if ( size(array) == dim1*dim2*dim3*dim4*dim5 ) then
1948        return
1949      else
1950        deallocate(array)
1951      end if
1952    end if
1953
1954    allocate(array(dim1,dim2,dim3,dim4,dim5))
1955    array(:,:,:,:,:) = 0.0d0
1956
1957  end subroutine utl_reAllocate_r8_5d
1958
1959
1960  subroutine utl_heapsort2d(array)
1961    !
1962    !:Purpose: Sort a real 2D array in ascending order according
1963    !          to the first column
1964    ! 
1965    implicit none
1966
1967    ! Arguments:
1968    real(4), intent(inout) :: array(:,:)
1969
1970    ! Locals:
1971    real(4) :: values(2) ! temporary value
1972    integer :: i,j,nsize
1973    integer :: ileft,iright
1974
1975    nsize = size(array,1)
1976    ileft=nsize/2+1
1977    iright=nsize
1978
1979    if (nsize == 1) return                  
1980
1981    do 
1982      if(ileft > 1)then
1983        ileft=ileft-1
1984        values(:) = array(ileft,:)
1985      else
1986        values(:) = array(iright,:)
1987        array(iright,:) = array(1,:)
1988        iright = iright-1
1989        if (iright == 1) then
1990          array(1,:) = values(:)
1991          return
1992        end if
1993      end if
1994      i = ileft
1995      j = 2*ileft
1996      do while (j <= iright) 
1997        if (j < iright) then
1998          if (array(j,1) < array(j+1,1)) j=j+1
1999        endif
2000        if (values(1) < array(j,1)) then
2001          array(i,:) = array(j,:)
2002          i = j
2003          j = j+j
2004        else
2005          j = iright+1
2006        end if
2007      end do
2008      array(i,:) = values(:)
2009    end do
2010
2011  end subroutine utl_heapsort2d
2012
2013
2014  subroutine utl_splitString(string,separator,stringArray)
2015    implicit none
2016
2017    ! Arguments:
2018    character(len=*),              intent(in) :: string
2019    character(len=*),              intent(in) :: separator
2020    character(len=*), allocatable, intent(inout) :: stringArray(:)
2021
2022    ! Locals:
2023    integer :: stringArraySize, stringIndex
2024
2025    stringArraySize = count(transfer(string, 'a', len(string)) == separator) + 1
2026
2027    allocate(stringArray(stringArraySize))
2028
2029    read(string, *) stringArray(1:stringArraySize)
2030
2031    write(*,*) 'utl_splitString: stringArraySize = ', stringArraySize
2032    write(*,*) 'utl_splitString: stringArray     = ', &
2033               (trim(stringArray(stringIndex))//' ', stringIndex=1,stringArraySize)
2034    
2035  end subroutine utl_splitString
2036
2037
2038  subroutine utl_combineString(string,separator,stringArray)
2039    implicit none
2040
2041    ! Arguments:
2042    character(len=*), intent(out) :: string
2043    character(len=*), intent(in)  :: separator
2044    character(len=*), intent(in)  :: stringArray(:)
2045
2046    ! Locals:
2047    integer :: stringArraySize, stringIndex, stringCount, stringCountTotal
2048
2049    stringArraySize = size(stringArray)
2050
2051    stringCountTotal = 0
2052    do stringIndex = 1, size(stringArray)
2053      if (trim(stringArray(stringIndex)) == '') cycle
2054      stringCountTotal = stringCountTotal + 1
2055    end do
2056
2057    string = ''
2058    stringCount = 0
2059    do stringIndex = 1, size(stringArray)
2060      if (trim(stringArray(stringIndex)) == '') cycle
2061      stringCount = stringCount + 1
2062      if (stringCount < stringCountTotal) then
2063        string = trim(string) // ' ' // trim(stringArray(stringIndex)) // trim(separator)
2064      else
2065        string = trim(string) // ' ' // trim(stringArray(stringIndex))
2066      end if
2067    end do
2068
2069    write(*,*) 'utl_combineString: stringCountTotal = ', stringCountTotal
2070    write(*,*) 'utl_combineString: string     = ', trim(string)
2071    
2072  end subroutine utl_combineString
2073
2074
2075  subroutine utl_removeEmptyStrings(stringArray)
2076    implicit none
2077
2078    ! Arguments:
2079    character(len=*), allocatable, intent(inout) :: stringArray(:)
2080
2081    ! Locals:
2082    integer :: stringArraySize, stringIndex, stringCount, stringCountTotal
2083    character(len=len(stringArray(1))), allocatable :: newStringArray(:)
2084
2085    stringArraySize = size(stringArray)
2086
2087    stringCountTotal = 0
2088    do stringIndex = 1, size(stringArray)
2089      if (trim(stringArray(stringIndex)) == '') cycle
2090      stringCountTotal = stringCountTotal + 1
2091    end do
2092
2093    allocate(newStringArray(stringCountTotal))
2094
2095    stringCount = 0
2096    do stringIndex = 1, size(stringArray)
2097      if (trim(stringArray(stringIndex)) == '') cycle
2098      stringCount = stringCount + 1
2099      newStringArray(stringCount) = stringArray(stringIndex)
2100    end do
2101
2102    deallocate(stringArray)
2103    allocate(stringArray(stringCountTotal))
2104    stringArray(:) = newStringArray(:)
2105    deallocate(newStringArray)
2106
2107    write(*,*) 'utl_removeEmptyStrings: stringCountTotal = ', stringCountTotal
2108    write(*,*) 'utl_removeEmptyStrings: stringArray      = ', &
2109               (trim(stringArray(stringIndex))//' ', stringIndex = 1, size(stringArray))
2110    
2111  end subroutine utl_removeEmptyStrings
2112
2113
2114  subroutine utl_stringArrayToIntegerArray(stringArray,integerArray)
2115    implicit none
2116
2117    ! Arguments:
2118    character(len=256),   intent(in)  :: stringArray(:)
2119    integer, allocatable, intent(out) :: integerArray(:)
2120
2121    ! Locals:
2122    integer :: arraySize, arrayIndex
2123
2124    arraySize = size(stringArray)
2125
2126    allocate(integerArray(arraySize))
2127
2128    do arrayIndex = 1, arraySize
2129      read(stringArray(arrayIndex),'(i5)')  integerArray(arrayIndex)
2130    end do
2131
2132    write(*,*) 'utl_stringArrayToIntegerArray: integerArray = ', integerArray(:)
2133
2134  end subroutine utl_stringArrayToIntegerArray
2135
2136  !--------------------------------------------------------------------------
2137  ! utl_isNamelistPresent
2138  !--------------------------------------------------------------------------
2139  function utl_isNamelistPresent(namelistSectionName, namelistFileName) result(found)
2140    !
2141    !:Purpose: To find if a namelist name tag is present in a namelist file
2142    ! 
2143    implicit none
2144
2145    ! Arguments:
2146    character(len=*), intent(in) :: namelistSectionName
2147    character(len=*), intent(in) :: namelistFileName
2148    ! Result:
2149    logical :: found
2150
2151    ! Locals:
2152    integer :: unit, fnom, fclos, ierr
2153    character (len=1000) :: text
2154    character (len=100)  :: word, namelistSectionNameUpper
2155    logical :: namelistExist
2156
2157    ! Check if namelistFileName is present
2158    inquire(file=namelistFileName,exist=namelistExist)
2159    if (.not. namelistExist) then
2160      call utl_abort('utl_isNamelistPresent: namelist file is missing : '// namelistFileName)
2161    end if
2162
2163    ! Open the namelist file
2164    unit=0
2165    ierr=fnom(unit,namelistFileName,'FTN+SEQ+R/O',0)
2166
2167    ! Search for namelistSectionName
2168    found = .false.
2169    namelistSectionNameUpper = namelistSectionName
2170    ierr = clib_toUpper(namelistSectionNameUpper)
2171    namelistLoop : do
2172      read (unit,"(a)",iostat=ierr) text ! read line into character variable
2173      if (ierr /= 0) exit
2174      if (trim(text) == "") cycle ! skip empty lines
2175      read (text,*) word ! read first word of line
2176      ierr = clib_toUpper(word)
2177      if (trim(word) == '&'//trim(namelistSectionNameUpper)) then ! case insensitive 
2178        ! found search string at beginning of line
2179        found = .true.
2180        exit
2181      end if
2182    end do namelistLoop
2183
2184    ! Close the namelist file
2185    ierr=fclos(unit)
2186
2187  end function utl_isNamelistPresent
2188
2189  !-----------------------------------------------------------------
2190  ! utl_parseColumns
2191  !-----------------------------------------------------------------
2192  subroutine utl_parseColumns(line, numColumns, stringArray_opt)
2193    !
2194    !:Purpose: To return column values in array of strings and
2195    !          the number of space-delimited columns in a string
2196    ! 
2197    implicit none
2198
2199    ! Arguments:
2200    character(len=*),           intent(in)  :: line
2201    integer,                    intent(out) :: numColumns
2202    character(len=*), optional, intent(out) :: stringArray_opt(:)
2203
2204    ! Locals:
2205    integer :: linePosition, wordPosition, lineLength
2206
2207    linePosition = 1
2208    lineLength = len_trim(line)
2209    numColumns = 0
2210    
2211    do while(linePosition <= lineLength)
2212
2213      do while(line(linePosition:linePosition) == ' ') 
2214        linePosition = linePosition + 1
2215        if (lineLength < linePosition) return
2216      end do
2217
2218      numColumns = numColumns + 1
2219      wordPosition = 0
2220      if (present(stringArray_opt)) then
2221        stringArray_opt(numColumns) = ''
2222      end if
2223
2224      do
2225        if (linePosition > lineLength) return
2226        if (line(linePosition:linePosition) == ' ') exit
2227        if (present(stringArray_opt)) then
2228          wordPosition = wordPosition + 1
2229          stringArray_opt(numColumns)(wordPosition:wordPosition) = line(linePosition:linePosition)
2230        end if
2231        linePosition = linePosition + 1
2232      end do
2233
2234    end do
2235    
2236  end subroutine utl_parseColumns
2237
2238  !--------------------------------------------------------------------------
2239  ! utl_copyFile
2240  !--------------------------------------------------------------------------
2241  function utl_copyFile(filein, fileout) result(status)
2242    !
2243    !:Purpose: Copy the specified file to the new location and/or name
2244    !          This function is very general, but was initially written to
2245    !          copy files from the disk to the ram disk
2246    !
2247    !
2248    implicit none
2249
2250    ! Arguments:
2251    character(len=*), intent(in) :: filein
2252    character(len=*), intent(in) :: fileout
2253    ! Result:
2254    integer :: status
2255
2256    ! Locals:
2257    integer :: ierr, unitin, unitout
2258    integer(8) :: numChar
2259    character :: bufferB
2260    integer, parameter :: bufferSizeKB = 1024
2261    character :: bufferKB(bufferSizeKB)
2262    integer, parameter :: bufferSizeMB = 1024*1024
2263    character :: bufferMB(bufferSizeMB)
2264
2265    write(*,*) 'utl_copyFile: copy from ', trim(filein), ' to ', trim(fileout)
2266
2267    call utl_tmg_start(175,'low-level--utl_copyFile')
2268
2269    unitin=10
2270    open(unit=unitin, file=trim(filein), status='OLD', form='UNFORMATTED', &
2271         action='READ', access='STREAM')
2272
2273    unitout=11
2274    open(unit=unitout, file=trim(fileout), status='REPLACE', form='UNFORMATTED', &
2275         action='WRITE', access='STREAM')
2276
2277    numChar = 0
2278    do 
2279      read(unitin,iostat=ierr) bufferMB
2280      if (ierr < 0) exit
2281      numChar = numChar + bufferSizeMB
2282      write(unitout) bufferMB
2283    end do
2284
2285    do 
2286      read(unitin,iostat=ierr,pos=numChar+1) bufferKB
2287      if (ierr < 0) exit
2288      numChar = numChar + bufferSizeKB
2289      write(unitout) bufferKB
2290    end do
2291
2292    do 
2293      read(unitin,iostat=ierr,pos=numChar+1) bufferB
2294      if (ierr < 0) exit
2295      numChar = numChar + 1
2296      write(unitout) bufferB
2297    end do
2298
2299    write(*,*) 'utl_copyFile: copied ', numChar, ' bytes'
2300
2301    close(unit=unitin)
2302    close(unit=unitout)
2303
2304    if (numChar > 0) then
2305      status = 0
2306    else
2307      status = -1
2308      if (numChar == 0) then
2309        call utl_abort('utl_copyFile: ERROR, zero bytes copied')
2310      else
2311        ! Note: If 'numChar' becomes negative then it means it got bigger
2312        !       than the maximum integer the 'integer' type and so the
2313        !       variable 'numChar' wraps around and becomes negative.
2314        call utl_abort('utl_copyFile: ERROR, overflow detected since number of bytes copied is negative!')
2315      end if
2316    end if
2317
2318    call utl_tmg_stop(175)
2319
2320  end function utl_copyFile
2321
2322  !--------------------------------------------------------------------------
2323  ! utl_allReduce
2324  !--------------------------------------------------------------------------
2325  subroutine utl_allReduce(localGlobalValue)
2326    !
2327    !:Purpose: Perform mpi_allReduce to sum integer values over all
2328    !          mpi tasks and copy result back to same variable.
2329    !
2330    implicit none
2331
2332    ! Arguments:
2333    integer, intent(inout) :: localGlobalValue
2334
2335    ! Locals:
2336    integer :: localValue, globalValue, ierr
2337
2338    localValue = localGlobalValue
2339    call rpn_comm_allReduce(localValue, globalValue, 1, 'mpi_integer', &
2340                            'mpi_sum', 'grid', ierr)
2341    localGlobalValue = globalValue
2342    
2343  end subroutine utl_allReduce
2344
2345  !--------------------------------------------------------------------------
2346  ! utl_findloc_char
2347  !--------------------------------------------------------------------------
2348  function utl_findloc_char(charArray, value) result(location)
2349    !
2350    !:Purpose: A modified version of the fortran function `findloc`.
2351    !          If multiple matches are found in the array, a warning
2352    !          message is printed to the listing.
2353    !
2354    implicit none
2355
2356    ! Arguments:
2357    character(len=*), intent(in) :: charArray(:)
2358    character(len=*), intent(in) :: value
2359    ! Result:
2360    integer                      :: location
2361
2362    ! Locals:
2363    integer :: numFound, arrayIndex
2364
2365    numFound = 0
2366    LOOP: do arrayIndex = 1, size(charArray)
2367      if (trim(charArray(arrayIndex)) == trim(value)) then
2368        numFound = numFound + 1
2369        ! return the first location found
2370        if (numFound == 1) location = arrayIndex
2371      end if
2372    end do LOOP
2373
2374    ! give warning if more than 1 found
2375    if (numFound > 1) then
2376      write(*,*) 'utl_findloc_char: found multiple locations of ', trim(value)
2377      write(*,*) 'utl_findloc_char: number locations found =  ', numFound    
2378    end if
2379
2380    ! return zero if not found
2381    if (numFound == 0) then
2382      location = 0
2383    end if
2384
2385  end function utl_findloc_char
2386
2387  !--------------------------------------------------------------------------
2388  ! utl_findloc_int
2389  !--------------------------------------------------------------------------
2390  function utl_findloc_int(intArray, value) result(location)
2391    !
2392    !:Purpose: A modified version of the fortran function `findloc`.
2393    !          If multiple matches are found in the array, a warning
2394    !          message is printed to the listing.
2395    !
2396    implicit none
2397
2398    ! Arguments:
2399    integer, intent(in) :: intArray(:)
2400    integer, intent(in) :: value
2401    ! Result:
2402    integer             :: location
2403
2404    ! Locals:
2405    integer :: numFound, arrayIndex
2406
2407    numFound = 0
2408    LOOP: do arrayIndex = 1, size(intArray)
2409      if (intArray(arrayIndex) == value) then
2410        numFound = numFound + 1
2411        ! return the first location found
2412        if (numFound == 1) location = arrayIndex
2413      end if
2414    end do LOOP
2415
2416    ! give warning if more than 1 found
2417    if (numFound > 1) then
2418      write(*,*) 'utl_findloc_int: found multiple locations of ', value
2419      write(*,*) 'utl_findloc_int: number locations found =  ', numFound    
2420    end if
2421
2422    ! return zero if not found
2423    if (numFound == 0) then
2424      location = 0
2425    end if
2426
2427  end function utl_findloc_int
2428
2429  !--------------------------------------------------------------------------
2430  ! utl_findlocs_char
2431  !--------------------------------------------------------------------------
2432  function utl_findlocs_char(charArray, value) result(locations)
2433    !
2434    !:Purpose: A modified version of the fortran function `findloc`.
2435    !          Returns an array of all matches found in the array.
2436    !
2437    implicit none
2438
2439    ! Arguments:
2440    character(len=*), intent(in) :: charArray(:)
2441    character(len=*), intent(in) :: value
2442    ! Result:
2443    integer, allocatable         :: locations(:)
2444
2445    ! Locals:
2446    integer :: numFound, arrayIndex
2447
2448    if (allocated(locations)) deallocate(locations)
2449
2450    ! count number of matches found
2451    numFound = 0
2452    do arrayIndex = 1, size(charArray)
2453      if (trim(charArray(arrayIndex)) == trim(value)) numFound = numFound + 1
2454    end do
2455
2456    if (numFound > 0) then
2457
2458      ! return all found locations
2459      allocate(locations(numFound))
2460      numFound = 0
2461      do arrayIndex = 1, size(charArray)
2462        if (trim(charArray(arrayIndex)) == trim(value)) then
2463          numFound = numFound + 1
2464          locations(numFound) = arrayIndex
2465        end if
2466      end do
2467
2468    else
2469
2470      ! return zero if not found
2471      allocate(locations(1))
2472      locations(1) = 0
2473
2474    end if
2475
2476  end function utl_findlocs_char
2477
2478  !--------------------------------------------------------------------------
2479  ! utl_randomOrderInt
2480  !--------------------------------------------------------------------------
2481  subroutine utl_randomOrderInt(intArray,randomSeed)
2482    !
2483    !:Purpose: Randomly shuffle the order of the integer array elements.
2484    !
2485    implicit none
2486
2487    ! Arguments:
2488    integer, intent(inout) :: intArray(:)
2489    integer, intent(in)    :: randomSeed
2490
2491    ! Locals:
2492    integer              :: arraySize, arrayIndex, arrayIndexMin
2493    integer, allocatable :: intArrayOut(:)
2494    real(8), allocatable :: realRandomArray(:)
2495
2496    arraySize = size(intArray)
2497    allocate(realRandomArray(arraySize))
2498    allocate(intArrayOut(arraySize))
2499
2500    call rng_setup(randomSeed)
2501    do arrayIndex = 1, arraySize
2502      realRandomArray(arrayIndex) = rng_uniform()
2503    end do
2504
2505    do arrayIndex = 1, arraySize
2506      arrayIndexMin = minloc(realRandomArray,dim=1)
2507      realRandomArray(arrayIndexMin) = huge(1.0D0)
2508      intArrayOut(arrayIndex) = intArray(arrayIndexMin)
2509    end do
2510
2511    intArray(:) = intArrayOut(:)
2512
2513    deallocate(intArrayOut)
2514    deallocate(realRandomArray)
2515
2516  end subroutine utl_randomOrderInt
2517
2518  !--------------------------------------------------------------------------
2519  ! utl_tmg_start
2520  !--------------------------------------------------------------------------
2521  subroutine utl_tmg_start(blockIndex, blockLabel)
2522    !
2523    !:Purpose: Wrapper for rpnlib subroutine tmg_start
2524    !
2525    implicit none
2526
2527    ! Arguments:
2528    integer,          intent(in) :: blockIndex
2529    character(len=*), intent(in) :: blockLabel
2530
2531    ! Locals:
2532    integer            :: labelLength, omp_get_thread_num
2533    integer, parameter :: labelPaddedLength = 40
2534    character(len=labelPaddedLength) :: blockLabelPadded
2535
2536    ! only the first thread does the timing
2537    if (omp_get_thread_num() > 0) return
2538
2539    blockLabelPadded = '........................................'
2540    labelLength = min(len_trim(blockLabel), labelPaddedLength)
2541    blockLabelPadded(1:labelLength) = blockLabel(1:labelLength)
2542
2543    call tmg_start(blockIndex, blockLabelPadded)
2544
2545  end subroutine utl_tmg_start
2546  
2547  !--------------------------------------------------------------------------
2548  ! utl_tmg_stop
2549  !--------------------------------------------------------------------------
2550  subroutine utl_tmg_stop(blockIndex)
2551    !
2552    !:Purpose: Wrapper for rpnlib subroutine tmg_stop
2553    !
2554    implicit none
2555
2556    ! Arguments:
2557    integer,          intent(in) :: blockIndex
2558
2559    ! Locals:
2560    integer            :: omp_get_thread_num
2561
2562    ! only the first thread does the timing
2563    if (omp_get_thread_num() > 0) return
2564
2565    call tmg_stop(blockIndex)
2566
2567  end subroutine utl_tmg_stop  
2568  
2569  !--------------------------------------------------------------------------
2570  ! utl_medianIndex
2571  !--------------------------------------------------------------------------
2572  function utl_medianIndex(inputVector) result(medianIndex)
2573    ! 
2574    !:Purpose: to find the median index of an input vector
2575    !
2576    implicit none
2577    
2578    ! Arguments:
2579    real(4), intent(in) :: inputVector(:)
2580    ! Result:
2581    integer             :: medianIndex
2582
2583    ! Locals:
2584    integer :: vectorIndex, vectorDim
2585    logical :: maskVector(size(inputVector))
2586    real(4) :: sortedArray(size(inputVector))
2587    real(4) :: median
2588
2589    vectorDim = size(inputVector)
2590
2591    ! sorting array:
2592    maskVector(:) = .true.
2593    do vectorIndex = 1, vectorDim 
2594      sortedArray(vectorIndex) = minval(inputVector, maskVector)
2595      maskVector(minloc(inputVector, maskVector)) = .false.
2596    end do
2597  
2598    if (mod(size(inputVector), 2) == 0) then
2599      median = sortedArray(vectorDim / 2)
2600    else
2601      median = sortedArray((vectorDim + 1) / 2)
2602    end if
2603
2604    do vectorIndex = 1, vectorDim
2605      if (inputVector(vectorIndex) == median) then
2606        medianIndex = vectorIndex
2607        exit
2608      end if
2609    end do
2610
2611  end function utl_medianIndex
2612
2613end module utilities_mod