!copyright (C) 2001 MSC-RPN COMM %%%RPNPHY%%%
***S/P DIFUVDFj
*
#include "phy_macros_f.h"
SUBROUTINE DIFUVDFj(TU, U, KU, GU, R, ALFA, BETA, S, SK, 8,2
% TAU, TYPE, F, A, B, C, D, NU, NR, N, NK)
*
#include "impnone.cdk"
INTEGER NU, NR, N, NK
REAL TU(NU, NK), U(NU, NK), KU(NR, NK), GU(NR, NK), R(NR,NK)
REAL ALFA(N), BETA(N), S(n,NK), SK(n,NK), TAU, F
INTEGER TYPE
REAL A(N, NK), B(N, NK), C(N, NK), D(N, NK)
*
*Author
* R. Benoit (Mar 89)
*
*Revisions
* 001 R. Benoit (Aug 93) -Local sigma: s and sk become 2D
* 002 B. Bilodeau (Dec 94) - "IF" tests on integer
* instead of character.
* 003 J. Mailhot (Sept 00) - Add type 4='EB'
* 004 A. PLante (June 2003) - IBM conversion
* - calls to vrec routine (from massvp4 library)
*
*Object
* to solve a vertical diffusion equation by finite
* differences
*
*Arguments
*
* - Output -
* TU U tendency (D/DT U) due to the vertical diffusion and to
* term R
*
* - Input -
* U variable to diffuse (U,V,T,Q,E)
* KU diffusion coefficient
* GU optional countergradient term
* R optional inhomogeneous term
* ALFA inhomogeneous term for the surface flux (for type 1='U' or 2='UT')
* surface boundary condition (for type 4='EB')
* BETA homogeneous term for the surface flux
* S sigma coordinates of full levels
* SK sigma coordinates of diffusion coefficient levels
* TAU length of timestep
* TYPE type of variable to diffuse (1='U',2='UT', 3='E' or 4='EB')
* F waiting factor for time 'N+1'
* A work space (N,NK)
* B work space (N,NK)
* C work space (N,NK)
* D work space (N,NK)
* NU 1st dimension of TU and U
* NR 1st dimension of KU, GU and R
* N number of columns to process
* NK vertical dimension
*
*Notes
* D/DT U = D(U) + R
* D(U) = D/DS J(U)
* J(U) = KU*(D/DS U + GU)
* Limiting Conditions where S=ST: J=0(for 'U'), D=0(for 'UT'
* and ST=1)
* U=0(for 'E' and 'EB')
* Limiting Conditions where S=SB: J=ALFA+BETA*U(S(NK))(for
* 'U'/'UT'), J=0(for 'E'), U=ALFA(for 'EB')
* ST = S(1)-1/2 (S(2)-S(1)) (except for 'TU')
* SB = SK(NK)
*
**
*
INTEGER I, K, NKX
REAL SC, SCTU, ST, SB, HM, HP, HD, KUM, KUP, SCK1
EXTERNAL DIFUVD1, DIFUVD2
LOGICAL DEBUG
SAVE DEBUG
DATA DEBUG /.FALSE./
************************************************************************
* AUTOMATIC ARRAYS
************************************************************************
*
AUTOMATIC (VHM ,REAL ,(N,NK) )
AUTOMATIC (VHP ,REAL ,(N,NK) )
AUTOMATIC (RHD ,REAL*8 ,(N,NK) )
AUTOMATIC (RHMD ,REAL*8 ,(N,NK) )
AUTOMATIC (RHPD ,REAL*8 ,(N,NK) )
*
st(i)=s(i,1)
sb(i)=sk(i,nk)
*
IF (DEBUG) THEN
PRINT *,' S/R DIFUVDFj..TAU,TYPE,F,NU,N,NK=',
% TAU,TYPE,F,NU,N,NK
* PRINT *,' S=',S
* PRINT *,' SK=',SK
I=1
PRINT *,' ALFA(',I,')=',ALFA(I)
PRINT *,' BETA(',I,')=',BETA(I)
PRINT *,' U(',I,',*)=',(K,U(I,K),K=1,NK)
PRINT *,' KU(',I,',*)=',(K,KU(I,K),K=1,NK)
PRINT *,' GU(',I,',*)=',(K,GU(I,K),K=1,NK)
PRINT *,' R(',I,',*)=',(K,R(I,K),K=1,NK)
ENDIF
* these definitions hereunder for st, sb are always true.
* so use as statement functions in code directly
* to handle local sigma.
* ST=S(1)
* SB=SK(NK)
*
IF (TYPE.LE.2) THEN
SC=1
SCTU=0
NKX=NK
SCK1=1
IF (TYPE.EQ.2) THEN
SCK1=0
* line below for st=s(1) is commented out. see note above on
* definition of st and sb
* ST=S(1)
ENDIF
ELSE IF (TYPE.EQ.3 .OR. TYPE.EQ.4) THEN
SC=1-F
SCTU=1
NKX=NK-1
ELSE
PRINT *,' S/R DIFUVDFj. TYPE INCONNU= ',TYPE,' STOP...'
CALL QQEXIT(1)
ENDIF
*
* (1) CONSTRUIRE L'OPERATEUR TRIDIAGONAL DE DIFFUSION N=(A,B,C)
* ET LE TERME CONTRE-GRADIENT (DANS D)
*
IF (TYPE.LE.2) THEN
*
* K=1
*
HM=0
DO 10 I=1,N
HP=S(i,2)-S(i,1)
HD=0.5*(S(i,1)+S(i,2))-ST(i)
A(I,1)=0
B(I,1)=-SCK1*KU(I,1)/(HP*HD)
C(I,1)=-SCK1*B(I,1)
10 D(I,1)=SCK1*KU(I,1)*GU(I,1)/HD
*
* K=2...NK-1
*
DO K=2,NK-1,1
DO I=1,N
C THE FOLLOWING LHS ARE IN REAL
VHM(I,K)=S(I,K)-S(I,K-1)
VHP(I,K)=S(I,K+1)-S(I,K)
HD=0.5*(VHM(I,K)+VHP(I,K))
C THE FOLLOWING LHS ARE IN REAL*8
RHD(I,K)=HD
RHMD(I,K)=VHM(I,K)*HD
RHPD(I,K)=VHP(I,K)*HD
ENDDO
ENDDO
CALL VREC(RHD (1,2), RHD(1,2),N*(NK-2))
CALL VREC(RHMD(1,2),RHMD(1,2),N*(NK-2))
CALL VREC(RHPD(1,2),RHPD(1,2),N*(NK-2))
DO K=2,NK-1,1
DO I=1,N
A(I,K)=KU(I,K-1)*RHMD(I,K)
B(I,K)=-(KU(I,K-1)/VHM(I,K) +KU(I,K)/VHP(I,K))*RHD(I,K)
C(I,K)=KU(I,K)*RHPD(I,K)
D(I,K)=(KU(I,K)*GU(I,K)-KU(I,K-1)*GU(I,K-1))*
$ RHD(I,K)
ENDDO
ENDDO
*
* K=NK
*
HP=0
DO 12 I=1,N
HM=S(i,NK)-S(i,NK-1)
HD=SB(i)-0.5*(S(i,NK-1)+S(i,NK))
A(I,NK)=KU(I,NK-1)/(HM*HD)
B(I,NK)=-(KU(I,NK-1)/HM + 0)/HD
C(I,NK)=0
12 D(I,NK)=(0-KU(I,NK-1)*GU(I,NK-1))/HD
*
ELSE IF (TYPE.EQ.3 .OR. TYPE.EQ.4) THEN
*
* TYPE='E' or 'EB'
*
* K=1
*
DO 13 I=1,N
HM=S(i,2)-S(i,1)
HP=SK(i,2)-SK(i,1)
HD=0.5*(SK(i,2)+SK(i,1)) -S(i,1)
KUM=0.5*KU(I,1)
KUP=0.5*(KU(I,1)+KU(I,2))
A(I,1)=KUM/(HM*HD)
B(I,1)=-(KUM/HM+KUP/HP)/HD
C(I,1)=KUP/(HP*HD)
13 D(I,1)=(KUP*(GU(I,1)+GU(I,2))-KUM*GU(I,1))/(2*HD)
*
* K=2...NK-2
*
DO K=2,NK-2,1
DO I=1,N
C THE FOLLOWING LHS ARE IN REAL
VHM(I,K)=SK(I,K)-SK(I,K-1)
VHP(I,K)=SK(I,K+1)-SK(I,K)
HD=0.5*(VHM(I,K)+VHP(I,K))
C THE FOLLOWING LHS ARE IN REAL*8
RHD(I,K)=HD
RHMD(I,K)=VHM(I,K)*HD
RHPD(I,K)=VHP(I,K)*HD
ENDDO
ENDDO
CALL VREC( RHD(1,2), RHD(1,2),N*(NK-3))
CALL VREC(RHMD(1,2),RHMD(1,2),N*(NK-3))
CALL VREC(RHPD(1,2),RHPD(1,2),N*(NK-3))
DO K=2,NK-2,1
DO I=1,N
KUM=0.5*(KU(I,K-1)+KU(I,K))
KUP=0.5*(KU(I,K+1)+KU(I,K))
A(I,K)=KUM*RHMD(I,K)
B(I,K)=-(KUM/VHM(I,K) +KUP/VHP(I,K))*RHD(I,K)
C(I,K)=KUP*RHPD(I,K)
D(I,K)=.5*(KUP*(GU(I,K)+GU(I,K+1))
% -KUM*(GU(I,K-1)+GU(I,K)))*RHD(I,K)
ENDDO
ENDDO
*
* K=NK-1=NKX
*
IF (TYPE.EQ.3) THEN
*
* TYPE='E'
*
HP=0
DO 15 I=1,N
HM=SK(i,NK-1)-SK(i,NK-2)
HD=SB(i)-0.5*(SK(i,NK-1)+SK(i,NK-2))
KUM=0.5*(KU(I,NK-1)+KU(I,NK-2))
KUP=0
A(I,NKX)=KUM/(HM*HD)
B(I,NKX)=-(KUM/HM + 0)/HD
C(I,NKX)=0
15 D(I,NKX)=(0-KUM*(GU(I,NK-1)+GU(I,NK-2)))/(2*HD)
*
ELSE IF (TYPE.EQ.4) THEN
*
* TYPE='EB'
*
DO I=1,N
HM=SK(i,NK-1)-SK(i,NK-2)
HP=SB(i)-SK(i,NK-1)
HD=S(i,NK)-S(i,NK-1)
KUM=0.5*(KU(I,NK-1)+KU(I,NK-2))
KUP=0.5*(KU(I,NK)+KU(I,NK-1))
A(I,NKX)=KUM/(HM*HD)
B(I,NKX)=-(KUM/HM + KUP/HP)/HD
C(I,NKX)=0
D(I,NKX)=(KUP*(GU(I,NK)+GU(I,NK-1))
% -KUM*(GU(I,NK-1)+GU(I,NK-2)))/(2*HD)
% +KUP*ALFA(I)/(HD*HP)
ENDDO
*
ENDIF
*
ENDIF
*
*
* (2) CALCULER LE COTE DROIT D=TAU*(SC*N(U)+R+D/DS(KU*GU))
*
IF (DEBUG) THEN
PRINT *,' ETAPE(2) DE DIFUVDFj'
I=1
PRINT *,' A(',I,',*)=',(K,A(I,K),K=1,NK)
PRINT *,' B(',I,',*)=',(K,B(I,K),K=1,NK)
PRINT *,' C(',I,',*)=',(K,C(I,K),K=1,NK)
PRINT *,' D(',I,',*)=',(K,D(I,K),K=1,NK)
ENDIF
CALL DIFUVD1 (D, SC, A, B, C, U, D, N, NU, NKX)
DO 20 K=1,NKX
DO 20 I=1,N
20 D(I,K)=TAU*(D(I,K)+R(I,K)) +SCTU*U(I,K)
IF (DEBUG) THEN
PRINT *,' ETAPE(2) DE DIFUVDFj. D FINAL'
I=1
PRINT *,' D(',I,',*)=',(K,D(I,K),K=1,NK)
ENDIF
*
* (3) CALCULER OPERATEUR DU COTE GAUCHE
*
DO 30 K=1,NKX
DO 30 I=1,N
A(I,K)= -F*TAU*A(I,K)
B(I,K)=1-F*TAU*B(I,K)
30 C(I,K)= -F*TAU*C(I,K)
*
* (4) AJOUTER TERME DE FLUX DE SURFACE POUR TYPE='U'/'UT'
*
IF (TYPE.LE.2) THEN
DO 40 I=1,N
HD=SB(i)-0.5*(S(i,NK-1)+S(i,NK))
B(I,NKX)=B(I,NKX)-F*TAU*BETA(I)/HD
40 D(I,NKX)=D(I,NKX)+(ALFA(I)+BETA(I)*U(I,NKX))*TAU/HD
ENDIF
*
* (5) RESOUDRE SYSTEME TRIDIAGONAL [A,B,C] X = D. METTRE X DANS TU.
*
CALL DIFUVD2 (TU, A, B, C, D, D, NU, N, NKX)
*
* (6) OBTENIR TENDANCE
*
DO 60 K=1,NKX
DO 60 I=1,N
60 TU(I,K)=(TU(I,K)-SCTU*U(I,K))/TAU
* K=NKX+1..NK
DO 70 K=NKX+1,NK
DO 70 I=1,N
70 TU(I,K)=0
*
RETURN
END