!/===========================================================================/ ! Copyright (c) 2007, The University of Massachusetts Dartmouth ! Produced at the School of Marine Science & Technology ! Marine Ecosystem Dynamics Modeling group ! All rights reserved. ! ! FVCOM has been developed by the joint UMASSD-WHOI research team. For ! details of authorship and attribution of credit please see the FVCOM ! technical manual or contact the MEDM group. ! ! ! This file is part of FVCOM. For details, see http://fvcom.smast.umassd.edu ! The full copyright notice is contained in the file COPYRIGHT located in the ! root directory of the FVCOM code. This original header must be maintained ! in all distributed versions. ! ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" ! AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, ! THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR ! PURPOSE ARE DISCLAIMED. ! !/---------------------------------------------------------------------------/ ! CVS VERSION INFORMATION ! $Id$ ! $Name$ ! $Revision$ !/===========================================================================/ !==============================================================================| ! Calculate the Turbulent Kinetic Energy and Mixing Length Based on | ! The Mellor-Yamada Level 2.5 Turbulent Closure Model | !==============================================================================| SUBROUTINE ADV_Q(Q,QF) !------------------------------------------------------------------------------| USE MOD_UTILS USE ALL_VARS USE MOD_PAR USE MOD_WD USE MOD_SPHERICAL USE MOD_NORTHPOLE # if defined (SEMI_IMPLICIT) USE MOD_SEMI_IMPLICIT # endif # if defined (PLBC) USE MOD_PERIODIC_LBC # endif IMPLICIT NONE REAL(SP), DIMENSION(0:MT,KB) :: Q,QF,XFLUX REAL(SP), DIMENSION(0:MT) :: PUPX,PUPY,PVPX,PVPY REAL(SP), DIMENSION(0:MT) :: PQPX,PQPY,PQPXD,PQPYD,VISCOFF REAL(SP), DIMENSION(3*(NT),KBM1) :: DTIJ REAL(SP), DIMENSION(3*(NT),KBM1) :: UVN REAL(SP) :: UTMP,VTMP,SITAI,FFD,FF1 !,X11,Y11,X22,Y22,X33,Y33,TMP1,TMP2,XI,YI REAL(SP) :: DXA,DYA,DXB,DYB,FIJ1,FIJ2,UN REAL(SP) :: TXX,TYY,FXX,FYY,VISCOF,EXFLUX,TEMP,STPOINT REAL(SP) :: FACT,FM1 INTEGER :: I,I1,I2,IA,IB,J,J1,J2,K,JTMP,JJ,II REAL(SP) :: Q1MIN, Q1MAX, Q2MIN, Q2MAX !!$# if defined (SPHERICAL) !!$ REAL(DP) :: TY,TXPI,TYPI !!$ REAL(DP) :: XTMP1,XTMP !!$ REAL(DP) :: X1_DP,Y1_DP,X2_DP,Y2_DP,XII,YII !!$ REAL(DP) :: X11_TMP,Y11_TMP,X33_TMP,Y33_TMP !!$ REAL(DP) :: VX1_TMP,VY1_TMP,VX2_TMP,VY2_TMP !!$ REAL(DP) :: TXPI_TMP,TYPI_TMP !!$# endif REAL(SP) :: QMEAN1 REAL(SP), DIMENSION(0:NT,KB) :: UQ,VQ REAL(SP), ALLOCATABLE :: UQ1(:,:),VQ1(:,:) # if defined (SEMI_IMPLICIT) REAL(SP) :: UN1 REAL(SP), DIMENSION(3*(NT),KBM1) :: UVN1 REAL(SP), DIMENSION(3*(NT),KBM1) :: DTIJ1 # endif IF(DBG_SET(DBG_SBR)) WRITE(IPT,*) "Start: adv_q" !------------------------------------------------------------------------------! QMEAN1 = 1.E-8 # if defined (SEMI_IMPLICIT) ALLOCATE(UQ1(0:NT,KB)); UQ1 = 0.0_SP ALLOCATE(VQ1(0:NT,KB)); VQ1 = 0.0_SP # else ALLOCATE(UQ1(0,0)) ALLOCATE(VQ1(0,0)) # endif SELECT CASE(HORIZONTAL_MIXING_TYPE) CASE ('closure') FACT = 1.0_SP FM1 = 0.0_SP CASE('constant') FACT = 0.0_SP FM1 = 1.0_SP CASE DEFAULT CALL FATAL_ERROR("UNKNOW HORIZONTAL MIXING TYPE:",& & TRIM(HORIZONTAL_MIXING_TYPE) ) END SELECT ! !--Initialize Fluxes-----------------------------------------------------------! ! QF = 0.0_SP XFLUX = 0.0_SP UQ = 0.0_SP VQ = 0.0_SP UVN = 0.0_SP # if defined (SEMI_IMPLICIT) UVN1 = 0.0_SP # endif DO K=2,KBM1 DO I=1,NT UQ(I,K) = (U(I,K)*DZ1(I,K-1)+U(I,K-1)*DZ1(I,K))/(DZ1(I,K)+DZ1(I,K-1)) VQ(I,K) = (V(I,K)*DZ1(I,K-1)+V(I,K-1)*DZ1(I,K))/(DZ1(I,K)+DZ1(I,K-1)) # if defined (SEMI_IMPLICIT) UQ1(I,K) = (UF(I,K)*DZ1(I,K-1)+UF(I,K-1)*DZ1(I,K))/(DZ1(I,K)+DZ1(I,K-1)) VQ1(I,K) = (VF(I,K)*DZ1(I,K-1)+VF(I,K-1)*DZ1(I,K))/(DZ1(I,K)+DZ1(I,K-1)) # endif END DO END DO ! !--Loop Over Control Volume Sub-Edges And Calculate Normal Velocity------------! ! DO I=1,NCV I1=NTRG(I) ! DTIJ(I)=DT1(I1) DO K=2,KBM1 DTIJ(I,K)=DT1(I1)*DZZ1(I1,K-1) UVN(I,K) = VQ(I1,K)*DLTXE(I) - UQ(I1,K)*DLTYE(I) # if defined (PLBC) UVN(I,K) = 0.0_SP*DLTXE(I) - UQ(I1,K)*DLTYE(I) # endif # if defined (SEMI_IMPLICIT) DTIJ1(I,K)=D1(I1)*DZZ1(I1,K-1) UVN1(I,K) = VQ1(I1,K)*DLTXE(I) - UQ1(I1,K)*DLTYE(I) # endif END DO END DO ! !--Calculate the Advection and Horizontal Diffusion Terms----------------------! ! DO K=2,KBM1 PQPX = 0.0_SP PQPY = 0.0_SP PQPXD = 0.0_SP PQPYD = 0.0_SP DO I=1,M DO J=1,NTSN(I)-1 I1=NBSN(I,J) I2=NBSN(I,J+1) ! FFD=0.5_SP*(Q(I1,K)+Q(I2,K)-QMEAN1(I1,K)-QMEAN1(I2,K)) FFD=0.5_SP*(Q(I1,K)+Q(I2,K)-QMEAN1-QMEAN1) FF1=0.5_SP*(Q(I1,K)+Q(I2,K)) !!$# if defined (SPHERICAL) !!$ XTMP = VX(I2)*TPI-VX(I1)*TPI !!$ XTMP1 = VX(I2)-VX(I1) !!$ IF(XTMP1 > 180.0_SP)THEN !!$ XTMP = -360.0_SP*TPI+XTMP !!$ ELSE IF(XTMP1 < -180.0_SP)THEN !!$ XTMP = 360.0_SP*TPI+XTMP !!$ END IF !!$ TXPI=XTMP*COS(DEG2RAD*VY(I)) !!$ TYPI=(VY(I1)-VY(I2))*TPI !!$ !!$ ! ERROR HERE !!$ IF(NODE_NORTHAREA(I) == 1)THEN !!$ VX1_TMP = REARTH * COS(VY(I1)*DEG2RAD) * COS(VX(I1)*DEG2RAD) & !!$ * 2._SP /(1._SP+sin(VY(I1)*DEG2RAD)) !!$ VY1_TMP = REARTH * COS(VY(I1)*DEG2RAD) * SIN(VX(I1)*DEG2RAD) & !!$ * 2._SP /(1._SP+sin(VY(I1)*DEG2RAD)) !!$ !!$ VX2_TMP = REARTH * COS(VY(I2)*DEG2RAD) * COS(VX(I2)*DEG2RAD) & !!$ * 2._SP /(1._SP+sin(VY(I2)*DEG2RAD)) !!$ VY2_TMP = REARTH * COS(VY(I2)*DEG2RAD) * SIN(VX(I2)*DEG2RAD) & !!$ * 2._SP /(1._SP+sin(VY(I2)*DEG2RAD)) !!$ !!$ TXPI = (VX2_TMP-VX1_TMP)/(2._SP /(1._SP+sin(VY(I)*DEG2RAD))) !!$ TYPI = (VY1_TMP-VY2_TMP)/(2._SP /(1._SP+sin(VY(I)*DEG2RAD))) !!$ IF(I /= NODE_NORTHPOLE)THEN !!$ TXPI_TMP = TYPI*COS(VX(I)*DEG2RAD)-TXPI*SIN(VX(I)*DEG2RAD) !!$ TYPI_TMP = TXPI*COS(VX(I)*DEG2RAD)+TYPI*SIN(VX(I)*DEG2RAD) !!$ TYPI_TMP = -TYPI_TMP !!$ !!$ TXPI = TXPI_TMP !!$ TYPI = TYPI_TMP !!$ END IF !!$ END IF !!$ ! END ERROR !!$ !!$ !!$ PQPX(I)=PQPX(I)+FF1*TYPI !!$ PQPY(I)=PQPY(I)+FF1*TXPI !!$ PQPXD(I)=PQPXD(I)+FFD*TYPI !!$ PQPYD(I)=PQPYD(I)+FFD*TXPI !!$# else !!$ PQPX(I)=PQPX(I)+FF1*(VY(I1)-VY(I2)) !!$ PQPY(I)=PQPY(I)+FF1*(VX(I2)-VX(I1)) !!$ PQPXD(I)=PQPXD(I)+FFD*(VY(I1)-VY(I2)) !!$ PQPYD(I)=PQPYD(I)+FFD*(VX(I2)-VX(I1)) !!$# endif PQPX(I)=PQPX(I)+FF1*DLTYTRIE(i,j) PQPY(I)=PQPY(I)+FF1*DLTXTRIE(i,j) PQPXD(I)=PQPXD(I)+FFD*DLTYTRIE(i,j) PQPYD(I)=PQPYD(I)+FFD*DltXTRIE(i,j) END DO PQPX(I)=PQPX(I)/ART2(I) PQPY(I)=PQPY(I)/ART2(I) PQPXD(I)=PQPXD(I)/ART2(I) PQPYD(I)=PQPYD(I)/ART2(I) END DO # if defined (PLBC) CALL replace_node(pqpx) CALL replace_node(pqpxd) CALL replace_node(pqpy) CALL replace_node(pqpyd) # endif DO I=1,M VISCOFF(I) = (VISCOFH(I,K)*DZ(I,K-1)+VISCOFH(I,K-1)*DZ(I,K))/ & (DZ(I,K)+DZ(I,K-1)) END DO DO I=1,NCV_I IA=NIEC(I,1) IB=NIEC(I,2) !!$ XI=0.5_SP*(XIJE(I,1)+XIJE(I,2)) !!$ YI=0.5_SP*(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) !!$ CALL ARCC(X2_DP,Y2_DP,X1_DP,Y1_DP,XII,YII) !!$ XI=XII !!$ XTMP = XI*TPI-VX(IA)*TPI !!$ XTMP1 = XI-VX(IA) !!$ IF(XTMP1 > 180.0_SP)THEN !!$ XTMP = -360.0_SP*TPI+XTMP !!$ ELSE IF(XTMP1 < -180.0_SP)THEN !!$ XTMP = 360.0_SP*TPI+XTMP !!$ END IF !!$ !!$ DXA=XTMP*COS(DEG2RAD*VY(IA)) !!$ DYA=(YI-VY(IA))*TPI !!$ XTMP = XI*TPI-VX(IB)*TPI !!$ XTMP1 = XI-VX(IB) !!$ IF(XTMP1 > 180.0_SP)THEN !!$ XTMP = -360.0_SP*TPI+XTMP !!$ ELSE IF(XTMP1 < -180.0_SP)THEN !!$ XTMP = 360.0_SP*TPI+XTMP !!$ END IF !!$ !!$ DXB=XTMP*COS(DEG2RAD*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=Q(IA,K)+DXA*PQPX(IA)+DYA*PQPY(IA) !!$ FIJ2=Q(IB,K)+DXB*PQPX(IB)+DYB*PQPY(IB) FIJ1=Q(IA,K)+DLTXNCVE(I,1)*PQPX(IA)+DLTYNCVE(I,1)*PQPY(IA) FIJ2=Q(IB,K)+DLTXNCVE(I,2)*PQPX(IB)+DLTYNCVE(I,2)*PQPY(IB) Q1MIN=MINVAL(Q(NBSN(IA,1:NTSN(IA)-1),K)) Q1MIN=MIN(Q1MIN, Q(IA,K)) Q1MAX=MAXVAL(Q(NBSN(IA,1:NTSN(IA)-1),K)) Q1MAX=MAX(Q1MAX, Q(IA,K)) Q2MIN=MINVAL(Q(NBSN(IB,1:NTSN(IB)-1),K)) Q2MIN=MIN(Q2MIN, Q(IB,K)) Q2MAX=MAXVAL(Q(NBSN(IB,1:NTSN(IB)-1),K)) Q2MAX=MAX(Q2MAX, Q(IB,K)) IF(FIJ1 < Q1MIN) FIJ1=Q1MIN IF(FIJ1 > Q1MAX) FIJ1=Q1MAX IF(FIJ2 < Q2MIN) FIJ2=Q2MIN IF(FIJ2 > Q2MAX) FIJ2=Q2MAX UN=UVN(I,K) # if defined (SEMI_IMPLICIT) UN1=UVN1(I,K) # endif ! VISCOF=HORCON*(FACT*0.5_SP*(VISCOFF(IA)+VISCOFF(IB))/HPRNU + FM1) ! VISCOF=HORCON*(FACT*0.5_SP*(VISCOFF(IA)+VISCOFF(IB)) + FM1) ! David moved HPRNU and added HVC VISCOF=(FACT*0.5_SP*(VISCOFF(IA)*NN_HVC(IA)+VISCOFF(IB)*NN_HVC(IB)) + FM1*0.5_SP*(NN_HVC(IA)+NN_HVC(IB)))/HPRNU !JQI NOV2021 # if defined (WET_DRY) ! Added by Adi Nugraha !Allow no diffusive flux between a wet node and a dry node (Wen Long and Tarang Khangaonkar) IF(ISWETN(IA)==1.AND.ISWETN(IB)==1)THEN # endif TXX=0.5_SP*(PQPXD(IA)+PQPXD(IB))*VISCOF TYY=0.5_SP*(PQPYD(IA)+PQPYD(IB))*VISCOF # if defined (WET_DRY) ELSE TXX = 0.0_SP TYY = 0.0_SP ENDIF # endif !JQI NOV2021 FXX=-DTIJ(I,K)*TXX*DLTYE(I) FYY= DTIJ(I,K)*TYY*DLTXE(I) # if defined (PLBC) FYY= 0.0_SP # endif # if !defined (SEMI_IMPLICIT) EXFLUX=-UN*DTIJ(I,K)* & ((1.0_SP+SIGN(1.0_SP,UN))*FIJ2+(1.0_SP-SIGN(1.0_SP,UN))*FIJ1)*0.5_SP+FXX+FYY # else EXFLUX=-UN*DTIJ(I,K)* & ((1.0_SP+SIGN(1.0_SP,UN))*FIJ2+(1.0_SP-SIGN(1.0_SP,UN))*FIJ1)*0.5_SP EXFLUX=(1.0_SP-IFCETA)*EXFLUX+IFCETA*(-UN1*DTIJ1(I,K)*((1.0_SP+SIGN(1.0_SP,UN1))*FIJ2+(1.0_SP-SIGN(1.0_SP,UN1))*FIJ1)*0.5_SP)+FXX+FYY # endif XFLUX(IA,K)=XFLUX(IA,K)+EXFLUX XFLUX(IB,K)=XFLUX(IB,K)-EXFLUX END DO # if defined (SPHERICAL) # if !defined (SEMI_IMPLICIT) CALL ADV_Q_XY(XFLUX,PQPX,PQPY,PQPXD,PQPYD,VISCOFF,Q,UQ,VQ,K,UQ1,VQ1,0.0_SP) # else CALL ADV_Q_XY(XFLUX,PQPX,PQPY,PQPXD,PQPYD,VISCOFF,Q,UQ,VQ,K,UQ1,VQ1,IFCETA) # endif # endif END DO !!SIGMA LOOP ! !-Accumulate Fluxes at Boundary Nodes ! # if defined (MULTIPROCESSOR) IF(PAR)CALL NODE_MATCH(0,NBN,BN_MLT,BN_LOC,BNC,MT,KB,MYID,NPROCS,XFLUX) # endif !-------------------------------------------------------------------- ! The central difference scheme in vertical advection !-------------------------------------------------------------------- DO K=2,KBM1 DO I=1,M # if defined (WET_DRY) IF(ISWETN(I)*ISWETNT(I) == 1) THEN # endif TEMP=WTS(I,K-1)*Q(I,K-1)-WTS(I,K+1)*Q(I,K+1) XFLUX(I,K)=XFLUX(I,K)+TEMP*ART1(I)*DZZ(I,K-1)/(DZ(I,K-1)+DZ(I,K)) # if defined (WET_DRY) END IF # endif END DO END DO !! SIGMA LOOP ! !--Update Q or QL-------------------------------------------------------------! ! # if defined (WET_DRY) DO I=1,M IF(ISWETN(I)*ISWETNT(I) == 1 )THEN DO K=2,KBM1 # if !defined (SEMI_IMPLICIT) QF(I,K)=(Q(I,K)-XFLUX(I,K)/ART1(I)*(DTI/(DT(I)*DZZ(I,k-1))))*(DT(I)/D(I)) # else QF(I,K)=(TE_TMP(I,K)-RK_TE(STAGE)*XFLUX(I,K)/ART1(I)*(DTI/(DT(I)*DZZ(I,k-1))))*(DT(I)/D(I)) # endif END DO ELSE DO K=2,KBM1 # if !defined (SEMI_IMPLICIT) QF(I,K)=Q(I,K) # else QF(I,K)=TE_TMP(I,K) # endif END DO END IF END DO # else DO I=1,M DO K=2,KBM1 # if !defined (SEMI_IMPLICIT) QF(I,K)=(Q(I,K)-XFLUX(I,K)/ART1(I)*(DTI/(DT(I)*DZZ(I,k-1))))*(DT(I)/D(I)) # else QF(I,K)=(TE_TMP(I,K)-RK_TE(STAGE)*XFLUX(I,K)/ART1(I)*(DTI/(DT(I)*DZZ(I,k-1))))*(DT(I)/D(I)) # endif END DO END DO # endif IF(DBG_SET(DBG_SBR)) WRITE(IPT,*) "End: adv_q" END SUBROUTINE ADV_Q !==============================================================================|