mod_action_ex Module Reference

Functions/Subroutines
subroutine	action_allo
subroutine	action_deallo
subroutine	algorithm_fct (CAS, IG, DTW, IDCMIN, IDCMAX)
subroutine	algorithm_crank_nicolson (CAD, IG, DTW, IDCMIN, IDCMAX, DD)
subroutine	tridag (A, B, C, R, U, N)
subroutine	adv_n (DTW)
subroutine	adv_hl (FF1, I, I1, I2)
subroutine	adv_hl_fij (N32_TMPP3, N32_TMPP2, I, IA, IB)
subroutine	actt2 (IA, IB, RATIO)

Variables
real, dimension(:,:), allocatable	n31
real, dimension(:,:,:), allocatable	n32

Function/Subroutine Documentation

action_allo()

subroutine mod_action_ex::action_allo

(

)

Definition at line 27 of file mod_action_ex.f90.

   
  USE all_vars, ONLY : mt
  USE swcomm3, ONLY : mdc,msc
!  USE MOD_USGRID, ONLY : MDC,MSC
  IMPLICIT NONE
   
  ALLOCATE(n31(mdc,msc))   ;  n31 = 0.0
  ALLOCATE(n32(mdc,msc,0:mt));  n32 = 0.0
   
  RETURN

action_deallo()

subroutine mod_action_ex::action_deallo

(

)

Definition at line 44 of file mod_action_ex.f90.

   
  IMPLICIT NONE
   
  DEALLOCATE(n31)
  DEALLOCATE(n32)
   
  RETURN

actt2()

subroutine mod_action_ex::actt2	(	integer	IA,
		integer	IB,
		real, dimension(msc)	RATIO
	)

Definition at line 979 of file mod_action_ex.f90.

! This subroutine is for Ice induced wave attenuation 
! The parameter 'RATIO' is the dimentional wave attenuation coefficient
! yzhang
!
!****************************************************************
!
  USE m_genarr
  USE swcomm2
  USE swcomm3
  USE swcomm4
  USE ocpcomm4
  USE m_parall
  USE all_vars
  USE vars_wave
  USE ice_state

!  USE ice_state  , only :aice , vice!, floe_diam
!  USE ICE_THERM_VERTICAL ,only:hicen_old
!  USE ice_model_size ,only: ncat

  IMPLICIT NONE 
  REAL    :: C_ICE,COES(7,9),COES2(7,9),ICETHI,SPCSIGL,DEPLOC,TWAVE
  REAL    :: ICEATT,ICEATTS(7),ICEATT1,ICEATT2,RATE,RATIO(MSC)
  REAL    :: mean_diam,mean_diam_u,mean_diam_d,floe_diam
  REAL ,parameter   :: ee=2,ff=0.9,dmin=20
  REAL(SP):: DEP3(MT),DIS,DTMP_DP,X1_DP,X2_DP,Y1_DP,Y2_DP
  INTEGER :: I,ISS,ID,ni,mm,IA,IB
 
  dep3(1:mt) = compda(1:mt,jdp2)
  coes2(1,:)=(/0.0,0.0,0.0,0.0,-0.003509,0.1473,-2.121,11.24,-22.42/)
  coes2(2,:)=(/0.0,0.0,0.0,0.0,-0.003723,0.173,-2.833,18.33,-43.89/)
  coes2(3,:)=(/0.0,0.0,0.0,0.0,-0.001264,0.07156,-1.364,9.616,-25.34/)
  coes2(4,:)=(/0.0,0.0,0.0,0.0,0.001401,-0.05012,0.5915,-3.329,5.372/)
  coes2(5,:)=(/0.0,0.0,0.0,0.0,0.0009987,-0.0403,0.5441,-3.471,6.658/)
  coes2(6,:)=(/0.0,0.0,0.0,0.0,0.0002636,-0.01267,0.1853,-1.403,2.459/)
  coes2(7,:)=(/0.0,0.0,0.0,0.0,1.551e-5,-0.002103,0.03268,-0.4304,0.3751/)

  c_ice=(aice(ia,1)+aice(ib,1))/2
  floe_diam=200
  dis=((vy(ia)-vy(ib))**2+(vx(ia)-vx(ib))**2)**0.5*2.0
  if(c_ice>0.00001)then
    mean_diam_u=0.0
    mean_diam_d=0.0
    mean_diam=0.0
    DO mm=1,3
      mean_diam_u=mean_diam_u+(ee**2*ff)**mm*ee**(-mm)*(floe_diam+floe_diam)/2
      mean_diam_d=mean_diam_d+(ee**2*ff)**mm
    END DO
    mean_diam=mean_diam_u/mean_diam_d
    if(c_ice.eq.1)then
      icethi=3.0
    else
      icethi=0.3*(1.0+4*c_ice)
    endif
    DO iss=1,msc
      spcsigl=spcsig(iss)
      deploc=dep3(i)
!      CALL KSCIP1(1,SPCSIGL,DEPLOC,KWAVEL,CGOL,NN,ND)
      twave=1/(spcsigl/pi2_w)
      IF(twave<=16.0)THEN
        IF(twave<=6.0)twave=6.0
          iceatts(1)=coes2(1,1)*twave**8+coes2(1,2)*twave**7+coes2(1,3)*twave**6+coes2(1,4)*twave**5+coes2(1,5)*twave**4+&
          coes2(1,6)*twave**3+coes2(1,7)*twave**2+coes2(1,8)*twave+coes2(1,9)

          iceatts(2)=coes2(2,1)*twave**8+coes2(2,2)*twave**7+coes2(2,3)*twave**6+coes2(2,4)*twave**5+coes2(2,5)*twave**4+&
          coes2(2,6)*twave**3+coes2(2,7)*twave**2+coes2(2,8)*twave+coes2(2,9)

          iceatts(3)=coes2(3,1)*twave**8+coes2(3,2)*twave**7+coes2(3,3)*twave**6+coes2(3,4)*twave**5+coes2(3,5)*twave**4+&
          coes2(3,6)*twave**3+coes2(3,7)*twave**2+coes2(3,8)*twave+coes2(3,9)

          iceatts(4)=coes2(4,1)*twave**8+coes2(4,2)*twave**7+coes2(4,3)*twave**6+coes2(4,4)*twave**5+coes2(4,5)*twave**4+&
          coes2(4,6)*twave**3+coes2(4,7)*twave**2+coes2(4,8)*twave+coes2(4,9)

          iceatts(5)=coes2(5,1)*twave**8+coes2(5,2)*twave**7+coes2(5,3)*twave**6+coes2(5,4)*twave**5+coes2(5,5)*twave**4+&
          coes2(5,6)*twave**3+coes2(5,7)*twave**2+coes2(5,8)*twave+coes2(5,9)

          iceatts(6)=coes2(6,1)*twave**8+coes2(6,2)*twave**7+coes2(6,3)*twave**6+coes2(6,4)*twave**5+coes2(6,5)*twave**4+&
          coes2(6,6)*twave**3+coes2(6,7)*twave**2+coes2(6,8)*twave+coes2(6,9)

          iceatts(7)=coes2(7,1)*twave**8+coes2(7,2)*twave**7+coes2(7,3)*twave**6+coes2(7,4)*twave**5+coes2(7,5)*twave**4+&
          coes2(7,6)*twave**3+coes2(7,7)*twave**2+coes2(7,8)*twave+coes2(7,9)
          IF (icethi<=0.4)THEN
            iceatt=iceatts(1)
          ELSE IF(icethi==0.6)THEN
            iceatt=iceatts(2)
          ELSE IF(icethi==0.8)THEN
            iceatt=iceatts(3)
          ELSE IF(icethi==1.2)THEN
            iceatt=iceatts(4)
          ELSE IF(icethi==1.6)THEN
            iceatt=iceatts(5)
          ELSE IF(icethi==2.4)THEN
            iceatt=iceatts(6)
          ELSE IF(icethi>=3.2)THEN
            iceatt=iceatts(7)
          ELSE IF(icethi>0.4.AND.icethi<0.6)THEN
            iceatt1=iceatts(1)
            iceatt2=iceatts(2)
            rate=(icethi-0.4)/0.2
            iceatt=iceatt1*(1-rate)+iceatt2*rate
          ELSE IF(icethi>0.6.AND.icethi<0.8)THEN
            iceatt1=iceatts(2)
            iceatt2=iceatts(3)
            rate=(icethi-0.6)/0.2
            iceatt=iceatt1*(1-rate)+iceatt2*rate
          ELSEIF(icethi>0.8.AND.icethi<1.2)THEN
            iceatt1=iceatts(3)
            iceatt2=iceatts(4)
            rate=(icethi-0.8)/0.4
            iceatt=iceatt1*(1-rate)+iceatt2*rate
          ELSE IF(icethi>1.2.AND.icethi<1.6)THEN
            iceatt1=iceatts(4)
            iceatt2=iceatts(5)
            rate=(icethi-1.2)/0.4
            iceatt=iceatt1*(1-rate)+iceatt2*rate
          ELSE IF(icethi>1.6.AND.icethi<2.4)THEN
            iceatt1=iceatts(5)
            iceatt2=iceatts(6)
            rate=(icethi-1.6)/0.8
            iceatt=iceatt1*(1-rate)+iceatt2*rate
          ELSE IF(icethi>2.4.AND.icethi<3.2)THEN
            iceatt1=iceatts(6)
            iceatt2=iceatts(7)
            rate=(icethi-2.4)/0.8
            iceatt=iceatt1*(1-rate)+iceatt2*rate
          ENDIF
          iceatt=exp(iceatt)
          iceatt=iceatt*c_ice/mean_diam
          ratio(iss)=exp(-iceatt*dis)
        ELSE
          iceatt=0.02*exp(-0.386*twave)*c_ice
          ratio(iss)=exp(-iceatt*dis)
        ENDIF
      END DO
    else
      DO iss=1,msc
      ratio(iss)=1.0
    end do
  endif

  RETURN

adv_hl()

subroutine mod_action_ex::adv_hl	(	real(sp), dimension(mdc,msc)	FF1,
		integer	I,
		integer	I1,
		integer	I2
	)

Definition at line 661 of file mod_action_ex.f90.

! This subroutine is for wave advection at high latitude, 
! in order to solve the invalid scalar assumption problem.
! yzhang
!==========YZHANG_NORTHPOLE_TESTING================================
  USE vars_wave
  USE swcomm3
  USE m_genarr
  USE mod_northpole
  
  IMPLICIT NONE
  REAL(SP) :: FF1(MDC,MSC),N32_TMPP(MDC,MSC),N32_TMPP4(MDC,MSC)
  INTEGER  :: ISS, ID,I1,I2,IDT,IDD,III,I,IDD1,IDD2
  REAL     :: UL_DEGREE,CENTER_DEGREE,DL_DEGREE,DIFLAT1,DIFLAT2,RATELAT1,RATELAT2


  DO iss = 1,msc  
    DO id = 1,mdc
      IF(i1==node_northpole)THEN
        ul_degree=22.5
        center_degree=0
        dl_degree=337.5

        IF (vx(i)<ul_degree .OR. vx(i)>dl_degree)THEN
          idt=(mdc*2/8)
          idd=id+idt
!          IDD=ID
          IF(idd>mdc)THEN
            idd=idd-mdc
          END IF
          n32_tmpp(id,iss)=n32(idd,iss,i1)
        END IF
        DO iii=1,7
          center_degree=iii*45.0
          ul_degree=center_degree+22.5
          dl_degree=center_degree-22.5
          IF(vx(i)<ul_degree .AND. vx(i)>dl_degree) THEN
            idt=(mdc*(iii+2)/8)
!            IDT=(MDC*III/8)
            IF(idt>=mdc)THEN
              idt=idt-mdc
            END IF
            idd=id+idt
            IF(idd>mdc)THEN
              idd=idd-mdc
            END IF
            n32_tmpp(id,iss)=n32(idd,iss,i1)
          END IF
        END DO
        ff1(id,iss)=0.5*(n32_tmpp(id,iss)+n32(id,iss,i2))
      ELSE IF(i2==node_northpole)THEN
        ul_degree=22.5
        center_degree=0
        dl_degree=337.5
        IF (vx(i)<ul_degree .OR. vx(i)>dl_degree)THEN
          idt=(mdc*2/8)
          idd=id+idt
!          IDD=ID
          IF(idd>mdc)THEN
            idd=idd-mdc
          END IF
          n32_tmpp(id,iss)=n32(idd,iss,i2)
        END IF
        DO iii=1,7
          center_degree=iii*45.0
          ul_degree=center_degree+22.5
          dl_degree=center_degree-22.5
          IF(vx(i)<ul_degree .AND. vx(i)>dl_degree) THEN
            idt=(mdc*(iii+2)/8)
!            IDT=(MDC*III/8)
            IF(idt>=mdc)THEN
              idt=idt-mdc
            END IF
            idd=id+idt
            IF(idd>mdc)THEN
              idd=idd-mdc
            END IF
            n32_tmpp(id,iss)=n32(idd,iss,i2)
          END IF
        END DO
        ff1(id,iss)=0.5*(n32_tmpp(id,iss)+n32(id,iss,i1))
      ELSE IF(vy(i)>80.AND.i1/=node_northpole.AND.i2/=node_northpole)THEN
        diflat1=vx(i)-vx(i1)
        diflat2=vx(i)-vx(i2)
        IF(diflat1 >  180.0)THEN
          diflat1 = -360.0+diflat1
        ELSE IF(diflat1 < -180.0)THEN
          diflat1 =  360.0+diflat1
        END IF
        IF(diflat2 >  180.0)THEN
          diflat2 = -360.0+diflat2
        ELSE IF(diflat2 < -180.0)THEN
          diflat2 =  360.0+diflat2
        END IF
        ratelat1=diflat1/360.
        ratelat2=diflat2/360.
        idd1=anint(ratelat1*mdc)+id
        idd2=anint(ratelat2*mdc)+id
        IF(idd1<=0)THEN
          idd1=idd1+mdc
        END IF

        IF(idd1>mdc)THEN
          idd1=idd1-mdc
        END IF

        IF(idd2<=0)THEN
          idd2=idd2+mdc
        END IF

        IF(idd2>mdc)THEN
          idd2=idd2-mdc
        END IF
        n32_tmpp(id,iss)=n32(idd1,iss,i1)
        n32_tmpp4(id,iss)=n32(idd2,iss,i2)
        ff1(id,iss)=0.5*(n32_tmpp(id,iss)+n32_tmpp4(id,iss))
!      ELSE
!        FF1(ID,ISS,I)=0.5*(N32(ID,ISS,I1)+N32(ID,ISS,I2))
      END IF
    END DO 
  END DO

  RETURN

adv_hl_fij()

subroutine mod_action_ex::adv_hl_fij	(	real, dimension(mdc,msc)	N32_TMPP3,
		real, dimension(mdc,msc)	N32_TMPP2,
		integer	I,
		integer	IA,
		integer	IB
	)

Definition at line 787 of file mod_action_ex.f90.

! This subroutine is for wave energy trasport at high latitudes
! yzhang
!=========================================================================
  USE vars_wave
  USE swcomm3
  USE m_genarr
  USE mod_northpole

  IMPLICIT NONE 
  REAL     ::FF1(MDC,MSC),N32_TMPP2(MDC,MSC),N32_TMPP3(MDC,MSC)
  INTEGER  :: ISS, ID,IA,IB,IDT,IDD,III,I,IDD1,IDD2
  REAL     :: UL_DEGREE,CENTER_DEGREE,DL_DEGREE,DIFLAT1,DIFLAT2,RATELAT1,RATELAT2

 
  n32_tmpp2=n32(:,:,ia)
  n32_tmpp3=n32(:,:,ib)

  DO iss=1,msc
    DO id=1,mdc
      IF(ia==node_northpole)THEN
        ul_degree=22.5
        center_degree=0
        dl_degree=337.5
        IF (vx(ib)<ul_degree .OR. vx(ib)>dl_degree)THEN
          idt=(mdc*2/8)
          idd=id+idt
!          IDD=ID
          IF(idd>mdc)THEN
            idd=idd-mdc
          END IF
          n32_tmpp2(id,iss)=n32(idd,iss,ia)
        END IF
        DO iii=1,7
          center_degree=iii*45.0
          ul_degree=center_degree+22.5
          dl_degree=center_degree-22.5
          IF(vx(ib)<ul_degree .AND. vx(ib)>dl_degree) THEN
            idt=(mdc*(iii+2)/8)
!            IDT=(MDC*III/8)
            IF(idt>=mdc)THEN
              idt=idt-mdc
            END IF
            idd=id+idt
            IF(idd>mdc)THEN
              idd=idd-mdc
            END IF
            n32_tmpp2(id,iss)=n32(idd,iss,ia)
          END IF
        END DO
      ELSE IF(ib==node_northpole)THEN
        ul_degree=22.5
        center_degree=0
        dl_degree=337.5
        IF (vx(ia)<ul_degree .OR. vx(ia)>dl_degree)THEN
          idt=(mdc*2/8)
          idd=id+idt
!          IDD=ID
          IF(idd>mdc)THEN
            idd=idd-mdc
          END IF
          n32_tmpp2(id,iss)=n32(idd,iss,ib)
        END IF
        DO iii=1,7
          center_degree=iii*45.0
          ul_degree=center_degree+22.5
          dl_degree=center_degree-22.5
          IF(vx(ia)<ul_degree .AND. vx(ia)>dl_degree) THEN
            idt=(mdc*(iii+2)/8)
!            IDT=(MDC*III/8)
            IF(idt>=mdc)THEN
              idt=idt-mdc
            END IF
            idd=id-idt
            IF(idd<1)THEN
              idd=idd+mdc
            END IF
            n32_tmpp2(id,iss)=n32(idd,iss,ia)
          END IF
        END DO
      ELSE IF((vy(ia)>80.OR.vy(ib)>80).AND.ib/=node_northpole.AND.ia/=node_northpole)THEN
        diflat1=vx(ib)-vx(ia)
        IF(diflat1 >  180.0)THEN
          diflat1 = -360.0+diflat1
        ELSE IF(diflat1 < -180.0)THEN
          diflat1 =  360.0+diflat1
        END IF
        ratelat1=diflat1/360.
        idd1=anint(ratelat1*mdc)+id
        IF(idd1<=0)THEN
          idd1=idd1+mdc
        END IF

        IF(idd1>mdc)THEN
          idd1=idd1-mdc
        END IF

        n32_tmpp2(id,iss)=n32(idd1,iss,ia)
!       ELSE
!        N32_TMPP2(ID,ISS,IA)=N32(ID,ISS,IA)
      END IF
      IF(ib==node_northpole)THEN
        ul_degree=22.5
        center_degree=0
        dl_degree=337.5
        IF (vx(ia)<ul_degree .OR. vx(ia)>dl_degree)THEN
          idt=(mdc*2/8)
          idd=id+idt
!          IDD=ID
          IF(idd>mdc)THEN
            idd=idd-mdc
          END IF
          n32_tmpp3(id,iss)=n32(idd,iss,ib)
        END IF
        DO iii=1,7
          center_degree=iii*45.0
          ul_degree=center_degree+22.5
          dl_degree=center_degree-22.5
          IF(vx(ia)<ul_degree .AND. vx(ia)>dl_degree) THEN
            idt=(mdc*(iii+2)/8)
!            IDT=(MDC*III/8)
            IF(idt>=mdc)THEN
              idt=idt-mdc
            END IF
            idd=id+idt
            IF(idd>mdc)THEN
              idd=idd-mdc
            END IF
            n32_tmpp3(id,iss)=n32(idd,iss,ib)
          END IF
        END DO
      ELSE IF(ia==node_northpole)THEN
        ul_degree=22.5
        center_degree=0
        dl_degree=337.5
        IF (vx(ib)<ul_degree .OR. vx(ib)>dl_degree)THEN
          idt=(mdc*2/8)
          idd=id+idt
!          IDD=ID
          IF(idd>mdc)THEN
            idd=idd-mdc
          END IF
          n32_tmpp3(id,iss)=n32(idd,iss,ia)
        END IF
        DO iii=1,7
          center_degree=iii*45.0
          ul_degree=center_degree+22.5
          dl_degree=center_degree-22.5
          IF(vx(ib)<ul_degree .AND. vx(ib)>dl_degree) THEN
            idt=(mdc*(iii+2)/8)
!            IDT=(MDC*III/8)
            IF(idt>=mdc)THEN
              idt=idt-mdc
            END IF
            idd=id-idt
            IF(idd<1)THEN
              idd=idd+mdc
            END IF
            n32_tmpp3(id,iss)=n32(idd,iss,ib)
          END IF
        END DO
      ELSE IF((vy(ia)>80.OR.vy(ib)>80).AND.ib/=node_northpole.AND.ia/=node_northpole)THEN
        diflat1=vx(ia)-vx(ib)
        IF(diflat1 >  180.0)THEN
          diflat1 = -360.0+diflat1
        ELSE IF(diflat1 < -180.0)THEN
          diflat1 =  360.0+diflat1
        END IF
        ratelat1=diflat1/360.
        idd1=anint(ratelat1*mdc)+id
        IF(idd1<=0)THEN
          idd1=idd1+mdc
        END IF

        IF(idd1>mdc)THEN
          idd1=idd1-mdc
        END IF

        n32_tmpp3(id,iss)=n32(idd1,iss,ib)
!       ELSE
!       N32_TMPP3(ID,ISS,IB)=N32(ID,ISS,IB)
      END IF
!      FIJ1A(ID,ISS,IA)=N32(ID,ISS,IA)+DLTXNCVE(I,1)*PSPX(ID,ISS,IA)+DLTYNCVE(I,1)*PSPY(ID,ISS,IA)
!      FIJ2A(ID,ISS,IB)=N32_TMPP2(ID,ISS,IB)+DLTXNCVE(I,2)*PSPX(ID,ISS,IB)+DLTYNCVE(I,2)*PSPY(ID,ISS,IB)
    END DO
  END DO

  RETURN

adv_n()

subroutine mod_action_ex::adv_n

(

real

DTW

)

Definition at line 327 of file mod_action_ex.f90.

!------------------------------------------------------------------------------|

  USE vars_wave
!  USE ALL_VARS
  USE swcomm3
  USE m_genarr

  IMPLICIT NONE
  
  REAL(SP),ALLOCATABLE :: EXFLUX(:,:),UN(:,:),PSPX(:,:,:),PSPY(:,:,:),FF1(:,:),FIJ1(:,:,:)
  REAL(SP),ALLOCATABLE :: XFLUX(:,:,:)
  REAL(SP),ALLOCATABLE :: XFLUX_TMP(:)
  REAL(SP) :: UTMP,VTMP,UN1,X11,Y11,X22,Y22,X33,Y33,TMP1,TMP2,XI,YI
  REAL(SP) :: DXA,DYA,DXB,DYB,FIJ1_TMP
  REAL(SP) :: TXX,TYY,FXX,FYY,VISCOF,TEMP,STPOINT
  REAL(SP) :: FACT,FM1
  INTEGER  :: I,I1,I2,IA,IB,J,J1,J2,K,JTMP,JJ,II
  REAL(SP) :: AC1MIN, AC1MAX, AC2MIN, AC2MAX
  INTEGER  :: ID,ISS,IG,IG2 
  REAL(SP) :: XIN,YIN,XIC,YIC,CANX,CANY
  REAL(SP) :: CAX(MDC,MSC),CAY(MDC,MSC)
!  REAL(SP) :: DTW,RF(MDC,MSC),DF(MDC,MSC)
  REAL(SP) :: RF(MDC,MSC),DF(MDC,MSC)
  REAL     :: DTW
  REAL(SP) :: SPCDIR2(MDC),SPCDIR3(MDC),CANX_2D(MDC,MSC),CANY_2D(MDC,MSC),CGO_1D(MSC),CGO_2D(MDC,MSC)
  INTEGER  :: NTSN_T
  REAL(SP) :: N32_T,XFLUX_T
  REAL(SP),ALLOCATABLE :: DEP2(:),AC2LOC(:)
  REAL(SP) :: DEPLOC,KWAVELOC,CGLOC,NN,ND,SPCSIGL
  REAL(SP) :: UA_NODE,VA_NODE
  INTEGER  :: CNT

  REAL,ALLOCATABLE ::EXFLUXA(:,:),EXFLUXB(:,:),UNA(:,:),UNBB(:,:),FIJ2(:,:,:),&
                     FIJ1A(:,:,:),FIJ1B(:,:,:),FIJ2A(:,:,:),FIJ2B(:,:,:),     &
             N32_TMPP2(:,:),N32_TMPP3(:,:)  
  REAL :: DLTXEA,DLTXEB,DLTXETMP,DLTYETMP,C_ICE
  REAL     :: CANXAB(MDC,MSC),CANYA(MDC,MSC),CANYB(MDC,MSC)
  REAL(SP)   :: XFLUX_1D(0:MT),PSPX_1D(M),PSPY_1D(M)
  REAL   :: RATIO(MSC)
  REAL   :: EXFLUX_TMPA(MDC,MSC),EXFLUX_TMPB(MDC,MSC)
  REAL   :: EXFLUX_TMP(MDC,MSC) 

!------------------------------------------------------------------------------!
  ALLOCATE(dep2(mt),ac2loc(0:mt),xflux_tmp(0:mt))
  ALLOCATE(exflux(mdc,msc),ff1(mdc,msc),fij1(mdc,msc,2),un(mdc,msc))
  ALLOCATE(xflux(mdc,msc,0:mt));xflux=0.0
  ALLOCATE(pspx(mdc,msc,m));pspx=0.0
  ALLOCATE(pspy(mdc,msc,m));pspy=0.0
  IF(high_latitude_wave)THEN 
    ALLOCATE(exfluxa(mdc,msc),exfluxb(mdc,msc),fij2(mdc,msc,0:m),&
             fij1a(mdc,msc,0:m),fij1b(mdc,msc,0:m),fij2a(mdc,msc,0:m),fij2b(mdc,msc,0:m),una(mdc,msc),unbb(mdc,msc))
    ALLOCATE(n32_tmpp2(mdc,msc),n32_tmpp3(mdc,msc))
  END IF
  dep2(1:mt) = compda(1:mt,jdp2)

  cax = 0.0
  cay = 0.0
    
!  DO ISS = 1,MSC
!    DO ID = 1,MDC
!!
!!--Initialize Fluxes-----------------------------------------------------------!
!!
!      XFLUX = 0.0
!      PSPX  = 0.0 
!      PSPY  = 0.0 
!# if defined(PLBC)
!   call replace_N32(N32,ID,ISS)
!# endif
!
  IF(.NOT. high_latitude_wave)THEN
  
  DO i = 1,m
    DO j=1,ntsn(i)-1
      i1=nbsn(i,j)
      i2=nbsn(i,j+1)
!!$          IF(DEP2(I1) <= DEPMIN .AND. DEP2(I2) > DEPMIN)THEN
!!$           FF1=0.5*(N32(ID,ISS,I)+N32(ID,ISS,I2))
!!$          ELSE IF(DEP2(I1) > DEPMIN .AND. DEP2(I2) <= DEPMIN)THEN
!!$           FF1=0.5*(N32(ID,ISS,I1)+N32(ID,ISS,I))
!!$          ELSE IF(DEP2(I1) <= DEPMIN .AND. DEP2(I2) <= DEPMIN)THEN
!!$           FF1=N32(ID,ISS,I)
!!$          ELSE
!!$           FF1=0.5*(N32(ID,ISS,I1)+N32(ID,ISS,I2))
!!$          END IF
      ff1(:,:)=0.5*(n32(:,:,i1)+n32(:,:,i2))
!!$#         if defined (SPHERICAL)
!!$          XTMP  = VX(I2)*TPI-VX(I1)*TPI
!!$       XTMP1 = VX(I2)-VX(I1)
!!$       IF(XTMP1 >  180.0)THEN
!!$         XTMP = -360.0*TPI+XTMP
!!$       ELSE IF(XTMP1 < -180.0)THEN
!!$         XTMP =  360.0*TPI+XTMP
!!$       END IF  
!!$          TXPI=XTMP*COS(DEG2RAD*VY(I))
!!$          TYPI=(VY(I1)-VY(I2))*TPI
!!$          PSPX(I)=PSPX(I)+FF1*TYPI
!!$          PSPY(I)=PSPY(I)+FF1*TXPI
!!$#         else
!!$          PSPX(I)=PSPX(I)+FF1*(VY(I1)-VY(I2))
!!$          PSPY(I)=PSPY(I)+FF1*(VX(I2)-VX(I1))
!!$#         endif
      pspx(:,:,i)=pspx(:,:,i)+ff1(:,:)*dltytrie(i,j)
      pspy(:,:,i)=pspy(:,:,i)+ff1(:,:)*dltxtrie(i,j)
    END DO
    pspx(:,:,i)=pspx(:,:,i)/art2(i)
    pspy(:,:,i)=pspy(:,:,i)/art2(i)
  END DO
!# if defined (PLBC)
!  PSPY=0.0_SP
!# endif
  DO i=1,ncv_i
    i1 = ntrg(i)
    ia = niec(i,1)
    ib = niec(i,2)
!!$     XI = 0.5*(XIJE(I,1)+XIJE(I,2))
!!$     YI = 0.5*(YIJE(I,1)+YIJE(I,2))
!!$      
!!$#       if defined (SPHERICAL)
!!$        X1_DP=XIJE(I,1)
!!$        Y1_DP=YIJE(I,1)
!!$        X2_DP=XIJE(I,2)
!!$        Y2_DP=YIJE(I,2)
!!$        XII = XCG2(I)
!!$        YII = YCG2(I)
!!$        XI=XII               
!!$        XTMP  = XI*TPI-VX(IA)*TPI
!!$        XTMP1 = XI-VX(IA)
!!$        IF(XTMP1 >  180.0)THEN
!!$          XTMP = -360.0*TPI+XTMP
!!$        ELSE IF(XTMP1 < -180.0)THEN
!!$          XTMP =  360.0*TPI+XTMP
!!$        END IF        
!!$
!!$        DXA=XTMP*VAL_COS_VY(IA)    
!!$        DYA=(YI-VY(IA))*TPI
!!$        XTMP  = XI*TPI-VX(IB)*TPI
!!$        XTMP1 = XI-VX(IB)
!!$        IF(XTMP1 >  180.0)THEN
!!$          XTMP = -360.0*TPI+XTMP
!!$        ELSE IF(XTMP1 < -180.0)THEN
!!$          XTMP =  360.0*TPI+XTMP
!!$        END IF        
!!$
!!$        DXB=XTMP*VAL_COS_VY(IB) 
!!$        DYB=(YI-VY(IB))*TPI
!!$#       else
!!$        DXA = XI - VX(IA)
!!$        DYA = YI - VY(IA)
!!$     DXB = XI - VX(IB)
!!$     DYB = YI - VY(IB)
!!$#       endif
!!$      
!!$        FIJ1=N32(ID,ISS,IA)+DXA*PSPX(IA)+DYA*PSPY(IA)
!!$        FIJ2=N32(ID,ISS,IB)+DXB*PSPX(IB)+DYB*PSPY(IB)

! DEVELOPMENT TESTING - FIRST ORDER WAVE ACTION ADVECTION IN GEOG. SPACE
    fij1(:,:,1)=n32(:,:,ia)+dltxncve(i,1)*pspx(:,:,ia)+dltyncve(i,1)*pspy(:,:,ia)
    fij1(:,:,2)=n32(:,:,ib)+dltxncve(i,2)*pspx(:,:,ib)+dltyncve(i,2)*pspy(:,:,ib)
     
!DO ISS=1,MSC
!DO ID=1,MDC
!        AC1MIN=MINVAL(N32(ID,ISS,NBSN(IA,1:NTSN(IA)-1)))
!        AC1MIN=MIN(AC1MIN, N32(ID,ISS,IA))
!        AC1MAX=MAXVAL(N32(ID,ISS,NBSN(IA,1:NTSN(IA)-1)))
!        AC1MAX=MAX(AC1MAX, N32(ID,ISS,IA))
!        AC2MIN=MINVAL(N32(ID,ISS,NBSN(IB,1:NTSN(IB)-1)))
!        AC2MIN=MIN(AC2MIN, N32(ID,ISS,IB))
!        AC2MAX=MAXVAL(N32(ID,ISS,NBSN(IB,1:NTSN(IB)-1)))
!        AC2MAX=MAX(AC2MAX, N32(ID,ISS,IB))
!        IF(FIJ1(ID,ISS,1) < AC1MIN) FIJ1(ID,ISS,1) = AC1MIN
!        IF(FIJ1(ID,ISS,1) > AC1MAX) FIJ1(ID,ISS,1) = AC1MAX
!        IF(FIJ1(ID,ISS,2) < AC2MIN) FIJ1(ID,ISS,2) = AC2MIN
!        IF(FIJ1(ID,ISS,2) > AC2MAX) FIJ1(ID,ISS,2) = AC2MAX
!END DO
!END DO

    deploc = (dep2(nv(i1,1))+dep2(nv(i1,2))+dep2(nv(i1,3)))/3.0
    IF(deploc <= depmin)THEN
!   *** depth is negative ***
!     KWAVEL = -1.
     cgo_1d   = 0.
    ELSE
!   *** call KSCIP1 to compute KWAVE and CGO ***
     CALL kscip2(spcsig,deploc,cgo_1d)
    END IF

!        CALL SWAPAR1(I1,ISS,ID,DEP2(1),KWAVELOC,CGLOC)
      
    utmp = (compda(nv(i1,1),jvx2)+compda(nv(i1,2),jvx2)+        &
            compda(nv(i1,3),jvx2))/3.0
    vtmp = (compda(nv(i1,1),jvy2)+compda(nv(i1,2),jvy2)+        &
            compda(nv(i1,3),jvy2))/3.0
      
    DO id=1,mdc
      cgo_2d(id,:)=cgo_1d
    END DO  
    spcdir2=spcdir(:,2)
    spcdir3=spcdir(:,3)
    CALL sproxy2(canx_2d ,cany_2d ,   &
                 cgo_2d  ,spcdir2,spcdir3,utmp ,vtmp )

    un = cany_2d*dltxe(i) - canx_2d*dltye(i)  
!# if defined (PLBC)
!    UN =  - CANX*DLTYE(I)
!# endif
    exflux=-un*((1.0+sign(1.0,un))*fij1(:,:,2)+(1.0-sign(1.0,un))*fij1(:,:,1))*0.5

    xflux(:,:,ia)=xflux(:,:,ia)+exflux
    xflux(:,:,ib)=xflux(:,:,ib)-exflux
 
  END DO

  ELSE
   
  
  END IF  
     
  DO iss=1,msc
    DO id=1,mdc
      DO i = 1,m
        IF(isonb_w(i) /= 0)THEN
!#          if !defined(PLBC)
           deploc = dep2(i)  
!#          else        
!           CALL replace_node_wave(I,IG2)
!           !print*,'I,IG2===',I,IG2
!           DEPLOC = DEP2(IG2)
!#          endif           
          IF(deploc <= depmin)THEN
!         *** depth is negative ***
            kwaveloc = -1.                                           
            cgloc   = 0.                                            
          ELSE
!         *** call KSCIP1 to compute KWAVE and CGO ***
            spcsigl = spcsig(iss)
            CALL kscip1(1,spcsigl,deploc,kwaveloc,cgloc,nn,nd)                                 
          ENDIF
          cax(id,iss) = cgloc * spcdir(id,2)
          cay(id,iss) = cgloc * spcdir(id,3)
!
!         --- adapt the velocities in case of diffraction
!
          IF(idiffr == 1 .AND. pdiffr(3) /= 0.)THEN 
!JQI            CAX(ID,ISS) = CAX(ID,ISS)*DIFPARAM(I)      
!JQI            CAY(ID,ISS) = CAY(ID,ISS)*DIFPARAM(I)      
          END IF
!
!         --- ambient currents added
!
          IF(icur == 1)THEN 
            cax(id,iss) = cax(id,iss) + compda(i,jvx2)
            cay(id,iss) = cay(id,iss) + compda(i,jvy2)
          END IF
        
          canx = cax(id,iss)
      cany = cay(id,iss)
!          write(100,*)I,ID,ISS,CANX,CANY
     
      fij1_tmp = n32(id,iss,i)
     
          IF(nbsn(i,ntsn(i)-1) > m)THEN

        un1 = cany*(vx(i)-vx(nbsn(i,2)))      &
            -canx*(vy(i)-vy(nbsn(i,2)))
!# if defined (PLBC)
!            UN1 =-CANX*(VY(I)-VY(NBSN(I,2)))
!# endif
          ELSE IF(nbsn(i,2) > m)THEN
        un1 = cany*(vx(nbsn(i,ntsn(i)-1))-vx(i))      &
            -canx*(vy(nbsn(i,ntsn(i)-1))-vy(i))
!# if defined (PLBC)
!            UN1 = -CANX*(VY(NBSN(I,NTSN(I)-1))-VY(I))
!# endif
          ELSE
        un1 = cany*(vx(nbsn(i,ntsn(i)-1))-vx(nbsn(i,2)))      &
            -canx*(vy(nbsn(i,ntsn(i)-1))-vy(nbsn(i,2)))
!# if defined (PLBC)
!            UN1 =-CANX*(VY(NBSN(I,NTSN(I)-1))-VY(NBSN(I,2)))
!# endif
          END IF
         
      un1 = 0.5*un1
         
!          EXFLUX   = MAX(0.0,-UN1*FIJ1_TMP)
          xflux(id,iss,i) = xflux(id,iss,i)+max(0.0,-un1*fij1_tmp)
        END IF   
      END DO  
    END DO
  END DO
!
!-Accumulate Fluxes at Boundary Nodes
!
  DO iss=1,msc
    DO id=1,mdc
    END DO
  END DO

  DO i = 1,m
!--Update Action Density ------------------------------------------------------!
!
    IF(dep2(i) > depmin)THEN
      ac2(:,:,i) = n32(:,:,i)-xflux(:,:,i)/art1(i)*dtw
      ac2(:,:,i) = max(0.0,ac2(:,:,i))
    ELSE
      ac2(:,:,i) = 0.0_sp
    END IF  
  END DO

!# if defined(PLBC)
!   CALL replace_ac2(ID,ISS)
!# endif

  DO iss=1,msc
    DO id=1,mdc
    END DO
  END DO

!----yzhang----
!----yzhang----
  DEALLOCATE(dep2,ac2loc,xflux_tmp,exflux,ff1,fij1,un)
  DEALLOCATE(xflux,pspx,pspy)
  IF(high_latitude_wave)THEN
    DEALLOCATE(exfluxa,exfluxb,fij2)
    DEALLOCATE(fij1a,fij1b,fij2a,fij2b,una,unbb)
    DEALLOCATE(n32_tmpp2,n32_tmpp3)
  END IF
   
  RETURN

Here is the call graph for this function:

algorithm_crank_nicolson()

subroutine mod_action_ex::algorithm_crank_nicolson	(	real, dimension(:,:,:)	CAD,
		integer	IG,
		real	DTW,
		integer, dimension(msc)	IDCMIN,
		integer, dimension(msc)	IDCMAX,
		real	DD
	)

Definition at line 239 of file mod_action_ex.f90.

  USE swcomm3, ONLY : mdc,msc
!  USE MOD_USGRID, ONLY : MDC,MSC
  IMPLICIT NONE
  INTEGER :: ISS,ID,IDM1,IDM2,IDP1,MDCM,IG,II
  INTEGER :: IDDUM
  INTEGER, DIMENSION(MSC) :: IDCMIN,IDCMAX
  REAL, PARAMETER :: ZETA = 0.5
  REAL :: CAD(:,:,:)
  REAL :: DTW,DD
  REAL :: N32M,N32P
   
  REAL,DIMENSION(MDC) :: A,B,C,R,U
   
  DO iss = 1,msc
    DO iddum = idcmin(iss),idcmax(iss)
      id = mod(iddum-1+mdc,mdc)+1
      idp1 = mod(iddum+mdc,mdc)+1
      idm1 = mod(iddum-2+mdc,mdc)+1
 
      b(id) = 1.0
      IF(id == 1)THEN
        a(id) = 0.0
      ELSE
        a(id) = -0.5*zeta*dtw*cad(idm1,iss,1)/dd
      END IF
      IF(id == mdc)THEN
        c(id) = 0.0
      ELSE
        c(id) = 0.5*zeta*dtw*cad(idp1,iss,1)/dd
      END IF
 
      IF(id == 1)THEN
        r(id) = cad(idp1,iss,1)*n31(idp1,iss) 
        r(id) = (1.0-zeta)*0.5*dtw*r(id)/dd
        r(id) = n31(id,iss)-r(id)
      ELSE IF(id == mdc)THEN
        r(id) = -cad(idm1,iss,1)*n31(idm1,iss)
        r(id) = (1.0-zeta)*0.5*dtw*r(id)/dd
        r(id) = n31(id,iss)-r(id)       
      ELSE
        r(id) = cad(idp1,iss,1)*n31(idp1,iss)-cad(idm1,iss,1)*n31(idm1,iss) 
        r(id) = (1.0-zeta)*0.5*dtw*r(id)/dd
        r(id) = n31(id,iss)-r(id)     
      END IF

    END DO

    CALL tridag(a,b,c,r,u,mdc)
       
    DO iddum = idcmin(iss),idcmax(iss)
      id = mod(iddum-1+mdc,mdc)+1
      n32(id,iss,ig) = u(id)
    END DO   
  END DO
   
  RETURN

Here is the call graph for this function:

algorithm_fct()

subroutine mod_action_ex::algorithm_fct	(	real, dimension(mdc,msc,10)	CAS,
		integer	IG,
		real	DTW,
		integer, dimension(msc)	IDCMIN,
		integer, dimension(msc)	IDCMAX
	)

Definition at line 58 of file mod_action_ex.f90.

  USE swcomm3, ONLY : mdc,msc
  USE m_genarr, ONLY : spcsig
!  USE MOD_USGRID, ONLY : SPCSIG,MDC,MSC
  USE vars_wave, ONLY : mt,ac2
!  USE ALL_VARS, ONLY : MT,AC2
  IMPLICIT NONE

  INTEGER :: ISS,ID,IG
  REAL :: CASR,CASL,FLUX1,FLUX2,FLUXLP,FLUXLM,FLUXHP,FLUXHM
  REAL :: MIN11,MIN22,MIN33,ADLP,ADLM
  REAL, DIMENSION(MDC,MSC) :: ADP,ADM,NL
  REAL :: CAS(MDC,MSC,10)
  REAL :: DTW,DS
   
  INTEGER, DIMENSION(MSC) :: IDCMIN,IDCMAX
   
  n31 = 0.0
   
  DO iss = 1,msc
    IF(iss == 1)THEN
      ds   = spcsig(iss+1) - spcsig(iss)                              
      DO id = 1,mdc
        casr = 0.5*(cas(id,iss,1)+cas(id,iss+1,1))
        flux1 = 0.5*(casr+abs(casr))*ac2(id,iss,ig)
        flux2 = 0.5*(casr-abs(casr))*ac2(id,iss+1,ig)
        fluxlp = flux1+flux2

        fluxlm = 0.0
         
        fluxhp = casr*0.5*(ac2(id,iss,ig)+ac2(id,iss+1,ig))
        fluxhm = 0.0
   
        nl(id,iss) = ac2(id,iss,ig)-(fluxlp-fluxlm)*dtw/ds

        adp(id,iss) = fluxhp-fluxlp
        adm(id,iss) = fluxhm-fluxlm
      END DO     
    ELSE IF(iss == msc)THEN
      ds   = spcsig(iss) - spcsig(iss-1)                              
      DO id = 1,mdc
        casr = cas(id,iss,1)
        flux1 = 0.5*(casr+abs(casr))*ac2(id,iss,ig)
        flux2 = 0.0
        fluxlp = flux1+flux2

        casl = cas(id,iss-1,1)
        flux1 = 0.5*(casl+abs(casl))*ac2(id,iss-1,ig)
        flux2 = 0.5*(casl-abs(casl))*ac2(id,iss,ig)
        fluxlm = flux1+flux2
         
        fluxhp = casr*ac2(id,iss,ig)
        fluxhm = casl*ac2(id,iss-1,ig)
   
        nl(id,iss) = ac2(id,iss,ig)-(fluxlp-fluxlm)*dtw/ds

        adp(id,iss) = fluxhp-fluxlp
        adm(id,iss) = fluxhm-fluxlm
      END DO     
    ELSE
      ds   = spcsig(iss) - spcsig(iss-1)                              
      DO id = 1,mdc
        casr = 0.5*(cas(id,iss,1)+cas(id,iss+1,1))
        flux1 = 0.5*(casr+abs(casr))*ac2(id,iss,ig)
        flux2 = 0.5*(casr-abs(casr))*ac2(id,iss+1,ig)
        fluxlp = flux1+flux2

        casl = 0.5*(cas(id,iss,1)+cas(id,iss-1,1))
    flux1 = 0.5*(casl+abs(casl))*ac2(id,iss-1,ig)
        flux2 = 0.5*(casl-abs(casl))*ac2(id,iss,ig)
        fluxlm = flux1+flux2
         
        fluxhp = casr*0.5*(ac2(id,iss,ig)+ac2(id,iss+1,ig))
        fluxhm = casl*0.5*(ac2(id,iss,ig)+ac2(id,iss-1,ig))
   
        nl(id,iss) = ac2(id,iss,ig)-(fluxlp-fluxlm)*dtw/ds

        adp(id,iss) = fluxhp-fluxlp
        adm(id,iss) = fluxhm-fluxlm
      END DO     
    END IF
  END DO
   
  DO iss = 1,msc
    IF(iss == 1)THEN
      ds = spcsig(iss+1) - spcsig(iss)                              
      DO id = 1,mdc
        min11 = abs(adp(id,iss))
        min22 = sign(1.,adp(id,iss))*(nl(id,iss+2)-nl(id,iss+1))*ds/dtw
        adlp = min(min11,min22)
        adlp = max(0.,adlp)
        adlp = sign(1.,adp(id,iss))*adlp

        min11 = abs(adm(id,iss))
        min22 = sign(1.,adm(id,iss))*(nl(id,iss+1)-nl(id,iss))*ds/dtw
        adlm = min(min11,min22)
        adlm = max(0.,adlm)
        adlm = sign(1.,adm(id,iss))*adlm

        n31(id,iss) = nl(id,iss)-(adlp-adlm)*dtw/ds
      END DO
    ELSE IF(iss == 2)THEN
      ds   = spcsig(iss) - spcsig(iss-1)                              
      DO id = 1,mdc
        min11 = abs(adp(id,iss))
        min22 = sign(1.,adp(id,iss))*(nl(id,iss+2)-nl(id,iss+1))*ds/dtw
        min33 = sign(1.,adp(id,iss))*(nl(id,iss)-nl(id,iss-1))*ds/dtw
        adlp = min(min11,min22,min33)
        adlp = max(0.,adlp)
        adlp = sign(1.,adp(id,iss))*adlp

        min11 = abs(adm(id,iss))
        min22 = sign(1.,adm(id,iss))*(nl(id,iss+1)-nl(id,iss))*ds/dtw
        adlm = min(min11,min22)
        adlm = max(0.,adlm)
        adlm = sign(1.,adm(id,iss))*adlm
    
        n31(id,iss) = nl(id,iss)-(adlp-adlm)*dtw/ds
      END DO
    ELSE IF(iss == msc-1)THEN   
      ds   = spcsig(iss) - spcsig(iss-1)                              
      DO id = 1,mdc
        min11 = abs(adp(id,iss))
        min33 = sign(1.,adp(id,iss))*(nl(id,iss)-nl(id,iss-1))*ds/dtw
        adlp = min(min11,min33)
        adlp = max(0.,adlp)
        adlp = sign(1.,adp(id,iss))*adlp

        min11 = abs(adm(id,iss))
        min22 = sign(1.,adm(id,iss))*(nl(id,iss+1)-nl(id,iss))*ds/dtw
        min33 = sign(1.,adm(id,iss))*(nl(id,iss-1)-nl(id,iss-2))*ds/dtw
        adlm = min(min11,min22,min33)
        adlm = max(0.,adlm)
        adlm = sign(1.,adm(id,iss))*adlm
    
        n31(id,iss) = nl(id,iss)-(adlp-adlm)*dtw/ds
      END DO
    ELSE IF(iss == msc)THEN
      ds   = spcsig(iss) - spcsig(iss-1)                              
      DO id = 1,mdc
        min11 = abs(adp(id,iss))
        min33 = sign(1.,adp(id,iss))*(nl(id,iss)-nl(id,iss-1))*ds/dtw
        adlp = min(min11,min33)
        adlp = max(0.,adlp)
        adlp = sign(1.,adp(id,iss))*adlp

        min11 = abs(adm(id,iss))
        min33 = sign(1.,adm(id,iss))*(nl(id,iss-1)-nl(id,iss-2))*ds/dtw
        adlm = min(min11,min33)
        adlm = max(0.,adlm)
        adlm = sign(1.,adm(id,iss))*adlm
    
        n31(id,iss) = nl(id,iss)-(adlp-adlm)*dtw/ds
      END DO
    ELSE    
      ds   = spcsig(iss) - spcsig(iss-1)                              
      DO id = 1,mdc
        min11 = abs(adp(id,iss))
        min22 = sign(1.,adp(id,iss))*(nl(id,iss+2)-nl(id,iss+1))*ds/dtw
        min33 = sign(1.,adp(id,iss))*(nl(id,iss)-nl(id,iss-1))*ds/dtw
        adlp = min(min11,min22,min33)
        adlp = max(0.,adlp)
        adlp = sign(1.,adp(id,iss))*adlp

        min11 = abs(adm(id,iss))
        min22 = sign(1.,adm(id,iss))*(nl(id,iss+1)-nl(id,iss))*ds/dtw
        min33 = sign(1.,adm(id,iss))*(nl(id,iss-1)-nl(id,iss-2))*ds/dtw
        adlm = min(min11,min22,min33)
        adlm = max(0.,adlm)
        adlm = sign(1.,adm(id,iss))*adlm
    
        n31(id,iss) = nl(id,iss)-(adlp-adlm)*dtw/ds
      END DO
    END IF  
  END DO    

  RETURN

tridag()

subroutine mod_action_ex::tridag	(	real, dimension(n)	A,
		real, dimension(n)	B,
		real, dimension(n)	C,
		real, dimension(n)	R,
		real, dimension(n)	U,
		integer	N
	)

Definition at line 301 of file mod_action_ex.f90.

  IMPLICIT NONE
  INTEGER  :: N,J
  REAL,DIMENSION(N) :: A,B,C,R,U
  INTEGER, PARAMETER :: NMAX = 500
  REAL BET,GAM(NMAX)
    
  IF(b(1) == 0.)pause 'TRIDAG: REWRITE EQUATIONS'
  bet = b(1)
  u(1) = r(1)/bet
  DO j=2,n
    gam(j) = c(j-1)/bet
    bet = b(j)-a(j)*gam(j)
    IF(bet == 0.)pause 'TRIDAG FAILED'
    u(j) = (r(j)-a(j)*u(j-1))/bet
  END DO
  DO j=n-1,1,-1
    u(j) = u(j)-gam(j+1)*u(j+1)
  END DO
    
  RETURN

Here is the caller graph for this function:

Variable Documentation

n31

real, dimension(:,:), allocatable mod_action_ex::n31

Definition at line 17 of file mod_action_ex.f90.

  REAL, ALLOCATABLE :: N31(:,:)

n32

real, dimension(:,:,:), allocatable mod_action_ex::n32

Definition at line 18 of file mod_action_ex.f90.

  REAL, ALLOCATABLE :: N32(:,:,:)