C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& SUBROUTINE HDIFFS C C ****************************************************************** C$$$ SUBPROGRAM DOCUMENTATION BLOCK C . . . C SUBPROGRAM: HDIFF HORIZONTAL DIFFUSION C PRGRMMR: JANJIC ORG: W/NP22 DATE: 93-11-17 C C ABSTRACT: C HDIFF CALCULATES THE CONTRIBUTION OF THE HORIZONTAL DIFFUSION C TO THE TENDENCIES OF TEMPERATURE, SPECIFIC HUMIDITY, WIND C COMPONENTS, AND TURBULENT KINETIC ENERGY AND THEN UPDATES THOSE C VARIABLES. A SECOND-ORDER NONLINEAR SCHEME SIMILAR TO C SMAGORINSKYS IS USED WHERE THE DIFFUSION COEFFICIENT IS C A FUNCTION OF THE DEFORMATION FIELD AND OF THE TURBULENT C KINETIC ENERGY. C C PROGRAM HISTORY LOG: C 87-06-?? JANJIC - ORIGINATOR C 95-03-25 BLACK - CONVERSION FROM 1-D TO 2-D IN HORIZONTAL C 96-03-28 BLACK - ADDED EXTERNAL EDGE C 98-10-30 BLACK - MODIFIED FOR DISTRIBUTED MEMORY C 13-09-09 MESINGER, SUEIRO, GOMES, MORELLI C Code changed to account for slopes (Lines marked by Cfm C and CGCM) C Changes due to diffusion limited not to change the sign C of the Laplacian, in order to prevent instability due C to large values of DEF C C USAGE: CALL HDIFF FROM MAIN PROGRAM EBU C C INPUT ARGUMENT LIST: C NONE C C OUTPUT ARGUMENT LIST: C NONE C C OUTPUT FILES: C NONE C C SUBPROGRAMS CALLED: C C UNIQUE: NONE C C LIBRARY: NONE C C COMMON BLOCKS: CTLBLK C MASKS C PHYS C VRBLS C PVRBLS C INDX C C ATTRIBUTES: C LANGUAGE: FORTRAN 90 C MACHINE : IBM SP C$$$ C*********************************************************************** Cfm Calculaton of hdiff v at points that have a neighboring "blocked" v C-- switched on (loops 410 and 420), using velocities at points on C-- slopes, but only half of the diffusion coefficient (note Fig. 2 of C-- the "upgraded Eta" paper) C-- fm and Sandra Morelli, June-July 2013 C----------------------------------------------------------------------- C ****************************************************************** P A R A M E T E R & (DEFC=8.0,DEFM=32.0,SCQ2=50.0 &, EPSQ2=0.2,FCDIF=1.0,RFCP=.25/1004.6) C---------------------------------------------------------------------- INCLUDE "parmeta" INCLUDE "parm.tbl" INCLUDE "mpp.h" #include "sp.h" C---------------------------------------------------------------------- P A R A M E T E R & (IMJM=IM*JM-JM/2,LP1=LM+1,KSMUD=1) P A R A M E T E R &(JAM=6+2*(JM-10),JAMD=(JAM*2-10)*3) C----------------------------------------------------------------------- L O G I C A L & RUN,FIRST,RESTRT,SIGMA &,SECOND,HEAT,STTDF C---------------------------------------------------------------------- INCLUDE "CTLBLK.comm" C----------------------------------------------------------------------- INCLUDE "MASKS.comm" C----------------------------------------------------------------------- INCLUDE "PHYS.comm" C----------------------------------------------------------------------- INCLUDE "VRBLS.comm" C----------------------------------------------------------------------- INCLUDE "PVRBLS.comm" C----------------------------------------------------------------------- INCLUDE "CLDWTR.comm" C----------------------------------------------------------------------- INCLUDE "INDX.comm" CGSM----------------------------------------------------------------------- INCLUDE "SLOPES.comm" C----------------------------------------------------------------------- D I M E N S I O N & Q2L (IDIM1:IDIM2,JDIM1:JDIM2),UT (IDIM1:IDIM2,JDIM1:JDIM2) &,HKNE (IDIM1:IDIM2,JDIM1:JDIM2),HKSE (IDIM1:IDIM2,JDIM1:JDIM2) &,VKNE (IDIM1:IDIM2,JDIM1:JDIM2),VKSE (IDIM1:IDIM2,JDIM1:JDIM2) &,HMASK(IDIM1:IDIM2,JDIM1:JDIM2),HMSKL(IDIM1:IDIM2,JDIM1:JDIM2) C D I M E N S I O N & TNE (IDIM1:IDIM2,JDIM1:JDIM2),TSE (IDIM1:IDIM2,JDIM1:JDIM2) &,QNE (IDIM1:IDIM2,JDIM1:JDIM2),QSE (IDIM1:IDIM2,JDIM1:JDIM2) &,Q2NE (IDIM1:IDIM2,JDIM1:JDIM2),Q2SE (IDIM1:IDIM2,JDIM1:JDIM2) &,UNE (IDIM1:IDIM2,JDIM1:JDIM2),USE (IDIM1:IDIM2,JDIM1:JDIM2) &,VNE (IDIM1:IDIM2,JDIM1:JDIM2),VSE (IDIM1:IDIM2,JDIM1:JDIM2) &,TDIF (IDIM1:IDIM2,JDIM1:JDIM2),QDIF (IDIM1:IDIM2,JDIM1:JDIM2) &,UDIF (IDIM1:IDIM2,JDIM1:JDIM2),VDIF (IDIM1:IDIM2,JDIM1:JDIM2) &,Q2DIF(IDIM1:IDIM2,JDIM1:JDIM2) &,DEF (IDIM1:IDIM2,JDIM1:JDIM2),CKE (IDIM1:IDIM2,JDIM1:JDIM2) C----------------------------------------------------------------------- C*** C*** DIFFUSING Q2 AT GROUND LEVEL DOESNT MATTER, USTAR2 IS RECALCULATED C*** C----------------------------------------------------------------------- SECOND=.TRUE. HEAT=.FALSE. LUL=UL(1) C----------------------------------------------------------------------- C---------------------MAIN VERTICAL INTEGRATION LOOP-------------------- C----------------------------------------------------------------------- !$omp parallel do !$omp& private(cke,def,defsk,deftk,hkne,hkse,hmskl,q2dif,q2l,q2ne,q2se, !$omp& qdif,qne,qse,tdif,tne,tse,udif,une,use,utk,vdif,vkne, !$omp& vkse,vne,vse,vtk) C----------------------------------------------------------------------- C Cfm Fill the v points at slopes with the wind above -------------------- IF (L.GT.1) THEN DO 205 J=MYJS_P1,MYJE_P1 DO 205 I=MYIS_P1,MYIE_P1 IF (VTMS(I,J,L).EQ.1) THEN U(I,J,L)=U(I,J,L-1) V(I,J,L)=V(I,J,L-1) ENDIF 205 CONTINUE END IF C----------------------------------------------------------------------- DO 500 L=1,LM C----------------------------------------------------------------------- CALL ZERO2(DEF) CALL ZERO2(Q2NE) CALL ZERO2(Q2SE) CALL ZERO2(QNE) CALL ZERO2(QSE) CALL ZERO2(TNE) CALL ZERO2(TSE) CALL ZERO2(UNE) CALL ZERO2(USE) CALL ZERO2(VSE) CALL ZERO2(VNE) CGSM CALL ZERO2(CKE) CGSM CALL ZERO2(TDIF) CALL ZERO2(QDIF) CALL ZERO2(UDIF) CALL ZERO2(VDIF) CALL ZERO2(Q2DIF) C----------------------------------------------------------------------- DO 210 J=MYJS_P1,MYJE_P1 DO 210 I=MYIS_P1,MYIE_P1 Q2L(I,J)=AMAX1(Q2(I,J,L),EPSQ2) 210 CONTINUE C--------------DEFORMATIONS--------------------------------------------- DO 220 J=MYJS1_P1,MYJE1_P1 DO 220 I=MYIS_P1,MYIE1_P1 C DEFTK =U(I+IHE(J),J,L)-U(I+IHW(J),J,L)-V(I,J+1,L)+V(I,J-1,L) DEFSK =U(I,J+1,L)-U(I,J-1,L)+V(I+IHE(J),J,L)-V(I+IHW(J),J,L) DEF (I,J)=DEFTK*DEFTK+DEFSK*DEFSK DEF (I,J)=SQRT(DEF(I,J)+DEF(I,J))*HBM2(I,J) Cfm C Maximazing DEF so as not to be less than DEFC, the line below, C that someone uncommented, led to the crash of a 1 km run at CPTEC. C Commenting the line back, as it was some years ago, removed the C problem. Apparently, the line led to the diffusion coefficient C greater that it is allowed for stability. C DEF(I,J)=AMAX1(DEF(I,J),DEFC) c DEF(I,J)=AMIN1(DEF(I,J),DEFM) CGSM DEF(I,J)=DEF(I,J)*HMSKL(I,J) 220 CONTINUE C--------------T, Q, Q2 DIAGONAL CONTRIBUTIONS-------------------------- DO 250 J=MYJS_P1,MYJE1_P1 DO 250 I=MYIS_P1,MYIE1_P1 HKNE(I,J)=(DEF(I,J)+DEF(I+IHE(J),J+1)) 1 *HTM(I,J,L)*HTM(I+IHE(J),J+1,L) CGSM 2 *HMSKL(I,J)*HMSKL(I+IHE(J),J+1) TNE (I,J)=(T (I+IHE(J),J+1,L)-T (I,J,L))*HKNE(I,J) QNE (I,J)=(Q (I+IHE(J),J+1,L)-Q (I,J,L))*HKNE(I,J) Q2NE(I,J)=(Q2(I+IHE(J),J+1,L)-Q2(I,J,L))*HKNE(I,J) 250 CONTINUE C DO 260 J=MYJS1_P1,MYJE_P1 DO 260 I=MYIS_P1,MYIE1_P1 HKSE(I,J)=(DEF(I+IHE(J),J-1)+DEF(I,J)) 1 *HTM(I+IHE(J),J-1,L)*HTM(I,J,L) CGSM 2 *HMSKL(I+IHE(J),J-1)*HMSKL(I,J) TSE (I,J)=(T (I+IHE(J),J-1,L)-T (I,J,L))*HKSE(I,J) QSE (I,J)=(Q (I+IHE(J),J-1,L)-Q (I,J,L))*HKSE(I,J) Q2SE(I,J)=(Q2(I+IHE(J),J-1,L)-Q2(I,J,L))*HKSE(I,J) 260 CONTINUE C----------------------------------------------------------------------- DO 270 J=MYJS1,MYJE1 DO 270 I=MYIS1,MYIE TDIF (I,J)=(TNE (I,J)-TNE (I+IHW(J),J-1) 1 +TSE (I,J)-TSE (I+IHW(J),J+1))*HDAC(I,J) QDIF (I,J)=(QNE (I,J)-QNE (I+IHW(J),J-1) 1 +QSE (I,J)-QSE (I+IHW(J),J+1))*HDAC(I,J)*FCDIF Q2DIF(I,J)=(Q2NE(I,J)-Q2NE(I+IHW(J),J-1) 1 +Q2SE(I,J)-Q2SE(I+IHW(J),J+1))*HDAC(I,J) 270 CONTINUE C--------------2-ND ORDER DIFFUSION------------------------------------- IF(SECOND)THEN DO 280 J=MYJS2,MYJE2 DO 280 I=MYIS1,MYIE1 C Cfm Do not permit sign change due to horizontal diffusion. This will C prevent instability due to large DEF as well. TFPAV=.25*(T(I+IHE(J),J+1,L)+T(I+IHW(J),J+1,L) 1 +T(I+IHW(J),J-1,L)+T(I+IHE(J),J-1,L)) IF (ABS(TDIF(I,J)).GT.ABS(TFPAV-T(I,J,L))) THEN TDIF(I,J)=TFPAV-T(I,J,L) ENDIF C T (I,J,L)=T (I,J,L)+TDIF (I,J) C QFPAV=.25*(Q(I+IHE(J),J+1,L)+Q(I+IHW(J),J+1,L) 1 +Q(I+IHW(J),J-1,L)+Q(I+IHE(J),J-1,L)) IF (ABS(QDIF(I,J)).GT.ABS(QFPAV-Q(I,J,L))) THEN QDIF(I,J)=QFPAV-Q(I,J,L) ENDIF C Q (I,J,L)=Q (I,J,L)+QDIF (I,J) C 280 CONTINUE C C----------------------------------------------------------------------- IF(L.NE.LM)THEN DO 290 J=MYJS2,MYJE2 DO 290 I=MYIS1,MYIE1 Q2(I,J,L)=Q2(I,J,L)+Q2DIF(I,J)*HTM(I,J,L+1) 290 CONTINUE ENDIF C GO TO 360 ENDIF C--------------4-TH ORDER DIAGONAL CONTRIBUTIONS------------------------ DO 310 J=MYJS,MYJE1 DO 310 I=MYIS,MYIE1 TNE (I,J)=(TDIF (I+IHE(J),J+1)-TDIF (I,J))*HKNE(I,J) QNE (I,J)=(QDIF (I+IHE(J),J+1)-QDIF (I,J))*HKNE(I,J) Q2NE(I,J)=(Q2DIF(I+IHE(J),J+1)-Q2DIF(I,J))*HKNE(I,J) 310 CONTINUE C DO 320 J=MYJS1,MYJE DO 320 I=MYIS,MYIE1 TSE (I,J)=(TDIF (I+IHE(J),J-1)-TDIF (I,J))*HKSE(I,J) QSE (I,J)=(QDIF (I+IHE(J),J-1)-QDIF (I,J))*HKSE(I,J) Q2SE(I,J)=(Q2DIF(I+IHE(J),J-1)-Q2DIF(I,J))*HKSE(I,J) 320 CONTINUE C----------------------------------------------------------------------- DO 330 J=MYJS2,MYJE2 DO 330 I=MYIS1,MYIE1 T(I,J,L)=T(I,J,L)-(TNE (I,J)-TNE (I+IHW(J),J-1) 1 +TSE (I,J)-TSE (I+IHW(J),J+1))*HDAC(I,J) Q(I,J,L)=Q(I,J,L)-(QNE (I,J)-QNE (I+IHW(J),J-1) 1 +QSE (I,J)-QSE (I+IHW(J),J+1))*HDAC(I,J) 2 *FCDIF 330 CONTINUE C C----------------------------------------------------------------------- IF(L.NE.LM)THEN DO 340 J=MYJS2,MYJE2 DO 340 I=MYIS1,MYIE1 Q2(I,J,L)=Q2(I,J,L)-(Q2NE(I,J)-Q2NE(I+IHW(J),J-1) 1 +Q2SE(I,J)-Q2SE(I+IHW(J),J+1))*HDAC(I,J) 2 *HTM(I,J,L+1) 340 CONTINUE ENDIF C--------------U,V, DIAGONAL CONTRIBUTIONS------------------------------ 360 DO 410 J=MYJS_P1,MYJE1_P1 DO 410 I=MYIS_P1,MYIE1_P1 VKNE(I,J)=(DEF(I+IVE(J),J)+DEF(I,J+1)) CGSM 1 *VTM(I,J,L)*VTM(I+IVE(J),J+1,L) CGSM 2 *HMASK(I+IVE(J),J)*HMASK(I,J+1) 1 *MAX(VTM(I+IVE(J),J+1,L),VTMS(I+IVE(J),J+1,L)) UNE(I,J)=(U(I+IVE(J),J+1,L)-U(I,J,L))*VKNE(I,J) VNE(I,J)=(V(I+IVE(J),J+1,L)-V(I,J,L))*VKNE(I,J) 410 CONTINUE C DO 420 J=MYJS1_P1,MYJE_P1 DO 420 I=MYIS_P1,MYIE1_P1 VKSE(I,J)=(DEF(I,J-1)+DEF(I+IVE(J),J)) CGSM 1 *VTM(I+IVE(J),J-1,L)*VTM(I,J,L) CGSM 2 *HMASK(I,J-1)*HMASK(I+IVE(J),J) 1 *MAX(VTM(I+IVE(J),J-1,L),VTMS(I+IVE(J),J-1,L)) USE(I,J)=(U(I+IVE(J),J-1,L)-U(I,J,L))*VKSE(I,J) VSE(I,J)=(V(I+IVE(J),J-1,L)-V(I,J,L))*VKSE(I,J) 420 CONTINUE C----------------------------------------------------------------------- DO 430 J=MYJS1,MYJE1 DO 430 I=MYIS,MYIE1 UDIF(I,J)=(UNE(I,J)-UNE(I+IVW(J),J-1) 1 +USE(I,J)-USE(I+IVW(J),J+1))*HDACV(I,J) VDIF(I,J)=(VNE(I,J)-VNE(I+IVW(J),J-1) 1 +VSE(I,J)-VSE(I+IVW(J),J+1))*HDACV(I,J) 430 CONTINUE C--------------2-ND ORDER DIFFUSION------------------------------------- IF(SECOND)THEN DO 440 J=MYJS2,MYJE2 DO 440 I=MYIS1,MYIE1 CGSM U(I,J,L)=U(I,J,L)+UDIF(I,J) CGSM V(I,J,L)=V(I,J,L)+VDIF(I,J) C Cfm Do not permit sign change due to horizontal diffusion. This will C prevent instability due to large DEF as well. UFPAV=.25*(U(I+IVE(J),J+1,L)+U(I+IVW(J),J+1,L) 1 +U(I+IVW(J),J-1,L)+U(I+IVE(J),J-1,L)) IF (ABS(UDIF(I,J)).GT.ABS(UFPAV-U(I,J,L))) THEN UDIF(I,J)=UFPAV-U(I,J,L) ENDIF VFPAV=.25*(V(I+IVE(J),J+1,L)+V(I+IVW(J),J+1,L) 1 +V(I+IVW(J),J-1,L)+V(I+IVE(J),J-1,L)) IF (ABS(VDIF(I,J)).GT.ABS(VFPAV-V(I,J,L))) THEN VDIF(I,J)=VFPAV-V(I,J,L) ENDIF C UTK=U(I,J,L) VTK=V(I,J,L) NNTMP= VTM(I+IVW(J),J+1,L)*VTM(I+IVE(J),J+1,L) 1 *VTM(I+IVW(J),J-1,L)*VTM(I+IVE(J),J-1,L) DCMD=NNTMP+0.5*MOD(NNTMP+1,2) U(I,J,L)=U(I,J,L)+UDIF(I,J)*VTM(I,J,L)*DCMD V(I,J,L)=V(I,J,L)+VDIF(I,J)*VTM(I,J,L)*DCMD CKE(I,J)=0.5*(U(I,J,L)*U(I,J,L)-UTK*UTK 1 +V(I,J,L)*V(I,J,L)-VTK*VTK) Cfm C Let individual CKE be positive, but not the four point average used C to increase T C IF(CKE(I,J).GT.0)THEN C CKE(I,J)=0. Cfm ENDIF C 440 CONTINUE ELSE c GO TO 500 c ENDIF C--------------4-TH ORDER DIAGONAL CONTRIBUTIONS------------------------ DO 450 J=MYJS,MYJE1 DO 450 I=MYIS,MYIE1 UNE(I,J)=(UDIF(I+IVE(J),J+1)-UDIF(I,J))*VKNE(I,J) VNE(I,J)=(VDIF(I+IVE(J),J+1)-VDIF(I,J))*VKNE(I,J) 450 CONTINUE C DO 460 J=MYJS1,MYJE DO 460 I=MYIS,MYIE1 USE(I,J)=(UDIF(I+IVE(J),J-1)-UDIF(I,J))*VKSE(I,J) VSE(I,J)=(VDIF(I+IVE(J),J-1)-VDIF(I,J))*VKSE(I,J) 460 CONTINUE C----------------------------------------------------------------------- DO 470 J=MYJS2,MYJE2 DO 470 I=MYIS1,MYIE1 UTK=U(I,J,L) VTK=V(I,J,L) U(I,J,L)=U(I,J,L)-(UNE(I,J)-UNE(I+IVW(J),J-1) 1 +USE(I,J)-USE(I+IVW(J),J+1))*HDACV(I,J) V(I,J,L)=V(I,J,L)-(VNE(I,J)-VNE(I+IVW(J),J-1) 1 +VSE(I,J)-VSE(I+IVW(J),J+1))*HDACV(I,J) CKE(I,J)=0.5*(U(I,J,L)*U(I,J,L)-UTK*UTK 1 +V(I,J,L)*V(I,J,L)-VTK*VTK) 470 CONTINUE ENDIF C----------------------------------------------------------------------- IF(HEAT)THEN DO 480 J=MYJS2,MYJE2 DO 480 I=MYIS1,MYIE1 CKET=CKE(I+IHE(J),J)+CKE(I,J+1)+CKE(I+IHW(J),J)+CKE(I,J-1) IF (CKET.GT.0.) CKET=0. T(I,J,L)=-RFCP*CKET*HBM2(I,J) +T(I,J,L) 480 CONTINUE ENDIF C----------------------------------------------------------------------- Cfm Fill the v points at slopes back with zeros ------------------------ IF (L.GT.1) THEN DO 485 J=MYJS_P1,MYJE_P1 DO 485 I=MYIS_P1,MYIE_P1 IF (VTMS(I,J,L).EQ.1) THEN U(I,J,L)=0. V(I,J,L)=0. ENDIF 485 CONTINUE ENDIF C----------------------------------------------------------------------- 500 CONTINUE C----------------------------------------------------------------------- RETURN END