borice.f

Go to the documentation of this file.

!                     *****************
                      SUBROUTINE borice
!                     *****************
!
     &(h,u,v,f,at,lt,dt,tra05,tra06,liubor,nptfr,numliq,klog,mask,mesh
     & ,s)
!
!***********************************************************************
! KHIONE   V7P3
!***********************************************************************
!
!brief    INCLUDE ICE PROCESSES TO THE BOUNDARIES SETTINGS.
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| AT             |-->| TIME OF THE DATASET
!| DT             |-->| TIME STEP
!| F              |-->| TRACER VALUES
!| H              |-->| WATER DEPTH
!| KLOG           |-->| CONVENTION FOR SOLID BOUNDARY
!| LIUBOR         |-->| TYPE OF BOUNDARY CONDITION FOR THE VELOCITY
!| LT             |-->| NUMBER OF TIME STEPS
!| MASK           |-->| BLOCK OF MASKS FOR DIFFERENT BOUNDARY CONDITIONS
!| MESH           |-->| MESH STRUCTURE
!| NPTFR          |-->| NUMBER OF BOUNDARY POINT
!| NUMLIQ         |-->| LIQUID BOUNDARY NUMBER
!| S              |-->| VOID STRUCTURE
!| TRA05          |-->| WORK ARRAY IN A BIEF_OBJ STRUCTURE
!| TRA06          |-->| WORK ARRAY IN A BIEF_OBJ STRUCTURE
!| U              |-->| X-COMPONENT OF THE VELOCITY
!| V              |-->| Y-COMPONENT OF THE VELOCITY
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE bief
      USE declarations_khione, ONLY:
     &         titicecas,ice_files,clgrfo,leoprd,
     &         lines,clogging,bd_ice,ind_fra,t1,t2,
     &         nfrclog,numclog,clog_ef,clog_theta,
     &         clog_vdist,clog_vdiam,clog_vlgth,
     &         clog_tdist,clog_tdiam,clog_tlgth,
     &         icetype,clog_twdth,clog_vwdth,
     &         rho_ice,clog_ef,clog_volum,iceloc,t3,nseclog,ntotclog
      USE freezeup_khione, ONLY: clogged_on_bar
      USE declarations_special
      USE interface_parallel, ONLY : p_dsum
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER, INTENT(IN) :: KLOG, NPTFR
!
      INTEGER,          INTENT(IN)    :: NUMLIQ(nptfr)
      INTEGER,          INTENT(IN)    :: LIUBOR(nptfr)
      TYPE(bief_mesh),  INTENT(INOUT) :: MESH
      TYPE(bief_obj),   INTENT(IN)    :: H,U,V,F,S
      TYPE(bief_obj),   INTENT(INOUT) :: MASK, TRA05,TRA06
      INTEGER,          INTENT(IN)    :: LT
      DOUBLE PRECISION, INTENT(IN)    :: AT,DT
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER               K,N,I,IFRLIQ,IELEB,IFILE
      LOGICAL               OUTPUT
      DOUBLE PRECISION PI, ANG,RV0,RV1,RT0,RT1, RATE, QCV,CV,CA,
     &      mfbar,mfvb,mftb, mfgrd,mfvt,mftt,mass0, dw,dv,taw, nv,nt
      DOUBLE PRECISION, ALLOCATABLE :: QF(:),QA(:),QV(:)
      CHARACTER(LEN=250)    CELL
!
!-----------------------------------------------------------------------
!
      INTRINSIC adjustl,trim
!
!=======================================================================
!
!     5 - CLOGGING
!
!     NO IMPACT ON THE HYDRODYNAMICS FOR NOW
!     WARNING: ONLY WORKS FOR MONO-CLASS FRAZIL ICE
!
      IF( clogging ) THEN
!
        IF(nseclog.GT.0) THEN
          ALLOCATE(qf(nseclog))
          ALLOCATE(qa(nseclog))
          ALLOCATE(qv(nseclog))
          CALL section_khione(mesh,mesh%IKLE%I,mesh%NELMAX,
     &                        mesh%IFABOR%I,f,u,v,qf,qa,qv,h,s,lt)
        ENDIF
!
        IF( lt.LE.1 .AND. ntotclog.GT.0 ) THEN
!
          IF(nfrclog.GT.0) THEN
!
!-----------------------------------------------------------------------
!         COMPUTES THE LENGTH OF THE OPEN BOUNDARY (CLOG_TLGTH) ONCE
            CALL os('X=C     ', x=t3, c=1.d0)
            DO ifrliq = 1,nfrclog
              n = numclog(ifrliq)
              CALL os( 'X=0     ', x=tra05 )
              DO ieleb = 1,mesh%NELEB
                k = mesh%IKLBOR%I(ieleb)
                IF( numliq(k).EQ.n )
     &            tra05%R(ieleb) = mask%ADR(8)%P%R(ieleb)
              ENDDO
              CALL vector(tra06,'=','MASVEC          ',
     &          ielbor(11,1),1.d0,t3,t3,t3,t3,t3,t3,mesh,.true.,tra05)
              clog_tlgth(ifrliq) = bief_sum(tra06)
              IF( ncsize.GT.1 )
     &          clog_tlgth(ifrliq) = p_dsum(clog_tlgth(ifrliq))
            ENDDO
          ENDIF
!
!-----------------------------------------------------------------------
!         HEADER OF ASCII CLOGGING FILE
!
          IF( ice_files(clgrfo)%NAME.NE.' ' ) THEN
            IF(ncsize.EQ.1.OR.ipid.EQ.0) THEN
              ifile = ice_files(clgrfo)%LU
              WRITE(ifile,301) ''
              WRITE(ifile,301) titicecas
              WRITE(ifile,301) ''
 301          FORMAT('# ',a)
              lines%HTXT = 'T'
              WRITE(ifile,300,advance='NO') trim(adjustl(lines%HTXT))
              lines%HUNT = 'S'
              DO ifrliq = 1,ntotclog
                IF( lines%CELLS(ifrliq)%NVAL.EQ.0 ) THEN
                  lines%CELLS(ifrliq)%NVAL = 8
                  ALLOCATE(lines%CELLS(ifrliq)%TXT(8))
                  ALLOCATE(lines%CELLS(ifrliq)%UNT(8))
                  ALLOCATE(lines%CELLS(ifrliq)%VAL(8))
                ENDIF
                WRITE(cell,'(I12)') ifrliq
                lines%CELLS(ifrliq)%TXT(1) = 'FRAZIL('
     &            // trim(adjustl(cell)) // ')'
                lines%CELLS(ifrliq)%UNT(1) = 'SI'
                lines%CELLS(ifrliq)%TXT(2) = 'WIDTH_V('
     &            // trim(adjustl(cell)) // ')'
                lines%CELLS(ifrliq)%UNT(2) = 'M'
                lines%CELLS(ifrliq)%TXT(3) = 'WIDTH_H('
     &            // trim(adjustl(cell)) // ')'
                lines%CELLS(ifrliq)%UNT(3) = 'M'
                lines%CELLS(ifrliq)%TXT(4) = 'CLG_AREA('
     &            // trim(adjustl(cell)) // ')'
                lines%CELLS(ifrliq)%UNT(4) = 'M^2'
                lines%CELLS(ifrliq)%TXT(5) = 'IMASS_V('
     &            // trim(adjustl(cell)) // ')'
                lines%CELLS(ifrliq)%UNT(5) = 'G'
                lines%CELLS(ifrliq)%TXT(6) = 'IMASS_H('
     &            // trim(adjustl(cell)) // ')'
                lines%CELLS(ifrliq)%UNT(6) = 'G'
                lines%CELLS(ifrliq)%TXT(7) = 'CLG_VOL('
     &            // trim(adjustl(cell)) // ')'
                lines%CELLS(ifrliq)%UNT(7) = 'M^3'
                lines%CELLS(ifrliq)%TXT(8) = 'CLG_MAS('
     &            // trim(adjustl(cell)) // ')'
                lines%CELLS(ifrliq)%UNT(8) = 'SI'
                WRITE(ifile,302,advance='NO') ( ' ' //
     &           trim(adjustl(lines%CELLS(ifrliq)%TXT(k))), k=1,8)
              ENDDO
              WRITE(ifile,300) ''
              WRITE(ifile,300,advance='NO') trim(adjustl(lines%HUNT))
              DO i = 1,ntotclog
                WRITE(ifile,302,advance='NO')
     &            ( ' ' // trim(adjustl(lines%CELLS(i)%UNT(k)))
     &            , k=1,8 )
              ENDDO
              WRITE(ifile,300) ''
 300          FORMAT(a)
 302          FORMAT(8a)
            ENDIF
          ENDIF
        ENDIF
!
!-----------------------------------------------------------------------
!       PREPARE COMPUTATION OF THE FLUXES AND THE AVERAGE BAR HEIGHTS
!
!
        CALL os('X=Y     ', x=t2, y=h )
        CALL os('X=YZ    ', x=t1, y=h, z=f%ADR(ind_fra)%P )

!
        pi  = 4.d0*atan(1.d0)
        output = .false.
!
!       ACUMULATION ANGLE
        ang = pi * clog_theta / 180.d0
!       INITIAL FRAZIL RADIUS
        rv0 = clog_vdiam*cos(ang)
        rt0 = clog_tdiam*cos(ang)
!
!       CLOGGING RATES ( AF = 1.D0 )
        rate = 2.d0*sin(ang)*sin(ang)*(1.d0)/( 1.d0-clog_ef )/ang
!
!-----------------------------------------------------------------------
!
        IF(nfrclog.GT.0) THEN
          DO ifrliq = 1,nfrclog
            n = numclog(ifrliq)
!
!           MASKING ONE LIQUID BOUNDARY AT A TIME
            CALL os( 'X=0     ', x=tra05 )
            DO ieleb = 1,mesh%NELEB
              k = mesh%IKLBOR%I(ieleb)
              IF( numliq(k).EQ.n )
     &          tra05%R(ieleb) = mask%ADR(8)%P%R(ieleb)
            ENDDO
!
!           FRAZIL FLUX = INTEGRAL( F * H * UV ) ACROSS BOUNDARY N
            CALL vector(tra06,'=','FLUBDF          ',
     &        ielbor(11,1),1.d0,t1,t1,t1,u,v,v,mesh,.true.,tra05)
            qcv = bief_sum(tra06)
            IF( ncsize.GT.1 ) qcv = p_dsum(qcv)
!           FRAZIL FLUX BY UNIT OF WIDTH
            qcv = abs(qcv) / clog_tlgth(ifrliq)
!
!           CROSS SECTIONAL AREA = INTEGRAL( H ) ACROSS BOUNDARY N
            CALL vector(tra06,'=','MASVEC          ',
     &        ielbor(11,1),1.d0,t2,t2,t2,t2,t2,t2,mesh,.true.,tra05)
            ca = bief_sum(tra06)
            IF( ncsize.GT.1 ) ca = p_dsum(ca)
!           AVERAGE BAR HEIGHT
            IF( clog_tlgth(ifrliq).GT.0.d0 ) THEN
              clog_vlgth(ifrliq) = ca/clog_tlgth(ifrliq)
            ELSE
              clog_vlgth(ifrliq) = 0.d0
            ENDIF
!
!           AVERAGE CONCENTRATION = INTEGRAL( F * H ) ACROSS BOUNDARY N
            CALL vector(tra06,'=','MASVEC          ',
     &        ielbor(11,1),1.d0,t1,t1,t1,t1,t1,t1,mesh,.true.,tra05)
            cv = bief_sum(tra06)
            IF( ncsize.GT.1 ) cv = p_dsum(cv)
!           ... DIVIDED BY THE AREA OF THE CROSS SECTION
            IF( ca.GT.0.d0 ) THEN
              cv = cv/ca
            ELSE
              cv = 0.d0
            ENDIF
!
!           REPRESENTATIVE NUMBERS OF BARS
            IF( clog_vdist*clog_vdiam.GT.0.d0 ) THEN
!             THIS COULD BE AN ODD NUMBER
              nv = clog_tlgth(ifrliq)/clog_vdist
            ELSE
              nv = 0.d0
            ENDIF
            IF( clog_tdist*clog_tdiam.GT.0.d0 ) THEN
!             COUNTIN ONLY THE FULLY SUBMERGED
              nt = 1.d0*int( clog_vlgth(ifrliq)/clog_tdist )
            ELSE
              nt = 0.d0
            ENDIF
!
!           VERTICAL BARS
            IF( clog_vlgth(ifrliq)*nv.GT.0.d0 ) THEN
              clog_vwdth(ifrliq) = clog_vwdth(ifrliq) +
     &                             dt*rate*qcv/clog_vlgth(ifrliq)
              rv1 = clog_vwdth(ifrliq)/2.d0/sin(ang)
              CALL clogged_on_bar( rv0,rv1,
     &          clog_vdiam,clog_vlgth(ifrliq),nv,ang,mfvb,mfvt)
            ELSE
              rv1 = 0.d0
              mfvb = 0.d0
              mfvt = 0.d0
            ENDIF
!
!           TRANSVERSE BARS
            IF( clog_tlgth(ifrliq)*nt.GT.0.d0 ) THEN
              clog_twdth(ifrliq) = clog_twdth(ifrliq) +
     &                             dt*rate*qcv/clog_vlgth(ifrliq)
              rt1 = clog_twdth(ifrliq)/2.d0/sin(ang)
              CALL clogged_on_bar( rt0,rt1,
     &          clog_tdiam,clog_tlgth(ifrliq),nt,ang,mftb,mftt)
            ELSE
              rt1 = 0.d0
              mftb = 0.d0
              mftt = 0.d0
            ENDIF
!
            IF( ice_files(clgrfo)%NAME.EQ.' ' ) cycle
            IF( lt.NE.leoprd*int(lt/leoprd) .AND. lt.GT.1 ) cycle
!
!           ADDITIONAL PRINTOUTS
            output = .true.
!
            IF(ncsize.EQ.1.OR.ipid.EQ.0) THEN
              WRITE(lines%CELLS(ifrliq)%VAL(1),303) cv
              WRITE(lines%CELLS(ifrliq)%VAL(2),303) clog_vwdth(ifrliq)
              WRITE(lines%CELLS(ifrliq)%VAL(3),303) clog_twdth(ifrliq)
!             WIDTH OF THE VERTICAL BARS EXCLUDED FROM THE TRANSVERSE BARS
              dw = clog_vdist - clog_vwdth(ifrliq)
              dv = max( clog_vwdth(ifrliq),clog_vdiam )
!             AVAILABLE AREA FOR THE RACK
              taw = clog_vlgth(ifrliq)*dw*(nv-1.d0)
     &          - clog_twdth(ifrliq) * ( clog_tlgth(ifrliq)-nv*dv ) * nt
              taw = max(taw,0.d0)
              WRITE(lines%CELLS(ifrliq)%VAL(4),303) taw
!             FRAZIL VOLUME FOR THE GRID
              clog_volum(ifrliq) = clog_volum(ifrliq) +
     &          qcv * ( 1.d0-clog_ef ) * dt *
     &          ( clog_vwdth(ifrliq)*nv + clog_twdth(ifrliq)*nt )
!             FRAZIL MASS FOR THE GRID
              mass0 = clog_volum(ifrliq)*rho_ice
!
              mfbar = mfvb + mftb
              WRITE(lines%CELLS(ifrliq)%VAL(5),303) mfbar
              mfgrd = mfvt + mftt
              WRITE(lines%CELLS(ifrliq)%VAL(6),303) mfgrd
              WRITE(lines%CELLS(ifrliq)%VAL(7),303) clog_volum(ifrliq)
              WRITE(lines%CELLS(ifrliq)%VAL(8),303) mass0
            ENDIF
!
          ENDDO
        ENDIF
!
        IF(nseclog.GT.0) THEN
          DO i=nfrclog+1,ntotclog
            clog_vlgth(i) = qa(i-nfrclog)/clog_tlgth(i)
            qcv = qf(i-nfrclog)
            IF(qa(i-nfrclog).GT.0.d0) THEN
              cv = qv(i-nfrclog) / qa(i-nfrclog)
            ELSE
              cv=0.d0
            ENDIF
!           REPRESENTATIVE NUMBERS OF BARS
            IF( clog_vdist*clog_vdiam.GT.0.d0 ) THEN
!             THIS COULD BE AN ODD NUMBER
              nv = clog_tlgth(i)/clog_vdist
            ELSE
              nv = 0.d0
            ENDIF
            IF( clog_tdist*clog_tdiam.GT.0.d0 ) THEN
!             COUNTIN ONLY THE FULLY SUBMERGED
              nt = 1.d0*int( clog_vlgth(i)/clog_tdist )
            ELSE
              nt = 0.d0
            ENDIF
!
!           VERTICAL BARS
            IF( clog_vlgth(i)*nv.GT.0.d0 ) THEN
              clog_vwdth(i) = clog_vwdth(i) +
     &                             dt*rate*qcv/clog_vlgth(i)
              rv1 = clog_vwdth(i)/2.d0/sin(ang)
              CALL clogged_on_bar( rv0,rv1,
     &          clog_vdiam,clog_vlgth(i),nv,ang,mfvb,mfvt)
            ELSE
              rv1 = 0.d0
              mfvb = 0.d0
              mfvt = 0.d0
            ENDIF
!
!           TRANSVERSE BARS
            IF( clog_tlgth(i)*nt.GT.0.d0 ) THEN
              clog_twdth(i) = clog_twdth(i) +
     &                             dt*rate*qcv/clog_vlgth(i)
              rt1 = clog_twdth(i)/2.d0/sin(ang)
              CALL clogged_on_bar( rt0,rt1,
     &          clog_tdiam,clog_tlgth(i),nt,ang,mftb,mftt)
            ELSE
              rt1 = 0.d0
              mftb = 0.d0
              mftt = 0.d0
            ENDIF
!
            IF( ice_files(clgrfo)%NAME.EQ.' ' ) cycle
            IF( lt.NE.leoprd*int(lt/leoprd) .AND. lt.GT.1 ) cycle
!
!           ADDITIONAL PRINTOUTS
            output = .true.
!
            IF(ncsize.EQ.1.OR.ipid.EQ.0) THEN
              WRITE(lines%CELLS(i)%VAL(1),303) cv
              WRITE(lines%CELLS(i)%VAL(2),303) clog_vwdth(i)
              WRITE(lines%CELLS(i)%VAL(3),303) clog_twdth(i)
!             WIDTH OF THE VERTICAL BARS EXCLUDED FROM THE TRANSVERSE BARS
              dw = clog_vdist - clog_vwdth(i)
              dv = max( clog_vwdth(i),clog_vdiam )
!             AVAILABLE AREA FOR THE RACK
              taw = clog_vlgth(i)*dw*(nv-1.d0)
     &          - clog_twdth(i) * ( clog_tlgth(i)-nv*dv ) * nt
              taw = max(taw,0.d0)
              WRITE(lines%CELLS(i)%VAL(4),303) taw
!             FRAZIL VOLUME FOR THE GRID
              clog_volum(i) = clog_volum(i) +
     &          qcv * ( 1.d0-clog_ef ) * dt *
     &          ( clog_vwdth(i)*nv + clog_twdth(i)*nt )
!             FRAZIL MASS FOR THE GRID
              mass0 = clog_volum(i)*rho_ice
!
              mfbar = mfvb + mftb
              WRITE(lines%CELLS(i)%VAL(5),303) mfbar
              mfgrd = mfvt + mftt
              WRITE(lines%CELLS(i)%VAL(6),303) mfgrd
              WRITE(lines%CELLS(i)%VAL(7),303) clog_volum(i)
              WRITE(lines%CELLS(i)%VAL(8),303) mass0
            ENDIF

          ENDDO
        ENDIF
!
!-----------------------------------------------------------------------
!
 303    FORMAT(g16.9e3)
        IF( ice_files(clgrfo)%NAME.NE.' ' .AND. output ) THEN
          IF(ncsize.EQ.1.OR.ipid.EQ.0) THEN
            WRITE(lines%HVAL,303) at
            WRITE(ice_files(clgrfo)%LU,302,advance='NO')
     &        trim(adjustl(lines%HVAL))
            DO i=1,ntotclog
              WRITE(ice_files(clgrfo)%LU,302,advance='NO')
     &          ( ' ' // trim(adjustl(lines%CELLS(i)%VAL(k)))
     &          , k=1,8 )
            ENDDO
            WRITE(ice_files(clgrfo)%LU,300) ''
          ENDIF
        ENDIF
!
      ENDIF
!
!=======================================================================
!
!     7 - STATIC BORDER ICE TRACKING
!
      IF( bd_ice ) THEN
!
        IF( lt.LE.1 ) THEN
!
          DO k = 1,nptfr
            IF( liubor(k).EQ.klog ) THEN
              n = mesh%NBOR%I(k)
              IF( icetype%I(n).EQ.1 ) iceloc%I(n) = 2
            ENDIF
          ENDDO
!
        ENDIF
!
      ENDIF
!
!=======================================================================
!
!     ? - OTHER PROCESSES
!
      IF(nseclog.GT.0) THEN
        DEALLOCATE(qf)
        DEALLOCATE(qa)
        DEALLOCATE(qv)
      ENDIF
!
!-----------------------------------------------------------------------
!
      RETURN
      END SUBROUTINE