text:
subroutine wrcomt(comfil    ,lundia    ,error     ,itcur     ,ntcur     , &
                & itimc     ,mmax      ,nmax      ,kmax      ,nmaxus    , &
                & nsrc      ,mnksrc    ,lstsci    ,lsal      ,ltem      , &
                & lsecfl    ,kfu       ,kfv       ,ibuff     ,s1        , &
                & u1        ,v1        ,qu        ,qv        ,taubmx    , &
                & r1        ,dicuv     ,dicww     ,discum    ,rbuff     , &
                & windu     ,windv     ,dzu1      ,dzv1      ,kmaxz     , &
                & hu        ,hv        ,thick     ,gdp       )
!----- GPL ---------------------------------------------------------------------
!
!  Copyright (C)  Stichting Deltares, 2011-2016.
!
!  This program is free software: you can redistribute it and/or modify
!  it under the terms of the GNU General Public License as published by
!  the Free Software Foundation version 3.
!
!  This program is distributed in the hope that it will be useful,
!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!  GNU General Public License for more details.
!
!  You should have received a copy of the GNU General Public License
!  along with this program.  If not, see <http://www.gnu.org/licenses/>.
!
!  contact: delft3d.support@deltares.nl
!  Stichting Deltares
!  P.O. Box 177
!  2600 MH Delft, The Netherlands
!
!  All indications and logos of, and references to, "Delft3D" and "Deltares"
!  are registered trademarks of Stichting Deltares, and remain the property of
!  Stichting Deltares. All rights reserved.
!
!-------------------------------------------------------------------------------
!  $Id: wrcomt.f90 5717 2016-01-12 11:35:24Z mourits $
!  $HeadURL: https://svn.oss.deltares.nl/repos/delft3d/tags/6686/src/engines_gpl/flow2d3d/packages/io/src/output/wrcomt.f90 $
!!--description-----------------------------------------------------------------
! NONE
!!--pseudo code and references--------------------------------------------------
! NONE
!!--declarations----------------------------------------------------------------
    use precision
    !
    use globaldata
    !
    implicit none
    !
    type(globdat),target :: gdp
    !
    ! The following list of pointer parameters is used to point inside the gdp structure
    !
    logical , pointer :: wind
!
! Global variables
!
    integer                                                                     :: itcur  !!  Current time counter for the com-
                                                                                !!  munication file, where starting
                                                                                !!  point depend on CYCLIC
    integer                                                                     :: itimc  !!  Current time step counter for 2D
                                                                                !!  system
    integer                                                                     :: kmax   !  Description and declaration in esm_alloc_int.f90
    integer                                                                     :: kmaxz  !!  = KMAX for Z-model, = 0 for sigma-model
                                                                                !!  Needed for correct dimensioning of DZU1 and DZV1
    integer                                                                     :: lsal   !  Description and declaration in dimens.igs
    integer                                                                     :: lsecfl !  Description and declaration in dimens.igs
    integer                                                                     :: lstsci !  Description and declaration in esm_alloc_int.f90
    integer                                                                     :: ltem   !  Description and declaration in dimens.igs
    integer                                                                     :: lundia !  Description and declaration in inout.igs
    integer                                                                     :: mmax   !  Description and declaration in esm_alloc_int.f90
    integer                                                                     :: nmax   !  Description and declaration in esm_alloc_int.f90
    integer                                                                     :: nmaxus !  Description and declaration in esm_alloc_int.f90
    integer                                                                     :: nsrc   !  Description and declaration in esm_alloc_int.f90
    integer                                                                     :: ntcur  !!  Total number of timesteps on com-
                                                                                !!  munication file (to write to)
    integer , dimension(7, nsrc)                                                :: mnksrc !  Description and declaration in esm_alloc_int.f90
    integer , dimension(gdp%d%nlb:gdp%d%nub, gdp%d%mlb:gdp%d%mub)               :: kfu    !  Description and declaration in esm_alloc_int.f90
    integer , dimension(gdp%d%nlb:gdp%d%nub, gdp%d%mlb:gdp%d%mub)               :: kfv    !  Description and declaration in esm_alloc_int.f90
    integer , dimension(nmaxus, mmax)                                           :: ibuff  !  Description and declaration in esm_alloc_int.f90
    real(fp), dimension(gdp%d%nlb:gdp%d%nub, gdp%d%mlb:gdp%d%mub)               :: s1     !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nlb:gdp%d%nub, gdp%d%mlb:gdp%d%mub)               :: taubmx !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nlb:gdp%d%nub, gdp%d%mlb:gdp%d%mub, 0:kmax)       :: dicww  !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nlb:gdp%d%nub, gdp%d%mlb:gdp%d%mub, kmax + 2)     :: dicuv  !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nlb:gdp%d%nub, gdp%d%mlb:gdp%d%mub, kmaxz)        :: dzu1   !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nlb:gdp%d%nub, gdp%d%mlb:gdp%d%mub, kmaxz)        :: dzv1   !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nlb:gdp%d%nub, gdp%d%mlb:gdp%d%mub)               :: hu     !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nlb:gdp%d%nub, gdp%d%mlb:gdp%d%mub)               :: hv     !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nlb:gdp%d%nub, gdp%d%mlb:gdp%d%mub, kmax)         :: qu     !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nlb:gdp%d%nub, gdp%d%mlb:gdp%d%mub, kmax)         :: qv     !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nlb:gdp%d%nub, gdp%d%mlb:gdp%d%mub, kmax)         :: u1     !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nlb:gdp%d%nub, gdp%d%mlb:gdp%d%mub, kmax)         :: v1     !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nlb:gdp%d%nub, gdp%d%mlb:gdp%d%mub)               :: windu  !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nlb:gdp%d%nub, gdp%d%mlb:gdp%d%mub)               :: windv  !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nlb:gdp%d%nub, gdp%d%mlb:gdp%d%mub, kmax, lstsci) :: r1     !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(nmaxus, mmax, kmax)                                     :: rbuff  !  Description and declaration in r-i-ch.igs
    real(fp), dimension(nsrc)                                                   :: discum !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(kmax)                                                   :: thick  !  Description and declaration in esm_alloc_real.f90
    logical                                                                     :: error  !!  Flag=TRUE if an error is encountered
    character(*)                                                                :: comfil !!  Name for communication file
                                                                                !!  com-<case><label>
!
! Local variables
!
!
!! executable statements -------------------------------------------------------
!
    wind       => gdp%gdprocs%wind
    !
    ! Write groups KENMNT and KENMTIM
    !
    call wrkent(comfil    ,lundia    ,error     ,itcur     ,ntcur     , &
              & itimc     ,mmax      ,nmax      ,nmaxus    ,kfu       , &
              & kfv       ,ibuff     ,gdp       )
    if (error) goto 9999
    !
    ! Write groups CURNT and CURTIM
    !
    call wrcurt(comfil    ,lundia    ,error     ,itcur     ,ntcur     , &
              & itimc     ,mmax      ,nmax      ,kmax      ,nmaxus    , &
              & lstsci    ,lsecfl    ,s1        ,u1        ,v1        , &
              & r1        ,qu        ,qv        ,dzu1      ,dzv1      , &
              & rbuff     ,kmaxz     ,hu        ,hv        ,thick     , &
              & gdp       )
    !
    if (error) goto 9999
    !
    ! Write group DWQTIM
    !
    call wrdwqt(comfil    ,lundia    ,error     ,itcur     ,itimc     , &
              & nsrc      ,mnksrc    ,mmax      ,nmax      ,kmax      , &
              & nmaxus    ,lstsci    ,lsal      ,ltem      ,r1        , &
              & dicuv     ,dicww     ,discum    ,taubmx    ,rbuff     , &
              & gdp       )
    if (error) goto 9999
    !
    !
    ! Write group WIND
    !
    if (wind) then
       call wrcomwind(error     ,comfil    ,itcur     ,itimc     , &
                    & mmax      ,nmaxus    , &
                    & windu     ,windv     ,gdp       )
       if (error) goto 9999
    endif
 9999 continue
end subroutine wrcomt

max_stars_repo_name: liujiamingustc/phd

max_issues_repo_name: liujiamingustc/phd

max_forks_repo_name: liujiamingustc/phd


text:
SUBROUTINE idot(g)

!-------------------------------------------------------------
!
!  Integral discrete ordinates transport
!
!  Add sv to kmat*psii to get the new phi
!  Compute psio with jpsi*(phi+src) + kpsi*psii
!
!-------------------------------------------------------------

USE invar
USE solvar
IMPLICIT NONE
INTEGER, INTENT(IN) :: g
INTEGER :: i, j, ky, ly, kx, lx, ix1, ix2, ix3, ix4
REAL*8, DIMENSION(ordsq*nx*ny) :: q

! Solve for the new phi with the sum of sv with the product of kmat and psi
phi(:,g) = sv + MATMUL(kmat(:,:,1),psi1) + MATMUL(kmat(:,:,2),psi2) &
              + MATMUL(kmat(:,:,3),psi3) + MATMUL(kmat(:,:,4),psi4)
cnvf(g) = 1

! Compute the new angular flux values at the boundaries from jpsi and kpsi
q = phi(:,g) + src

! Set up the loop independent indices
ix1 = 1
ix2 = nx*order
ix3 = ix2 + 1
ix4 = ix2 + ny*order
DO j = 1, 4
   psio(:,j,g) = MATMUL(jpsi(:,:,j),q)
   DO i = 1, apo
      ky = (i-1)*nx*order + 1
      ly = i*nx*order
      kx = apo*nx*order + (i-1)*ny*order + 1
      lx = apo*nx*order + i*ny*order
      IF (j == 1) THEN
         psio(ky:ly,j,g) = psio(ky:ly,j,g) + MATMUL(kpsi(ix1:ix2,ix1:ix2,i,j),psi1(ky:ly)) &
                                           + MATMUL(kpsi(ix1:ix2,ix3:ix4,i,j),psi1(kx:lx))
         psio(kx:lx,j,g) = psio(kx:lx,j,g) + MATMUL(kpsi(ix3:ix4,ix1:ix2,i,j),psi1(ky:ly)) &
                                           + MATMUL(kpsi(ix3:ix4,ix3:ix4,i,j),psi1(kx:lx))
      ELSE IF (j == 2) THEN
         psio(ky:ly,j,g) = psio(ky:ly,j,g) + MATMUL(kpsi(ix1:ix2,ix1:ix2,i,j),psi2(ky:ly)) &
                                           + MATMUL(kpsi(ix1:ix2,ix3:ix4,i,j),psi2(kx:lx))
         psio(kx:lx,j,g) = psio(kx:lx,j,g) + MATMUL(kpsi(ix3:ix4,ix1:ix2,i,j),psi2(ky:ly)) &
                                           + MATMUL(kpsi(ix3:ix4,ix3:ix4,i,j),psi2(kx:lx))
      ELSE IF (j == 3) THEN
         psio(ky:ly,j,g) = psio(ky:ly,j,g) + MATMUL(kpsi(ix1:ix2,ix1:ix2,i,j),psi3(ky:ly)) &
                                           + MATMUL(kpsi(ix1:ix2,ix3:ix4,i,j),psi3(kx:lx))
         psio(kx:lx,j,g) = psio(kx:lx,j,g) + MATMUL(kpsi(ix3:ix4,ix1:ix2,i,j),psi3(ky:ly)) &
                                           + MATMUL(kpsi(ix3:ix4,ix3:ix4,i,j),psi3(kx:lx))
      ELSE
         psio(ky:ly,j,g) = psio(ky:ly,j,g) + MATMUL(kpsi(ix1:ix2,ix1:ix2,i,j),psi4(ky:ly)) &
                                           + MATMUL(kpsi(ix1:ix2,ix3:ix4,i,j),psi4(kx:lx))
         psio(kx:lx,j,g) = psio(kx:lx,j,g) + MATMUL(kpsi(ix3:ix4,ix1:ix2,i,j),psi4(ky:ly)) &
                                           + MATMUL(kpsi(ix3:ix4,ix3:ix4,i,j),psi4(kx:lx))
      END IF
   END DO
END DO

RETURN
END SUBROUTINE idot

max_stars_repo_name: lanl/pidots

max_issues_repo_name: lanl/pidots

max_forks_repo_name: lanl/pidots


text:
!!  Copyright (C)  Stichting Deltares, 2012-2016.
!!
!!  This program is free software: you can redistribute it and/or modify
!!  it under the terms of the GNU General Public License version 3,
!!  as published by the Free Software Foundation.
!!
!!  This program is distributed in the hope that it will be useful,
!!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!!  GNU General Public License for more details.
!!
!!  You should have received a copy of the GNU General Public License
!!  along with this program. If not, see <http://www.gnu.org/licenses/>.
!!
!!  contact: delft3d.support@deltares.nl
!!  Stichting Deltares
!!  P.O. Box 177
!!  2600 MH Delft, The Netherlands
!!
!!  All indications and logos of, and references to registered trademarks
!!  of Stichting Deltares remain the property of Stichting Deltares. All
!!  rights reserved.

      subroutine pdfch2 ( lurep , procesdef, procha)

      ! Fill in the names of the charon coupling, modify stochiometric factors

      use processet
      implicit none

      ! declaration of arguments

      integer               :: lurep               ! unit number report file
      type(procespropcoll)  :: procesdef           ! all processes
      type(ProcesProp)      :: procha              ! charon process definition

!
!     Local declarations
!
      integer                   :: nproc           ! number of processes
      integer                   :: iproc           ! loop counter processes
      type(procesprop), pointer :: proc            ! process description
      integer                   :: i_input         ! index input item
      integer                   :: indx            ! index in io list
      integer                   :: iret            ! return code
      integer                   :: i_output        ! index output item
      integer                   :: iflux           ! index flux
      integer                   :: istochi         ! index stochi
      integer                   :: naij2           ! naij2
      type(ioitempropcoll)      :: input_item      ! one collection of items
      type(ioitempropcoll)      :: output_item     ! one collection of items
      type(ioitempropcoll)      :: fluxoutput      ! one collection of items
      type(stochipropcoll)      :: fluxstochi      ! one collection of stochis

      ! Loop over the processes

      nproc = procesdef%cursize
      do iproc = 1, nproc
         proc => procesdef%procesprops(iproc)

         if ( proc%routine .eq. 'D40CHA' ) then

            ! add extra input

            input_item%cursize=proc%no_input
            input_item%maxsize=proc%no_input
            input_item%ioitemprops => proc%input_item
            do i_input = 1 , procha%no_input
               indx = proc%no_input + i_input
               procha%input_item(i_input)%indx = indx
               iret = ioitempropcolladdindx( input_item, procha%input_item(i_input), indx )
            enddo
            proc%no_input     = input_item%cursize
            proc%input_item   =>input_item%ioitemprops

            ! add extra output

            output_item%cursize=proc%no_output
            output_item%maxsize=proc%no_output
            output_item%ioitemprops => proc%output_item
            do i_output = 1 , procha%no_output
               indx = proc%no_output + i_output
               procha%output_item(i_output)%indx = indx
               iret = ioitempropcolladdindx( output_item, procha%output_item(i_output), indx )
            enddo
            proc%no_output     = output_item%cursize
            proc%output_item   =>output_item%ioitemprops

            ! add fluxes

            fluxoutput%cursize=proc%no_fluxoutput
            fluxoutput%maxsize=proc%no_fluxoutput
            fluxoutput%ioitemprops => proc%fluxoutput
            do iflux = 1 , procha%no_fluxoutput
               indx = proc%no_fluxoutput + iflux
               procha%fluxoutput(iflux)%indx = indx
               iret = ioitempropcolladdindx( fluxoutput, procha%fluxoutput(iflux), indx )
            enddo
            proc%no_fluxoutput     = fluxoutput%cursize
            proc%fluxoutput        =>fluxoutput%ioitemprops

            ! add stochi

            fluxstochi%cursize=proc%no_fluxstochi
            fluxstochi%maxsize=proc%no_fluxstochi
            fluxstochi%stochiprops => proc%fluxstochi
            do istochi = 1 , procha%no_fluxstochi
               iret = StochiPropCollAdd( fluxstochi, procha%fluxstochi(istochi) )
            enddo
            proc%no_fluxstochi     = fluxstochi%cursize
            proc%fluxstochi        =>fluxstochi%stochiprops

         endif

      enddo

      ! modify fluxes

      call chphas ( naij2 )
      call chsto2 ( lurep , naij2 , procesdef)

      return
 2000 format('ERROR: allocating workarray in PDFCH2',i8)
      end

ext: f


text:
!==========================================================================
elemental function gsw_pot_enthalpy_ice_freezing_poly (sa, p)
!==========================================================================
!
!  Calculates the potential enthalpy of ice at which seawater freezes.
!  The error of this fit ranges between -2.5 and 1 J/kg with an rms of
!  1.07, between SA of 0 and 120 g/kg and p between 0 and 10,000 dbar (the
!  error in the fit is between -0.7 and 0.7 with an rms of
!  0.3, between SA of 0 and 120 g/kg and p between 0 and 5,000 dbar) when
!  compared with the potential enthalpy calculated from the exact in-situ
!  freezing temperature which is found by a Newton-Raphson iteration of the
!  equality of the chemical potentials of water in seawater and in ice.
!
!  SA  =  Absolute Salinity                                        [ g/kg ]
!  p   =  sea pressure                                             [ dbar ]
!         ( i.e. absolute pressure - 10.1325 dbar )
!
!  pot_enthalpy_ice_freezing = potential enthalpy of ice at freezing
!                              of seawater                         [ J/kg ]
!--------------------------------------------------------------------------

use gsw_mod_kinds

implicit none

real (r8), intent(in) :: sa, p

real (r8) :: gsw_pot_enthalpy_ice_freezing_poly

real (r8) :: p_r, sa_r, x

real (r8), parameter :: c0  = -3.333548730778702e5_r8
real (r8), parameter :: c1  = -1.249490228128056e4_r8
real (r8), parameter :: c2  =  0.891189273859881e4_r8
real (r8), parameter :: c3  = -2.405994758887321e4_r8
real (r8), parameter :: c4  =  3.217945710496395e4_r8
real (r8), parameter :: c5  = -2.374817375023954e4_r8
real (r8), parameter :: c6  =  0.651630522289954e4_r8
real (r8), parameter :: c7  = -2.034535061416256e4_r8
real (r8), parameter :: c8  = -0.252580687014574e4_r8
real (r8), parameter :: c9  =  0.021290274388826e4_r8
real (r8), parameter :: c10 =  0.315423710959628e3_r8
real (r8), parameter :: c11 = -0.239518382138314e3_r8
real (r8), parameter :: c12 =  0.379377450285737e3_r8
real (r8), parameter :: c13 =  0.822414256564615e3_r8
real (r8), parameter :: c14 = -1.781443326566310e3_r8
real (r8), parameter :: c15 = -0.160245473297112e3_r8
real (r8), parameter :: c16 = -1.923856387576336e3_r8
real (r8), parameter :: c17 =  2.522158744711316e3_r8
real (r8), parameter :: c18 =  0.268604113069031e3_r8
real (r8), parameter :: c19 =  0.967023925992424e3_r8
real (r8), parameter :: c20 = -1.052684746354551e3_r8
real (r8), parameter :: c21 = -0.184147500983788e3_r8
real (r8), parameter :: c22 = -0.263384562367307e3_r8

sa_r = sa*1e-2_r8
x = sqrt(sa_r)
p_r = p*1e-4_r8

gsw_pot_enthalpy_ice_freezing_poly = &
         c0 + sa_r*(c1 + x*(c2 + x*(c3 + x*(c4 + x*(c5 + c6*x))))) &
            + p_r*(c7 + p_r*(c8 + c9*p_r)) + sa_r*p_r*(c10 + p_r*(c12 &
            + p_r*(c15 + c21*sa_r)) + sa_r*(c13 + c17*p_r + c19*sa_r) &
            + x*(c11 + p_r*(c14 + c18*p_r) + sa_r*(c16 + c20*p_r + c22*sa_r)))

return
end function

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

max_stars_repo_name: minsukji/ci-debug

max_issues_repo_name: minsukji/ci-debug

max_forks_repo_name: minsukji/ci-debug


text:
      SUBROUTINE DFBS
C
C     FBS   L,U,B/X/V,Y,ISYM=0/V,Y,KSIGN=1/V,Y,IPREC=0/V,Y,ITYPE=0 $
C
C     ISYM  =  1  USE FBS
C           = -1  USE GFBS
C           =  0  CHOOSE WHICH BASED ON SUPPLIED INPUT
C     KSIGN =  1, SOLVE LUX= B
C             -1,       LUX=-B
C     IPREC = REQUESTED PRECISION - DEFAULT BASED ON INPUT OR SYSTEM(55)
C     ITYPE = REQUESTED TYPE OF X - DEFAULT IS LOGICAL CHOICE ON INPUT
C
C     REVISED  12/91 BY G.CHAN/UNISYS
C     FATAL ERROR IN FBS (NOT GFBS) IF INPUT MATRIX IS NOT A LOWER
C     TRIANGULAR FACTOR
C
      INTEGER         L,U,B,X,SBNM(2),DOSI(3),REFUS(3),OUTPT,SCR
      DIMENSION       ZZ(1)
      CHARACTER       UFM*23,UWM*25,UIM*29,SFM*25,SWM*27
      COMMON /XMSSG / UFM,UWM,UIM,SFM,SWM
      COMMON /BLANK / ISYM, KSIGN, IPREC, ITYPE
      COMMON /SYSTEM/ KSYSTM(65)
      COMMON /FBSX  / IL(7),IU(7),IB(7),IX(7),INX,IP1,IS1,ISCR
      COMMON /GFBSX / JL(7),JU(7),JB(7),JX(7),JNX,JP1,JS1
CZZ   COMMON /ZZDFB1/ Z(1)
      COMMON /ZZZZZZ/ Z(20000)
CZZ   COMMON /ZZDFB2/ ZZ(1)
      EQUIVALENCE     (ZZ(1),Z(1))
      EQUIVALENCE     (KSYSTM(55),KPREC),(KSYSTM(2),OUTPT)
      DATA    L, U, B, X, SCR   / 101,102,103,201,301 /
      DATA    SBNM  / 4HDFBS,1H /
      DATA    DOSI  / 4HSING, 4HDOUB, 4HMLTP/,  REFUS / 2*3H   ,3HREF/
C
C
      JU(1) = U
      CALL RDTRL (JU)
   10 IF (ISYM) 150,20,30
   20 ISYM  = -1
      IF (JU(1) .LT. 0) ISYM = 1
      GO TO 10
C
C     SET UP CALL TO FBS
C
   30 NOGO  = 0
      IL(1) = L
      CALL RDTRL (IL)
      IF (IL(1) .GT. 0) GO TO 40
      CALL MESAGE (30,198,L)
      NOGO  = 1
   40 CONTINUE
      IF (IL(4) .NE. 4) GO TO 100
      N     = IL(2)
      IB(1) = B
      CALL RDTRL (IB)
      IF (NOGO .EQ. 0) GO TO 50
      CALL MESAGE (-30,199,SBNM)
   50 CONTINUE
      INX   = KORSZ(Z)
      IPREC1= MAX0(IL(5),IB(5),IU(5))
      IF (IPREC1 .GT. 2) IPREC1 = IPREC1 - 2
      IF (IPREC1.LT.1 .OR. IPREC1.GT.2) IPREC1 = KPREC
      IF (IPREC.EQ.IPREC1 .OR. IPREC.EQ.0) GO TO 70
      IF (IPREC.LT.1 .OR. IPREC.GT.2) IPREC = 3
      WRITE  (OUTPT,60) SWM,DOSI(IPREC),REFUS(IPREC),SBNM,DOSI(IPREC1)
   60 FORMAT (A27,' 2163, REQUESTED ',A4,'LE PRECISION ',A3,' USED BY ',
     1        2A4,2H. ,A4,'LE PRECISION IS LOGICAL CHOICE')
      IF (IPREC .NE. 3) IPREC1 = IPREC
   70 IPREC = IPREC1
      IP1   = IPREC1
      IS1   = KSIGN
      LTYPE = IPREC1
      IF (IL(5).EQ.3 .OR. IL(5).EQ.4 .OR. IU(5).EQ.3 .OR. IU(5).EQ.4 .OR
CWKBR spr 93014  1   .IL(5).EQ.3 .OR. IL(5).EQ.4)  LTYPE = IPREC1 + 2
     1   .IB(5).EQ.3 .OR. IB(5).EQ.4)  LTYPE = IPREC1 + 2
      IF (ITYPE.EQ.0 .OR. ITYPE.EQ.LTYPE) GO TO 90
      JJ    = 1
      IF (ITYPE.LT.1 .OR. ITYPE.GT.4) JJ = 3
      WRITE  (OUTPT,80) SWM,ITYPE,REFUS(JJ),SBNM,LTYPE
   80 FORMAT (A27,' 2164, REQUESTED TYPE ',I4,2H, ,A3,' USED BY ',2A4,
     1       '. TYPE ',I4,' IS LOGICAL CHOICE.')
      IF (JJ .NE. 3) LTYPE = ITYPE
   90 ITYPE = LTYPE
      IX(5) = ITYPE
      IX(1) = X
      ISCR  = SCR
      CALL FBS (Z,Z)
      IX(3) = N
      IX(4) = 2
      IF (IX(3) .EQ. IX(2)) IX(4) = 1
      CALL WRTTRL (IX)
      GO TO 200
C
  100 CALL FNAME (IL(1),IL(2))
      WRITE  (OUTPT,110) IL(2),IL(3),IL(4)
  110 FORMAT ('0*** INPUT MATRIX ',2A4,' TO FBS MODULE IS NOT A LOWER ',
     1        'TRIANGULAR FACTOR.  FORM =',I4)
      CALL ERRTRC ('DFBS    ',110)
      GO TO 200
C
C     SET UP CALL TO GFBS
C
  150 JL(1) = L
      CALL RDTRL (JL)
      N     = JL(2)
      JB(1) = B
      CALL RDTRL (JB)
      JNX   = KORSZ(ZZ)
      IPREC1= MAX0(JL(5),JB(5),JU(5))
      IF (IPREC1 .GT. 2) IPREC1 = IPREC1 - 2
      IF (IPREC1.LT.1 .OR. IPREC1.GT.2) IPREC1 = KPREC
      IF (IPREC.EQ.IPREC1 .OR. IPREC.EQ.0) GO TO 160
      IF (IPREC.LT.1 .OR. IPREC.GT.2) IPREC = 3
      WRITE (OUTPT,60) SWM,DOSI(IPREC),REFUS(IPREC),SBNM,DOSI(IPREC1)
      IF (IPREC .NE. 3) IPREC1 = IPREC
  160 IPREC = IPREC1
      JP1   = IPREC1
      JS1   = KSIGN
      JX(1) = X
      LTYPE = IPREC1
      IF (JL(5).EQ.3 .OR. JL(5).EQ.4 .OR. JU(5).EQ.3 .OR. JU(5).EQ.4 .OR
     1   .JL(5).EQ.3 .OR. JL(5).EQ.4) LTYPE = IPREC1 + 2
      IF (ITYPE.EQ.0 .OR. ITYPE.EQ.LTYPE) GO TO 170
      JJ    = 1
      IF (ITYPE.LT.1 .OR. ITYPE.GT.4) JJ = 3
      WRITE (OUTPT,80) SWM,ITYPE,REFUS(JJ),SBNM,LTYPE
      IF (JJ .NE. 3) LTYPE = ITYPE
  170 ITYPE = LTYPE
      JX(5) = ITYPE
      CALL GFBS (ZZ,ZZ)
      JX(3) = N
      JX(4) = 2
      IF (JX(3) .EQ. JX(2)) JX(4) =  1
      CALL WRTTRL (JX)
C
  200 RETURN
      END


text:
*----------------------------------------------------------------------*
      subroutine del_zero_terms(fl_tgt,mode,op_info,thrsh)
*----------------------------------------------------------------------*
*     look for terms with factor below threshold and delete them
*----------------------------------------------------------------------*
      implicit none

      integer, parameter ::
     &     ntest = 000

      include 'stdunit.h'
      include 'opdim.h'
      include 'def_contraction.h'
      include 'mdef_operator_info.h'
      include 'def_formula_item.h'
      include 'def_formula_item_array.h'
      include 'def_formula_item_list.h'

      type(formula_item), target, intent(inout) ::
     &     fl_tgt
      ! only for debug output:
      type(operator_info), intent(in) ::
     &     op_info
      real(8) ::
     &     thrsh
      character(len=*), intent(in) ::
     &     mode

      integer ::
     &     idxop_tgt, iterm, iterms
      type(formula_item), pointer ::
     &     fl_tgt_pnt, fl_tgt_pnt_next, fl_tgt_current

      ! if first entry is an [END]: do nothing
      if (fl_tgt%command.eq.command_end_of_formula) return

      if (ntest.ge.100) then
        write(lulog,*) '========================='
        write(lulog,*) ' info from del_zero_terms'
        write(lulog,*) '========================='
        write(lulog,*) ' mode = ',trim(mode)
      end if

      if (fl_tgt%command.ne.command_set_target_init) then
        if (fl_tgt%command.ne.command_add_contribution)
     &       call quit(1,'del_zero_terms',
     &         'must start with [INIT] or [ADD]')
        idxop_tgt = fl_tgt%target
      end if

      ! first sum identical terms?
      if (trim(mode).eq.'sum'.or.trim(mode).eq.'SUM')
     &   call sum_terms(fl_tgt,iterms,op_info)

      iterm = 0
      fl_tgt_current => fl_tgt
      ! loop over target items
      tgt_loop: do

        ! new operator target ?
        if (fl_tgt_current%command.eq.command_set_target_init) then
          idxop_tgt = fl_tgt_current%target
          if (.not.associated(fl_tgt_current%next))
     &         call quit(1,'del_zero_terms',
     &         'unexpected end of list (target)')
          if (ntest.ge.100) then
            write(lulog,'(70("="))')
            write(lulog,*) 'New operator target: ',idxop_tgt
            write(lulog,'(70("="))')
          end if
          fl_tgt_current => fl_tgt_current%next
          if (fl_tgt_current%command.eq.command_end_of_formula)
     &                                             exit tgt_loop
        end if

        iterm = iterm+1
        if (ntest.ge.100) then
          write(lulog,*) 'current term: # ',iterm
          call prt_contr2(lulog,fl_tgt_current%contr,op_info)
        end if

        fl_tgt_pnt => fl_tgt_current%next

        ! delete term with factor below threshold
        if (abs(fl_tgt_current%contr%fac).lt.thrsh) then
          call delete_fl_node(fl_tgt_current)
          deallocate(fl_tgt_current)
        end if

        ! advance to next item
        fl_tgt_current => fl_tgt_pnt
        if (.not.associated(fl_tgt_current))
     &         call quit(1,'del_zero_terms',
     &         'unexpected end of list (target)')
        if (fl_tgt_current%command.eq.command_end_of_formula)
     &                                             exit tgt_loop

      end do tgt_loop
c dbg
c      print *,' del_zero_terms yields',iterm,' terms.'
c dbgend

      return
      end

max_stars_repo_name: ak-ustutt/GeCCo-public

max_issues_repo_name: ak-ustutt/GeCCo-public

max_forks_repo_name: ak-ustutt/GeCCo-public


module GCM_io_module
! contains subroutines handle to GCM outputs (ie, annual climatologies)
!   1. read the paths of configuration variables in the config file (config/GCM_io_condition)
!   2. read the continental variables (area, landfrac, temperature and runoff) from GCM outputs

    implicit none
    integer, parameter, private:: fmaxlen=500, vmaxlen=200
    double precision, parameter, private:: axis_relat_accuracy=1d-5
    ! maximum allowed relative difference between axis (lon, lat) from different files

    contains

    subroutine load_climatology(fID)

        use netcdf

        include 'combine_foam.inc'

        integer, intent(in):: fID
        character(len=fmaxlen), dimension(2):: area_file
        character(len=fmaxlen), dimension(nclimber):: climo_file, gmst_climo_file
        character(len=vmaxlen), dimension(2):: area_x_varname, area_y_varname, area_varname
        character(len=vmaxlen):: x_varname, y_varname, temp_varname, runf_varname, var_units, vname
        character(len=vmaxlen):: gmst_x_varname, gmst_y_varname, glob_temp_varname
        double precision:: fillvalue

        double precision, dimension(2+2*nclimber, nlon):: all_x
        double precision, dimension(2+2*nclimber, nlat):: all_y
        double precision, dimension(nlon, nlat):: totarea2D, landarea2D
        double precision, dimension(nclimber, nlon, nlat):: temperature2D, runoff2D, glob_temperature
        integer:: ierr, nax

        print *
        print *
        print *, 'Read GCM conditions file and load land inputs from GCM climatology files'
        print *, '------------------------------------------------------------------------'

        ! get netCDF file info
        call read_GCM_condition(fID, &
        !                       <><><><><>
                                co2climber, &
        !                       <><><><><>
                                area_file, area_x_varname, area_y_varname, area_varname, &
                                climo_file, x_varname, y_varname, temp_varname, runf_varname, &
                                gmst_climo_file, gmst_x_varname, gmst_y_varname, glob_temp_varname)

        ! get total area
        call load_variable('area', area_varname(1), area_x_varname(1), area_y_varname(1), &
                           single_input_file=area_file(1), varout2D=totarea2D, x=all_x(1,:), y=all_y(1,:), &
                           fill_missing=.true.)

        ! get land fraction/area
        call load_variable('landarea', area_varname(2), area_x_varname(2), area_y_varname(2), &
                           single_input_file=area_file(2), varout2D=landarea2D, x=all_x(2,:), y=all_y(2,:), &
                           totarea=totarea2D, fill_missing=.true.)

        ! get temperature
        call load_variable('temperature', temp_varname, x_varname, y_varname, &
                           multiple_input_file=climo_file, varout3D=temperature2D, &
                           xvec=all_x(3:2+nclimber,:), yvec=all_y(3:2+nclimber,:), &
                           landarea=landarea2D, fill_missing=.false., fillval_handling_option=ERROR_HANDLING_OPTION(2))
        ! note: ERROR_HANDLING_OPTION(2) tells the code what to do if it found missing-value where area > 0

        ! get runoff
        call load_variable('runoff', runf_varname, x_varname, y_varname, &
                           multiple_input_file=climo_file, varout3D=runoff2D, &
                           landarea=landarea2D, fill_missing=.false., fillval_handling_option=ERROR_HANDLING_OPTION(2))

        ! get optional global temperature
        if (glob_temp_varname /= '-') then
            call load_variable('temperature', glob_temp_varname, gmst_x_varname, gmst_y_varname, &
                            multiple_input_file=gmst_climo_file, varout3D=glob_temperature, &
                            xvec=all_x(3+nclimber:2+2*nclimber, :), yvec=all_y(3+nclimber:2+2*nclimber, :), &
                            fill_missing=.true.)
            ! number of axis to check
            nax = 2+2*nclimber
        else
            ! number of axis to check
            nax = 2+nclimber
        end if

        call check_axis(all_x(1:nax,:), all_y(1:nax,:), ERROR_HANDLING_OPTION(1))
        ! note: ERROR_HANDLING_OPTION(1) tells the code what to do if it find axis mismatch

        !<><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><>!
        !       Put the loaded variables in GEOCLIM global variables       !
        !
        ref_x_axis = all_x(1,:)
        ref_y_axis = all_y(1,:)
        ! unravel 2D-arrays in 1D-array, the first dimension (x) begin the most rapidly varying => "natural" order
        areaEarth   = reshape(totarea2D,  shape=(/npixel/))
        areaclimber = reshape(landarea2D, shape=(/npixel/))
        !    Transpose CO2 axis, from 1st dimension to last dimension --vvvvvvvvvvvvv
        Tairclimber = reshape(temperature2D, shape=(/npixel,nclimber/), order=(/2,1/))
        Runclimber  = reshape(runoff2D,      shape=(/npixel,nclimber/), order=(/2,1/))
        ! Optional: Compute GMST
        if (glob_temp_varname /= '-') then
            totarea = sum(totarea2D)
            do k=1,nclimber
                GMSTclimber(k) = sum(glob_temperature(k,:,:)*totarea2D) / totarea
            end do
        else
            GMSTclimber = -1d99 ! "internal" missing-value -> will be replaced by "output" missing-value
        end if
        !                                                                  !
        !<><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><>!

    end subroutine

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

    subroutine read_GCM_condition(fID, &
                                  CO2_levels, &
                                  area_file, area_x_varname, area_y_varname, area_varname, &
                                  climo_file, x_varname, y_varname, temp_varname, runf_varname, &
                                  gmst_climo_file, gmst_x_varname, gmst_y_varname, glob_temp_varname)

        use utils, only: add_path, read_comment

        integer, intent(in):: fID
        double precision, dimension(:), intent(out):: CO2_levels
        character(len=fmaxlen), dimension(2), intent(out):: area_file
        character(len=vmaxlen), dimension(2), intent(out):: area_x_varname, area_y_varname, area_varname
        character(len=fmaxlen), dimension(:), intent(out):: climo_file, gmst_climo_file
        character(len=vmaxlen), intent(out):: x_varname, y_varname, temp_varname, runf_varname
        character(len=vmaxlen), intent(out):: gmst_x_varname, gmst_y_varname, glob_temp_varname
        character(len=50):: tag

        integer:: nlevels
        integer:: k

        nlevels = size(CO2_levels)

        ! Get CO2 values
        call read_comment(fID)
        read(unit=fID, fmt=*) CO2_levels

        ! Get area and landfrac info
        do k = 1,2
            call read_comment(fID)
            read(unit=fID, fmt=*) area_file(k), area_x_varname(k), area_y_varname(k), area_varname(k)
            call add_path(area_file(k))
        end do

        ! Get climatology variables names
        call read_comment(fID)
        read(unit=fID, fmt='(A)') x_varname
        call read_comment(fID)
        read(unit=fID, fmt='(A)') y_varname
        call read_comment(fID)
        read(unit=fID, fmt='(A)') temp_varname
        call read_comment(fID)
        read(unit=fID, fmt='(A)') runf_varname

        ! Get climatology file paths
        call read_comment(fID)
        read(unit=fID, fmt=*) tag
        if (tag /= '<<--START-->>') then
            print *
            print *, 'Error: Expected "START" tag not found while reading GCM output file names'
            stop
        end if
        !
        do k = 1,nlevels
            call read_comment(fID)
            read(unit=fID, fmt=*) climo_file(k)
            call add_path(climo_file(k))
        end do
        !
        call read_comment(fID)
        read(unit=fID, fmt=*) tag
        if (tag /= '<<--STOP-->>') then
            print *
            print *, 'Error: Expected "STOP" tag not found while reading GCM output file names'
            print *, '       => number of input GCM files greater than size of CO2 axis'
            stop
        end if

        ! Get global temperature climatology variables names
        call read_comment(fID)
        read(unit=fID, fmt='(A)') gmst_x_varname
        call read_comment(fID)
        read(unit=fID, fmt='(A)') gmst_y_varname
        call read_comment(fID)
        read(unit=fID, fmt='(A)') glob_temp_varname

        if (glob_temp_varname /= '-') then

            ! Get global temperature climatology file paths
            call read_comment(fID)
            read(unit=fID, fmt=*) tag
            if (tag /= '<<--START-->>') then
                print *
                print *, 'Error: Expected "START" tag not found while reading GCM output file names for GMST'
                stop
            end if
            call read_comment(fID)
            read(unit=fID, fmt=*) tag
            if (tag == '<<--STOP-->>') then ! No file specified => keep same as "main" climatology
                gmst_x_varname = x_varname
                gmst_y_varname = y_varname
                gmst_climo_file = climo_file
            else
                backspace(unit=fID)
                !
                do k = 1,nlevels
                    call read_comment(fID)
                    read(unit=fID, fmt=*) gmst_climo_file(k)
                    call add_path(gmst_climo_file(k))
                end do
                !
                call read_comment(fID)
                read(unit=fID, fmt=*) tag
                if (tag /= '<<--STOP-->>') then
                    print *
                    print *, 'Error: Expected "STOP" tag not found while reading GCM output file names for GMST'
                    print *, '       => number of input GCM files greater than size of CO2 axis'
                    stop
                end if
            end if

        end if

    end subroutine

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

    subroutine load_variable(internal_varname, varname, x_varname, y_varname,         &
                             single_input_file, multiple_input_file,                  &
                             varout2D, varout3D, x, y, xvec, yvec, landarea, totarea, &
                             fill_missing, fillval_handling_option                    )
        use netcdf
        include 'shape.inc'

        character(len=*), intent(in):: internal_varname
        character(len=vmaxlen), intent(in):: varname, x_varname, y_varname
        character(len=fmaxlen), intent(in), optional:: single_input_file
        character(len=fmaxlen), dimension(nclimber), intent(in), optional:: multiple_input_file
        double precision, dimension(nlon,nlat), intent(inout), optional:: landarea
        double precision, dimension(nlon,nlat), intent(in), optional:: totarea
        double precision, dimension(nlon,nlat), intent(out), optional:: varout2D
        double precision, dimension(nclimber, nlon,nlat), intent(out), optional:: varout3D
        double precision, intent(out), optional:: x(nlon), y(nlat), xvec(nclimber,nlon), yvec(nclimber,nlat)
        integer, intent(in), optional:: fillval_handling_option
        logical, intent(in), optional:: fill_missing

        double precision, dimension(nlon,nlat):: dummyvar2D
        double precision, dimension(nclimber,nlon,nlat):: dummyvar3D
        double precision:: fillvalue
        character(len=vmaxlen):: var_units, vname
        character(len=1):: oper
        logical:: loc_fill_missing
        integer:: ierr, check, loc_fillval_handling_option
        integer:: k, n, noperations

        ! How to handle missing-value
        if (present(fill_missing)) then
            loc_fill_missing = fill_missing
        else
            loc_fill_missing = .false. ! default behaviour
        end if
        !
        ! Note: checking missing-points is done only if fill_missing==.false. and "landarea" argument is given
        if (present(fillval_handling_option)) then
            loc_fillval_handling_option = fillval_handling_option
        else
            loc_fillval_handling_option = -1 ! default value => ask user
        end if

        ! Check consistency between 2D variable or 3D (ie, multiple CO2 levels) variable
        ! => the right "input_file" and "var" arguments must be given.
        check = -11
        if (present(varout2D)) then
            check = check + 1
        end if
        if (present(single_input_file)) then
            check = check + 10
        end if
        if (present(varout3D)) then
            check = check + 2
        end if
        if (present(multiple_input_file)) then
            check = check + 9
        end if
        if (check /= 0) then
            print *
            print *, 'INTERNAL ERROR: inconsistent set of optional variable passed to the subroutine'
            print *, '"load_variable", in module "GCM_io_module"'
            stop
        end if

        ! Instead of reading just 1 variable name, they may be several variables,
        ! with arithmetic operation between them (addition or substraction)
        noperations = get_arithmetic_operations(varname)
        ! => write the scratch files:
        !        unit=334: list of individual variables
        !        unit=335: list of mathematics operators (+ or -)

        if (present(varout2D)) then

            ! Initialization
            varout2D = 0d0

            ! Loop to get variables and perform arithmetic operations:
            do n = 1,noperations

                ! load current variable
                read(unit=334, fmt=*) vname
                if (loc_fill_missing) then
                    ! load variable
                    if (present(x) .and. present(y)) then
                        call load_netcdf_dble2D(single_input_file, x_varname, y_varname, vname, dummyvar2D, &
                                                x=x, y=y, varunits=var_units, fillval=fillvalue, fillval_iostat=ierr)
                    else
                        call load_netcdf_dble2D(single_input_file, x_varname, y_varname, vname, dummyvar2D, &
                                                varunits=var_units, fillval=fillvalue, fillval_iostat=ierr)
                    end if
                    ! set var=0 on "missing" cells
                    if (ierr==NF90_NOERR) where (dummyvar2D==fillvalue) dummyvar2D = 0d0
                    ! check variable units
                    if (present(totarea)) then
                        call check_units(internal_varname, vname, dummyvar2D, var_units, totarea=totarea)
                    else
                        call check_units(internal_varname, vname, dummyvar2D, var_units)
                    end if
                else
                    ! load variable
                    if (present(x) .and. present(y)) then
                        call load_netcdf_dble2D(single_input_file, x_varname, y_varname, vname, dummyvar2D, &
                                                x=x, y=y, varunits=var_units, fillval=fillvalue)
                    else
                        call load_netcdf_dble2D(single_input_file, x_varname, y_varname, vname, dummyvar2D, &
                                                varunits=var_units, fillval=fillvalue)
                    end if
                    ! check missingpoint
                    if (present(landarea)) then
                        call check_missingpoints(vname, single_input_file, dummyvar2D, fillvalue, landarea, &
                                                 loc_fillval_handling_option)
                    end if
                    ! check variable units
                    if (present(totarea)) then
                        call check_units(internal_varname, vname, dummyvar2D, var_units, fillvalue=fillvalue, totarea=totarea)
                    else
                        call check_units(internal_varname, vname, dummyvar2D, var_units, fillvalue=fillvalue)
                    end if
                end if

                ! perform arithmetic operation
                read(unit=335, fmt=*) oper
                select case (oper)
                    case ("+")
                        varout2D = varout2D + dummyvar2D
                    case ("-")
                        varout2D = varout2D - dummyvar2D
                end select

            end do

        else ! varout3D case => load var for each CO2 level

            ! Initialization
            varout3D = 0d0

            ! Loop to get variables and perform arithmetic operations:
            do n = 1,noperations

                ! load current variable
                read(unit=334, fmt=*) vname
                do k = 1,nclimber
                    if (loc_fill_missing) then
                        ! load variable
                        if (present(xvec) .and. present(yvec)) then
                            call load_netcdf_dble2D(multiple_input_file(k), x_varname, y_varname, vname, dummyvar3D(k,:,:), &
                                               x=xvec(k,:), y=yvec(k,:), varunits=var_units, fillval=fillvalue, fillval_iostat=ierr)
                        else
                            call load_netcdf_dble2D(multiple_input_file(k), x_varname, y_varname, vname, dummyvar3D(k,:,:), &
                                                    varunits=var_units, fillval=fillvalue, fillval_iostat=ierr)
                        end if
                        ! set var=0 on "missing" cells
                        if (ierr==NF90_NOERR) where (dummyvar3D(k,:,:)==fillvalue) dummyvar3D(k,:,:) = 0d0
                        ! check variable units
                        if (present(totarea)) then
                            call check_units(internal_varname, vname, dummyvar3D(k,:,:), var_units, totarea=totarea)
                        else
                            call check_units(internal_varname, vname, dummyvar3D(k,:,:), var_units)
                        end if
                    else
                        ! load variable
                        if (present(xvec) .and. present(yvec)) then
                            call load_netcdf_dble2D(multiple_input_file(k), x_varname, y_varname, vname, dummyvar3D(k,:,:), &
                                                    x=xvec(k,:), y=yvec(k,:), varunits=var_units, fillval=fillvalue)
                        else
                            call load_netcdf_dble2D(multiple_input_file(k), x_varname, y_varname, vname, dummyvar3D(k,:,:), &
                                                    varunits=var_units, fillval=fillvalue)
                        end if
                        ! check missingpoint
                        if (present(landarea)) then
                            call check_missingpoints(vname, multiple_input_file(k), dummyvar3D(k,:,:), fillvalue, landarea, &
                                                     loc_fillval_handling_option)
                        end if
                        ! check variable units
                        if (present(totarea)) then
                            call check_units(internal_varname, vname, dummyvar3D(k,:,:), var_units, fillvalue=fillvalue, &
                                             totarea=totarea)
                        else
                            call check_units(internal_varname, vname, dummyvar3D(k,:,:), var_units, fillvalue=fillvalue)
                        end if
                    end if
                end do

                ! perform arithmetic operation
                read(unit=335, fmt=*) oper
                select case (oper)
                    case ("+")
                        varout3D = varout3D + dummyvar3D
                    case ("-")
                        varout3D = varout3D - dummyvar3D
                end select

            end do

        end if

        ! close scratch files
        close(unit=334)
        close(unit=335)

    end subroutine

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

    function get_arithmetic_operations(varstring)
    ! => create scratch files 334 and 335 containing (respectively) the lists of variables and operators
    ! Note: the first operators (corresponding to the first variable) is automatically '+'
    ! return the number of variables/operations
        character(len=*), intent(in):: varstring
        integer:: get_arithmetic_operations
        character(len=1), dimension(2), parameter:: operators = (/'+', '-'/) ! Legal arithmetic operators
        integer:: i0, i1
        open(unit=334, status='scratch') ! individual variables
        open(unit=335, status='scratch') ! operators
        i0 = 1
        i1 = 1
        get_arithmetic_operations = 1
        write(unit=335, fmt='(A)') '+'
        do i1 = 1,len(varstring)
            if (any(operators == varstring(i1:i1))) then
                write(unit=334, fmt='(A)') trim(adjustl(varstring(i0:i1-1)))
                write(unit=335, fmt='(A)') varstring(i1:i1)
                i0 = i1+1
                get_arithmetic_operations = get_arithmetic_operations + 1
            end if
        end do
        write(unit=334, fmt='(A)') trim(adjustl(varstring(i0:)))
        rewind(unit=334)
        rewind(unit=335)
    end function

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

    subroutine load_netcdf_dble2D(fname, x_varname, y_varname, varname, var, x, y, varunits, fillval, fillval_iostat)

        use netcdf

        character(len=*), intent(in):: fname, x_varname, y_varname, varname
        double precision, dimension(:,:), intent(out):: var
        double precision, dimension(:), intent(out), optional:: x, y
        character(len=*), intent(out), optional:: varunits
        double precision, intent(out), optional:: fillval
        integer, intent(out), optional:: fillval_iostat
        integer:: nx, ny, ndim, k, n
        integer:: ierr, fid, varid, xdimid, ydimid, shp_ix, shp_iy
        integer, dimension(:), allocatable:: dimids, shp
        logical:: transp
        double precision, dimension(:,:), allocatable:: loc_var

        nx = size(var, 1)
        ny = size(var, 2)
        ! must also match the size of "x" and "y"

        ! open netCDF file
        ierr = nf90_open(fname, NF90_NOWRITE, fid)
        call nf90_check(ierr, 'Error while openning file '//trim(fname))

        ! Get axis data
        if (present(x)) then
            call load_axis(fname, fid, x_varname, xdimid, x)
        else
            call load_axis(fname, fid, x_varname, xdimid)
        end if
        if (present(y)) then
            call load_axis(fname, fid, y_varname, ydimid, y)
        else
            call load_axis(fname, fid, y_varname, ydimid)
        end if

        ! Get variable ID
        ierr = nf90_inq_varid(fid, varname, varid)
        call nf90_check(ierr, 'Error while inquiring variable "'//trim(varname)//'" ID in file "'//trim(fname)//'"')

        ! Get variable number of dimension
        ierr = nf90_inquire_variable(fid, varid, ndims=ndim)
        call nf90_check(ierr, 'Error while inquiring number of dim. of variable "'//trim(varname)//'" in file "'//trim(fname)//'"')

        allocate(dimids(ndim))
        allocate(shp(ndim))

        ! Get variable dimension IDs
        ierr = nf90_inquire_variable(fid, varid, dimids=dimids)
        call nf90_check(ierr, 'Error while inquiring dimensions IDs of variable "'//trim(varname)//'" in file "'//trim(fname)//'"')

        ! Get variable shape (inquire length of all dimensions)
        do k = 1,ndim
            ierr = nf90_inquire_dimension(fid, dimids(k), len=shp(k))
            call nf90_check(ierr, 'Error while inquiring length of dimension in file "'//trim(fname)//'"')
        end do

        ! Check variable shape (must have exactly 2 non degenerated dimensions)
        if (count(shp > 1) /= 2) then
            print *
            print *, 'Error: variable "'//trim(varname)//'" of file "'//trim(fname)// &
                     '" must have exactly 2 non-degenerated (ie, size-1) dimensions'
            stop
        else
            k = 1
            do while (shp(k)==1)
                k = k + 1
            end do
            shp_ix = k
            k = k + 1
            do while (shp(k)==1)
                k = k + 1
            end do
            shp_iy = k
        end if

        ! Check that variable is defined on the given dimensions
        if (dimids(shp_ix)==xdimid .and. dimids(shp_iy)==ydimid) then
            transp = .false.
            allocate(loc_var(nx,ny))
        elseif (dimids(shp_ix)==ydimid .and. dimids(shp_iy)==xdimid) then
            transp = .true.
            allocate(loc_var(ny,nx))
        else
            print *
            print *
            print *, 'Error: variable "'//trim(varname)//'" of file "'//trim(fname)// &
            '" is not defined on the given dimensions "'//trim(x_varname)//'" and "'//trim(y_varname)//'"'
            stop
        end if

        ! load variable
        if (ndim>2) then
            print *
            print *,                 'Variable "'//trim(varname)//'" of file "'//trim(fname)//'"'
            write(*, fmt='(A,I2,A)') '    Ignore ', ndim-2, ' degenerated (size-1) dimension(s).'
        end if
        ! + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + !
        ierr = nf90_get_var(fid, varid, loc_var)
        call nf90_check(ierr, 'Error while getting variable "'//trim(varname)//'" of file "'//trim(fname)//'"')
        ! + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + !
        if (transp) then
            var = transpose(loc_var)
        else
            var = loc_var
        end if

        ! Get variable units (if asked)
        if (present(varunits)) then
            varunits = ''
            ierr = nf90_get_att(fid, varid, 'units', varunits)
            call nf90_check(ierr, 'Warning: unable to get attribute "units" of variable "'//trim(varname)//'" in file "' &
                                   //trim(fname)//'". Variable assumed to be dimensionless.', kill=.false.)
        end if

        ! Get variale fill-value (if asked)
        if (present(fillval)) then
            ierr = nf90_get_att(fid, varid, '_FillValue', fillval)
            if (present(fillval_iostat)) then
                fillval_iostat = ierr
            else
                call nf90_check(ierr, 'Error while getting attribute "_FillValue" of variable "'//trim(varname)//'" in file "' &
                                                                                //trim(fname)//'"')
            end if
        end if

        ! close netCDF file
        ierr = nf90_close(fid)
        call nf90_check(ierr, 'Error while closing file '//trim(fname), kill=.false.)

        deallocate(dimids)
        deallocate(shp)
        deallocate(loc_var)

    end subroutine

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

    subroutine load_axis(fname, fid, axname, axdimid, ax)
        use netcdf
        character(len=*), intent(in):: fname, axname
        integer, intent(in):: fid
        integer, intent(out):: axdimid
        double precision, dimension(:), intent(out), optional:: ax
        integer:: varid, ierr

        ! Inquire dim ID
        ierr = nf90_inq_dimid(fid, axname, axdimid)
        call nf90_check(ierr, 'Error while inquiring dimension "'//trim(axname)//'" ID in file "'//trim(fname)//'"')
        ! Inquire var ID
        ierr = nf90_inq_varid(fid, axname, varid)
        call nf90_check(ierr, 'Error while inquiring variable "'//trim(axname)//'" ID in file "'//trim(fname)//'"')
        ! Get variable
        if (present(ax)) then
            ierr = nf90_get_var(fid, varid, ax)
            call nf90_check(ierr, 'Error while getting variable "'//trim(axname)//'" in file "'//trim(fname)//'"')
        end if

    end subroutine

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

    subroutine compare_units(units_string, known_units, conversion_factor, conversion_offset, passed)

        use physical_units, only: units

        character(len=*), intent(in):: units_string
        type(units), intent(in):: known_units
        double precision, intent(out):: conversion_factor, conversion_offset
        logical, intent(out):: passed
        integer:: k

        conversion_factor = 1d0
        conversion_offset = 0d0

        passed = .false.

        if (units_string == known_units%reference) then
            passed = .true.
        else
            do k = 1,known_units%naccepted
                if (units_string == known_units%accepted(k)%string) then
                    passed = .true.
                    conversion_factor = known_units%accepted(k)%conversion(1)
                    conversion_offset = known_units%accepted(k)%conversion(2)
                end if
            end do
        end if

    end subroutine

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

    subroutine check_units(which_variable, varname, var, varunits, fillvalue, totarea, error_handling)

        use physical_units, only: units, area_units, fraction_units, temperature_units, runoff_units

        character(len=*), intent(in):: which_variable, varname, varunits
        double precision, dimension(:,:), intent(inout):: var
        double precision, dimension(:,:), intent(in), optional:: totarea
        double precision, intent(in), optional:: fillvalue
        integer, intent(in), optional:: error_handling
        type(units):: known_units
        double precision:: factor, offset
        logical:: multiply_by_area, passed, loop
        integer:: k, loc_error_handling, answer

        if (present(error_handling)) then
            loc_error_handling = error_handling
        else
            loc_error_handling = -1 ! default behaviour: ask user
        end if

        multiply_by_area = .false.

        ! case-dependent statement
        select case (which_variable)

            case ("area")
                known_units = area_units()

            case ("landarea")
                known_units = area_units()

            case ("temperature")
                known_units = temperature_units()

            case ("runoff")
                known_units = runoff_units()

            case default
                print *
                print *, 'INTERNAL ERROR in function "check_units" of module "GCM_io_module": unkown variable case "' &
                                                                                //trim(which_variable)//'"'
                stop

        end select

        ! Compare given units to reference units, and get conversion values if not equals
        call compare_units(varunits, known_units, factor, offset, passed)
        ! in the case "landarea", also try "fraction" units
        if (which_variable == "landarea" .and. (.not. passed)) then
            call compare_units(varunits, fraction_units(), factor, offset, passed)
            multiply_by_area = passed
        end if

        ! Print message
        if (passed) then
            if (factor/=1 .or. offset/=0 .or. multiply_by_area) then
                print *
                print *, 'Automatic conversion of variable "'//trim(varname)//'":'
                print *, '    "'//trim(varunits)//'" => "'//trim(known_units%reference)//'"'
            end if
        else
            print *
            print *, 'WARNING: unkown units of variable "'//trim(varname)//'".'
            print *, '    got units:      "'//trim(varunits)//'"'
            print *, '    expected units: "'//trim(known_units%reference)//'"'
            ! Error handling:
            select case (loc_error_handling)
                case (-1) ! ask user interactively
                    print *
                    loop = .true.
                    do while (loop)
                        loop = .false.
                        print *, 'Enter one of the following options:'
                        print *, '    0: abort the program'
                        print *, '    2: ignore the issue and continue the execution'
                        read *, answer
                        select case (answer)
                            case (0); stop
                            case (2) ! do nothing
                            case default; loop=.true.
                        end select
                    end do

                case (0) ! abort the program
                    stop

                case (2) ! do nothing

                case default
                    print *
                    print *, 'INTERNAL ERROR: illegal error handling option:', error_handling
                    stop
            end select
        end if

        ! Conversion
        if (present(fillvalue)) then
            where (var/=fillvalue)  var = factor*(var + offset)
        else
            var = factor*(var + offset)
        end if

        ! Multiply value by cell area
        if (multiply_by_area) then
            if (present(totarea)) then
                print *, 'Variable "'//trim(varname)//'" multiplied by total area'
                if (present(fillvalue)) then
                    where (var/=fillvalue) var = var * totarea
                else
                    var = var * totarea
                end if
            else
                print *
                print *, 'INTERNAL ERROR: total area not passed to the function "check_units" in module "GCM_io_module".'
                print *, 'You have the right to be mad at the developer.'
                stop
            end if
        end if

    end subroutine

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

    subroutine check_missingpoints(varname, filename, var, fillvalue, landarea, error_handling)

        character(len=*), intent(in):: varname, filename
        double precision, dimension(:,:), intent(in):: var
        double precision, intent(in):: fillvalue
        double precision, dimension(:,:), intent(inout):: landarea
        integer, intent(in):: error_handling
        logical, dimension(:,:), allocatable:: errormask
        double precision:: area_err, tot_landarea
        integer:: nerr, answer
        logical:: loop

        tot_landarea = sum(landarea)

        allocate(errormask(size(landarea,1), size(landarea,2)))

        errormask = (landarea>0 .and. var==fillvalue)

        nerr = count(errormask)
        area_err = sum(landarea, mask=errormask)

        ! If error found:
        if (nerr > 0) then

            print *
            write(*,'(A)')         ' WARNING: found missing-value on continental points (land area > 0) in variable "' &
                                    //trim(varname)//'" of file "'//trim(filename)//'"'
            write(*,'(A,I0,A,I0)') '     Number of problematic continent cells:  ', nerr, ' / ', count(landarea>0)
            write(*,'(A,E14.7)')   '     Total area of those cells (km2):        ', 1d-6*area_err
            write(*,'(A,E14.7)')   '     Which is a fraction of total land area: ', area_err/tot_landarea

            ! Error handling:
            select case (error_handling)
                case (-1) ! ask user interactively
                    print *
                    loop = .true.
                    do while (loop)
                        loop = .false.
                        print *, 'Enter one of the following options:'
                        print *, '    0: abort the program'
                        print *, '    1: remove erratic points (set area=0)'
                        print *, '    2: ignore the issue and continue the execution'
                        read *, answer
                        select case (answer)
                            case (0); stop
                            case (1); where (errormask) landarea=0
                            case (2) ! do nothing
                            case default; loop=.true.
                        end select
                    end do

                case (0) ! abort the program
                    stop

                case (1) ! automatic correction
                    print *
                    print *, 'Automatic correction => remove erratic points (set area=0)'
                    where (errormask) landarea=0

                case (2) ! do nothing

                case default
                    print *
                    print *, 'INTERNAL ERROR: illegal error handling option:', error_handling
                    stop

            end select

        end if

        deallocate(errormask)

    end subroutine

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

    subroutine raise_axis_error(which_axis, nerr, axis_len, max_mismatch, error_handling)
        character(len=1), intent(in):: which_axis
        integer, intent(in):: nerr, axis_len
        double precision, intent(in):: max_mismatch
        integer, intent(in):: error_handling
        integer:: answer
        logical:: loop

        print *
        write(*,'(A)')         ' WARNING: found mismatch of '//which_axis//' axis between GCM input files.'
        write(*,'(A,I0,A,I0)') '     Number of points concerned / axis length:  ', nerr, ' / ', axis_len
        write(*,'(A,E14.7)')   '     Maximum mismatch found (in axis units):    ', max_mismatch

        ! Error handling:
        select case (error_handling)
            case (-1) ! ask user interactively
                print *
                loop = .true.
                do while (loop)
                    loop = .false.
                    print *, 'Enter one of the following options:'
                    print *, '    0: abort the program'
                    print *, '    2: ignore the issue and continue the execution'
                    read *, answer
                    select case (answer)
                        case (0); stop
                        case (2) ! do nothing
                        case default; loop=.true.
                    end select
                end do

            case (0) ! abort the program
                stop

            case (1) ! asked for automatic correction
                print *
                print *, 'Bad error handling option "1": cannot remove points with axis mismatch.'
                stop

            case (2) ! do nothing

            case default
                print *
                print *, 'INTERNAL ERROR: illegal error handling option:', error_handling
                stop

        end select

    end subroutine

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

    subroutine check_axis(all_x, all_y, error_handling)
        double precision, dimension(:,:), intent(in):: all_x, all_y ! dim #1: files (temp, runoff, ...), dim #2: axis
        integer, intent(in):: error_handling
        double precision:: daxis, max_mismatch
        logical, dimension(:), allocatable:: errormask
        integer:: n, k, axlen, nerr

        n = size(all_x, 1) ! also = size(all_y, 1)

        ! x axis
        ! ------

        axlen = size(all_x, 2)
        allocate(errormask(axlen))
        errormask = .false.

        daxis = abs(all_x(1,2) - all_x(1,1))
        max_mismatch = 0d0

        do k = 2,n
            errormask = (errormask  .or.  maxval(abs(all_x(k,:) - all_x(1,:)))/daxis > axis_relat_accuracy)
            max_mismatch = max(max_mismatch, maxval(abs(all_x(k,:) - all_x(1,:))))
        end do
        nerr = count(errormask)
        if (nerr > 0) then
            call raise_axis_error('x', nerr, axlen, max_mismatch, error_handling)
        end if

        deallocate(errormask)

        ! y axis
        ! ------

        axlen = size(all_y, 2)
        allocate(errormask(axlen))
        errormask = .false.

        daxis = abs(all_y(1,2) - all_y(1,1))
        max_mismatch = 0d0

        do k = 2,n
            errormask = (errormask  .or.  maxval(abs(all_y(k,:) - all_y(1,:)))/daxis > axis_relat_accuracy)
            max_mismatch = max(max_mismatch, maxval(abs(all_y(k,:) - all_y(1,:))))
        end do
        nerr = count(errormask)
        if (nerr > 0) then
            call raise_axis_error('y', nerr, axlen, max_mismatch, error_handling)
        end if

        deallocate(errormask)

    end subroutine

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

    subroutine nf90_check(ierr, message, kill)

        use netcdf
        integer, intent(in):: ierr
        character(len=*), intent(in):: message
        logical, optional, intent(in):: kill
        logical:: loc_kill

        if (present(kill)) then
            loc_kill = kill
        else
            loc_kill = .true.
        end if

        if (ierr/=NF90_NOERR) then
            print *
            print *, message
            print *, nf90_strerror(ierr)
            if (loc_kill) stop
        end if

    end subroutine

end module


text:
      subroutine dblas()
      d = dasum(n,dx,incx)
      call daxpy(n,da,dx,incx,dy,incy)
      call  dcopy(n,dx,incx,dy,incy)
      d = ddot(n,dx,incx,dy,incy)
      d = dmach(job)
      d = dnrm2 ( n, dx, incx)
      call  drot (n,dx,incx,dy,incy,c,s)
      call drotg(da,db,c,s)
      call  dscal(n,da,dx,incx)
      call  dswap (n,dx,incx,dy,incy)
      i = idamax(n,dx,incx)
      stop
      end

max_stars_repo_name: NetlibRepository/lapack-www

max_issues_repo_name: NetlibRepository/lapack-www

max_forks_repo_name: NetlibRepository/lapack-www


text:
*deck s0000
      subroutine s0000( etai,           etaj,          etak,
     &  etal,           g4,             g3,            g2,
     &  g1,             gout,           zet )
      implicit real*8 (a-h,o-z)
      logical esfc, esfcij, esfckl, igueq1, jgueq1, kgueq1, lgueq1
      parameter (a1s2=0.5d0, asrtpi=1.1283791670955126d0,
     1  pie4=7.85398163397448d-01)
      common /parmr/ cutoff, tol
c
      integer
     & mblu,  mccu, mconu, mcu,  mcxu, mpru, nbft, nnbft,
     & mrcru, mstu, msu,   nblu, ng,   ns,   nst
      common /parmi/
     & mblu,  mccu, mconu, mcu,  mcxu, mpru, nbft, nnbft,
     & mrcru, mstu, msu,   nblu, ng,   ns,   nst
c
      common /shlinf/ xi,yi,zi,rij,xk,yk,zk,rkl,igu,jgu,kgu,lgu,nc4,
     1  nc3,nc2,nc1,igueq1,jgueq1,kgueq1,lgueq1
      common/misc/xij,xijm,yij,yijm,zij,zijm,xkl,xklm,ykl,yklm,zkl,zklm,
     1icons,ircru,jcons,jrcru,kcons,krcru,lcons,lrcru,esfc,esfcij,esfckl
      dimension etai(mrcru,*),etaj(mrcru,*),etak(mrcru,*),
     1  etal(mrcru,*),g4(*),g3(*),g2(*),g1(*),gout(*),zet(mconu,*)
c
      abx=xij-xkl
      aby=yij-ykl
      abz=zij-zkl
      if(igueq1) then
        fctr1=asrtpi*etai(1,1)
      else
        fctr1=asrtpi
      endif
      if(jgueq1) fctr1=fctr1*etaj(1,1)
      if(kgueq1) fctr1=fctr1*etak(1,1)
      if(lgueq1) fctr1=fctr1*etal(1,1)
      kgl=1
c
c  i primitive
c
      do 550 ig=1,igu
      if(.not.igueq1) call wzero(nc3,g3,1)
      fctr2=fctr1
      ai=zet(ig,icons)
      arij=ai*rij
      if(esfc) kgl=ig
c
c  j primitive
c
      if(esfcij) jgu=ig
      do 350 jg=1,jgu
      if(.not.jgueq1) call wzero(nc2,g2,1)
      if((esfcij.and..not.igueq1).and.jg.eq.ig) fctr2=a1s2*fctr2
      aj=zet(jg,jcons)
      aa=ai+aj
      aarij=aj*arij/aa
      if(aarij.gt.tol) go to 350
      fctr3=fctr2/aa
      aaa=(ai-aj)/aa
      px=abx+aaa*xijm
      py=aby+aaa*yijm
      pz=abz+aaa*zijm
c
c  k primitive
c
      do 250 kg=kgl,kgu
      if(.not.kgueq1) call wzero(nc1,g1,1)
      bk=zet(kg,kcons)
      brkl=bk*rkl
c
c  l primitive
c
      if(esfc.and.(kg.eq.ig)) then
        lgl=jg
      else
        lgl=1
      endif
      if(esfckl) lgu=kg
      do 150 lg=lgl,lgu
      bl=zet(lg,lcons)
      bb=bk+bl
      dum=aarij+bl*brkl/bb
      if(dum.gt.tol) go to 150
      expe=fctr3*exp(-dum)/bb
      if(((esfc.and..not.(igueq1.and.jgueq1)).and.(kg.eq.ig)).and.
     1  lg.eq.jg) expe=a1s2*expe
      if((esfckl.and..not.kgueq1).and.lg.eq.kg) expe=a1s2*expe
      bbb=(bk-bl)/bb
      apb=aa+bb
      x=((px-bbb*xklm)**2+(py-bbb*yklm)**2+(pz-bbb*zklm)**2)*aa*bb/apb
      if(x.gt.5.0d+00) go to 50
      if(x.gt.1.0d+00) go to 30
      if(x.gt.3.0d-7) go to 20
      ww1=1.0d+00-x/3.0d+00
      go to 100
   20 f1=          ((((((((-8.36313918003957d-08*x+1.21222603512827d-06
     1)*x-1.15662609053481d-05 )*x+9.25197374512647d-05
     2)*x-6.40994113129432d-04 )*x+3.78787044215009d-03
     3)*x-1.85185172458485d-02 )*x+7.14285713298222d-02
     4)*x-1.99999999997023d-01 )*x+3.33333333333318d-01
      ww1=(x+x)*f1+exp(-x)
      go to 100
   30 if(x.gt.3.0d+00) go to 40
      y=x-2.0d+00
      f1=        ((((((((((-1.61702782425558d-10*y+1.96215250865776d-09
     1)*y-2.14234468198419d-08 )*y+2.17216556336318d-07
     2)*y-1.98850171329371d-06 )*y+1.62429321438911d-05
     3)*y-1.16740298039895d-04 )*y+7.24888732052332d-04
     4)*y-3.79490003707156d-03 )*y+1.61723488664661d-02
     5)*y-5.29428148329736d-02 )*y+1.15702180856167d-01
      ww1=(x+x)*f1+exp(-x)
      go to 100
   40 y=x-4.0d+00
      f1=        ((((((((((-2.62453564772299d-11*y+3.24031041623823d-10
     1)*y-3.614965656163d-09)*y+3.760256799971d-08)*y-3.553558319675d-07
     2)*y+3.022556449731d-06)*y-2.290098979647d-05)*y+1.526537461148d-04
     3)*y-8.81947375894379d-04 )*y+4.33207949514611d-03
     4)*y-1.75257821619926d-02 )*y+5.28406320615584d-02
      ww1=(x+x)*f1+exp(-x)
      go to 100
   50 if(x.gt.15.0d+00) go to 70
      e=exp(-x)
      if(x.gt.10.0d+00) go to 60
      ww1=    (((((( 4.6897511375022d-01/x-6.9955602298985d-01)/x
     1+5.3689283271887d-01)/x-3.2883030418398d-01)/x
     2+2.4645596956002d-01)/x-4.9984072848436d-01)/x
     3-3.1501078774085d-06)*e + sqrt(pie4/x)
      go to 100
   60 ww1=       (((-1.8784686463512d-01/x+2.2991849164985d-01)/x
     1-4.9893752514047d-01)/x-2.1916512131607d-05)*e + sqrt(pie4/x)
      go to 100
   70 if(x.gt.33.0d+00) go to 90
      e=exp(-x)
      ww1=        (( 1.9623264149430d-01/x-4.9695241464490d-01)/x
     1-6.0156581186481d-05)*e + sqrt(pie4/x)
      go to 100
   90 ww1=sqrt(pie4/x)
  100 gout(1)=ww1*expe/sqrt(apb)
c
c  l transformation
c
      if(lgueq1) go to 200
c
      do 140 lrcr=1,lrcru
        g1(lrcr)=g1(lrcr)+gout(1)*etal(lrcr,lg)
  140 continue
  150 continue
c
c  k transformation
c
  200 if(kgueq1) go to 280
c
      kl2=0
      do 240 krcr=1,krcru
        if(esfckl) then
          do 210 lrcr=1,krcr
            g2(kl2+lrcr)=g2(kl2+lrcr)+g1(lrcr)*etak(krcr,kg)+
     1                                g1(krcr)*etak(lrcr,kg)
  210     continue
          kl2=kl2+krcr
        else
          do 230 lrcr=1,lrcru
            g2(kl2+lrcr)=g2(kl2+lrcr)+g1(lrcr)*etak(krcr,kg)
  230     continue
          kl2=kl2+lrcru
        endif
  240 continue
  250 continue
c
c  j transformation
c
  280 if(jgueq1) go to 410
c
      jkl3=0
      if(esfc.and.igueq1) then
        do 310 jrcr=1,jrcru
          do 300 lrcr=1,jrcr
            g3(jkl3+lrcr)=g3(jkl3+lrcr)+g2(lrcr)*etaj(jrcr,jg)+
     1                                  g2(jrcr)*etaj(lrcr,jg)
  300     continue
          jkl3=jkl3+jrcr
  310   continue
      else
        do 340 jrcr=1,jrcru
          kl2=0
          do 330 krcr=1,krcru
            if(esfckl) lrcru=krcr
            do 320 lrcr=1,lrcru
              g3(jkl3+lrcr)=g3(jkl3+lrcr)+g2(kl2+lrcr)*etaj(jrcr,jg)
  320       continue
            jkl3=jkl3+lrcru
            kl2=kl2+lrcru
  330     continue
  340   continue
      endif
  350 continue
c
c  i transformation
c
  410 if(igueq1) return
c
      if(esfc) ij3=-nc2+1
      if(esfcij) ikl3i=0
      ijkl4=0
      do 550 ircr=1,ircru
        jkl3j=0
        if(esfcij) jrcru=ircr
        do 540 jrcr=1,jrcru
          jkl3=jkl3j
          if(esfcij) ikl3=ikl3i
          if(esfc) krcru=ircr
          do 520 krcr=1,krcru
            if(esfc) then
              if(esfckl) then
                if(krcr.eq.ircr) then
                  lrcru=jrcr
                else
                  lrcru=krcr
                endif
                do 450 lrcr=1,lrcru
                  g4(ijkl4+lrcr)=g4(ijkl4+lrcr)+
     1                           g3(jkl3+lrcr)*etai(ircr,ig)+
     2                           g3(ikl3+lrcr)*etai(jrcr,ig)+
     3                           g3(ij3+nc2*lrcr)*etai(krcr,ig)+
     4                           g3(ij3+nc2*krcr)*etai(lrcr,ig)
  450           continue
                ijkl4=ijkl4+lrcru
                jkl3=jkl3+krcr
                ikl3=ikl3+krcr
              else
                if(krcr.eq.ircr) then
                  lrcru=jrcr
                else
                  lrcru=jrcru
                endif
                do 470 lrcr=1,lrcru
                  g4(ijkl4+lrcr)=g4(ijkl4+lrcr)+
     1                           g3(jkl3+lrcr)*etai(ircr,ig)+
     2                           g3(ij3+nc2*lrcr)*etai(krcr,ig)
  470           continue
                ijkl4=ijkl4+lrcru
                jkl3=jkl3+jrcru
              endif
            else
              if(esfckl) lrcru=krcr
              if(esfcij) then
                do 490 lrcr=1,lrcru
                  g4(ijkl4+lrcr)=g4(ijkl4+lrcr)+
     1                           g3(jkl3+lrcr)*etai(ircr,ig)+
     2                           g3(ikl3+lrcr)*etai(jrcr,ig)
  490           continue
                ijkl4=ijkl4+lrcru
                jkl3=jkl3+lrcru
                ikl3=ikl3+lrcru
              else
                do 510 lrcr=1,lrcru
                  g4(ijkl4+lrcr)=g4(ijkl4+lrcr)+
     1                           g3(jkl3+lrcr)*etai(ircr,ig)
  510           continue
                ijkl4=ijkl4+lrcru
                jkl3=jkl3+lrcru
              endif
            endif
  520     continue
          jkl3j=jkl3j+nc2
          if(esfc) ij3=ij3+1
  540   continue
        if(esfcij) ikl3i=ikl3i+nc2
  550 continue
      return
      end

max_stars_repo_name: MOLFDIR/MOLFDIR

max_issues_repo_name: yingxingcheng/MOLFDIR

max_forks_repo_name: yingxingcheng/MOLFDIR


text:
C
C**************************  VCOPYI ************************************
C
C  copy integer vector IVEC2 into integer vector IVEC1
C
C***********************************************************************
C
      SUBROUTINE VCOPYI(NUM,IVEC1,IVEC2)
C
      IMPLICIT  NONE
      INTEGER   I
      INTEGER   NUM
      INTEGER   IVEC1(*),IVEC2(*)
C
      DO I=1,NUM
         IVEC1(I)=IVEC2(I)
      END DO
C
      RETURN
      END


text:
      SUBROUTINE DLBPT2 (INPUT,W1JK,W2JK)
C
      INTEGER         W1JK,W2JK,SYSBUF,ECORE,TW1JK,TW2JK,NAME(2)
      DIMENSION       A(4),IZ(1)
      COMMON /PACKX / ITI,ITO,II,NN,INCR
      COMMON /AMGP2 / TW1JK(7),TW2JK(7)
      COMMON /AMGMN / MCB(7),NROW,ND,NE,REFC,FMACH,RFK
      COMMON /SYSTEM/ SYSBUF
      COMMON /ZZZZZZ/ NJ1,NK1,NP,NB,NTP,NBZ,NBY,NTZ,NTY,NT0,NTZS,NTYS,
     1                Z(1)
      EQUIVALENCE     (IZ(1),Z(1))
      DATA    NAME  / 4HDLBP,4HT2  /
C
C     GET CORE THEN SET POINTERS TO ACPT TABLE ARRAYS
C
      ECORE = KORSZ(IZ) - 4*SYSBUF
C
C     READ LENGTHS OF ARRAYS
C
      CALL FREAD (INPUT,NJ1,13,0)
C
C     COMPUTE POINTERS TO OPEN CORE
C
      IF (NTP .EQ. 0) CALL FREAD (INPUT,0,0,1)
      IF (NTP .EQ. 0) GO TO 50
      LNS  = IZ(1)
      INC  = 1
      INS  = INC
      INB  = INS + NP
      INAS = INB + NP
      IZIN = INAS
      IYIN = IZIN
      INBEA1 = IYIN   + NP
      INBEA2 = INBEA1 + NB
      INSBEA = INBEA2 + NB
      IZB  = INSBEA + NB
      IYB  = IZB  + NB
      IAVR = IYB  + NB
      IARB = IAVR + NB
      INFL = IARB + NB
      IXLE = INFL + NB
      IXTE = IXLE + NB
      INT121 = IXTE   + NB
      INT122 = INT121 + NB
      IZS = INT122 + NB
      N = 3*NP + 12*NB
C
C     READ FIXED ARRAYS
C
      IF (N .GT. ECORE) GO TO 180
      CALL FREAD (INPUT,IZ,N,0)
C
C     GET LENGTHS OF VARIABLE ARRAYS, PANELS THEN BODIES
C
      LNAS = 0
      IF (NP .EQ. 0) GO TO 20
      DO 10 I = 1,NP
   10 LNAS = LNAS + IZ(INAS+I-1)
   20 LNB  = 0
      LNSB = 0
      LNFL = 0
      LT1  = 0
      LT2  = 0
      DO 30 I = 1,NB
      K    = I - 1
      LNB  = LNB  + IZ(INBEA1+K)
      LNSB = LNSB + IZ(INSBEA+K)
      LNFL = LNFL + IZ(INFL+K)
      LT1  = LT1  + IZ(INT121+K)
   30 LT2  = LT2  + IZ(INT122+K)
C
C     READ VARIABLE  ARRAYS AND SET POINTERS TO CORE
C
      NEXT = N + 1
      N = 2*NB + 5*LNS + 4*NTP + 3*LNB + 4*LNSB + LNAS + 2*LNFL
     *  + LT1  + LT2
      IF (NEXT+N .GE. ECORE) GO TO 180
      CALL READ (*190,*190,INPUT,IZ(NEXT),N,1,NW)
      NEXT = NEXT+ N + 1
      IYS  = IZS + NB + LNS
      ICS  = IYS
      IEE  = ICS + NB + LNS
      ISG  = IEE + LNS
      ICG  = ISG + LNS
      IXIJ = ICG
      IX   = IXIJ+ LNS
      IDELX= IX  + NTP + LNB
C
C     COMPUTE TERMS AND PACK
C
      NN = II + 1
      DO 40 I = 1,NTP
      A(1) = 0.0
      A(2) = 1.0
      CALL PACK (A,W1JK,TW1JK)
      A(1) = -(2.0/REFC)
      A(2) = Z(IDELX+I-1)/(2.0*REFC)
      CALL PACK (A,W2JK,TW2JK)
C
C     BUMP PACK INDEXES
C
      II = II + 2
      IF (I .EQ. NTP) GO TO 40
      NN = NN + 2
   40 CONTINUE
   50 NTZY = NTZ + NTY
      IF (NTZY .EQ. 0) GO TO 70
      NN   = II + 1
      A(1) = 0.0
      A(2) = 0.0
      DO 60 I = 1,NTZY
      CALL PACK (A,W1JK,TW1JK)
      CALL PACK (A,W2JK,TW2JK)
   60 CONTINUE
   70 NTZY = NTZS + NTYS
      IF (NTZY .EQ. 0) GO TO 200
C
C     ANOTHER HARDER SHUFFLE
C
      III = II
      INBEA2 = INBEA2 - 1
      INSBEA = INSBEA - 1
      IFY = II
      IF (NBZ .EQ. 0) GO TO 120
      DO 110 I = 1,NBZ
      IBT = IZ(INBEA2+I)
      NBE = IZ(INSBEA+I)
      IF (IBT .EQ. 2) GO TO 90
      A(1) = 0.0
      A(2) = 1.0
      A(3) = -2.0/REFC
      A(4) = 0.0
      DO 80 J = 1,NBE
      NN = II + 1
      CALL PACK (A,W1JK,TW1JK)
      CALL PACK (A(3),W2JK,TW2JK)
      II  = II + 2
      IFY = II
   80 CONTINUE
      GO TO 110
   90 A(1) = 0.0
      A(4) = 0.0
      DO 100 J = 1,NBE
      NN   = II + 3
      A(2) = 0.0
      A(3) = 1.0
      CALL PACK (A,W1JK,TW1JK)
      A(2) = -2.0/REFC
      A(3) = 0.0
      CALL PACK (A,W2JK,TW2JK)
      II = II + 4
  100 CONTINUE
  110 CONTINUE
  120 IF (NBY .EQ. 0) GO TO 170
      II = IFY
      NBTD = NB - NBY + 1
      DO 160 I = NBTD,NB
      IBT = IZ(INBEA2+I)
      NBE = IZ(INSBEA+I)
      IF (IBT .EQ. 3) GO TO 140
      A(2) = 0.0
      A(3) = 0.0
      DO 130 J = 1,NBE
      NN   = II + 3
      A(1) = 0.0
      A(4) =-1.0
      CALL PACK (A,W1JK,TW1JK)
      A(1) = -2.0/REFC
      A(4) = 0.0
      CALL PACK (A,W2JK,TW2JK)
      II  = II + 4
  130 CONTINUE
      GO TO 160
  140 A(1) = 0.0
      A(2) =-1.0
      A(3) =-2.0/REFC
      A(4) = 0.0
      DO 150 J = 1,NBE
      NN = II + 1
      CALL PACK (A,W1JK,TW1JK)
      CALL PACK (A(3),W2JK,TW2JK)
      II = II + 2
  150 CONTINUE
  160 CONTINUE
  170 II = III + NTZY*2
      NN = II  - 1
      GO TO 200
C
C     ERROR MESSAGES
C
  180 CALL MESAGE (-8,0,NAME)
  190 CALL MESAGE (-7,0,NAME)
  200 RETURN
      END


! params.f90
!
!>  The model parameters module.
!
!> @copyright
!> 2015 Lesley De Cruz & Jonathan Demaeyer.
!> See LICENSE.txt for license information.
!
!---------------------------------------------------------------------------
!
!>  @remark
!>  Once the init_params() subroutine is called, the parameters are loaded
!>  globally in the main program and its subroutines and function
!
!---------------------------------------------------------------------------

MODULE params

  IMPLICIT NONE

  PUBLIC

  REAL(KIND=8) :: n         !< \f$n = 2 L_y / L_x\f$ - Aspect ratio
  REAL(KIND=8) :: phi0      !< Latitude in radian
  REAL(KIND=8) :: rra       !< Earth radius
  REAL(KIND=8) :: sig0      !< \f$\sigma_0\f$ - Non-dimensional static stability of the atmosphere.
  REAL(KIND=8) :: k         !< Bottom atmospheric friction coefficient.
  REAL(KIND=8) :: kp        !< \f$k'\f$ - Internal atmospheric friction coefficient.
  REAL(KIND=8) :: r         !< Frictional coefficient at the bottom of the ocean.
  REAL(KIND=8) :: d         !< Merchanical coupling parameter between the ocean and the atmosphere.
  REAL(KIND=8) :: f0        !< \f$f_0\f$ - Coriolis parameter
  REAL(KIND=8) :: gp        !< \f$g'\f$Reduced gravity
  REAL(KIND=8) :: H         !< Depth of the active water layer of the ocean.
  REAL(KIND=8) :: phi0_npi  !< Latitude exprimed in fraction of pi.

  REAL(KIND=8) :: lambda    !< \f$\lambda\f$ - Sensible + turbulent heat exchange between the ocean and the atmosphere.
  REAL(KIND=8) :: Co        !< \f$C_a\f$ - Constant short-wave radiation of the ocean.
  REAL(KIND=8) :: Go        !< \f$\gamma_o\f$ - Specific heat capacity of the ocean.
  REAL(KIND=8) :: Ca        !< \f$C_a\f$ - Constant short-wave radiation of the atmosphere.
  REAL(KIND=8) :: To0       !< \f$T_o^0\f$ -  Stationary solution for the 0-th order ocean temperature.
  REAL(KIND=8) :: Ta0       !< \f$T_a^0\f$ -  Stationary solution for the 0-th order atmospheric temperature.
  REAL(KIND=8) :: epsa      !< \f$\epsilon_a\f$ - Emissivity coefficient for the grey-body atmosphere.
  REAL(KIND=8) :: Ga        !< \f$\gamma_a\f$ - Specific heat capacity of the atmosphere.
  REAL(KIND=8) :: RR        !< \f$R\f$ - Gas constant of dry air

  REAL(KIND=8) :: scale     !< \f$L_y = L \, \pi\f$ - The characteristic space scale.
  REAL(KIND=8) :: pi        !< \f$\pi\f$
  REAL(KIND=8) :: LR        !< \f$L_R\f$ - Rossby deformation radius
  REAL(KIND=8) :: G         !< \f$\gamma\f$
  REAL(KIND=8) :: rp        !< \f$r'\f$ - Frictional coefficient at the bottom of the ocean.
  REAL(KIND=8) :: dp        !< \f$d'\f$ - Non-dimensional mechanical coupling parameter between the ocean and the atmosphere.
  REAL(KIND=8) :: kd        !< \f$k_d\f$ - Non-dimensional bottom atmospheric friction coefficient.
  REAL(KIND=8) :: kdp       !< \f$k'_d\f$ - Non-dimensional internal atmospheric friction coefficient.

  REAL(KIND=8) :: Cpo       !< \f$C'_a\f$ - Non-dimensional constant short-wave radiation of the ocean.
  REAL(KIND=8) :: Lpo       !< \f$\lambda'_o\f$ - Non-dimensional sensible + turbulent heat exchange from ocean to atmosphere.
  REAL(KIND=8) :: Cpa       !< \f$C'_a\f$ - Non-dimensional constant short-wave radiation of the atmosphere. @remark Cpa acts on psi1-psi3, not on theta.
  REAL(KIND=8) :: Lpa       !< \f$\lambda'_a\f$ - Non-dimensional sensible + turbulent heat exchange from atmosphere to ocean.
  REAL(KIND=8) :: sBpo      !< \f$\sigma'_{B,o}\f$ - Long wave radiation lost by ocean to atmosphere & space.
  REAL(KIND=8) :: sBpa      !< \f$\sigma'_{B,a}\f$ - Long wave radiation from atmosphere absorbed by ocean.
  REAL(KIND=8) :: LSBpo     !< \f$S'_{B,o}\f$ - Long wave radiation from ocean absorbed by atmosphere.
  REAL(KIND=8) :: LSBpa     !< \f$S'_{B,a}\f$ - Long wave radiation lost by atmosphere to space & ocean.
  REAL(KIND=8) :: L         !< \f$L\f$ - Domain length scale
  REAL(KIND=8) :: sc        !< Ratio of surface to atmosphere temperature.
  REAL(KIND=8) :: sB        !< Stefan–Boltzmann constant
  REAL(KIND=8) :: betp      !< \f$\beta'\f$ - Non-dimensional beta parameter

  REAL(KIND=8) :: nua=0.D0  !< Dissipation in the atmosphere
  REAL(KIND=8) :: nuo=0.D0  !< Dissipation in the ocean

  REAL(KIND=8) :: nuap      !< Non-dimensional dissipation in the atmosphere
  REAL(KIND=8) :: nuop      !< Non-dimensional dissipation in the ocean

  REAL(KIND=8) :: t_trans   !< Transient time period
  REAL(KIND=8) :: t_run     !< Effective intergration time (length of the generated trajectory)
  REAL(KIND=8) :: dt        !< Integration time step
  REAL(KIND=8) :: tw        !< Write all variables every tw time units
  REAL(KIND=8) :: tw_solo
  REAL(KIND=8) :: t_run_da
  REAL(KIND=8) :: tw_da
  REAL(KIND=8) :: infl
  REAL(KIND=8) :: ini_err
  INTEGER :: ens_num
  INTEGER :: nobs

  LOGICAL :: writeout       !< Write to file boolean

  INTEGER :: nboc   !< Number of atmospheric blocks
  INTEGER :: nbatm  !< Number of oceanic blocks
  INTEGER :: natm=0 !< Number of atmospheric basis functions
  INTEGER :: noc=0  !< Number of oceanic basis functions
  INTEGER :: ndim   !< Number of variables (dimension of the model)
  INTEGER, DIMENSION(:,:), ALLOCATABLE :: oms   !< Ocean mode selection array
  INTEGER, DIMENSION(:,:), ALLOCATABLE :: ams   !< Atmospheric mode selection array

  PRIVATE :: init_nml

CONTAINS

  !> Read the basic parameters and mode selection from the namelist.
  SUBROUTINE init_nml
    INTEGER :: AllocStat

    NAMELIST /aoscale/  scale,f0,n,rra,phi0_npi
    NAMELIST /oparams/  gp,r,H,d,nuo
    NAMELIST /aparams/  k,kp,sig0,nua
    NAMELIST /toparams/ Go,Co,To0
    NAMELIST /taparams/ Ga,Ca,epsa,Ta0
    NAMELIST /otparams/ sc,lambda,RR,sB

    NAMELIST /modeselection/ oms,ams
    NAMELIST /numblocs/ nboc,nbatm

    NAMELIST /int_params/ t_trans,t_run,dt,tw,writeout
    NAMELIST /int_params_solo/ tw_solo
    NAMELIST /da_params/ t_run_da, tw_da, nobs, infl
    NAMELIST /ens_params/ ens_num, ini_err

    OPEN(8, file="params.nml", status='OLD', recl=80, delim='APOSTROPHE')

    READ(8,nml=aoscale)
    READ(8,nml=oparams)
    READ(8,nml=aparams)
    READ(8,nml=toparams)
    READ(8,nml=taparams)
    READ(8,nml=otparams)

    CLOSE(8)

    OPEN(8, file="modeselection.nml", status='OLD', recl=80, delim='APOSTROPHE')
    READ(8,nml=numblocs)

    ALLOCATE(oms(nboc,2),ams(nbatm,2), STAT=AllocStat)
    IF (AllocStat /= 0) STOP "*** Not enough memory ! ***"

    READ(8,nml=modeselection)
    CLOSE(8)

    OPEN(8, file="int_params.nml", status='OLD', recl=80, delim='APOSTROPHE')
    READ(8,nml=int_params)
    READ(8,nml=int_params_solo)

    CLOSE(8)

    OPEN(8, file="da_params.nml", status='OLD', recl=80, delim='APOSTROPHE')
    READ(8,nml=da_params)
    READ(8,nml=ens_params)

    CLOSE(8)

  END SUBROUTINE init_nml

  !> Parameters initialisation routine
  SUBROUTINE init_params
    INTEGER, DIMENSION(2) :: s
    INTEGER :: i
    CALL init_nml

    !---------------------------------------------------------!
    !                                                         !
    ! Computation of the dimension of the atmospheric         !
    ! and oceanic components                                  !
    !                                                         !
    !---------------------------------------------------------!

    natm=0
    DO i=1,nbatm
       IF (ams(i,1)==1) THEN
          natm=natm+3
       ELSE
          natm=natm+2
       ENDIF
    ENDDO
    s=shape(oms)
    noc=s(1)

    ndim=2*natm+2*noc

    !---------------------------------------------------------!
    !                                                         !
    ! Some general parameters (Domain, beta, gamma, coupling) !
    !                                                         !
    !---------------------------------------------------------!

    pi=dacos(-1.D0)
    L=scale/pi
    phi0=phi0_npi*pi
    LR=sqrt(gp*H)/f0
    G=-L**2/LR**2
    betp=L/rra*cos(phi0)/sin(phi0)
    rp=r/f0
    dp=d/f0
    kd=k*2
    kdp=kp

    !-----------------------------------------------------!
    !                                                     !
    ! DERIVED QUANTITIES                                  !
    !                                                     !
    !-----------------------------------------------------!

    Cpo=Co/(Go*f0) * RR/(f0**2*L**2)
    Lpo=lambda/(Go*f0)
    Cpa=Ca/(Ga*f0) * RR/(f0**2*L**2)/2 ! Cpa acts on psi1-psi3, not on theta
    Lpa=lambda/(Ga*f0)
    sBpo=4*sB*To0**3/(Go*f0) ! long wave radiation lost by ocean to atmosphere space
    sBpa=8*epsa*sB*Ta0**3/(Go*f0) ! long wave radiation from atmosphere absorbed by ocean
    LSBpo=2*epsa*sB*To0**3/(Ga*f0) ! long wave radiation from ocean absorbed by atmosphere
    LSBpa=8*epsa*sB*Ta0**3/(Ga*f0) ! long wave radiation lost by atmosphere to space & ocea
    nuap=nua/(f0*L**2)
    nuop=nuo/(f0*L**2)

  END SUBROUTINE init_params
END MODULE params


text:
!! R832 if-construct
!    is  if-then-stmt
!              block
!          [ else-if-smt
!              block ] ...
!          [ else-stmt
!              block ]
!         end-if-stmt
!
! Also tests rules R841 (if-then-stmt), R842 (else-if-stmt), R843 (else-stmt),
! R844 (end-if-stmt), and if-construct-name.  Does not test
! scalar-logical-expr.

10 if(x .eq. y) then
   x = 1
11 else if(x .lt. y) then
   x = 0
12 else
   x = -1
13 end if

if(x .eq. y) then
   x = 1
else
   x = -1
end if

if(x .eq. y) then
   x = 1
end if

my_if: if(x .eq. y) then
   x = 0
else if(x .lt. y) then my_if
   if (.false.) x = 1
else if(x .gt. y) then
   x = 2
   y = 13
else my_if
end if my_if

end

max_stars_repo_name: ouankou/rose

max_issues_repo_name: WildeGeist/rose

max_forks_repo_name: WildeGeist/rose


text:
      INTEGER FUNCTION IMRKOPR(NEMO)

C$$$  SUBPROGRAM DOCUMENTATION BLOCK
C
C SUBPROGRAM:    IMRKOPR
C   PRGMMR: J. ATOR          ORG: NCEP       DATE: 2016-05-04
C
C ABSTRACT: THIS FUNCTION DETERMINES WHETHER THE GIVEN MNEMONIC
C   CONTAINS A TABLE C MARKER OPERATOR.
C
C PROGRAM HISTORY LOG:
C 2016-05-04  J. ATOR    -- ORIGINAL AUTHOR
C
C USAGE:    IMRKOPR (NEMO)
C   INPUT ARGUMENT LIST:
C     NEMO     - CHARACTER*(*): MNEMONIC
C
C   OUTPUT ARGUMENT LIST:
C     IMRKOPR  - INTEGER: RETURN CODE INDICATING WHETHER NEMO CONTAINS
C                A TABLE C MARKER OPERATOR
C                  0 - NO
C                  1 - YES
C
C REMARKS:
C    THIS ROUTINE CALLS:        None
C    THIS ROUTINE IS CALLED BY: IGETRFEL IOKOPER  STSEQ
C                               Normally not called by any application
C                               programs.
C
C ATTRIBUTES:
C   LANGUAGE: FORTRAN 77
C   MACHINE:  PORTABLE TO ALL PLATFORMS
C
C$$$

      INCLUDE 'bufrlib.prm'

      CHARACTER*(*)  NEMO

C-----------------------------------------------------------------------
C-----------------------------------------------------------------------

      IF (LEN(NEMO).LT.6) THEN
        IMRKOPR = 0
      ELSE IF ( ( NEMO(4:6).EQ.'255' )
     +			.AND.
     +	       ( ( NEMO(1:3).EQ.'223' ) .OR. ( NEMO(1:3).EQ.'224' ) .OR.
     +	         ( NEMO(1:3).EQ.'225' ) .OR. ( NEMO(1:3).EQ.'232' ) ) )
     +			THEN
        IMRKOPR = 1
      ELSE
        IMRKOPR = 0
      ENDIF

      RETURN
      END