lecdon_khione.f

Go to the documentation of this file.

!                   ************************
                    SUBROUTINE lecdon_khione
!                   ************************
!
     & (file_desc,path,ncar)
!
!***********************************************************************
! KHIONE   V7P2                                              02/11/2016
!***********************************************************************
!
!brief    READS THE STEERING FILE THROUGH A DAMOCLES CALL.
!
!history  F. SOUILLE (EDF)
!+        30/09/2019
!+        V8P0
!+        Updated variables and sections
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| FILE_DESC      |-->| STORES THE FILES 'SUBMIT' ATTRIBUTES
!|                |   | IN DICTIONARIES. IT IS FILLED BY DAMOCLES.
!| NCAR           |-->| LENGTH OF PATH
!| PATH           |-->| NAME OF CURRENT DIRECTORY
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE declarations_special
      USE declarations_khione
      USE meteo_telemac
      USE freezeup_khione, ONLY : buoyancy_velocity
!
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      CHARACTER(LEN=PATH_LEN), INTENT(INOUT) :: FILE_DESC(4,maxkeyword)
      INTEGER, INTENT(IN)                    :: NCAR
      CHARACTER(LEN=PATH_LEN), INTENT(IN)    :: PATH
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      CHARACTER(LEN=8) ::   MNEMO(maxvar)
      INTEGER          ::   K,I
      DOUBLE PRECISION ::   ANG
      CHARACTER(LEN=2) CHAR2
!
      CHARACTER(LEN=PATH_LEN) :: NOM_CAS
      CHARACTER(LEN=PATH_LEN) :: NOM_DIC
!
!-----------------------------------------------------------------------
!
! ARRAYS USED IN THE DAMOCLES CALL
!
      INTEGER            ADRESS(4,maxkeyword)
      INTEGER            DIMENS(4,maxkeyword)
      DOUBLE PRECISION   MOTREA(maxkeyword)
      INTEGER            MOTINT(maxkeyword)
      LOGICAL            MOTLOG(maxkeyword)
      CHARACTER(LEN=PATH_LEN) MOTCAR(maxkeyword)
      CHARACTER(LEN=72)  MOTCLE(4,maxkeyword,2)
      INTEGER            TROUVE(4,maxkeyword)
      LOGICAL            DOC
      INTEGER :: ID_DICO, ID_CAS
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      IF (lng.EQ.1) WRITE(lu,1)
      IF (lng.EQ.2) WRITE(lu,2)
1     FORMAT(1x,/,19x, '********************************************',/,
     &            19x, '*     SOUS-PROGRAMME LECDON_KHIONE         *',/,
     &            19x, '*           APPEL DE DAMOCLES              *',/,
     &            19x, '*     VERIFICATION DES DONNEES LUES        *',/,
     &            19x, '*           SUR LE FICHIER CAS             *',/,
     &            19x, '********************************************',/)
2     FORMAT(1x,/,19x, '********************************************',/,
     &            19x, '*        SUBROUTINE LECDON_KHIONE          *',/,
     &            19x, '*           CALL OF DAMOCLES               *',/,
     &            19x, '*        VERIFICATION OF READ DATA         *',/,
     &            19x, '*            ON STEERING FILE              *',/,
     &            19x, '********************************************',/)
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
!     DAMOCLE CALL - PARSING DICO AND CAS FILES
!
!     __________________________________________________________________
!     INITIALISES THE VARIABLES FOR DAMOCLES CALL
      DO k=1,maxkeyword
!       A FILENAME NOT GIVEN BY DAMOCLES WILL BE RECOGNIZED AS A WHITE SPACE
!       (IT MAY BE THAT NOT ALL COMPILERS WILL INITIALISE LIKE THAT)
        motcar(k)(1:1)=' '
!
        dimens(1,k) = 0
        dimens(2,k) = 0
        dimens(3,k) = 0
        dimens(4,k) = 0
      ENDDO
!
!     __________________________________________________________________
!     WRITES OUT INFO
      doc = .false.
!
!     __________________________________________________________________
!     OPENS DICTIONNARY AND STEERING FILES
!
      IF(ncar.GT.0) THEN
!
        nom_dic=path(1:ncar)//'ICEDICO'
        nom_cas=path(1:ncar)//'ICECAS'
!
      ELSE
!
        nom_dic='ICEDICO'
        nom_cas='ICECAS'
!
      ENDIF
!
      CALL get_free_id(id_dico)
      OPEN(id_dico,file=nom_dic,form='FORMATTED',action='READ')
      CALL get_free_id(id_cas)
      OPEN(id_cas,file=nom_cas,form='FORMATTED',action='READ')
!
!     __________________________________________________________________
!     CALL DAMOCLE
      CALL damocle
     &( adress, dimens, maxkeyword, doc, lng, lu, motint,
     &  motrea, motlog, motcar, motcle , trouve, id_dico, id_cas,
     &  .false.,file_desc)
!
!     __________________________________________________________________
!     CLOSES DICTIONNARY AND STEERING FILES
      CLOSE(id_dico)
      CLOSE(id_cas)
!
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
!     RETRIEVES FILE NUMBERS FROM KHIONE FORTRAN PARAMETERS
!
!     DECODES 'SUBMIT' CHAINS
      CALL read_submit(ice_files,maxlu_ice,file_desc,maxkeyword)
!
!     AT THIS LEVEL LOGICAL UNITS ARE EQUAL TO THE FILE NUMBER
      DO i=1,maxlu_ice
        IF    (ice_files(i)%TELNAME.EQ.'ICECLI') THEN
          icecli=i
        ELSEIF(ice_files(i)%TELNAME.EQ.'ICEGEO') THEN
          icegeo=i
        ELSEIF(ice_files(i)%TELNAME.EQ.'ICEREF') THEN
          iceref=i
        ELSEIF(ice_files(i)%TELNAME.EQ.'ICERES') THEN
          iceres=i
        ELSEIF(ice_files(i)%TELNAME.EQ.'ICECOV') THEN
          icecov=i
        ELSEIF(ice_files(i)%TELNAME.EQ.'ICEBLK') THEN
          iceblk=i
        ELSEIF(ice_files(i)%TELNAME.EQ.'CLGRFO') THEN
          clgrfo=i
        ENDIF
      ENDDO
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
!     ALGORITHMIC DIFFERENTIATION
!
      nadvar    = motint( adress(1,13) )
!
!     NAMES OF DIFFERENTIATED VARIABLES
      n_names_advar = dimens(4,13)
      nadvar = max( nadvar,n_names_advar ) ! WARNING HERE ?
      IF(nadvar.GT.0) THEN
        DO i=1,nadvar
          WRITE(char2,'(I2)') i
          names_advar(i) =  'DERIVATIVE '//adjustl(char2)//'   '
     &                   // '??              '
        ENDDO
      ENDIF
      IF(n_names_advar.GT.0) THEN
        DO k=1,n_names_advar
          names_advar(k) = motcar(adress(4,13)+k-1)(1:32)
        ENDDO
      ENDIF
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
!     LISTING AND GRAPHICAL OUTPUTS
!     __________________________________________________________________
!     TITLE OF THE STUDY
      titicecas = motcar( adress(4, 3) )(1:72)
!     __________________________________________________________________
!     GRAPHIC PRINTOUT PERIOD
      leoprd    = motint( adress(1,  1) )
!     __________________________________________________________________
!     LISTING PRINTOUT PERIOD
      lisprd    = motint( adress(1,  2) )
!     __________________________________________________________________
!     FILES IN THE STEERING FILE
!
      ice_files(icecli)%NAME=motcar( adress(4,4 ) )
      ice_files(icegeo)%NAME=motcar( adress(4,5 ) )
      ice_files(icegeo)%FMT=motcar( adress(4,35) )(1:8)
      ice_files(iceref)%NAME=motcar( adress(4,6) )
      ice_files(iceref)%FMT=motcar( adress(4,36) )(1:8)
      ice_files(iceres)%NAME=motcar( adress(4,7 ) )
      ice_files(iceres)%FMT=motcar( adress(4,37) )(1:8)
      ! PREVIOUS ICE COVER FILE
      ice_files(icecov)%NAME=motcar( adress(4,8) )
      ice_files(icecov)%FMT=motcar( adress(4,14) )(1:8)
      CALL majus(ice_files(icecov)%FMT)
      ! PREVIOUS ICE BLOCKS FILE
      ice_files(iceblk)%NAME=motcar( adress(4,15) )
      ice_files(iceblk)%FMT=motcar( adress(4,16) )(1:8)
      CALL majus(ice_files(icecov)%FMT)
      ! CLOGGING RESULTS FILE
      ice_files(clgrfo)%NAME=motcar( adress(4,25 ) )
!     __________________________________________________________________
!     UPDATE OF MNEMO AND NOMVARS
      DO i=1,maxvar
        mnemo(i) = '        '
      ENDDO
      CALL nomvar_khione( texte,textpr,mnemo,nadvar,names_advar )
!     __________________________________________________________________
!     GRAPHICAL OUTPUT VARIABLES
      vardes = motcar( adress(4, 10) )(1:72)
      CALL majus( vardes )
      CALL sortie( vardes , mnemo , maxvar , sorleo )
      varimp = motcar( adress(4, 11) )(1:72)
      CALL majus( varimp )
      CALL sortie( varimp , mnemo , maxvar , sorimp )
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
!     PHYSICAL PRPERTIES AND CONSTANTS RELATED KEYWORDS
!     __________________________________________________________________
!     WATER KINEMATIC VISCOSITY
      xnu = motrea( adress(2,8) )
!     __________________________________________________________________
!     DENSITIES
      ro0       = motrea( adress(2,  1) ) !WATER
      rho_air   = motrea( adress(2, 11) )
      rho_ice   = motrea( adress(2, 12) )
!     __________________________________________________________________
!     SPECIFIC AND LATENT HEATS
      cp_eau    = motrea( adress(2, 5) )
      cp_ice    = motrea( adress(2, 16) )
      lh_ice    = motrea( adress(2, 17) )
!     __________________________________________________________________
!     ICE ALBEDO
      albe = motrea( adress(2,34) )
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
!     ATMOSPHERIC PROPERTY RELATED KEYWORDS (IN METEO_TELEMAC)
!
      atmoexch  = motint( adress(1,17) )
!     __________________________________________________________________
!     VERTICAL DATUM
      windz     = motrea( adress(2,22) )
      modelz    = motrea( adress(2,23) )
!     __________________________________________________________________
!     SUN RISE AND SUN SET ANGLES
      alphsd    = motrea( adress(2,25) )     !? sun exit angle
      alphrd    = motrea( adress(2,26) )     !? sun emission angle
!     __________________________________________________________________
!     SOLAR CONSTANT
      sio       = motrea( adress(2,33) )     !? solar constant
!     __________________________________________________________________
!     HEAT EXCHANGE CALIBRATION COEFFICIENTS
!     (ATMOSPHERE-WATER EXCHANGE MODEL=3)
      lin_watair = motrea( adress(2, 18) )   ! replaces HWA
      cst_watair = motrea( adress(2, 20) )   ! replaces ALPW
      lin_iceair = motrea( adress(2, 19) )   ! replaces HIA
      cst_iceair = motrea( adress(2, 21) )   ! replaces ALP
!     __________________________________________________________________
!     HEAT EXCHANGE CALIBRATION COEFFICIENTS
      coef_phib = motrea( adress(2, 56) )
      coef_phie = motrea( adress(2, 57) )
      coef_phih = motrea( adress(2, 58) )
      coef_phip = motrea( adress(2, 59) )
!     __________________________________________________________________
!     METEOROLOGY, IN CASE ABSENT FROM METEO FILES
!     CST_TAIR  = MOTREA( ADRESS(2, 41) )
      cst_tdew  = motrea( adress(2, 42) )
!     CST_CLDC  = MOTREA( ADRESS(2, 43) )
      cst_visbi = motrea( adress(2, 45) )
!     CST_RAINFALL = MOTREA( ADRESS(2, 46) )
!     IF( DIMENS(2,47).EQ.2 ) THEN
!       CST_WINDS = MOTREA( ADRESS(2,47) + 0 )
!       CST_WINDD = MOTREA( ADRESS(2,47) + 1 )
!       DTRS4 = ATAN(1.D0) / 45.D0
!       CST_WINDX = - CST_WINDS * SIN( CST_WINDD*DTRS4 )
!       CST_WINDY = - CST_WINDS * COS( CST_WINDD*DTRS4 )
!     ENDIF
!     IF( DIMENS(2,48).EQ.2 ) THEN
!       CST_WINDX = MOTREA( ADRESS(2,48) + 0 )
!       CST_WINDY = MOTREA( ADRESS(2,48) + 1 )
!     ENDIF
!     CST_PATMOS  = MOTREA( ADRESS(2, 49) )
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
!     WATER PROPERTY RELATED KEYWORDS
!
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
!     FREEZEUP RELATED KEYWORDS
!
      de    = motrea( adress(2,28) ) ! DIAMETER THICKNESS RATIO
      nuss  = motrea( adress(2,24) ) ! CONSTANT NUSSELT NUMBER
      nussi = motint( adress(1, 9) ) ! MODEL FOR NUSSLET NUMBER COMPUTATION
      ibuoy = motint( adress(1,10) ) ! MODEL FOR BUOYANCY VELOCITY
      isnuc = motint( adress(1,12) ) ! MODEL FOR SECONDARY NUCLEATION
      snnmax= motrea( adress(2,13) ) ! NMAX PARAM FOR SECONDARY NUCLEATION
      ifloc = motint( adress(1,15) ) ! MODEL FOR FLOCULATION AND BREAKUP
      afloc = motrea( adress(2,14) ) ! AFLOC PARAM FOR FLOCULATION
      inrjb = motint( adress(1,18) ) ! ENERGY BALANCE VERSION
      iseed = motint( adress(1,19) ) ! SEEDING MODEL
      minnk = motint( adress(1,20) ) ! MINIMUM NUMBER OF FRAZIL PARTICLES PER UNIT VOLUME
      seedr = motrea( adress(2,40) ) ! FRAZIL SEEDING RATE
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
!     ICE COVER PROPERTY RELATED KEYWORDS
!
      ifrot = motint(adress(1,8))       ! ICE COVER FRICTION LAW
      ifice = motint(adress(1,16))      ! LAW FOR FRICTION COEF
      fice = motrea(adress(2,4))        ! ICE COVER FRICTION COEFFICIENT
      fice_max = motrea(adress(2,52))   ! ICE COVER MAX FRICTION COEFFICIENT
      thie = motrea( adress(2,75) )     ! EQUIVALENT SURFACE ICE THCIKNESS
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
!     STATIC BORDER ICE RELATED KEYWORDS
!
!
      bch = motrea( adress(2,98) )       ! (15) CHANNEL WIDTH FOR THE COMPUTATION OF SURFACE WATER TEMPERATURE
      tc = motrea( adress(2,71) )        ! (-1.1) CRITICAL WATER SURFACE TEMPERATURE FOR BORDER ICE FORMATION
      vcrbor = motrea(adress(2,50))      ! (0.07) THRESHOLD VELOCITY FOR STATIC BORDER ICE
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
!     DYNAMIC BORDER ICE RELATED KEYWORDS
!
      vcrbom = motrea(adress(2,73))      ! (0.4) THRESHOLD VELOCITY ABOVE WHICH DYNAMIC BORDER ICE WILL NOT FORM
!
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
!     ICE PROPERTY RELATED KEYWORDS
!
!*******************************
!     REAL KEYWORDS            *
!*******************************
!
      alsm = motrea( adress(2,30) )     ! standard longitude, in degrees
      allm = motrea( adress(2,31) )     ! local longitude, in degrees
      etadir = motrea( adress(2,32) )   ! - for west longitude, + for east longitude
      af   = motrea( adress(2,29) )      ! (1.) DEPOSITION COEFFICIENT OF FRAZIL ON THE BAR
      surf_ef  = motrea( adress(2,39) )  ! POROSITY OF SURFACE ICE
      tc_wt = motrea( adress(2,44) )     ! => XKWP: (0.56594) WATER-ICE THERMAL CONDUCTIVITY
!
      clog_ef    = motrea( adress(2,37) )! POROSITY OF FRAZIL ICE
      clog_theta = motrea( adress(2,36) )
      IF( dimens(2,35).EQ.4 ) THEN
        clog_tdist = motrea( adress(2,35) + 0 )
        clog_tdiam = motrea( adress(2,35) + 1 )
        clog_vdist = motrea( adress(2,35) + 2 )
        clog_vdiam = motrea( adress(2,35) + 3 )
      ENDIF
!
      nfrclog = dimens(1,14)
      IF(modulo(dimens(1,3),2).EQ.1) THEN
        IF ( motint( adress(1,3) ).EQ.0 ) THEN
          nseclog = 0
        ELSE
          WRITE(lu,*) ' '
          WRITE(lu,*) 'LECDON_KHIONE'
          WRITE(lu,*) 'CLOOGING SECTION WORKS BY PAIR'
          WRITE(lu,*) 'PLEASE GIVE A MULTIPLE OF 2 NUMBER OF NODES'
          CALL plante(1)
          stop
        ENDIF
      ELSE
        nseclog = dimens(1,3)/2
      ENDIF
      ntotclog = nfrclog + nseclog
      IF( motint( adress(1,14) ).EQ.0 ) nfrclog = 0
      IF( ntotclog.GT.0 ) THEN
!
        ALLOCATE( lines%CELLS(ntotclog) )
        lines%NVAL = ntotclog
        IF(nfrclog.GT.0) ALLOCATE(numclog(nfrclog))
        IF(nseclog.GT.0) ALLOCATE(seclog(dimens(1,3)))
        IF(nseclog.GT.0) ALLOCATE(un1(nseclog))
        IF(nseclog.GT.0) ALLOCATE(un2(nseclog))
        ALLOCATE(clog_tlgth(ntotclog))
        ALLOCATE(clog_twdth(ntotclog))
        ALLOCATE(clog_vlgth(ntotclog))
        ALLOCATE(clog_vwdth(ntotclog))
        ALLOCATE(clog_volum(ntotclog))
!
!       ACUMULATION ANGLE
        ang = 4.d0*atan(1.d0) * clog_theta / 180.0
        DO k = 1,nfrclog
          numclog(k) = motint( adress(1,14) + k-1 )
        ENDDO
        DO k = 1,nseclog*2
          seclog(k) = motint( adress(1,3) + k-1 )
        ENDDO
        DO k = 1,ntotclog
          lines%CELLS(k)%NVAL = 0
          clog_tlgth(k) = 0.d0
          clog_twdth(k) = 2.0*clog_tdiam*cos(ang)*sin(ang)
          clog_vlgth(k) = 0.d0
          clog_vwdth(k) = 2.0*clog_vdiam*cos(ang)*sin(ang)
          clog_volum(k) = 0.d0
        ENDDO
!
      ENDIF
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
!     ICE PROCESS
!
      cwi1 = motrea( adress(2,62) )       ! CONSTANT FOR HEAT TRANSFER BETWEEN TURBULENT WATER AND ICE
      ciw1 = motrea( adress(2,63) )       ! CONSTANT FOR HEAT TRANSFER FOR SUPERCOOLED TURBULENT FLOW
      ata = motrea( adress(2,64) )        ! NUSSELT NUMBER FOR HEAT TRANSFER BETWEEN LAMINAR WATER AND ICE
      iturb = motint( adress(1,4) )       ! TURBULENCE MODEL FOR K AND EPSILON APPROXIMATION
      tc_bi = motrea( adress(2,65) )      ! THERMAL CONDUCTIVITY OF BLACK ICE (W/M/OC) => XKI
!
      tc_s = motrea( adress(2,67) )       ! THERMAL CONDUCTIVITY OF SNOW(W/M/OC) => XKS
      sgma = motrea( adress(2,69) )       ! BOLTZMANN CONSTANT (WM-2K-4)
      cst_tmelt = motrea( adress(2,70) )  ! FREEZING POINT OF WATER
!
      anfem0 = motrea( adress(2,74) )     ! (1.) MAXIMUM CONCENTRATION FOR BORDER ICE FORMATION
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
!     FREEZEUP PARAMETERS INIT
!
      ALLOCATE(rk_frzl(nc_fra))
      ALLOCATE(ek_frzl(nc_fra))
      ALLOCATE(vk_frzl(nc_fra))
      ALLOCATE(nusst(nc_fra))
      ALLOCATE(vbk(nc_fra))
!
      IF(nc_fra.GT.0) THEN
        DO k=1,nc_fra
          IF(dimens(2, 27).GE.k) THEN
            rk_frzl(k) = motrea( adress(2, 27) + k-1 )
          ELSE
            IF(k.EQ.1) THEN
              rk_frzl(k) = 1.0e-4
            ELSE
              rk_frzl(k) = rk_frzl(1) + (k-1)*1.0e-4
            ENDIF
            WRITE(lu,*) 'WARNING: THE VALUE OF FRAZIL CRISTAL RADIUS ',k
            WRITE(lu,*) ' IS NOT GIVEN, RADIUS IS SET TO ', rk_frzl(k)
          ENDIF
          nusst(k) = nuss
          pi = 4.d0*atan(1.d0)
          ek_frzl(k) = 2.d0*rk_frzl(k)/de
          vk_frzl(k) = 2.d0*pi*rk_frzl(k)**3/de
          vbk(k) = buoyancy_velocity(rk_frzl(k), ek_frzl(k))
        ENDDO
      ENDIF
!
!-----------------------------------------------------------------------
!
      RETURN
      END SUBROUTINE