SUBROUTINE VTADV !>-------------------------------------------------------------------------------------------------- !> SUBROUTINE VTADV !> !> SUBPROGRAM: VTADV - VERTICAL ADVECTION !> PROGRAMMER: JANJIC !> ORG: W/NP22 !> DATE: 93-11-17 !> !> ABSTRACT: !> VTADV CALCULATES THE CONTRIBUTION OF THE VERTICAL ADVECTION TO THE TENDENCIES OF TEMPERATURE, !> SPECIFIC HUMIDITY, WIND COMPONENTS, AND TURBULENT KINETIC ENERGY AND THEN UPDATES THOSE !> VARIABLES. FOR ALL VARIABLES EXCEPT SPECIFIC HUMIDITY A SIMPLE CENTERED DIFFERENCE SCHEME IN !> SPACE IS USED IN CONJUNCTION WITH THE PURE EULER-BACKWARD TIME SCHEME. !> A PIECEWISE LINEAR SCHEME IS USED TO CALCULATE THE VERTICAL ADVECTION OF SPECIFIC HUMIDITY SO !> THAT NO FALSE MAXIMA OR MINIMA ARE PRODUCED. !> !> PIECEWISE LINEAR SCHEME (MESINGER AND JOVIC, NCEP OFFICE NOTE #439) HERE USED FOR ALL VARIABLES, !> AVOIDING FALSE ADVECTION FROM BELOW GROUND, AND AT THE SAME TIME MAKING THE ETA VERY NEARLY A !> FINITE VOLUME MODEL. !> !> PROGRAM HISTORY LOG: !> 87-06-?? JANJIC - ORIGINATOR !> 90-??-?? MESINGER - INSERTED PIECEWISE LINEAR SCHEME FOR SPECIFIC HUMIDITY !> 95-03-25 BLACK - CONVERSION FROM 1-D TO 2-D IN HORIZONTAL !> 95-11-20 ABELES - PARALLEL OPTIMIZATION !> 96-03-29 BLACK - ADDED EXTERNAL EDGE; REMOVED SCRCH COMMON !> 98-10-30 BLACK - MODIFIED FOR DISTRIBUTED MEMORY !> 00-02-04 BLACK - ADDED CLOUD WATER/ICE !> 01-12-11 BLACK - SMOOTHING FOR CFL VIOLATION !> 18-01-15 LUCCI - MODERNIZATION OF THE CODE, INCLUDING: !> * F77 TO F90/F95 !> * INDENTATION & UNIFORMIZATION CODE !> * REPLACEMENT OF COMMONS BLOCK FOR MODULES !> * DOCUMENTATION WITH DOXYGEN !> * OPENMP FUNCTIONALITY !> !> INPUT ARGUMENT LIST: !> NONE !> !> OUTPUT ARGUMENT LIST: !> NONE !> !> INPUT/OUTPUT ARGUMENT LIST: !> NONE !> !> OUTPUT FILES: !> NONE !> !> USE MODULES: CLDWTR !> CONTIN !> CTLBLK !> DYNAM !> F77KINDS !> GLB_TABLE !> INDX !> LOOPS !> MAPPINGS !> MASKS !> MPPCOM !> NHYDRO !> PARMETA !> PVRBLS !> TEMPCOM !> TOPO !> TEMPCOM !> TOPO !> VRBLS !v !> DRIVER : EBU !> NEWFLT !> !> CALLS : ----- !>-------------------------------------------------------------------------------------------------- ! !-------------------------------------------------------------------------------------------- ! THIS IS THE CODE AS E-MAILED APRIL 10, 07, BUT WITH BUGS CORRECTED: ! 1) FFD1, TWICE, 1 IS REMOVED, ! 2) F4D, IN LOOP 1320, THE SECOND PAIR CHANGED INTO F4Q; ! 3) ..TEND1.. 2ND TIME IN LOOP 1330 REPL. BY ..TEND2.. ! AND ALSO: ! A) DEAD PIECES OF THE CODE - ABANDONED FINITE DIFFERENCE VERTICAL ADVECTION OF T AND U, V ! (TWICE "GO TO ..." SEGMENTS) - DELETED; ! B) SWITCHING TO CALLS EVERY TIME STEP COMMENTED OUT, AS THIS TAKES NON-NEGLIGIBLE TIME, AND ! IN THE ETA HORIZONTAL ADVECTION IS CALLED EVERY SECOND TIME STEP ALSO (CONSISTENCY) !-------------------------------------------------------------------------------------------- USE CLDWTR USE CONTIN USE CTLBLK USE DYNAM USE F77KINDS USE GLB_TABLE USE INDX USE LOOPS , ONLY : JAM USE MAPPINGS USE MASKS USE MPPCOM USE NHYDRO USE PARMETA USE PVRBLS USE TEMPCOM USE TOPO USE VRBLS ! IMPLICIT NONE ! REAL (KIND=R4KIND), PARAMETER :: EDQMX = 2.E-5 REAL (KIND=R4KIND), PARAMETER :: EDQMN = -2.E-5 REAL (KIND=R4KIND), PARAMETER :: EPSQ = 1.E-12 REAL (KIND=R4KIND), PARAMETER :: EPSQ2 = 0.12 REAL (KIND=R4KIND), PARAMETER :: CFL_MAX = 0.97 ! INTEGER(KIND=I4KIND), PARAMETER :: KSMUD = 0 ! #include "sp.h" ! INTEGER(KIND=I4KIND), PARAMETER :: IMJM = IM * JM - JM / 2 ! LOGICAL(KIND=L4KIND) ::& & NOSLA , NOSLAW ! REAL (KIND=R4KIND), DIMENSION(LM1) ::& & WFA , WFB ! REAL (KIND=R4KIND), DIMENSION(IDIM1:IDIM2, JDIM1:JDIM2) ::& & ETADTL , & & TQ2B , & & DQTI , DQBI , & & RPDX , RPDY , & & QDEDB , QDEUB , & & EDBD , & & DQDEB , & & DUTI , DUBI , & & DVTI , DVBI , & & UDEDB , UDEUB , & & VDEDB , VDEUB , & & EDBDU , EDBDV , & & DUDEB , DVDEB , & & FNE , FSE ! REAL (KIND=R4KIND), DIMENSION(IDIM1:IDIM2, JDIM1:JDIM2, LM) ::& & SAM , & & QBI , & & Q2ST , & & SAMU , & & UBI , & & SAMV , & & VBI ! REAL (KIND=R4KIND), DIMENSION(IDIM1:IDIM2, JDIM1:JDIM2, LM1) ::& & ARRAY1 , & & ARRAY2 , & & ARRAYU1 , & & ARRAYU2 , & & ARRAYV1 , & & ARRAYV2 ! REAL (KIND=R4KIND), DIMENSION(:,:,:), ALLOCATABLE ::& & S ! REAL (KIND=R4KIND), DIMENSION(IDIM1:IDIM2, JDIM1:JDIM2, LM) ::& & VAD_TEND1 , VAD_TEND2 ! REAL (KIND=R4KIND), DIMENSION(LM) ::& & VAD_TNDX1 , VAD_TNDX2 ! INTEGER(KIND=I4KIND), DIMENSION(IDIM1:IDIM2, JDIM1:JDIM2) ::& & LBOT_CFL_T , LTOP_CFL_T , & & LBOT_CFL_U , LTOP_CFL_U , & & LBOT_CFL_V , LTOP_CFL_V !------------------------ ! IMPLICIT NONE VARIABLES !------------------------ INTEGER(KIND=I4KIND) ::& & K , NMSAP , NMSAPW , I , J , NSMUD , KS , NS , MSA , & & IER , LSTART , LSTOP ! REAL (KIND=R4KIND) ::& & DTAD , CFL , EXTREM , DQTIK , ASTIK , ASBIK , QDEDTK , QDEUTK , SEDBK , & & DQDEK , EDBFK , EDTDK , TQ2AK , EXTREMU , EXTREMV , DUTIK , DVTIK , ASTIKU , & & ASTIKV , ASBIKU , ASBIKV , VMIJ , UDEDTK , VDEDTK , VDEUTK , UDEUTK , SEDBKU , & & SEDBKV , DUDEK , DVDEK , EDBFKU , EDBFKV , EDTDKU , EDTDKV ! DTAD = IDTAD * DT !------------------------------------------ ! DEFINE ADDED UPSTREAM ADVECTION CONSTANTS !------------------------------------------ ! DO 25 K=1,LM1 ! WFA(K) = DETA(K ) / (DETA(K) + DETA(K+1)) ! WFB(K) = DETA(K+1) / (DETA(K) + DETA(K+1)) ! 25 END DO WFA(1:LM1) = DETA(1:LM1 ) / (DETA(1:LM1) + DETA(2:LM)) WFB(1:LM1) = DETA(2:LM) / (DETA(1:LM1) + DETA(2:LM)) !------------------------------------------- ! NO MOISTURE SLOPE ADJUSTMENT IF NOT WANTED !------------------------------------------- NOSLA = .FALSE. NOSLAW = .FALSE. !----------------------------------------------------- ! IF FALSE, NUMBER OF MOISTURE SLOPE ADJUSTMENT PASSES !----------------------------------------------------- NMSAP = 3 NMSAPW = 3 !---------------------------------------- ! SMOOTHING VERTICAL VELOCITY AT H POINTS !---------------------------------------- IF (KSMUD > 0) THEN ! !$omp parallel do ! DO 90 K=1,LM1 DO 50 J=MYJS_P4,MYJE_P4 DO 50 I=MYIS_P4,MYIE_P4 ETADT(I,J,K) = ETADT(I,J,K) * HBM2(I,J) 50 END DO ! NSMUD = KSMUD !------------------------------------------------------------------------- ! HE FNE, FSE, ETADTL, AND ETADT ARRAYS ARE ON OR ASSOCIATED WITH H POINTS !------------------------------------------------------------------------- DO 90 KS=1,NSMUD DO 80 J=MYJS_P3,MYJE1_P3 DO 80 I=MYIS_P3,MYIE_P3 FNE(I,J) = (ETADT(I+IHE(J),J+1,K ) - ETADT(I ,J ,K )) & & * HTM(I ,J ,K+1) * HTM(I+IHE(J),J+1,K+1) 80 END DO ! DO 82 J=MYJS1_P3,MYJE_P3 DO 82 I=MYIS_P3,MYIE_P3 FSE(I,J) = (ETADT(I+IHE(J),J-1,K ) - ETADT(I,J,K )) & & * HTM(I+IHE(J),J-1,K+1) * HTM(I,J,K+1) 82 END DO ! DO 84 J=MYJS2_P1,MYJE2_P1 DO 84 I=MYIS_P1,MYIE_P1 ETADTL(I,J) = (FNE(I,J) - FNE(I+IHW(J),J-1) & & + FSE(I,J) - FSE(I+IHW(J),J+1)) & & * HBM2(I,J) 84 END DO ! DO 86 J=MYJS2_P1,MYJE2_P1 DO 86 I=MYIS_P1,MYIE_P1 ETADT(I,J,K) = ETADTL(I,J) * 0.125 + ETADT(I,J,K) 86 END DO ! 90 END DO ! END IF !------------------------------------------------------------------------------------------ ! IF THE CFL CRITERION IS VIOLATED THEN LOCATE VERTICAL LIMITS BETWEEN WHICH TO SMOOTH THE ! TENDENCIES !------------------------------------------------------------------------------------------ ! !$omp parallel do ! DO J=MYJS,MYJE DO I=MYIS,MYIE LTOP_CFL_T(I,J) = 0 LBOT_CFL_T(I,J) = 0 LTOP_CFL_U(I,J) = 0 LBOT_CFL_U(I,J) = 0 LTOP_CFL_V(I,J) = 0 LBOT_CFL_V(I,J) = 0 END DO END DO ! DO K=1,LM1 ! !$omp parallel do private (CFL) ! DO J=MYJS2,MYJE2 DO I=MYIS,MYIE !------------ ! MASS POINTS !------------ CFL = ETADT(I,J,K) * DTAD * HBM2(I,J)/(0.5*(DETA(K)+DETA(K+1))) IF (ABS(CFL) > CFL_MAX) THEN IF (LTOP_CFL_T(I,J) == 0) LTOP_CFL_T(I,J) = MAX(K,2) IF (LBOT_CFL_T(I,J) < K ) LBOT_CFL_T(I,J) = MAX(K,2) END IF !------------ ! U COMPONENT !------------ CFL = (ETADT(I+IVW(J),J,K) + ETADT(I+IVE(J),J,K)) & & * DTAD * VBM2(I,J) / (DETA(K)+DETA(K+1)) ! IF (ABS(CFL) > CFL_MAX) THEN IF (LTOP_CFL_U(I,J) == 0) LTOP_CFL_U(I,J) = MAX(K,2) IF (LBOT_CFL_U(I,J) < K ) LBOT_CFL_U(I,J) = MAX(K,2) END IF !------------ ! V COMPONENT !------------ CFL = (ETADT(I,J-1,K) + ETADT(I,J+1,K)) * DTAD * VBM2(I,J) / (DETA(K) + DETA(K+1)) ! IF (ABS(CFL) > CFL_MAX) THEN IF (LTOP_CFL_V(I,J) == 0) LTOP_CFL_V(I,J) = MAX(K,2) IF (LBOT_CFL_V(I,J) < K ) LBOT_CFL_V(I,J) = MAX(K,2) END IF ! END DO END DO ! END DO !------------------------------------------------------------ ! PIECEWISE LINEAR UPSTREAM VERTICAL ADVECTION OF Q AND CLOUD !------------------------------------------------------------ ALLOCATE(S(IDIM1:IDIM2, JDIM1:JDIM2,LM), STAT=I) !------------------------------------------------------------------------------------------------- ! INTIALIZE Q AT THE BOTTOM INTERFACE AND THE SLOPE ADJUSTMENT MASK (SAM=1 FOR SA PERMITTED, 0 FOR ! NOT PERMITTED) !------------------------------------------------------------------------------------------------- ! !-------------------------------- ! LOOP OVER S VARIABLES (Q,CWM,T) !-------------------------------- DO 400 NS=1,3 ! IF (NS == 1) THEN ! !$omp parallel do ! DO K=1,LM DO J=JDIM1,JDIM2 DO I=IDIM1,IDIM2 S(I,J,K) = Q(I,J,K) END DO END DO END DO ELSE IF (NS == 2) THEN ! !$omp parallel do ! DO K=1,LM DO J=JDIM1,JDIM2 DO I=IDIM1,IDIM2 S(I,J,K) = CWM(I,J,K) END DO END DO END DO ELSE IF (NS == 3) THEN ! !$omp parallel do ! DO K=1,LM DO J=JDIM1,JDIM2 DO I=IDIM1,IDIM2 S(I,J,K) = T(I,J,K) END DO END DO END DO ELSE WRITE(0,*)'ERROR IN VTADV. WILL STOP NOW.' CALL EXIT(103) END IF ! !$omp parallel do ! DO 175 K=1,LM DO 175 J=MYJS2,MYJE2 DO 175 I=MYIS,MYIE QBI(I,J,K) = S(I,J,K) SAM(I,J,K) = 1. 175 END DO ! IF (NOSLA) GOTO 290 !------------------------------------------------------ ! THE SLOPE ADJUSTMENT CODE ! NO SLOPE PERMITTED AT THE TOP AND AT THE BOTTOM LAYER !------------------------------------------------------ ! !$omp parallel do ! DO 190 J=MYJS2,MYJE2 DO 190 I=MYIS,MYIE SAM(I,J, 1) = 0. SAM(I,J,LM) = 0. 190 END DO ! !$omp parallel do ! DO 200 K=1,LM1 DO 200 J=MYJS2,MYJE2 DO 200 I=MYIS,MYIE SAM(I,J,K) = SAM(I,J,K) * HTM(I,J,K+1) 200 END DO !---------------------------------------------------------------------------------------- ! NOW, SEARCH FOR THE MAXIMA AND MINIMA OF Q (AT THE FIRST PASS) AND FOR LAYERS WHICH HAD ! OVERADJUSTED (AT SUBSEQUENT PASSES) DUE TO ROUND-OFF ERRORS !---------------------------------------------------------------------------------------- ! !$omp parallel do private (DQBI, DQTI, EXTREM) ! DO 220 K=2,LM1 DO 220 J=MYJS2,MYJE2 DO 220 I=MYIS,MYIE DQTI(I,J) = S(I,J,K ) - S(I,J,K-1) DQBI(I,J) = S(I,J,K+1) - S(I,J,K ) ! EXTREM = DQTI(I,J) * DQBI(I,J) ! IF (EXTREM <= 0.) SAM(I,J,K) = 0. 220 END DO ! !$omp parallel do ! DO 230 K=2,LM1 DO 230 J=MYJS2,MYJE2 DO 230 I=MYIS,MYIE ARRAY1(I,J,K) = WFA(K-1) * (1.-SAM(I,J,K-1)) + WFB(K-1) ARRAY2(I,J,K) = WFA(K ) + WFB(K) * (1.-SAM(I,J,K+1)) 230 END DO ! DO 260 MSA=1,NMSAP !------------------------------------------------------------------------------------------------- ! CALCULATE DQ AT INTERFACES AND ADJUST THE SLOPES WHERE AND TO THE EXTENT PERMITTED OBSERVING THE ! MONOTONICITY CONDITION (E.G. VAN LEER, J. COMP. PHYS. 1977, 276-299) !------------------------------------------------------------------------------------------------- ! !$omp parallel do ! DO 240 J=MYJS2,MYJE2 DO 240 I=MYIS,MYIE DQBI(I,J) = 2. * S(I,J,2) - QBI(I,J,2) - QBI(I,J,1) 240 END DO ! DO 250 K=2,LM1 ! !$omp parallel do private (ASBIK, ASTIK, DQTIK) ! DO 250 J=MYJS2,MYJE2 DO 250 I=MYIS,MYIE DQTIK = DQBI(I,J) ASTIK = ARRAY1(I,J,K) * DQTIK DQBI(I,J) = 2. * S(I,J,K+1) - QBI(I,J,K+1) - QBI(I,J,K) ASBIK = ARRAY2(I,J,K) * DQBI(I,J) QBI(I,J,K) = QBI(I,J,K) & & + (ASTIK-SIGN(1.,ASTIK) * DIM(ABS(ASTIK), ABS(ASBIK))) & & * SAM(I,J,K) 250 END DO ! 260 END DO !------------------------------------------------------------------------------------------------ ! SLOPE ADJUSTMENT OF THE LAYERS ABOVE THAT NEXT TO THE SURFACE IS DONE; NOW ADJUST THE LOWERMOST ! LAYER !------------------------------------------------------------------------------------------------ DO 270 K=9,LM1 ! !$omp parallel do ! DO 270 J=MYJS2,MYJE2 DO 270 I=MYIS,MYIE IF (HTM(I,J,K+1) == 0.) QBI(I,J,K) = 2. * S(I,J,K) - QBI(I,J,K-1) 270 END DO ! !$omp parallel do ! DO 280 J=MYJS2,MYJE2 DO 280 I=MYIS,MYIE QBI(I,J,LM) = 2. * S(I,J,LM) - QBI(I,J,LM1) 280 END DO !--------------------------------- ! END OF THE SLOPE ADJUSTMENT CODE !--------------------------------- 290 CONTINUE ! !$omp parallel do ! DO 300 J=MYJS2,MYJE2 DO 300 I=MYIS,MYIE QDEDB(I,J) = 0. QDEUB(I,J) = 0. DQDEB(I,J) = 2. * (QBI(I,J,1) - S(I,J,1)) * RDETA(1) EDBD (I,J) = 0. 300 END DO ! DO 320 K=1,LM1 ! !$omp parallel do private (DQDEK, EDBFK, EDTDK, QDEDTK, QDEUTK, SEDBK) ! DO 320 J=MYJS2,MYJE2 DO 320 I=MYIS,MYIE QDEDTK = QDEDB(I,J) QDEUTK = QDEUB(I,J) SEDBK = SIGN(1.,ETADT(I,J,K)) DQDEK = DQDEB(I,J) DQDEB(I,J) = 2. * (QBI(I,J,K+1) - S(I,J,K+1)) * RDETA(K+1) EDBFK = ETADT(I,J,K) * F4D QDEDB(I,J) = (1.+SEDBK) * (QBI(I,J,K) + DQDEK * EDBFK) * (-EDBFK) QDEUB(I,J) = (1.-SEDBK) * (2.*S(I,J,K+1) - QBI(I,J,K+1) + DQDEB(I,J)*EDBFK) & & * EDBFK ! EDTDK = EDBD(I,J) EDBD (I,J) = ETADT(I,J,K) * (-F4Q) S(I,J,K) = S(I,J,K) + (QDEDTK - QDEUTK - QDEDB(I,J) + QDEUB(I,J) & & + S(I,J,K) * (EDBD(I,J)-EDTDK)) * RDETA(K) 320 END DO ! !$omp parallel do ! DO 330 J=MYJS2,MYJE2 DO 330 I=MYIS,MYIE S(I,J,LM) = S(I,J,LM) + (QDEDB(I,J) - QDEUB(I,J) & & + S(I,J,LM) * (-EDBD(I,J))) * RDETA(LM) 330 END DO !-------------------------------------------------------- ! NEGATIVE MOISTURE MAY OCCUR DUE TO VIOLATION OF THE CFL !-------------------------------------------------------- DO 350 K=1,LM1 ! !$omp parallel do ! DO 350 J=MYJS2,MYJE2 DO 350 I=MYIS,MYIE IF (S(I,J,K) < EPSQ) THEN DQBI(I,J) = S(I,J,K ) S(I,J,K ) = EPSQ S(I,J,K+1) = S(I,J,K+1) + DETA(K) * RDETA(K+1) * DQBI(I,J) END IF 350 END DO ! !$omp parallel do ! DO 360 J=MYJS2,MYJE2 DO 360 I=MYIS,MYIE IF (S(I,J,LM) < EPSQ) S(I,J,LM) = EPSQ 360 END DO ! IF (NS == 1) THEN ! !$omp parallel do ! DO K=1,LM DO J=JDIM1,JDIM2 DO I=IDIM1,IDIM2 Q(I,J,K) = S(I,J,K) END DO END DO END DO ELSE IF (NS == 2) THEN ! !$omp parallel do ! DO K=1,LM DO J=JDIM1,JDIM2 DO I=IDIM1,IDIM2 CWM(I,J,K) = S(I,J,K) END DO END DO END DO ELSE IF (NS == 3) THEN ! !$omp parallel do ! DO K=1,LM DO J=JDIM1,JDIM2 DO I=IDIM1,IDIM2 T(I,J,K) = S(I,J,K) END DO END DO END DO ELSE WRITE(0,*)'ERROR IN VTADV. WILL STOP NOW.' CALL EXIT(104) END IF ! 400 END DO ! DEALLOCATE(S,STAT=IER) !----------------------------------- ! VERTICAL (MATSUNO) ADVECTION OF Q2 !----------------------------------- ! !$omp parallel do ! DO 420 J=MYJS2,MYJE2 DO 420 I=MYIS,MYIE TQ2B(I,J) = Q2(I,J,1) * ETADT(I,J,1) * F4Q2(1) 420 END DO ! DO 425 K=1,LM2 ! !$omp parallel do private (TQ2AK) ! DO 425 J=MYJS2,MYJE2 DO 425 I=MYIS,MYIE TQ2AK = (Q2(I,J,K+1) - Q2(I,J,K)) * (ETADT(I,J,K) + ETADT(I,J,K+1)) * F4Q2(K+1) Q2ST(I,J,K) = TQ2AK + TQ2B(I,J) + Q2(I,J,K) TQ2B(I,J) = TQ2AK 425 END DO ! !$omp parallel do private (TQ2AK) ! DO 440 J=MYJS2,MYJE2 DO 440 I=MYIS,MYIE TQ2AK = (Q2(I,J,LM) - Q2(I,J,LM1)) * ETADT(I,J,LM1) * F4Q2(LM) Q2ST(I,J,LM1) = TQ2AK + TQ2B(I,J) + Q2(I,J,LM1) Q2ST(I,J,LM ) = Q2(I,J,LM) 440 END DO !------------------------------- ! SECOND (BACKWARD) MATSUNO STEP !------------------------------- ! !$omp parallel do ! DO 450 J=MYJS2,MYJE2 DO 450 I=MYIS,MYIE TQ2B(I,J) = Q2ST(I,J,1) * ETADT(I,J,1) * F4Q2(1) 450 END DO ! DO K=1,LM2 ! !$omp parallel do private (TQ2AK) ! DO J=MYJS2,MYJE2 DO I=MYIS,MYIE TQ2AK = (Q2ST(I,J,K+1) - Q2ST(I,J,K)) * (ETADT(I,J,K) + ETADT(I,J,K+1)) * F4Q2(K+1) VAD_TEND1(I,J,K) = TQ2AK + TQ2B(I,J) TQ2B(I,J) = TQ2AK END DO END DO ! END DO ! DO J=MYJS2,MYJE2 DO I=MYIS,MYIE TQ2AK = (Q2ST(I,J,LM) - Q2ST(I,J,LM1)) * ETADT(I,J,LM1) * F4Q2(LM) VAD_TEND1(I,J,LM1) = TQ2AK + TQ2B(I,J) END DO END DO !--------------------------------------------------------------------- ! IF THE CFL CRITERION IS VIOLATED THEN VERTICALLY SMOOTH THE TENDENCY !--------------------------------------------------------------------- ! !$omp parallel do private (CFL, LBOT_CFL, LSTART, LSTOP, LTOP_CFL, VAD_TEND1, VAD_TNDX1) ! DO J=MYJS2,MYJE2 DO I=MYIS,MYIE ! IF (LTOP_CFL_T(I,J) > 0) THEN LSTART = LTOP_CFL_T(I,J) LSTOP = MIN(LBOT_CFL_T(I,J),LM-2) ! DO K=LSTART,LSTOP VAD_TNDX1(K) = (VAD_TEND1(I,J,K-1) + VAD_TEND1(I,J,K+1) + 2. & & * VAD_TEND1(I,J,K )) * 0.25 END DO DO K=LSTART,LSTOP VAD_TEND1(I,J,K) = VAD_TNDX1(K) END DO END IF ! END DO END DO ! DO 470 K=1,LM2 ! !$omp parallel do ! DO 470 J=MYJS2,MYJE2 DO 470 I=MYIS,MYIE Q2(I,J,K) = VAD_TEND1(I,J,K) + Q2(I,J,K) Q2(I,J,K) = AMAX1(Q2(I,J,K), EPSQ2) 470 END DO ! !$omp parallel do ! DO 480 J=MYJS2,MYJE2 DO 480 I=MYIS,MYIE Q2(I,J,LM1) = VAD_TEND1(I,J,LM1) + Q2(I,J,LM1) Q2(I,J,LM1) = AMAX1(Q2(I,J,LM1), EPSQ2) 480 END DO !------------------------------------------- ! DEFINITION OF VARIABLES NEEDED AT V POINTS !------------------------------------------- ! !$omp parallel do ! DO 500 K=1,LM1 DO 500 J=MYJS_P1,MYJE_P1 DO 500 I=MYIS_P1,MYIE_P1 ETADT(I,J,K) = ETADT(I,J,K) * PDSL(I,J) * HBM2(I,J) 500 END DO ! !$omp parallel do ! DO 510 J=MYJS2,MYJE2 DO 510 I=MYIS,MYIE RPDX(I,J) = 1. / (PDSL(I+IVW(J),J ) + PDSL(I+IVE(J),J )) RPDY(I,J) = 1. / (PDSL(I ,J-1) + PDSL(I ,J+1)) 510 END DO !-------------------------------------------------------------------------------------------------- ! PIECEWISE LINEAR UPSTREAM VERTICAL ADVECTION OF U AND V ! INTIALIZE U AND V AT THE BOTTOM INTERFACE AND THE SLOPE ADJUSTMENT MASK (SAMU=1 AND SAMV=1 FOR SA ! PERMITTED, 0 FOR NOT PERMITTED) !-------------------------------------------------------------------------------------------------- ! !$omp parallel do ! DO 999 K=1,LM DO 999 J=MYJS2,MYJE2 DO 999 I=MYIS,MYIE UBI(I,J,K) = U(I,J,K) VBI(I,J,K) = V(I,J,K) SAMU(I,J,K) = 1. SAMV(I,J,K) = 1. 999 END DO ! IF (NOSLAW) GOTO 998 !------------------------------------------------------ ! THE SLOPE ADJUSTMENT CODE ! NO SLOPE PERMITTED AT THE TOP AND AT THE BOTTOM LAYER !------------------------------------------------------ ! !$omp parallel do ! DO 997 J=MYJS2,MYJE2 DO 997 I=MYIS,MYIE SAMU(I,J, 1) = 0. SAMU(I,J,LM) = 0. SAMV(I,J, 1) = 0. SAMV(I,J,LM) = 0. 997 END DO ! !$omp parallel do ! DO 996 K=1,LM1 DO 996 J=MYJS2,MYJE2 DO 996 I=MYIS,MYIE SAMU(I,J,K) = SAMU(I,J,K) * VTM(I,J,K+1) SAMV(I,J,K) = SAMV(I,J,K) * VTM(I,J,K+1) 996 END DO !-------------------------------------------------------------------------------------------- ! NOW, SEARCH FOR THE MAXIMA AND MINIMA OF U & V (AT THE FIRST PASS) AND FOR LAYERS WHICH HAD ! OVERADJUSTED (AT SUBSEQUENT PASSES) DUE TO ROUND-OFF ERRORS !-------------------------------------------------------------------------------------------- ! !$omp parallel do private (DUBI, DVBI, DUTI, DVTI, EXTREMU, EXTREMV) ! DO 995 K=2,LM1 DO 995 J=MYJS2,MYJE2 DO 995 I=MYIS,MYIE DUTI(I,J) = U(I,J,K ) - U(I,J,K-1) DVTI(I,J) = V(I,J,K ) - V(I,J,K-1) DUBI(I,J) = U(I,J,K+1) - U(I,J,K ) DVBI(I,J) = V(I,J,K+1) - V(I,J,K ) ! EXTREMU = DUTI(I,J) * DUBI(I,J) EXTREMV = DVTI(I,J) * DVBI(I,J) ! IF (EXTREMU < 0.) SAMU(I,J,K) = 0. IF (EXTREMV < 0.) SAMV(I,J,K) = 0. 995 END DO ! !$omp parallel do ! DO 994 K=2,LM1 DO 994 J=MYJS2,MYJE2 DO 994 I=MYIS,MYIE ARRAYU1(I,J,K) = WFA(K-1) * (1.-SAMU(I,J,K-1)) + WFB(K-1) ARRAYV1(I,J,K) = WFA(K-1) * (1.-SAMV(I,J,K-1)) + WFB(K-1) ARRAYU2(I,J,K) = WFA(K ) + WFB(K) * (1.-SAMU(I,J,K+1)) ARRAYV2(I,J,K) = WFA(K ) + WFB(K) * (1.-SAMV(I,J,K+1)) 994 END DO ! DO 993 MSA=1,NMSAPW !-------------------------------------------------------------------------------------------------- ! CALCULATE DU & DV AT INTERFACES AND ADJUST THE SLOPES WHERE AND TO THE EXTENT PERMITTED OBSERVING ! THE MONOTONICITY CONDITION (E.G. VAN LEER, J. COMP. PHYS. 1977, 276-299, SCHEME USED HERE OF ! MESINGER AND JOVIC, NCEP OFFICE NOTE #439) !-------------------------------------------------------------------------------------------------- ! !$omp parallel do ! DO 992 J=MYJS2,MYJE2 DO 992 I=MYIS,MYIE DUBI(I,J) = 2. * U(I,J,2) - UBI(I,J,2) - UBI(I,J,1) DVBI(I,J) = 2. * V(I,J,2) - VBI(I,J,2) - VBI(I,J,1) 992 END DO ! !$omp parallel do private (ASBIKU, ASBIKV, ASTIKU, ASTIKV, DUTIK, DVTIK) ! DO 991 K=2,LM1 DO 991 J=MYJS2,MYJE2 DO 991 I=MYIS,MYIE DUTIK = DUBI(I,J) DVTIK = DVBI(I,J) ASTIKU = ARRAYU1(I,J,K) * DUTIK ASTIKV = ARRAYV1(I,J,K) * DVTIK ! DUBI(I,J) = 2.0 * U(I,J,K+1) - UBI(I,J,K+1) - UBI(I,J,K) DVBI(I,J) = 2.0 * V(I,J,K+1) - VBI(I,J,K+1) - VBI(I,J,K) ! ASBIKU = ARRAYU2(I,J,K) * DUBI(I,J) ASBIKV = ARRAYV2(I,J,K) * DVBI(I,J) ! UBI(I,J,K) = UBI(I,J,K) + (ASTIKU-SIGN(1.,ASTIKU) & & * DIM(ABS(ASTIKU),ABS(ASBIKU))) * SAMU(I,J,K) ! VBI(I,J,K) = VBI(I,J,K) + (ASTIKV-SIGN(1.,ASTIKV) & & * DIM(ABS(ASTIKV),ABS(ASBIKV))) * SAMV(I,J,K) 991 END DO 993 END DO !------------------------------------------------------------------------------------------------ ! SLOPE ADJUSTMENT OF THE LAYERS ABOVE THAT NEXT TO THE SURFACE IS DONE; NOW ADJUST THE LOWERMOST ! LAYER !------------------------------------------------------------------------------------------------ DO 990 K=9,LM1 ! !$omp parallel do ! DO 990 J=MYJS2,MYJE2 DO 990 I=MYIS,MYIE IF (VTM(I,J,K+1) == 0.) UBI(I,J,K) = 2. * U(I,J,K) - UBI(I,J,K-1) IF (VTM(I,J,K+1) == 0.) VBI(I,J,K) = 2. * V(I,J,K) - VBI(I,J,K-1) 990 END DO ! !$omp parallel do ! DO 989 J=MYJS2,MYJE2 DO 989 I=MYIS,MYIE UBI(I,J,LM) = 2. * U(I,J,LM) - UBI(I,J,LM1) VBI(I,J,LM) = 2. * V(I,J,LM) - VBI(I,J,LM1) 989 END DO !--------------------------------- ! END OF THE SLOPE ADJUSTMENT CODE !--------------------------------- 998 CONTINUE ! !$omp parallel do ! DO 988 J=MYJS2,MYJE2 DO 988 I=MYIS,MYIE UDEDB(I,J) = 0. VDEDB(I,J) = 0. UDEUB(I,J) = 0. VDEUB(I,J) = 0. DUDEB(I,J) = 2. * (UBI(I,J,1) - U(I,J,1)) * RDETA(1) DVDEB(I,J) = 2. * (VBI(I,J,1) - V(I,J,1)) * RDETA(1) EDBDU(I,J) = 0. EDBDV(I,J) = 0. 988 END DO ! DO 987 K=1,LM1 ! !$omp parallel do private (DQDEK, EDBFK, EDTDK, QDEDTK, QDEUTK, SEDBK) ! DO 987 J=MYJS2,MYJE2 DO 987 I=MYIS,MYIE VMIJ = VTM(I,J,K) * VBM2(I,J) ! UDEDTK = UDEDB(I,J) VDEDTK = VDEDB(I,J) UDEUTK = UDEUB(I,J) VDEUTK = VDEUB(I,J) ! SEDBKU = SIGN(1.,(ETADT(I+IVW(J),J ,K) + ETADT(I+IVE(J),J ,K))) SEDBKV = SIGN(1.,(ETADT(I ,J-1,K) + ETADT(I ,J+1,K))) ! DUDEK = DUDEB(I,J) DVDEK = DVDEB(I,J) ! DUDEB(I,J) = 2. * (UBI(I,J,K+1) - U(I,J,K+1)) * RDETA(K+1) DVDEB(I,J) = 2. * (VBI(I,J,K+1) - V(I,J,K+1)) * RDETA(K+1) ! EDBFKU = (ETADT(I+IVW(J),J ,K) + ETADT(I+IVE(J),J ,K)) * RPDX(I,J) * F4D EDBFKV = (ETADT(I ,J-1,K) + ETADT(I ,J+1,K)) * RPDY(I,J) * F4D ! UDEDB(I,J) = (1.+SEDBKU) * (UBI(I,J,K) + DUDEK * EDBFKU) * (-EDBFKU) VDEDB(I,J) = (1.+SEDBKV) * (VBI(I,J,K) + DVDEK * EDBFKV) * (-EDBFKV) ! UDEUB(I,J) = (1.-SEDBKU) * (U(I,J,K+1) + U(I,J,K+1) & & - UBI(I,J,K+1) + DUDEB(I,J) & & * EDBFKU) * EDBFKU ! VDEUB(I,J) = (1.-SEDBKV) * (V(I,J,K+1) + V(I,J,K+1) & & - VBI(I,J,K+1) + DVDEB(I,J) & & * EDBFKV) * EDBFKV EDTDKU = EDBDU(I,J) EDTDKV = EDBDV(I,J) ! EDBDU(I,J) = (ETADT(I+IVW(J),J ,K) + ETADT(I+IVE(J),J ,K)) * RPDX(I,J) * (-F4Q) EDBDV(I,J) = (ETADT(I ,J-1,K) + ETADT(I ,J+1,K)) * RPDY(I,J) * (-F4Q) ! VAD_TEND1(I,J,K) = (UDEDTK - UDEUTK - UDEDB(I,J) + UDEUB(I,J) + U(I,J,K) & & * (EDBDU(I,J) - EDTDKU)) * RDETA(K) * VMIJ ! VAD_TEND2(I,J,K) = (VDEDTK - VDEUTK - VDEDB(I,J) + VDEUB(I,J) + V(I,J,K) & & * (EDBDV(I,J) - EDTDKV)) * RDETA(K) * VMIJ 987 END DO ! !$omp parallel do ! DO 986 J=MYJS2,MYJE2 DO 986 I=MYIS,MYIE VMIJ = VTM(I,J,LM) * VBM2(I,J) VAD_TEND1(I,J,LM) = (UDEDB(I,J) - UDEUB(I,J) + U(I,J,LM) * (-EDBDU(I,J))) & & * RDETA(LM) * VMIJ VAD_TEND2(I,J,LM) = (VDEDB(I,J) - VDEUB(I,J) + V(I,J,LM) * (-EDBDV(I,J))) & & * RDETA(LM) * VMIJ 986 END DO !----------------------------------------------------------------------- ! IF THE CFL CRITERION IS VIOLATED THEN VERTICALLY SMOOTH THE TENDENCIES !----------------------------------------------------------------------- ! !$omp parallel do private (LSTART, LSTOP, VAD_TNDX1, VAD_TNDX2) ! DO J=MYJS2,MYJE2 DO I=MYIS,MYIE !---------- ! COMPONENT !---------- IF (LTOP_CFL_U(I,J) > 0) THEN LSTART = LTOP_CFL_U(I,J) LSTOP = MIN(LBOT_CFL_U(I,J), LM-1) ! DO K=LSTART,LSTOP VAD_TNDX1(K) = (VAD_TEND1(I,J,K-1) + VAD_TEND1(I,J,K+1) + 2. & & * VAD_TEND1(I,J,K )) * 0.25 END DO ! DO K=LSTART,LSTOP VAD_TEND1(I,J,K) = VAD_TNDX1(K) END DO ! END IF !------------ ! V COMPONENT !------------ IF (LTOP_CFL_V(I,J) > 0) THEN LSTART = LTOP_CFL_V(I,J) LSTOP = MIN(LBOT_CFL_V(I,J),LM-1) ! DO K=LSTART,LSTOP VAD_TNDX2(K) = (VAD_TEND2(I,J,K-1) + VAD_TEND2(I,J,K+1) + 2. & & * VAD_TEND2(I,J,K )) * 0.25 END DO DO K=LSTART,LSTOP VAD_TEND2(I,J,K) = VAD_TNDX2(K) END DO ! END IF ! END DO END DO ! DO 580 K=1,LM ! !$omp parallel do ! DO 580 J=MYJS2,MYJE2 DO 580 I=MYIS,MYIE U(I,J,K) = U(I,J,K) + VAD_TEND1(I,J,K) V(I,J,K) = V(I,J,K) + VAD_TEND2(I,J,K) 580 END DO ! IF (.NOT. HYDRO) THEN ! !$omp parallel do ! DO K=1,LM1 DO J=MYJS_P1,MYJE_P1 DO I=MYIS_P1,MYIE_P1 ETADT(I,J,K) = ETADT(I,J,K) / PDSL(I,J) END DO END DO END DO ! END IF ! RETURN ! END SUBROUTINE VTADV