TELEMAC-MASCARET Forum: Operations on structures

Welcome, Guest

TOPIC: Operations on structures - tracer decay

Operations on structures - tracer decay 7 years 7 months ago #25798

kingja.x

Hi all,

I'm trying to modify the decay processes in TELEMAC but am having trouble with a formulation which relates decay to the depth of water. I've tried to base this on the example modification in difsou.f and although it compiles without error the model crashes on the first timestep and returns the error:

Program received signal SIGSEGV : Segmentation fault - invalid memory reference

I'm pretty sure this is due to the way in which I modify the implicit source term in difsou.f (extract below)

CALL OS('X=Y    ',X=TIMP%ADR(ITRAC)%P)
DO I=1,HPROP%DIM1
    DKi = ALFA*IRR*(1.D0-EXP(EXT*HPROP%R(I)))/(EXT*HPROP%R(I))
    DEK = DKi + DKt
    TIMP%ADR(ITRAC)%P%R(I)=-DEK*86400.D0*HPROP%R(I)
ENDDO

The error does not occur when lines 1 and 5 are commented out however I'm not sure of the correct way to allocate the implicit source term. I've looked through the guide for programming but can't work out what I need to do (I'm relatively new to Fortran). I want to do something like this (below) but this also does not work.

DO I=1,HPROP%DIM1
    DKi = ALFA*IRR*(1.D0-EXP(EXT*HPROP%R(I)))/(EXT*HPROP%R(I))
    DEK = DKi + DKt
    CALL OS('X=CY    ',X=TIMP%ADR(ITRAC)%P%R(I), Y=HPROP%R(I),
    &       C=-DEK*86400.D0)
ENDDO

Any help would be really appreciated.

Thanks
Jonathan

The administrator has disabled public write access.

Operations on structures - tracer decay 7 years 7 months ago #25802

pham OFFLINE Administrator Posts: 1559 Thank you received: 602	Hi Jonathan, Have you written YASMI(ITRAC) = .TRUE. somewhere? Chi-Tuan
	The administrator has disabled public write access.

Operations on structures - tracer decay 7 years 7 months ago #25809

kingja.x

Hi Chi-Tuan,

Yes, here is the full routine for implicit source terms:

!
!-----------------------------------------------------------------------
!
!     IMPLICIT SOURCE TERMS (DEPENDING ON THE LAW CHOSEN)
!
      DO ITRAC=1,NTRA
        IF(LOITRAC(ITRAC).EQ.0) THEN
          YASMI(ITRAC)=.FALSE.
        ELSEIF(LOITRAC(ITRAC).EQ.1) THEN
          YASMI(ITRAC)=.TRUE.
          CALL OS('X=CY    ',X=TIMP%ADR(ITRAC)%P,Y=HPROP,
     &            C=-2.3D0/COEF1TRAC(ITRAC)/3600.D0)
!     IMPLEMENT MANCINI (1978) / CHAPRA (1997)
        ELSEIF(LOITRAC(ITRAC).EQ.2) THEN
          YASMI(ITRAC)=.TRUE.
          IF ((FRACD > SUNRISE) .AND. (FRACD < SUNSET)) THEN
              IRR = 260.D0
!             CONVERT TO LANGLEY/HR
              IRR = 11.63D0*IRR
          ELSE
              IRR = 0.D0
          END IF
          EXT = 1.8D0 / SS
          DKs = 0.8D0 + 0.02D0*SAL
          DKt = DKs*(1.07D0**(TEMP-20D0))
	  CALL OS('X=Y    ',X=TIMP%ADR(ITRAC)%P)
	  DO I=1,HPROP%DIM1
	      DKi = ALFA*IRR*(1.D0-EXP(EXT*HPROP%R(I)))/(EXT*HPROP%R(I))
              DEK = DKi + DKt
	      TIMP%ADR(ITRAC)%P%R(I)=-DEK*86400.D0*HPROP%R(I)
          ENDDO
        ELSE
          IF(LNG.EQ.1) WRITE(LU,*) 'DIFSOU : LOI NON PROGRAMMEE'
          IF(LNG.EQ.2) WRITE(LU,*) 'DIFSOU : LAW NOT IMPLEMENTED'
        ENDIF
      ENDDO
!
!                                   N+1
!     EXAMPLE WHERE WE ADD -0.0001 F      IN THE RIGHT HAND-SIDEs
!     OF THE TRACER EQUATION THAT BEGINS WITH DF/DT=...
!     (T12=SMI WILL BE DIVIDED BY HPROP IN CVDFTR, THE EQUATION IS:
!     DT/DT=...+SMI*T(N+1)/H
!
!     HERE THIS IS DONE FOR TRACER 3 ONLY IN A RECTANGULAR ZONE
!
!     CALL OS('X=0     ',X=TIMP%ADR(3)%P)
!     DO I=1,HPROP%DIM1
!       IF(X(I).GE.263277.D0.AND.X(I).LE.265037.D0) THEN
!       IF(Y(I).GE.379007.D0.AND.Y(I).LE.380326.D0) THEN
!         TIMP%ADR(3)%P%R(I)=-0.00001D0*HPROP%R(I)
!       ENDIF
!       ENDIF
!     ENDDO
!     YASMI(3)=.TRUE.
!
!-----------------------------------------------------------------------

Thanks
Jonathan

The administrator has disabled public write access.

Operations on structures - tracer decay 7 years 7 months ago #25811

c.coulet OFFLINE Moderator Posts: 3722 Thank you received: 1031	Doesn't look in detail but the division by HPROP could be sometimes dangerous if your model has tidal flat! Regard
	Christophe The administrator has disabled public write access. The following user(s) said Thank You: kingja.x

Operations on structures - tracer decay 7 years 7 months ago #25958

kingja.x	Thank you, I think that was the problem as my model does have a large area of tidal flats. I've included a condition to prevent a division by any depth below a threshold and this has solved the issue.
	The administrator has disabled public write access.