SUBROUTINE EPS !>-------------------------------------------------------------------------------------------------- !> SUBROUTINE EPS !> !> SUBPROGRAM: EPS - ????? !> PROGRAMMER: JANJIC !> ORG: W/NP22 !> DATE: 9?-??-?? !> !> ABSTRACT: !> EPS COMPUTES THE VERTICAL AND HORIZONTAL ADVECTION OF DZ / DT. !> !> PROGRAM HISTORY LOG: !> 9?-??-?? JANJIC - ORIGINATOR !> 00-01-05 BLACK - DISTRIBUTED MEMORY AND THREADS !> 26-11-12 RISTIC IVAN - LOT OF CHANGES AND BUG FIX !> 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 !> !> USE MODULES: CONTIN !> CTLBLK !> DYNAM !> F77KINDS !> GLB_TABLE !> INDX !> LOOPS !> MAPPINGS !> MASKS !> MPPCOM !> NHYDRO !> PARMETA !> PVRBLS !> TEMPCOM !> TOPO !> VRBLS !> !> DRIVER : EBU !> !> CALLS : EXCH !> MPI_ALLREDUCE !-------------------------------------------------------------------------------------------------- USE CONTIN USE CTLBLK USE DYNAM USE EXCHM USE F77KINDS USE GLB_TABLE USE INDX USE LOOPS USE MAPOT USE MAPPINGS USE MASKS USE MPPCOM USE NHYDRO USE PARMETA USE PVRBLS USE TEMPCOM USE TOPO USE VRBLS ! IMPLICIT NONE ! INCLUDE "EXCHM.h" INCLUDE "mpif.h" ! INTEGER(KIND=I4KIND), PARAMETER :: IM1 = IM - 1 INTEGER(KIND=I4KIND), PARAMETER :: JAMD = (JAM * 2 - 10) * 3 INTEGER(KIND=I4KIND), PARAMETER :: NTSHY = 3 ! REAL (KIND=R4KIND), PARAMETER :: EPSQ = 1.E-12 REAL (KIND=R4KIND), PARAMETER :: EPSQ2 = 0.12 REAL (KIND=R4KIND), PARAMETER :: FF1 = 0.52500 REAL (KIND=R4KIND), PARAMETER :: FF2 = -0.64813 REAL (KIND=R4KIND), PARAMETER :: FF3 = 0.24520 REAL (KIND=R4KIND), PARAMETER :: FF4 = -0.12189 REAL (KIND=R4KIND), PARAMETER :: EPSFC = 1. / 1.05 REAL (KIND=R4KIND), PARAMETER :: EPSN = -EPSFC REAL (KIND=R4KIND), PARAMETER :: EPSP = EPSFC REAL (KIND=R4KIND), PARAMETER :: ZERO = 1.E-06 REAL (KIND=R4KIND), PARAMETER :: G = 9.8 REAL (KIND=R4KIND), PARAMETER :: CP = 1004.6 REAL (KIND=R4KIND), PARAMETER :: CAPA = 287.04 / CP REAL (KIND=R4KIND), PARAMETER :: GMA = -287.04 * (1. - CAPA) / 2. ! LOGICAL(KIND=L4KIND) ::& & TOP , BOT !------------------------ ! IMPLICIT NONE VARIABLES !------------------------ INTEGER(KIND=I4KIND) ::& & I , J , K , LMP , LAP , LLAP , JFP , JFQ , IRECV , & & JHL , JHH , IHL , IHH , JX , IX , KN , KP ! REAL (KIND=R4KIND) ::& & ADDT , RDT , RR , ARR , DWP , DETAL , RFACW , W4P , AFRP , & & D2PQW , WP , W00 , WP0 , ENH , HM , DVOLP , TTA , TTB , & & PP , QP , DWSTIJ , RFWIJ , W1IJ , WSTIJ , W0Q , WA , DWDTMX , & & DWDTMN , DWDTT , GDT , GDT2 , WGHT , FFC , PDP , RDP , DPPL , & & DPSTR , PP1 , TFC , TTFC , PSTRUP , RDPDN , RDPUP , RPD , PSTRDN , & & TMP , HBM2IJ , DPTU , FCT , DPTL , DELP !-------------------------------------- ! LOCAL VARIABLES REQUIRING DECLARATION !-------------------------------------- REAL (KIND=R4KIND), DIMENSION(LM+1) ::& & PONE , PSTR , PNP1 , COFF , DFRC , CHI , CHIN , WRES , CLIM ! REAL (KIND=R4KIND), DIMENSION(LM) ::& & RDPP , C0 , B1 , B2 , B3 , W3 , W4 ! REAL (KIND=R4KIND), DIMENSION(LM) ::& & ETADTL , DWL , AFR ! INTEGER(KIND=I4KIND), DIMENSION(LM) ::& & LA ! REAL (KIND=R4KIND), DIMENSION(IDIM1:IDIM2, JDIM1:JDIM2, LM) ::& & AFP , AFQ , W1 , DWST , DVOL ! REAL (KIND=R4KIND), DIMENSION(IDIM1:IDIM2, JDIM1:JDIM2) ::& & DARE , EMH , ANE , ASE ! INTEGER(KIND=I4KIND), DIMENSION(IDIM1:IDIM2, JDIM1:JDIM2, LM) ::& & IFPA , IFQA , & & IFPF , IFQF , & & JFPA , JFQA , & & JFPF , JFQF ! REAL (KIND=R4KIND), DIMENSION(2, LM) ::& & GSUMS , XSUMS ! REAL (KIND=R4KIND) ::& & SUMPW , SUMNW ! IF (NTSD <= NTSHY .OR. HYDRO) THEN !------- ! OPENMP !------- ! !$omp parallel do ! DO J=MYJS_P2,MYJE_P2 DO I=MYIS_P1,MYIE_P1 PINT(I,J,1) = PT END DO END DO DO K=1,LM !------- ! OPENMP !------- ! !$omp parallel do ! DO J=MYJS_P2,MYJE_P2 DO I=MYIS_P1,MYIE_P1 DWDT(I,J,K) = 1. PDWDT(I,J,K) = 1. PINT(I,J,K+1) = PDSL(I,J) * DETA(K) + PINT(I,J,K) END DO ENDDO END DO ! RETURN ! END IF ! ADDT = DT RDT = 1. / ADDT !-------------- ! TIME TENDENCY !-------------- !------- ! OPENMP !------- ! !$omp parallel do ! DO K=1,LM DO J=MYJS_P1,MYJE_P1 DO I=MYIS_P1,MYIE_P1 DWDT(I,J,K) = (W(I,J,K) - DWDT(I,J,K)) * HTM(I,J,K) * HBM2(I,J) * RDT END DO END DO END DO !------------------- ! VERTICAL ADVECTION !------------------- !------- ! OPENMP !------- ! !$omp parallel do !$omp private (AFR , AFRP , ARR , BOT , D2PQW , DETAL , DWL , DWP , & !$omp ETADTL , LA , LAP , LLAP , LMP , RFACW , RR , SUMNW , & !$omp SUMPW , TOP , W00 , W3 , W4 , W4P , WP , WP0 ) ! DO 200 J=MYJS_P1,MYJE_P1 DO 200 I=MYIS_P1,MYIE_P1 ! LMP=LMH(I,J) ! DO K=1,LMP W3(K) = W(I,J,K) W4(K) = W3(K) END DO ETADTL(1) = ETADT(I,J,1) * 0.5 DO K=2,LMP-1 ETADTL(K) = (ETADT(I,J,K-1) + ETADT(I,J,K)) * 0.5 END DO ! ETADTL(LMP) = ETADT(I,J,LMP-1) * 0.5 ! SUMPW = 0. SUMNW = 0. DO K=1,LMP RR = ETADTL(K) * (-ADDT) ARR = ABS(RR) IF (ARR > 0.) THEN LAP = RR / ARR ELSE LAP = 0 END IF ! LA(K) = LAP LLAP = K + LAP ! TOP = .FALSE. BOT = .FALSE. ! IF (LLAP > 0 .AND. LLAP < LMH(I,J)+1 .AND. LAP /= 0) THEN RR = ARR / ABS(AETA(LLAP) - AETA(K)) AFR(K) = (((FF4 * RR + FF3) * RR + FF2) * RR + FF1) * RR DWP = (W3(LLAP) - W3(K)) * RR DWL(K) = DWP ELSE TOP = LLAP == 0 BOT = LLAP == LMH(I,J)+1 ! RR = 0. AFR(K) = 0. DWL(K) = 0. END IF END DO ! IF (TOP) THEN IF (LA(2) < 0) THEN DWL(1) = -DWL(2) * DETA(2) / DETA(1) END IF END IF ! IF (BOT) THEN IF (LA(LMP-1) > 0)THEN DWL(LMP) = -DWL(LMP-1) * DETA(LMP-1) / DETA(LMP) END IF END IF ! DO K=1,LMP DETAL = DETA(K) DWP = DWL (K) * DETAL IF (DWP > 0.) THEN SUMPW = SUMPW + DWP ELSE SUMNW = SUMNW + DWP END IF END DO !------------------------------- ! FIRST MOMENT CONSERVING FACTOR !------------------------------- IF (SUMPW > 1.E-9) THEN RFACW = -SUMNW / SUMPW ELSE RFACW = 1. END IF ! IF(RFACW < 0.9 .OR. RFACW > 1.1) RFACW = 1. !--------------------------------- ! IMPOSE CONSERVATION ON ADVECTION !--------------------------------- IF (RFACW < 1.) THEN DO K=1,LMP DWP = DWL(K) IF (DWP < 0.) DWP = DWP / RFACW W4(K) = W3(K) + DWP END DO ELSE DO K=1,LMP DWP = DWL(K) IF (DWP >= 0.) DWP = DWP * RFACW W4(K) = W3(K) + DWP END DO END IF !-------------------- ! ANTI-FILTERING STEP !-------------------- SUMPW = 0. SUMNW = 0. !------------------------ ! ANTI-FILTERING LIMITERS !------------------------ DO 50 K=2,LMP-1 DETAL = DETA(K) ! W4P = W4(K) ! LAP = LA(K) ! IF (LAP /= 0) THEN AFRP = 2. * AFR(K) * (AETA(K+LAP) - AETA(K)) ** 2 / (AETA(K+LAP) - AETA(K-LAP)) D2PQW = ((W4(K+LAP) - W4P) / (AETA(K+LAP) - AETA(K)) & & - (W4P - W4(K-LAP)) / (AETA(K) - AETA(K-LAP))) & & * AFRP ELSE D2PQW = 0. END IF ! WP = W4P - D2PQW ! W00 = W3(K) WP0 = W3(K+LAP) ! WP = MAX(WP,MIN(W00,WP0)) WP = MIN(WP,MAX(W00,WP0)) ! DWP = WP - W4P ! DWL(K) = DWP ! DWP = DWP * DETAL ! IF (DWP > 0.) THEN SUMPW = SUMPW + DWP ELSE SUMNW = SUMNW + DWP END IF ! 50 END DO ! DWL(1) = 0. DWL(LMP) = 0. !------------------------------- ! FIRST MOMENT CONSERVING FACTOR !------------------------------- IF (SUMPW > 1.E-9) THEN RFACW = -SUMNW / SUMPW ELSE RFACW = 1. END IF ! IF (RFACW < 0.9 .OR. RFACW > 1.1) RFACW = 1. !-------------------------------------- ! IMPOSE CONSERVATION ON ANTI-FILTERING !-------------------------------------- IF (RFACW < 1.) THEN DO K=1,LMP DWP = DWL(K) IF (DWP >= 0.) DWP = DWP * RFACW DWDT(I,J,K) = DWDT(I,J,K) - (W4(K) + DWP - W3(K)) * RDT END DO ELSE DO K=1,LMP DWP = DWL(K) IF (DWP < 0.) DWP = DWP / RFACW DWDT(I,J,K) = DWDT(I,J,K) - (W4(K) + DWP - W3(K)) * RDT END DO END IF ! 200 END DO !-------------------------- ! END OF VERTICAL ADVECTION !-------------------------- ! ENH = ADDT / (08. * DY) !------- ! OPENMP !------- ! !$omp parallel do ! DO J=MYJS_P2,MYJE_P2 DO I=MYIS_P1,MYIE_P1 EMH(I,J) = ADDT / (08. * DX(I,J)) DARE(I,J) = HBM2(I,J) * DX(I,J) * DY END DO END DO !--------------------- ! HORIZONTAL ADVECTION !--------------------- !------- ! OPENMP !------- ! !$omp parallel do !$omp private (DVOLP , DWSTIJ , HM , JFP , JFQ , PP , QP , SUMNW , & !$omp SUMPW , TTA , TTB ) ! DO 300 K=1,LM ! DO J=MYJS_P2,MYJE_P2 DO I=MYIS_P1,MYIE_P1 DVOL(I,J,K) = DARE(I,J) * PDSL(I,J) * DETA(K) HM = HTM(I,J,K) W1(I,J,K) = W(I,J,K) * HM END DO END DO ! SUMPW = 0. SUMNW = 0. ! DO 225 J=MYJS2_P1,MYJE2_P1 DO 225 I=MYIS1_P1,MYIE1_P1 ! DVOLP = DVOL(I,J,K) * HBM3(I,J) TTA = (U(I,J-1,K) + U(I+IHW(J),J,K) + U(I+IHE(J),J,K) + U(I,J+1,K)) & & * HBM2(I,J) * EMH(I,J) ! TTB = (V(I,J-1,K) + V(I+IHW(J),J,K) + V(I+IHE(J),J,K) + V(I,J+1,K)) & & * HBM2(I,J) * ENH ! PP = -TTA - TTB QP = TTA - TTB ! JFP = INT(SIGN(1.,PP)) JFQ = INT(SIGN(1.,QP)) ! IFPA(I,J,K) = IHE(J) + I + ( JFP-1) / 2 IFQA(I,J,K) = IHE(J) + I + (-JFQ-1) / 2 ! JFPA(I,J,K) = J + JFP JFQA(I,J,K) = J + JFQ ! IFPF(I,J,K) = IHE(J) + I + (-JFP-1) / 2 IFQF(I,J,K) = IHE(J) + I + ( JFQ-1) / 2 ! JFPF(I,J,K) = J - JFP JFQF(I,J,K) = J - JFQ ! PP = ABS(PP) * HTM(I,J,K) * HTM(IFPA(I,J,K), JFPA(I,J,K),K) QP = ABS(QP) * HTM(I,J,K) * HTM(IFQA(I,J,K), JFQA(I,J,K),K) ! AFP(I,J,K) = (((FF4 * PP + FF3) * PP + FF2) * PP + FF1) * PP AFQ(I,J,K) = (((FF4 * QP + FF3) * QP + FF2) * QP + FF1) * QP ! DWSTIJ = (W(IFPA(I,J,K), JFPA(I,J,K),K) - W(I,J,K)) * PP & & + (W(IFQA(I,J,K), JFQA(I,J,K),K) - W(I,J,K)) * QP ! DWST(I,J,K) = DWSTIJ ! 225 END DO !---------------------------- ! GLOBAL SUM FOR CONSERVATION !---------------------------- DO 230 J=MYJS2,MYJE2 DO 230 I=MYIS1,MYIE1 ! DVOLP = DVOL(I,J,K) * HBM3(I,J) DWSTIJ = DWST(I,J,K) * DVOLP ! IF (DWSTIJ > 0.) THEN SUMPW = SUMPW + DWSTIJ ELSE SUMNW = SUMNW + DWSTIJ END IF ! 230 END DO ! XSUMS(1,K) = SUMPW XSUMS(2,K) = SUMNW CONTINUE !END OF VERTICAL LOOP ! 300 END DO !----------------- ! GLOBAL REDUCTION !----------------- CALL MPI_ALLREDUCE(XSUMS, GSUMS, 2*LM, MPI_REAL, MPI_SUM, MPI_COMM_COMP, IRECV) ! IF (IRECV /= 0) THEN PRINT*, ' RETURN CODE FROM 1ST ALLREDUCE IN EPS = ', IRECV CALL EXIT(101) END IF !------- ! OPENMP !------- ! !$omp parallel do !$omp private (D2PQW , DVOLP , DWSTIJ , RFACW , RFWIJ , SUMNW , SUMPW , W00 , & !$omp W0Q , W1IJ , WP0 , WSTIJ ) ! DO 400 K=1,LM ! SUMPW = GSUMS(1,K) SUMNW = GSUMS(2,K) !------------------------------- ! FIRST MOMENT CONSERVING FACTOR !------------------------------- IF (SUMPW > 1.) THEN RFACW = -SUMNW / SUMPW ELSE RFACW = 1. END IF ! IF (RFACW < 0.9 .OR. RFACW > 1.1) RFACW = 1. !--------------------------------- ! IMPOSE CONSERVATION ON ADVECTION !--------------------------------- IF (RFACW < 1.) THEN DO J=MYJS2_P1,MYJE2_P1 DO I=MYIS1_P1,MYIE1_P1 DWSTIJ = DWST(I,J,K) RFWIJ = HBM3(I,J) * (RFACW-1.) + 1. IF (DWSTIJ < 0.) DWSTIJ = DWSTIJ / RFWIJ W1(I,J,K) = W(I,J,K) + DWSTIJ END DO END DO ELSE DO J=MYJS2_P1,MYJE2_P1 DO I=MYIS1_P1,MYIE1_P1 DWSTIJ = DWST(I,J,K) RFWIJ = HBM3(I,J) * (RFACW-1.) + 1. IF (DWSTIJ >= 0.) DWSTIJ = DWSTIJ * RFWIJ W1(I,J,K) = W(I,J,K) + DWSTIJ END DO END DO END IF !-------------------- ! ANTI-FILTERING STEP !-------------------- SUMPW = 0. SUMNW = 0. !------------------------ ! ANTI-FILTERING LIMITERS !------------------------ DO 350 J=MYJS2,MYJE2 DO 350 I=MYIS1,MYIE1 ! DVOLP = DVOL(I,J,K) W1IJ = W1(I,J,K) ! D2PQW = (( W1(IFPA(I,J,K), JFPA(I,J,K), K) - W1IJ) & & - (W1IJ - W1(IFPF(I,J,K), JFPF(I,J,K), K)) & & * HTM(IFPF(I,J,K), JFPF(I,J,K), K)) & & * AFP(I,J,K) + ((W1(IFQA(I,J,K), JFQA(I,J,K), K) - W1IJ) & & - (W1IJ - W1(IFQF(I,J,K), JFQF(I,J,K), k)) & & * HTM(IFQF(I,J,K), JFQF(I,J,K), K)) * AFQ(I,J,K) ! WSTIJ = W1IJ - D2PQW ! W00 = W(I,J,K) WP0 = W(IFPA(I,J,K), JFPA(I,J,K), K) W0Q = W(IFQA(I,J,K), JFQA(I,J,K), K) ! WSTIJ = AMAX1(WSTIJ, AMIN1(W00,WP0,W0Q)) WSTIJ = AMIN1(WSTIJ, AMAX1(W00,WP0,W0Q)) ! DWSTIJ = WSTIJ - W1IJ ! DWST(I,J,K) = DWSTIJ ! DWSTIJ = DWSTIJ * DVOLP ! IF (DWSTIJ > 0.) THEN SUMPW = SUMPW + DWSTIJ ELSE SUMNW = SUMNW + DWSTIJ END IF ! 350 END DO ! XSUMS(1,K) = SUMPW XSUMS(2,K) = SUMNW ! CONTINUE !END OF VERTICAL LOOP 400 END DO !----------------- ! GLOBAL REDUCTION !----------------- CALL MPI_ALLREDUCE(XSUMS, GSUMS, 2*LM, MPI_REAL, MPI_SUM, MPI_COMM_COMP, IRECV) ! IF (IRECV /= 0) THEN PRINT*, ' RETURN CODE FROM 2ND ALLREDUCE IN EPS = ', IRECV CALL EXIT(102) END IF !------- ! OPENMP !------- ! !$omp parallel do !$omp private (DWSTIJ , RFACW , RFWIJ , SUMNW , SUMPW) ! DO 425 K=1,LM ! SUMPW = GSUMS(1,K) SUMNW = GSUMS(2,K) !------------------------------- ! FIRST MOMENT CONSERVING FACTOR !------------------------------- IF (SUMPW > 1.) THEN RFACW = -SUMNW / SUMPW ELSE RFACW = 1. END IF ! IF (RFACW < 0.9 .OR. RFACW > 1.1) RFACW = 1. !-------------------------------------- ! IMPOSE CONSERVATION ON ANTI-FILTERING !-------------------------------------- IF (RFACW < 1.) THEN DO J=MYJS2,MYJE2 DO I=MYIS1,MYIE1 DWSTIJ = DWST(I,J,K) RFWIJ = HBM2(I,J) * (RFACW-1.) + 1. IF (DWSTIJ >= 0.) DWSTIJ = DWSTIJ * RFWIJ DWDT(I,J,K) = (DWDT(I,J,K) - (W1(I,J,K) + DWSTIJ - W(I,J,K)) * RDT) * HBM3(I,J) END DO END DO ELSE DO J=MYJS2,MYJE2 DO I=MYIS1,MYIE1 DWSTIJ = DWST(I,J,K) RFWIJ = HBM2(I,J) * (RFACW-1.) + 1. IF (DWSTIJ < 0.) DWSTIJ = DWSTIJ / RFWIJ DWDT(I,J,K) = (DWDT(I,J,K) - (W1(I,J,K) + DWSTIJ - W(I,J,K)) * RDT) * HBM3(I,J) END DO END DO END IF ! 425 END DO !-------------------------------------------------- ! RESTRICTING THE ACCELERATION ALONG THE BOUNDARIES !-------------------------------------------------- JHL = 0 ! IF (JHL > 1) THEN JHH = JM - JHL + 1 ! IHL = JHL / 2 + 1 IHH = IM - IHL + MOD(J,2) !------- ! OPENMP !------- ! !$omp parallel do private (IX , JX) ! DO 450 K=1,LM ! DO J=1,JHL IF (J >= MY_JS_GLB-JBPAD2 .AND. J <= MY_JE_GLB+JTPAD2) THEN JX = J - MY_JS_GLB + 1 DO I=1,IM IF (I >= MY_IS_GLB-ILPAD2 .AND. I <= MY_IE_GLB+IRPAD2) THEN IX = I - MY_IS_GLB + 1 DWDT(IX,JX,K) = MAX(DWDT(IX,JX,k) * HBM3(IX,JX), -0.03) DWDT(IX,JX,K) = MIN(DWDT(IX,JX,K) * HBM3(IX,JX), 0.03) END IF END DO END IF END DO ! DO J=JHH,JM IF (J >= MY_JS_GLB-JBPAD2 .AND. J <= MY_JE_GLB+JTPAD2) THEN JX = J - MY_JS_GLB + 1 DO I=1,IM IF (I >= MY_IS_GLB-ILPAD2 .AND. I <= MY_IE_GLB+IRPAD2) THEN IX = I - MY_IS_GLB + 1 DWDT(IX,JX,K) = MAX(DWDT(IX,JX,K) * HBM3(IX,JX), -0.03) DWDT(IX,JX,K) = MIN(DWDT(IX,JX,K) * HBM3(IX,JX), 0.03) END IF END DO END IF END DO ! DO J=1,JM IF (J >= MY_JS_GLB-JBPAD2 .AND. J <= MY_JE_GLB+JTPAD2) THEN JX = J - MY_JS_GLB + 1 DO I=1,IHL IF (I >= MY_IS_GLB-ILPAD2 .AND. I <= MY_IE_GLB+IRPAD2) THEN IX = I - MY_IS_GLB + 1 DWDT(IX,JX,L) = MAX(DWDT(IX,JX,L) * HBM3(IX,JX), -0.03) DWDT(IX,JX,L) = MIN(DWDT(IX,JX,L) * HBM3(IX,JX), 0.03) END IF END DO END IF END DO ! DO J=1,JM IF (J >= MY_JS_GLB-JBPAD2 .AND. J <= MY_JE_GLB+JTPAD2) THEN JX = J - MY_JS_GLB + 1 DO I=IHH,IM IF (I >= MY_IS_GLB-ILPAD2 .AND. I <= MY_IE_GLB+IRPAD2) THEN IX = I - MY_IS_GLB + 1 DWDT(IX,JX,K) = MAX(DWDT(IX,JX,K) * HBM3(IX,JX), -0.03) DWDT(IX,JX,K) = MIN(DWDT(IX,JX,K) * HBM3(IX,JX), 0.03) END IF END DO END IF END DO ! 450 END DO ! END IF ! CALL EXCH(DWDT,LM,1,1) ! WA = 0.15 !------- ! OPENMP !------- ! !$omp parallel do private (ANE , ASE) ! DO 500 K=1,LM ! DO J=MYJS1_P1,MYJE2_P1 DO I=MYIS_P1,MYIE1_P1 ANE(I,J) = (DWDT(I+IHE(J),J+1,K) - DWDT(I,J,K)) * HTM(I,J,K) & & * HTM(I+IHE(J),J+1,K) END DO END DO ! DO J=MYJS2_P1,MYJE1_P1 DO I=MYIS_P1,MYIE1_P1 ASE(I,J) = (DWDT(I+IHE(J),J-1,K) - DWDT(I,J,K)) * HTM(I+IHE(J),J-1,K) & & * HTM(I,J,K) END DO END DO ! DO J=MYJS2,MYJE2 DO I=MYIS1,MYIE1 DWDT(I,J,K) = (ANE(I,J) - ANE(I+IHW(J),J-1) + ASE(I,J) - ASE(I+IHW(J),J+1)) & & * WA * HBM2(I,J) + DWDT(I,J,K) END DO END DO ! 500 END DO ! WP = 0.075 KN = 0 KP = 0 DWDTMX = 0. DWDTMN = 0. !------- ! OPENMP !------- ! !$omp parallel do !$omp private (DWDTMN , DWDTMX , DWDTT , HM , KN , KP) ! DO 525 K=1,LM ! DO J=MYJS,MYJE DO I=MYIS,MYIE HM = HTM(I,J,K) * HBM2(I,J) DWDTT = DWDT(I,J,K) * HM ! DWDTMX = AMAX1(DWDTT,DWDTMX) DWDTMN = AMIN1(DWDTT,DWDTMN) ! IF (DWDTT < EPSN) THEN DWDTT = EPSN KN = KN + 1 END IF ! IF (DWDTT > EPSP) THEN DWDTT = EPSP KP = KP + 1 END IF ! DWDT(I,J,K) = (DWDTT/G+1.) * (1.-WP) + PDWDT(I,J,K) * WP END DO END DO ! 525 END DO ! GDT = G * DT GDT2 = GDT * GDT WGHT = 0.35 FFC = -4. * R / GDT2 !------- ! OPENMP !------- ! !$omp parallel do !$omp private (CHI , CHIN , CHMOD , CLIM , COFF , DELP , DFRC , DP , & !$omp DPSTR , DPTL , DPTU , DWDTT , HBM2IJ , IMAX , INOT , ITER , & !$omp ITMX , JMAX , K , LMP , PDP , PNP1 , PONE , PP1 , & !$omp PSTR , RDP , RESDL , WIL , WRES ) ! DO 600 J=MYJS2,MYJE2 DO 600 I=MYIS1,MYIE1 ! LMP = LMH(I,J) PDP = PDSL(I,J) RPD = 1. / PDP ! PONE(1) = PT PSTR(1) = PT PNP1(1) = PT CHI(1) = 0. ! DO K=2,LMP+1 CHI(K) = 0. DPPL = DETA(K-1) * PDP RDPP(K-1) = 1. / DPPL PONE(K) = PINT(I,J,K) DPSTR = DWDT(I,J,K-1) * DPPL PSTR(K) = PSTR(K-1) + DPSTR PP1 = PNP1(K-1) + DPSTR PNP1(K) = (PP1-PONE(K)) * WGHT + PONE(K) TFC = Q(I,J,K-1) * 0.608 + 1. TTFC = -CAPA * TFC + 1. COFF(K-1) = T(I,J,K-1) * TTFC * TFC * DPPL * FFC / ((PNP1(K-1)+PNP1(K)) * & & (PNP1(K-1)+PNP1(K))) END DO ! PSTRUP = -(PSTR(1) + PSTR(2) - PONE(1) - PONE(2)) * COFF(1) ! DO K=2,LMP RDPDN = RDPP(K) RDPUP = RDPP(K-1) ! PSTRDN = -(PSTR(K) + PSTR(K+1) - PONE(K) - PONE(K+1)) * COFF(K) ! B1(K) = COFF(K-1) + RDPUP B2(K) = (COFF(K-1) + COFF(K)) - (RDPUP+RDPDN) B3(K) = COFF(K) + RDPDN C0(K) = PSTRUP + PSTRDN ! PSTRUP = PSTRDN END DO ! B1(2) = 0. B2(LMP) = B2(LMP) + B3(LMP) !------------ ! ELIMINATION !------------ DO K=3,LMP TMP = -B1(K) / B2(K-1) B2(K) = B3(K-1) * TMP + B2(K) C0(K) = C0(K-1) * TMP + C0(K) END DO ! CHI(1) = 0. !------------------ ! BACK SUBSTITUTION !------------------ CHI(LMP) = C0(LMP) / B2(LMP) CHI(LMP+1) = CHI(LMP) ! DO K=LMP-1,2,-1 CHI(K) = (-B3(k) * CHI(K+1) + C0(K)) / B2(K) END DO ! HBM2IJ = HBM2(I,J) DPTU = 0. FCT = 0.5 / CP * HBM2IJ ! DO K=1,LMP DPTL = (CHI(K+1) + PSTR(K+1) - PINT(I,J,K+1)) * HBM2IJ PINT(I,J,K+1) = PINT(I,J,K+1) + DPTL T (I,J,K) = (DPTU+DPTL) * RTOP(I,J,K) * FCT + T(I,J,K) DELP = (PINT(I,J,K+1) - PINT(I,J,K)) * RDPP(K) W(I,J,K) = ((DELP - DWDT(I,J,K)) * GDT + W(I,J,K)) * HBM2IJ DWDT(I,J,K) = (DELP-1.) * HBM2IJ + 1. ! DPTU = DPTL END DO ! DO K=LMP+1,LM PINT(I,J,K+1) = PDSL(I,J) * DETA(K) + PINT(I,J,K) END DO ! 600 END DO ! RETURN ! END SUBROUTINE EPS