utils.F90

Go to the documentation of this file.

module utilities

contains

 subroutine error_handler(string, rc)

 use mpi

 implicit none

 character(len=*), intent(in)    :: string
 
 integer,          intent(in)    :: rc

 integer :: ierr

 print*,"- FATAL ERROR: ", trim(string)
 print*,"- IOSTAT IS: ", rc
 call mpi_abort(mpi_comm_world, 999, ierr)

 end subroutine error_handler

 subroutine netcdf_err( err, string )

 use mpi
 use netcdf

 implicit none
 integer, intent(in) :: err
 character(len=*), intent(in) :: string
 character(len=256) :: errmsg
 integer :: iret

 if( err.EQ.nf90_noerr )return
 errmsg = nf90_strerror(err)
 print*,''
 print*,'FATAL ERROR: ', trim(string), ': ', trim(errmsg)
 print*,'STOP.'
 call mpi_abort(mpi_comm_world, 999, iret)

 return
 end subroutine netcdf_err
 
function to_upper(strIn) result(strOut)

     implicit none

     character(len=*), intent(in) :: strin
     character(len=len(strIn)) :: strout
     integer :: i,j

     do i = 1, len(strin)
          j = iachar(strin(i:i))
          if (j>= iachar("a") .and. j<=iachar("z") ) then
               strout(i:i) = achar(iachar(strin(i:i))-32)
          else
               strout(i:i) = strin(i:i)
          end if
     end do

end function to_upper

subroutine to_lower(strIn)

     implicit none

     character(len=*), intent(inout) :: strIn
     character(len=len(strIn)) :: strOut
     integer :: i,j

     do i = 1, len(strin)
          j = iachar(strin(i:i))
          if (j>= iachar("A") .and. j<=iachar("Z") ) then
               strout(i:i) = achar(iachar(strin(i:i))+32)
          else
               strout(i:i) = strin(i:i)
          end if
     end do
  strin(:) = strout(:)
end subroutine to_lower

subroutine handle_grib_error(vname,lev,method,value,varnum,read_from_input, iret,var,var8,var3d)

  use, intrinsic :: ieee_arithmetic
  use esmf

  implicit none
   
  real(esmf_kind_r4), intent(in)    :: value
  logical, intent(inout)            :: read_from_input(:)
  real(esmf_kind_r4), intent(inout), optional :: var(:,:)
  real(esmf_kind_r8), intent(inout), optional :: var8(:,:)
  real(esmf_kind_r8), intent(inout), optional  :: var3d(:,:,:)
  
  character(len=20), intent(in)     :: vname, lev, method
  character(len=200)                :: err_msg
  
  integer, intent(in)               :: varnum 
  integer, intent(inout)            :: iret
  
  iret = 0
  if (varnum == 9999) then
    print*, "WARNING: ", trim(vname), " NOT FOUND AT LEVEL ", lev, " IN EXTERNAL FILE ", &
            "AND NO ENTRY EXISTS IN VARMAP TABLE. VARIABLE WILL NOT BE USED."
    iret = 1

    return
  endif

  if (trim(method) == "skip" ) then
    print*, "WARNING: SKIPPING ", trim(vname), " IN FILE"
    read_from_input(varnum) = .false.
    iret = 1
  elseif (trim(method) == "set_to_fill") then
    print*, "WARNING: ,", trim(vname), " NOT AVAILABLE AT LEVEL ", trim(lev), &
           ". SETTING EQUAL TO FILL VALUE OF ", value
    if(present(var)) var(:,:) = value
    if(present(var8)) var8(:,:) = value
    if(present(var3d)) var3d(:,:,:) = value
  elseif (trim(method) == "set_to_NaN") then
    print*, "WARNING: ,", trim(vname), " NOT AVAILABLE AT LEVEL ", trim(lev), &
           ". SETTING EQUAL TO NaNs"
    if(present(var)) var(:,:) = ieee_value(var,ieee_quiet_nan)
    if(present(var8)) var8(:,:) = ieee_value(var8,ieee_quiet_nan)
    if(present(var3d)) var3d(:,:,:) = ieee_value(var3d,ieee_quiet_nan)
  elseif (trim(method) == "stop") then
    err_msg="READING " // trim(vname) // " at level " //lev// ". TO MAKE THIS NON-" // &
            "FATAL, CHANGE STOP TO SKIP FOR THIS VARIABLE IN YOUR VARMAP FILE."
    call error_handler(err_msg, iret)
  elseif (trim(method) == "intrp") then
    print*, "WARNING: ,"//trim(vname)//" NOT AVAILABLE AT LEVEL "//trim(lev)// &
          ". WILL INTERPOLATE INTERSPERSED MISSING LEVELS AND/OR FILL MISSING"//&
          " LEVELS AT EDGES."
  else
    err_msg="ERROR USING MISSING_VAR_METHOD. PLEASE SET VALUES IN" // &
            " VARMAP TABLE TO ONE OF: set_to_fill, set_to_NaN, intrp, skip, or stop."
    call error_handler(err_msg, 1)
  endif

end subroutine handle_grib_error

recursive subroutine quicksort(a, first, last)
  implicit none
  real*8  a(*), x, t
  integer first, last
  integer i, j

  x = a( (first+last) / 2 )
  i = first
  j = last
  do
     do while (a(i) < x)
        i=i+1
     end do
     do while (x < a(j))
        j=j-1
     end do
     if (i >= j) exit
     t = a(i);  a(i) = a(j);  a(j) = t
     i=i+1
     j=j-1
  end do
  if (first < i-1) call quicksort(a, first, i-1)
  if (j+1 < last)  call quicksort(a, j+1, last)
end subroutine quicksort


subroutine check_soilt(soilt, landmask, skint,ICET_DEFAULT,i_input,j_input,lsoil_input)
  use esmf
  implicit none
  integer, intent(in)               :: i_input, j_input, lsoil_input
  real(esmf_kind_r8), intent(inout) ::  soilt(i_input,j_input,lsoil_input)
  real(esmf_kind_r8), intent(in)    ::  skint(i_input,j_input)
  real,  intent(in) :: ICET_DEFAULT
  integer(esmf_kind_i4), intent(in)    ::  landmask(i_input,j_input)
  
  integer                           :: i, j, k

  do k=1,lsoil_input
    do j = 1, j_input
      do i = 1, i_input
        if (landmask(i,j) == 0_esmf_kind_i4 ) then 
          soilt(i,j,k) = skint(i,j)
        else if (landmask(i,j) == 1_esmf_kind_i4 .and. soilt(i,j,k) > 350.0_esmf_kind_r8) then 
          soilt(i,j,k) = skint(i,j)
        else if (landmask(i,j) == 2_esmf_kind_i4 ) then 
          soilt(i,j,k) = icet_default
        endif
      enddo
    enddo
  enddo
end subroutine check_soilt


subroutine check_cnwat(cnwat,i_input,j_input)
  use esmf
  implicit none 
  integer, intent(in)               :: i_input, j_input
  real(esmf_kind_r8), intent(inout) :: cnwat(i_input,j_input)
  
  real(esmf_kind_r8)                :: max_cnwat = 0.5
  
  integer :: i, j

  do i = 1,i_input
    do j = 1,j_input
      if (cnwat(i,j) .gt. max_cnwat) cnwat(i,j) = 0.0_esmf_kind_r8
    enddo
  enddo
end subroutine check_cnwat

SUBROUTINE dint2p(PPIN,XXIN,NPIN,PPOUT,XXOUT,NPOUT   &                   
   ,LINLOG,XMSG,IER)
      IMPLICIT NONE

! NCL code for pressure level interpolation
!
! This code was designed for one simple task. It has since
! been mangled and abused for assorted reasons. For example,
! early gfortran compilers had some issues with automatic arrays.
! Hence, the C-Wrapper was used to create 'work' arrays which
! were then passed to this code.  The original focused (non-NCL) 
! task was to handle PPIN & PPOUT that had the same 'monotonicity.' 
! Extra code was added to handle the more general case. 
! Blah-Blah:  Punch line: it is embarrassingly convoluted!!!
!
!                                                ! input types
      INTEGER NPIN,NPOUT,LINLOG,IER
      real*8 PPIN(NPIN),XXIN(NPIN),PPOUT(NPOUT),XMSG
                                                ! output
      real*8 XXOUT(NPOUT)
                                                ! work
      real*8 PIN(NPIN),XIN(NPIN),P(NPIN),X(NPIN)
      real*8 POUT(NPOUT),XOUT(NPOUT)

! local
      INTEGER NP,NL,NLMAX,NLSAVE,NP1,NO1,N1,N2,LOGLIN,   &
             NLSTRT
      real*8 SLOPE,PA,PB,PC

      loglin = abs(linlog)

! error check: enough points: pressures consistency?

      ier = 0
      IF (npout.GT.0) THEN
          DO np = 1,npout
              xxout(np) = xmsg
          END DO
      END IF
! Jili Dong input levels have to be the same as output levels:
! we only interpolate for levels with missing variables
!      IF (.not. all(PPIN .eq. PPOUT)) IER = IER+1

      IF (npin.LT.2 .OR. npout.LT.1) ier = ier + 1

      IF (ier.NE.0) THEN
!          PRINT *,'INT2P: error exit: ier=',IER
          RETURN
      END IF

! should *input arrays* be reordered? want p(1) > p(2) > p(3) etc
! so that it will match order for which code was originally designed
! copy to 'work'  arrays

      np1 = 0
      no1 = 0
      IF (ppin(1).LT.ppin(2)) THEN
          np1 = npin + 1
      END IF
      IF (ppout(1).LT.ppout(2)) THEN
          no1 = npout + 1
      END IF

      DO np = 1,npin
          pin(np) = ppin(abs(np1-np))
          xin(np) = xxin(abs(np1-np))
      END DO

      DO np = 1,npout
          pout(np) = ppout(abs(no1-np))
      END DO

! eliminate XIN levels with missing data. 
! .   This can happen with observational data.

      nl = 0
      DO np = 1,npin
          IF (xin(np).NE.xmsg .AND. pin(np).NE.xmsg) THEN
              nl = nl + 1
              p(nl) = pin(np)
              x(nl) = xin(np)
          END IF
      END DO
      nlmax = nl

                                                ! all missing data
      IF (nlmax.LT.2) THEN
          ier = ier + 1000
          print *,'INT2P: ier=',ier
          RETURN
      END IF

! ===============> pressure in decreasing order <================
! perform the interpolation  [pin(1)>pin(2)>...>pin(npin)]
!                                                      ( p ,x)
! ------------------------- p(nl+1), x(nl+1)   example (200,5)
! .
! ------------------------- pout(np), xout(np)         (250,?)
! .
! ------------------------- p(nl)  , x(nl)             (300,10)


! exact p-level matches
      nlstrt = 1
      nlsave = 1
      DO np = 1,npout
          xout(np) = xmsg
          DO nl = nlstrt,nlmax
              IF (pout(np).EQ.p(nl)) THEN
                  xout(np) = x(nl)
                  nlsave = nl + 1
                  GO TO 10
              END IF
          END DO
   10     nlstrt = nlsave
      END DO

      IF (loglin.EQ.1) THEN
          DO np = 1,npout
              DO nl = 1,nlmax - 1
                  IF (pout(np).LT.p(nl) .AND. pout(np).GT.p(nl+1)) THEN
                      slope = (x(nl)-x(nl+1))/ (p(nl)-p(nl+1))
                      xout(np) = x(nl+1) + slope* (pout(np)-p(nl+1))
                  END IF
              END DO
          END DO
      ELSE
          DO np = 1,npout
              DO nl = 1,nlmax - 1
                  IF (pout(np).LT.p(nl) .AND. pout(np).GT.p(nl+1)) THEN
                      pa = log(p(nl))
                      pb = log(pout(np))
! special case: In case someome inadvertently enter p=0.
                      if (p(nl+1).gt.0.d0) then
                          pc = log(p(nl+1))
                      else
                          pc = log(1.e-4)
                      end if

                      slope = (x(nl)-x(nl+1))/ (pa-pc)
                      xout(np) = x(nl+1) + slope* (pb-pc)
                  END IF
              END DO
          END DO
      END IF

! extrapolate?
! . use the 'last' valid slope for extrapolating

      IF (linlog.LT.0) THEN
          DO np = 1,npout
              DO nl = 1,nlmax
                  IF (pout(np).GT.p(1)) THEN
                      IF (loglin.EQ.1) THEN
                          slope = (x(2)-x(1))/ (p(2)-p(1))
                          xout(np) = x(1) + slope* (pout(np)-p(1))
                      ELSE
                          pa = log(p(2))
                          pb = log(pout(np))
                          pc = log(p(1))
                          slope = (x(2)-x(1))/ (pa-pc)
                          xout(np) = x(1) + slope* (pb-pc)
                      END IF
                  ELSE IF (pout(np).LT.p(nlmax)) THEN
                      n1 = nlmax
                      n2 = nlmax - 1
                      IF (loglin.EQ.1) THEN
                          slope = (x(n1)-x(n2))/ (p(n1)-p(n2))
                          xout(np) = x(n1) + slope* (pout(np)-p(n1))
                      ELSE
                          pa = log(p(n1))
                          pb = log(pout(np))
                          pc = log(p(n2))
                          slope = (x(n1)-x(n2))/ (pa-pc)
                          !XOUT(NP) = X(N1) + SLOPE* (PB-PC) !bug fixed below
                          xout(np) = x(n1) + slope* (pb-pa)
                      END IF
                  END IF
              END DO
          END DO
      END IF

! place results in the return array;
! .   possibly .... reverse to original order

      if (no1.GT.0) THEN
          DO np = 1,npout
             n1 = abs(no1-np)
             ppout(np) = pout(n1)
             xxout(np) = xout(n1)
          END DO
      ELSE
          DO np = 1,npout
             ppout(np) = pout(np)
             xxout(np) = xout(np)
          END DO
      END IF


      RETURN
      END SUBROUTINE dint2p
end module utilities