!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& SUBROUTINE DIGFLT !----------------------------------------------------------------- INCLUDE "EXCHM.h" INCLUDE "parmeta.f90" INCLUDE "parm.tbl.f90" INCLUDE "mpp.h" INCLUDE "mpif.h" #include "sp.h" !----------------------------------------------------------------------- PARAMETER & (PIQ=3.141592654,LP1=LM+1,JAM=6+2*(JM-10)) !----------------------------------------------------------------------- ! NTIM IS A HALF-SPAN IN TIME STEP UNITS ! ACTUAL TIME IS NTIM*DT PARAMETER & (NTIM=110) PARAMETER & (CM1=2937.4,CM2=4.9283,CM3=23.5518,EPS=0.622) !----------------------------------------------------------------------- LOGICAL :: & RUN,FIRST,RESTRT,SIGMA,NEST !----------------------------------------------------------------------- REAL :: & PDQ(IDIM1:IDIM2,JDIM1:JDIM2) & ,TQ(IDIM1:IDIM2,JDIM1:JDIM2,LM) & ,UQ(IDIM1:IDIM2,JDIM1:JDIM2,LM),VQ(IDIM1:IDIM2,JDIM1:JDIM2,LM) & ,RELHUM(IDIM1:IDIM2,JDIM1:JDIM2),PDSL0(IDIM1:IDIM2,JDIM1:JDIM2) !----------------------------------------------------------------------- INCLUDE "COMM_CTLBLK.f90" !----------------------------------------------------------------------- INCLUDE "COMM_MASKS.f90" !----------------------------------------------------------------------- INCLUDE "COMM_DYNAM.f90" !----------------------------------------------------------------------- INCLUDE "COMM_CONTIN.f90" !----------------------------------------------------------------------- INCLUDE "COMM_VRBLS.f90" !----------------------------------------------------------------------- INCLUDE "COMM_PVRBLS.f90" !----------------------------------------------------------------------- DATA KNT/0/,IUNRH/51/,IUNDF/52/ !----------------------------------------------------------------------- !&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& !----------------------------------------------------------------------- NBOCO=0 REWIND IUNRH REWIND IUNDF !*** !*** COMPUTE AND WRITE OUT THE RELATIVE HUMIDITY !*** DO J=MYJS,MYJE DO I=MYIS,MYIE PDSL0(I,J)=RES(I,J)*PD(I,J) ENDDO ENDDO DO L=1,LM CALL ZERO2(RELHUM) DO J=MYJS,MYJE DO I=MYIS,MYIE IF(HTM(I,J,L) > 0.5)THEN CLOGES=-CM1/T(I,J,L)-CM2*ALOG10(T(I,J,L))+CM3 ESE=10.**(CLOGES+2.) PRS=AETA(L)*PDSL0(I,J)+PT QIJ=Q(I,J,L) E=PRS*QIJ/(EPS*(1.-QIJ)+QIJ) RELHUM(I,J)=AMIN1(E/ESE,0.98) ELSE RELHUM(I,J)=0. ENDIF ENDDO ENDDO !*** SMOOTH THE INITIAL RH FIELD THEN SAVE IT ! ccccc CALL FILT25(RELHUM,HTM(IDIM1,JDIM1,L),5,l,0) DO J=MYJS,MYJE DO I=MYIS,MYIE RELHUM(I,J)=AMIN1(RELHUM(I,J),0.98) RELHUM(I,J)=AMAX1(RELHUM(I,J),0.02) ENDDO ENDDO !*** !*** ASSEMBLE GLOBAL RELHUM !*** CALL LOC2GLB(RELHUM,TEMP1) IF(MYPE == 0)WRITE(IUNRH)TEMP1 !*** SMOOTH THE INITIAL TEMPERATURE FIELD BEFORE EXECUTING !*** THE DIGITAL FILTER ! ccccc CALL FILT25(T(IDIM1,JDIM1,L),HTM(IDIM1,JDIM1,L),5,l,1) ENDDO !*** SAVE CURRENT PROG VARIABLES IN WORK FILE FOR LATER USE CALL LOC2GLB(PD,TEMP1) IF(MYPE == 0)WRITE(IUNDF)TEMP1 DO L=1,LM CALL LOC2GLB(T(IDIM1,JDIM1,L),TEMP1) IF(MYPE == 0)WRITE(IUNDF)TEMP1 CALL LOC2GLB(U(IDIM1,JDIM1,L),TEMP1) IF(MYPE == 0)WRITE(IUNDF)TEMP1 CALL LOC2GLB(V(IDIM1,JDIM1,L),TEMP1) IF(MYPE == 0)WRITE(IUNDF)TEMP1 ENDDO !*** SET UP ARRAYS TO HOLD SUMS FOR DIGITAL FILTER DO J=MYJS,MYJE DO I=MYIS,MYIE PDQ(I,J)=PD(I,J) ENDDO ENDDO DO L=1,LM DO J=MYJS,MYJE DO I=MYIS,MYIE TQ(I,J,L)=T(I,J,L) UQ(I,J,L)=U(I,J,L) VQ(I,J,L)=V(I,J,L) ENDDO ENDDO ENDDO HSPAN=FLOAT(NTIM)*DT/3600. IF(MYPE == 0)WRITE(6,100) HSPAN 100 FORMAT(' ','INITIALIZATION CALLED WITH HALF-SPAN',F5.1,' HOURS') CALL MPI_BARRIER(MPI_COMM_COMP,IRTN) !----------------------------------------------------------------------- ! RUN THE FORECAST MODEL FORWARD AND BACKWARD (DIGITAL FILTER) !***** FIRST, FORWARD INTEGRATION ************* !----------------------------------------------------------------------- ! CALCULATE NORM NEEDED FOR LANCZOS WINDOW ! 'TETC' DEFINES CUT-OFF FREQUENCY FOR THE WINDOW ! (FACTOR 2 APPEARS ONLY TO SHOW GENERAL FORMULA) ! (IT SHOULD BE TETC=PIQ/FLOAT(NTIM)) ! 'NTIM' IS A NUMBER OF TIME STEPS IN A HALF-SPAN ! 'TIME' CORRESPONDS TO NUMBER OF TIME STEPS OF A SPAN TIME=2.*FLOAT(NTIM) QNT=FLOAT(NTIM+1) TETC=2.*PIQ/TIME CSTT=TETC/PIQ SQ=CSTT DO NT=1,NTIM FNT=FLOAT(NT) AS1=PIQ*FNT/QNT AS2=FNT*TETC ASS1=SIN(AS1)/AS1 ASS2=SIN(AS2)/AS2 SQ=SQ+2.*CSTT*ASS1*ASS2 ENDDO !---------------------------------------------------------------- !*** NORMALIZATION OF THE WINDOW FUNCTION CSTT=CSTT/SQ DO J=MYJS,MYJE DO I=MYIS,MYIE PDQ(I,J)=PDQ(I,J)*CSTT ENDDO ENDDO DO L=1,LM DO J=MYJS,MYJE DO I=MYIS,MYIE TQ(I,J,L)=TQ(I,J,L)*CSTT UQ(I,J,L)=UQ(I,J,L)*CSTT VQ(I,J,L)=VQ(I,J,L)*CSTT ENDDO ENDDO ENDDO !---------------------------------------------------------------- !-------------------------------------------------------------------- NTSD=0 !----------------------------------------------------------------------- !********ENTRY INTO THE TIME LOOP*************************************** 2010 NTSD=NTSD+1 KNT=KNT+1 IF(MYPE == 0)WRITE(6,2015)NTSD,NTIM 2015 FORMAT(' NTSD=',I5,' NTSTM=',I4) !*********************************************************************** !*** DIVERGENCE AND HORIZONTAL PART OF THE OMEGA-ALPHA TERM IF(NTSD > 1)CALL EXCH(T,LM,U,LM,V,LM,2,2) !Exchange T, U, and V CALL DIVHOA !----------------------------------------------------------------------- !*** PRESS. TEND.,ETA DOT & VERTICAL OMEGA-ALPHA CALL PDTEDT !----------------------------------------------------------------------- !*** VERTICAL ADVECTION IF(MOD(NTSD-1,IDTAD) == 0)THEN CALL EXCH(ETADT,LM-1,1,1) CALL VTADVF CALL EXCH(T,LM,U,LM,V,LM,Q2,LM,1,1) ENDIF !----------------------------------------------------------------------- !*** UPDATE SURFACE PRESSURE (MINUS PTOP) CALL PDNEW !----------------------------------------------------------------------- !*** UPDATE H BOUNDARY POINTS IF(MOD(NTSD,IDTAD) == 0)THEN CALL EXCH(T,LM,Q2,LM,1,1) ENDIF CALL EXCH(PD,1,1,1) CALL BOCOHF !----------------------------------------------------------------------- !*** PRESSURE GRADIENT AND CORIOLIS FORCE TERMS CALL EXCH(PD,1,T,LM,2,2) !Exchange PD and T CALL PGCOR CALL EXCH(PDSL,1,5,5) !----------------------------------------------------------------------- !*** UPDATE V BOUNDARY POINTS CALL EXCH(U,LM,V,LM,1,1) !Exchange U and V CALL BOCOV !----------------------------------------------------------------------- !*** HORIZONTAL ADVECTION IF(MOD(NTSD,IDTAD) == 0)THEN CALL EXCH(T,LM,U,LM,V,LM,4,4) !Exchange T, U, and V CALL EXCH(Q2,LM,5,5) CALL HZADV CALL EXCH(U,LM,V,LM,2,2) !Exchange U and V ENDIF !----------------------------------------------------------------------- !*** CALCULATE WINDOW PARAMETER BNT=FLOAT(NTSD) BS1=BNT*PIQ/QNT BS2=BNT*TETC BSS1= SIN(BS1)/BS1 BSS2= SIN(BS2)/BS2 HNTSD=CSTT*BSS1*BSS2 !+++ HNTSD=HNTSD*1.042 DO L=1,LM DO J=MYJS,MYJE DO I=MYIS,MYIE TQ(I,J,L)=TQ(I,J,L)+HNTSD*T(I,J,L) UQ(I,J,L)=UQ(I,J,L)+HNTSD*U(I,J,L) VQ(I,J,L)=VQ(I,J,L)+HNTSD*V(I,J,L) ENDDO ENDDO ENDDO DO J=MYJS,MYJE DO I=MYIS,MYIE PDQ(I,J)=PDQ(I,J)+HNTSD*PD(I,J) ENDDO ENDDO !----------------------------------------------------------------------- IF(NTSD == NTIM) GO TO 2013 !*********************************************************************** GO TO 2010 !********EXIT FROM THE TIME LOOP**************************************** 2013 CONTINUE CALL MPI_BARRIER(MPI_COMM_COMP,IRTN) !----------------------------------------------------------------------- ! NOW BACKWARD INTEGRATION, STARTING FROM THE INITIAL TIME !----------------------------------------------------------------------- !*** CHANGE (SIGN ONLY OF) IMPORTANT TIME CONSTANTS DT =-DT CPGFV=-CPGFV EN =-EN ENT =-ENT F4D =-F4D F4Q =-F4Q EF4T =-EF4T DO JK=1,JAM EM (JK)=-EM (JK) EMT(JK)=-EMT(JK) ENDDO DO L=1,LM F4Q2(L)=-F4Q2(L) ENDDO DO J=MYJS,MYJE DO I=MYIS,MYIE WPDAR(I,J)=-WPDAR(I,J) CPGFU(I,J)=-CPGFU(I,J) CURV (I,J)=-CURV (I,J) FCP (I,J)=-FCP (I,J) FAD (I,J)=-FAD (I,J) F (I,J)=-F (I,J) ENDDO ENDDO CALL EXCH(WPDAR,1,CPGFU,1,2,2) CALL EXCH(CURV,1,FCP,1,2,2) CALL EXCH(FAD,1,F,1,2,2) !------------------- END OF CHANGE ------------------------- ! DEFINE INITIAL DATA FOR BACKWARD INTEGRATION ! (PDQ,TQ,UQ,VQ ARE SUMS NEEDED FOR DIGITAL FILTER) IF(MYPE == 0)THEN REWIND IUNDF READ(IUNDF)TEMP1 ENDIF CALL DSTRB(TEMP1,PD,1,1,1) DO L=1,LM IF(MYPE == 0)THEN READ(IUNDF)TEMP1 READ(IUNDF)TEMP2 READ(IUNDF)TEMP3 ENDIF CALL DSTRB(TEMP1,T,1,LM,L) CALL DSTRB(TEMP2,U,1,LM,L) CALL DSTRB(TEMP3,V,1,LM,L) ENDDO CALL EXCH(T,LM,U,LM,V,LM,5,5) !Exchange T, U, and V CALL EXCH(PD,1,5,5) !-------------------------------------------------------------------- NTSD=0 FIRST= .TRUE. TSPH=-3600./DT IF(MYPE == 0)THEN IF( .NOT. NEST)THEN REWIND NBC READ(NBC) READ(NBC)BCHR ELSE READ(NBC,REC=2)BCHR ENDIF ENDIF CALL MPI_BCAST(BCHR,1,MPI_REAL,0,MPI_COMM_COMP,IRTN) NBOCO=INT(BCHR*TSPH+0.5) !----------------------------------------------------------------------- !********ENTRY INTO THE TIME LOOP*************************************** 2020 NTSD=NTSD+1 KNT=KNT+1 IF(MYPE == 0)WRITE(6,2015)NTSD,NTIM !*********************************************************************** !*** DIVERGENCE AND HORIZONTAL PART OF THE OMEGA-ALPHA TERM IF(NTSD > 1)CALL EXCH(T,LM,U,LM,V,LM,2,2) !Exchange T, U, and V CALL DIVHOA !----------------------------------------------------------------------- !*** PRESS. TEND.,ETA DOT & VERTICAL OMEGA-ALPHA CALL PDTEDT !----------------------------------------------------------------------- !*** VERTICAL ADVECTION IF(MOD(NTSD-1,IDTAD) == 0)THEN CALL EXCH(ETADT,LM-1,1,1) CALL VTADVF CALL EXCH(T,LM,U,LM,V,LM,Q2,LM,1,1) !Exchange T, U, and V ENDIF !----------------------------------------------------------------------- !*** UPDATE SURFACE PRESSURE (MINUS PTOP) CALL PDNEW !----------------------------------------------------------------------- !*** UPDATE H BOUNDARY POINTS IF(MOD(NTSD,IDTAD) == 0)THEN CALL EXCH(T,LM,Q2,LM,1,1) ENDIF CALL EXCH(PD,1,1,1) CALL BOCOHF !----------------------------------------------------------------------- !*** PRESSURE GRADIENT AND CORIOLIS FORCE TERMS CALL EXCH(PD,1,T,LM,2,2) !Exchange PD and T CALL PGCOR CALL EXCH(PDSL,1,5,5) !----------------------------------------------------------------------- !*** UPDATE V BOUNDARY POINTS CALL EXCH(U,LM,V,LM,1,1) !Exchange U and V CALL BOCOV !----------------------------------------------------------------------- !*** HORIZONTAL ADVECTION IF(MOD(NTSD,IDTAD) == 0)THEN CALL EXCH(T,LM,U,LM,V,LM,4,4) !Exchange T, U, and V CALL EXCH(Q2,LM,5,5) CALL HZADV CALL EXCH(U,LM,V,LM,2,2) !Exchange U and V ENDIF !----------------------------------------------------------------------- !*** CALCULATE WINDOW PARAMETER BNT=FLOAT(NTSD) BS1=BNT*PIQ/QNT BS2=BNT*TETC BSS1= SIN(BS1)/BS1 BSS2= SIN(BS2)/BS2 HNTSD=CSTT*BSS1*BSS2 !+++ HNTSD=HNTSD*1.042 DO L=1,LM DO J=MYJS,MYJE DO I=MYIS,MYIE TQ(I,J,L)=TQ(I,J,L)+HNTSD*T(I,J,L) UQ(I,J,L)=UQ(I,J,L)+HNTSD*U(I,J,L) VQ(I,J,L)=VQ(I,J,L)+HNTSD*V(I,J,L) ENDDO ENDDO ENDDO DO J=MYJS,MYJE DO I=MYIS,MYIE PDQ(I,J)=PDQ(I,J)+HNTSD*PD(I,J) ENDDO ENDDO !----------------------------------------------------------------------- IF(NTSD == NTIM) GO TO 2022 !*********************************************************************** GO TO 2020 !********EXIT FROM THE TIME LOOP**************************************** 2022 CONTINUE !----------------------------------------------------------------------- !*** CHANGE BACK (SIGN ONLY) IMPORTANT TIME CONSTANTS DT =-DT CPGFV=-CPGFV EN =-EN ENT =-ENT F4D =-F4D F4Q =-F4Q EF4T =-EF4T DO JK=1,JAM EM (JK)=-EM (JK) EMT(JK)=-EMT(JK) ENDDO DO L=1,LM F4Q2(L)=-F4Q2(L) ENDDO DO J=MYJS,MYJE DO I=MYIS,MYIE WPDAR(I,J)=-WPDAR(I,J) CPGFU(I,J)=-CPGFU(I,J) CURV (I,J)=-CURV (I,J) FCP (I,J)=-FCP (I,J) FAD (I,J)=-FAD (I,J) F (I,J)=-F (I,J) ENDDO ENDDO CALL EXCH(WPDAR,1,CPGFU,1,2,2) CALL EXCH(CURV,1,FCP,1,2,2) CALL EXCH(FAD,1,F,1,2,2) !*** UPDATE INITIALIZED PROGNOSTIC VARIALBES DO J=MYJS,MYJE DO I=MYIS,MYIE PD(I,J)=PDQ(I,J) ENDDO ENDDO DO L=1,LM DO J=MYJS,MYJE DO I=MYIS,MYIE T(I,J,L)=TQ(I,J,L) U(I,J,L)=UQ(I,J,L) V(I,J,L)=VQ(I,J,L) ENDDO ENDDO ENDDO !*** !*** RETRIEVE THE INITIAL RELATIVE HUMIDITY AND COMPUTE Q !*** SO AS TO MAINTAIN THE RH GIVEN THE NEW TEMPERATURES !*** DO J=MYJS,MYJE DO I=MYIS,MYIE PDSL0(I,J)=RES(I,J)*PD(I,J) ENDDO ENDDO IF(MYPE == 0)REWIND IUNRH DO L=1,LM IF(MYPE == 0)READ(IUNRH)TEMP1 CALL DSTRB(TEMP1,RELHUM,1,1,1) DO J=MYJS,MYJE DO I=MYIS,MYIE IF(HTM(I,J,L) > 0.5)THEN CLOGES=-CM1/T(I,J,L)-CM2*ALOG10(T(I,J,L))+CM3 ESE=10.**(CLOGES+2.) PRS=AETA(L)*PDSL0(I,J)+PT E=RELHUM(I,J)*ESE Q(I,J,L)=EPS*E/(PRS+E*(EPS-1.)) ENDIF ENDDO ENDDO ENDDO CALL EXCH(T,LM,Q,LM,U,LM,V,LM,4,4) !Exchange T, Q, U, and V CALL EXCH(PD,1,5,5) ! RETURN BC FILE TO START FORECAST NTSD=0 IF(MYPE == 0)THEN IF( .NOT. NEST)THEN REWIND NBC READ(NBC) READ(NBC)BCHR ELSE READ(NBC,REC=2)BCHR ENDIF ENDIF CALL MPI_BCAST(BCHR,1,MPI_REAL,0,MPI_COMM_COMP,IRTN) NBOCO=INT(BCHR*TSPH+0.5) !----------------------------------------------------------- !------------------- END OF CHANGE ------------------------- !----------------------------------------------------------- RETURN END SUBROUTINE DIGFLT