C*********************************************************************** C Portions of Models-3/CMAQ software were developed or based on * C information from various groups: Federal Government employees, * C contractors working on a United States Government contract, and * C non-Federal sources (including research institutions). These * C research institutions have given the Government permission to * C use, prepare derivative works, and distribute copies of their * C work in Models-3/CMAQ to the public and to permit others to do * C so. EPA therefore grants similar permissions for use of the * C Models-3/CMAQ software, but users are requested to provide copies * C of derivative works to the Government without restrictions as to * C use by others. Users are responsible for acquiring their own * C copies of commercial software associated with Models-3/CMAQ and * C for complying with vendor requirements. Software copyrights by * C the MCNC Environmental Modeling Center are used with their * C permissions subject to the above restrictions. * C*********************************************************************** C RCS file, release, date & time of last delta, author, state, [and locker] C $Header: /project/work/rep/CCTM/src/hdiff/multiscale/hdiff.F,v 1.10 2005/05/27 15:53:36 yoj Exp $ C what(1) key, module and SID; SCCS file; date and time of last delta: C %W% %P% %G% %U% C::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: SUBROUTINE HDIFF ( CGRID, JDATE, JTIME, TSTEP ) C----------------------------------------------------------------------- C Function: C Horizontal diffusion with constant eddy diffusivity - gen. coord. C The process time step is set equal to TSTEP(2). Boundary concentrations C are set using a Dirichlet (no flux) condition C Preconditions: C Dates and times represented YYYYDDD:HHMMSS. C No "skipped" dates and times. All boundary input variables have the C same boundary perimeter structure with a thickness of 1 C CGRID in ppm units or micro-g/m**3, #/m**3 for aerosols C Subroutines and functions called: C INIT3, TIME2SEC, SEC2TIME, CGRID_MAP, NEXTIME, RHO_J, C HCDIFF3D C Revision history: C Jeff - 5 Nov 97, 1 Jan 98 C DWB - 1 Feb 98, use simple B/C (no conc gradient at domain boundary) C David Wong Sep. 1998 C -- parallelized the code C -- removed the intermediate constant CRHOJ_Q and placed the answer of C the calculation directly into CGRID. Removed the next immediate C loop completely. C David Wong 1/19/99 C -- add a loop_index call C -- change loop index ending point to avoid accessing invalid region. C (reason to do this is to prevent using boundary data from PINTERP, C which sets pseudo-boundary data to 0) C Daewon Byun 10/10/2000 C -- generalized 3d horizontal diffusivity module C -- accomdates 3d hdiff values C 15 Dec 00 J.Young: PE_COMM3 -> Dave Wong's f90 stenex COMM C 6 Aug 01 J.Young: Use HGRD_DEFN C 25 Mar 04 G.Hammond: RK11/RK22 ghost cell updates moved outside main loop; C use explicit boundary arrays for CGRID ghost cells; C use SNL's "swap3d". C 31 Jan 05 J.Young: dyn alloc - establish both horizontal & vertical C domain specifications in one module C----------------------------------------------------------------------- USE GRID_CONF ! horizontal & vertical domain specifications USE CGRID_SPCS ! CGRID species number and offsets USE SUBST_MODULES ! stenex ! USE SUBST_UTIL_MODULE ! stenex ! USE SUBST_COMM_MODULE ! stenex USE SWAP_SANDIA #ifdef snl_timing use timing #endif IMPLICIT NONE C Includes: ! INCLUDE SUBST_HGRD_ID ! horizontal dimensioning parameters ! INCLUDE SUBST_VGRD_ID ! vertical dimensioning parameters INCLUDE SUBST_GC_SPC ! gas chemistry species table INCLUDE SUBST_AE_SPC ! aerosol species table INCLUDE SUBST_NR_SPC ! non-reactive species table INCLUDE SUBST_TR_SPC ! tracer species table INCLUDE SUBST_GC_DIFF ! gas chem diffusion species and map table INCLUDE SUBST_AE_DIFF ! aerosol diffusion species and map table INCLUDE SUBST_NR_DIFF ! non-react diffusion species and map table INCLUDE SUBST_TR_DIFF ! tracer diffusion species and map table INCLUDE SUBST_CONST ! constants INCLUDE SUBST_IOPARMS ! I/O parameters definitions INCLUDE SUBST_IOFDESC ! file header data structure INCLUDE SUBST_IODECL ! I/O definitions and declarations ! INCLUDE SUBST_COORD_ID ! coordinate & domain definitions (req IOPARMS) C Arguments: ! REAL CGRID( NCOLS,NROWS,NLAYS,* ) ! concentrations ! REAL :: CGRID( :,:,:,: ) REAL, POINTER :: CGRID( :,:,:,: ) INTEGER JDATE ! current model date, coded YYYYDDD INTEGER JTIME ! current model time, coded HHMMSS INTEGER TSTEP( 2 ) ! time step vector (HHMMSS) ! TSTEP(1) = local output step ! TSTEP(2) = sciproc sync. step (chem) C External Functions not declared in IODECL3.EXT: INTEGER, EXTERNAL :: SEC2TIME, TIME2SEC, GETEFILE C Parameters: ! REAL, PARAMETER :: CMIN = 0.0 ! REAL, PARAMETER :: TH1 = 1.0 / 3.0 ! REAL, PARAMETER :: TH4 = 4.0 * TH1 C Advected species dimension INTEGER, PARAMETER :: N_SPC_DIFF = N_GC_DIFF & + N_AE_DIFF & + N_NR_DIFF & + N_TR_DIFF ! & + 1 ! diff RHO_J C File Variables: REAL CONC ( 0:NCOLS+1,0:NROWS+1 ) ! conc working array REAL RHOJ ( 0:NCOLS+1,0:NROWS+1,NLAYS ) ! density X Jacobian CHARACTER( 8 ), SAVE :: COMMSTR ! for both CONC and RHOJ ! REAL RRHOJ2( 0:NCOLS+1,0:NROWS+1,NLAYS ) ! reciprocal density ! X Jacobian 2nd time REAL RK11 ( NCOLS+1,NROWS+1,NLAYS ) ! RHOJ at x1 cell face ! reused as 11 eddy diff. factor REAL RK22 ( NCOLS+1,NROWS+1,NLAYS ) ! RHOJ at x2 cell face ! reused as 22 eddy diff. factor REAL K11BAR3D ( NCOLS+1,NROWS+1,NLAYS ) ! ave. Cx11 eddy diff REAL K22BAR3D ( NCOLS+1,NROWS+1,NLAYS ) ! ave. Cx22 eddy diff REAL DT ! diffusion time step REAL CRHOJ_Q ! intermediate, coupled conc. C New File Variables: integer m real bound_south( ncols,nlays,n_spc_diff ) real bound_north( ncols,nlays,n_spc_diff ) real bound_west ( nrows,nlays,n_spc_diff ) real bound_east ( nrows,nlays,n_spc_diff ) real buf_ns ( ncols*nlays*n_spc_diff ) real buf_ew ( nrows*nlays*n_spc_diff ) C Local Variables: CHARACTER( 16 ) :: PNAME = 'HDIFF' LOGICAL, SAVE :: FIRSTIME = .TRUE. REAL DX1 ! dx1 (meters) REAL DX2 ! dx2 (meters) REAL, SAVE :: RDX1S ! reciprocal dx1*dx1 REAL, SAVE :: RDX2S ! reciprocal dx2*dx2 REAL DTDX1S ! dt/dx1**2 REAL DTDX2S ! dt/dx2**2 REAL DTSEC ! model time step in seconds INTEGER NSTEPS ! diffusion time steps INTEGER STEP ! FIX dt INTEGER FDATE ! interpolation date INTEGER FTIME ! interpolation time ! INTEGER NDATE ! interpolation date, 2nd step ! INTEGER NTIME ! interpolation time, 2nd step INTEGER, SAVE :: DIFF_MAP( N_SPC_DIFF ) ! global diff map to CGRID INTEGER, SAVE :: LOGDEV ! CHARACTER( 16 ) :: OUTNAME = 'HDIFF_NEG_STATS' ! INTEGER, SAVE :: OUTDEV ! LOGICAL :: RDONLY = .FALSE. ! LOGICAL :: FMTTED = .TRUE. ! LOGICAL EXCNT ! write header flag INTEGER C, R, L, S, V, N ! loop counters INTEGER D2C CHARACTER( 96 ) :: XMSG = ' ' INTEGER MY_TEMP INTEGER, SAVE :: STARTROW, ENDROW INTEGER, SAVE :: STARTCOL, ENDCOL INTERFACE SUBROUTINE RHO_J ( JDATE, JTIME, TSTEP, RHOJ ) IMPLICIT NONE INTEGER, INTENT( IN ) :: JDATE, JTIME, TSTEP( 2 ) REAL, INTENT( OUT ) :: RHOJ( :,:,: ) END SUBROUTINE RHO_J SUBROUTINE HCDIFF3D ( JDATE, JTIME, K11BAR, K22BAR, DT ) IMPLICIT NONE INTEGER, INTENT( IN ) :: JDATE, JTIME REAL, INTENT( OUT ) :: K11BAR( :,:,: ), K22BAR( :,:,: ) REAL, INTENT( OUT ) :: DT END SUBROUTINE HCDIFF3D END INTERFACE C----------------------------------------------------------------------- IF ( FIRSTIME ) THEN FIRSTIME = .FALSE. LOGDEV = INIT3() C Get dx1 from COORD include file IF ( GDTYP_GD .EQ. LATGRD3 ) THEN DX1 = DG2M * XCELL_GD & * COS( PI180*( YORIG_GD + YCELL_GD*FLOAT( GL_NROWS/2 ))) ! in m. DX2 = DG2M * YCELL_GD ! in m. ELSE DX1 = XCELL_GD ! in m. DX2 = YCELL_GD ! in m. END IF RDX1S = 1.0 / ( DX1 * DX1 ) RDX2S = 1.0 / ( DX2 * DX2 ) C Get CGRID offsets CALL CGRID_MAP( NSPCSD, GC_STRT, AE_STRT, NR_STRT, TR_STRT ) C Create global map to CGRID S = 0 DO V = 1, N_GC_DIFF S = S + 1 DIFF_MAP( S ) = GC_STRT - 1 + GC_DIFF_MAP( V ) END DO DO V = 1, N_AE_DIFF S = S + 1 DIFF_MAP( S ) = AE_STRT - 1 + AE_DIFF_MAP( V ) END DO DO V = 1, N_NR_DIFF S = S + 1 DIFF_MAP( S ) = NR_STRT - 1 + NR_DIFF_MAP( V ) END DO DO V = 1, N_TR_DIFF S = S + 1 DIFF_MAP( S ) = TR_STRT - 1 + TR_DIFF_MAP( V ) END DO ! DIFF_MAP( N_SPC_DIFF ) = N_GC_SPCD ! RHO_J C open output ASCII file ! OUTDEV = GETEFILE ( OUTNAME, RDONLY, FMTTED, PNAME ) CALL SUBST_LOOP_INDEX ( 'R', 1, NROWS, 1, MY_TEMP, STARTROW, ENDROW ) CALL SUBST_LOOP_INDEX ( 'C', 1, NCOLS, 1, MY_TEMP, STARTCOL, ENDCOL ) END IF ! if firstime DTSEC = FLOAT( TIME2SEC( TSTEP( 2 ) ) ) FDATE = JDATE FTIME = JTIME C Get the computational grid ( rho X Jacobian ) for this step CALL RHO_J ( FDATE, FTIME, TSTEP, RHOJ ) C swap3d interface ... C swap3d( 1st value of array to be sent, C 1st value of array to be received, C number of values to send in 1st dimension of array, C number of values to send in 2nd dimension of array, C number of values to send in 3rd dimension of array, C size of 1st dimension in arrays, C size of 2nd dimension in arrays, C direction to receive from: NSEW) #ifdef snl_timing call start_timing( hdiff_int, comm_int, 2 ) #endif call swap3d( rhoj( 1,my_nrows,1 ), rhoj( 1,0,1 ), my_ncols, 1, nlays, & my_ncols+2, my_nrows+2, SOUTH ) call swap3d( rhoj( my_ncols,1,1 ), rhoj( 0,1,1 ), 1, my_nrows, nlays, & my_ncols+2, my_nrows+2, WEST ) #ifdef snl_timing call stop_timing( hdiff_int, comm_int ) #endif C get face values for RHOJ (assumes dx1 = dx2) DO L = 1, NLAYS DO R = STARTROW, ENDROW ! DO R = 1, NROWS + 1 DO C = STARTCOL, ENDCOL ! DO C = 1, NCOLS + 1 RK11( C,R,L ) = 0.5 * ( RHOJ( C,R,L ) + RHOJ( C-1,R, L ) ) RK22( C,R,L ) = 0.5 * ( RHOJ( C,R,L ) + RHOJ( C, R-1,L ) ) END DO END DO END DO C Get the computational grid for the next time step C*dwb Do not need to read next time step C*dwb CGRID includes RHOJ2 implicitly !!! (since coupled) ! NDATE = FDATE ! NTIME = FTIME ! CALL NEXTIME( NDATE, NTIME, TSTEP( 3 ) ) ! CALL RHO_J ( NDATE, NTIME, TSTEP( 3 ), RRHOJ2 ) ! DO L = 1, NLAYS ! DO R = 1, MY_NROWS + 1 ! DO C = 1, MY_NCOLS + 1 ! RRHOJ2( C,R,L ) = 1.0 / RRHOJ2( C,R,L ) ! END DO ! END DO ! END DO C Do the gridded computation for horizontal diffusion C Get the contravariant eddy diffusivities CALL HCDIFF3D ( FDATE, FTIME, K11BAR3D, K22BAR3D, DT ) C get number of steps based on eddy time NSTEPS = INT ( DTSEC / DT ) + 1 DT = DTSEC / FLOAT( NSTEPS ) WRITE( LOGDEV,1005 ) DT, NSTEPS DTDX1S = DT * RDX1S DTDX2S = DT * RDX2S DO L = 1, NLAYS DO R = STARTROW, ENDROW ! DO R = 1, NROWS + 1 DO C = STARTCOL, ENDCOL ! DO C = 1, NCOLS + 1 RK11( C,R,L ) = RK11( C,R,L ) * K11BAR3D( C,R,L ) RK22( C,R,L ) = RK22( C,R,L ) * K22BAR3D( C,R,L ) END DO END DO END DO #ifdef snl_timing call start_timing( hdiff_int, comm_int, 2 ) #endif call swap3d( rk11( 1,1,1 ), rk11( my_ncols+1,1,1 ), 1, my_nrows, nlays, & my_ncols+1, my_nrows+1, EAST ) call swap3d( rk22( 1,1,1 ), rk22( 1,my_nrows+1,1 ), my_ncols, 1, nlays, & my_ncols+1, my_nrows+1, NORTH ) #ifdef snl_timing call stop_timing( hdiff_int, comm_int ) #endif C Create and fill four 3d boundary arrays with cgrid/rhoj C send to SOUTH, recv from NORTH, fill bound_north C send to NORTH, recv from SOUTH, fill bound_south C send to EAST, recv from WEST, fill bound_west C send to WEST, recv from EAST, fill bound_east m = 0 do s = 1, n_spc_diff d2c = diff_map( s ) do l = 1, nlays do c = 1, my_ncols bound_south( c,l,s ) = cgrid( c,1,l,d2c ) / rhoj( c,1,l ) bound_north( c,l,s ) = cgrid( c,my_nrows,l,d2c ) / rhoj( c,my_nrows,l ) m = m + 1 buf_ns( m ) = bound_north( c,l,s ) end do end do end do #ifdef snl_timing call start_timing( hdiff_int, comm_int, 2 ) #endif call swap3d( bound_south, bound_north, my_ncols, nlays, n_spc_diff, & my_ncols, nlays, NORTH ) call swap3d( buf_ns, bound_south, my_ncols, nlays, n_spc_diff, & my_ncols, nlays, SOUTH ) #ifdef snl_timing call stop_timing( hdiff_int, comm_int ) #endif m = 0 do s = 1, n_spc_diff d2c = diff_map( s ) do l = 1, nlays do r = 1, my_nrows bound_west( r,l,s ) = cgrid( 1,r,l,d2c ) / rhoj( 1,r,l ) bound_east( r,l,s ) = cgrid( my_ncols,r,l,d2c ) / rhoj( my_ncols,r,l ) m = m + 1 buf_ew( m ) = bound_east( r,l,s ) end do end do end do #ifdef snl_timing call start_timing( hdiff_int, comm_int, 2 ) #endif call swap3d( bound_west, bound_east, my_nrows, nlays, n_spc_diff, & my_nrows, nlays, EAST ) call swap3d( buf_ew, bound_west, my_nrows, nlays, n_spc_diff, & my_nrows, nlays, WEST ) #ifdef snl_timing call stop_timing( hdiff_int, comm_int ) #endif C Loop over species, layers, nsteps DO 366 S = 1, N_SPC_DIFF D2C = DIFF_MAP( S ) DO 355 L = 1, NLAYS C Load working array (CGRID is coupled, CONC is mixing ratio) DO R = 1, MY_NROWS DO C = 1, MY_NCOLS CONC( C,R ) = CGRID( C,R,L,D2C ) / RHOJ( C,R,L ) END DO END DO DO 344 N = 1, NSTEPS C Fill 4 boundaries of working array CONC: south, north, west, east do c = 1, my_ncols conc( c,0 ) = bound_south( c,l,s ) conc( c,my_nrows+1 ) = bound_north( c,l,s ) end do do r = 1, my_nrows conc( 0,r ) = bound_west( r,l,s ) conc( my_ncols+1,r ) = bound_east( r,l,s ) end do C Update CGRID DO R = 1, MY_NROWS DO C = 1, MY_NCOLS CGRID( C,R,L,D2C ) = RHOJ( C,R,L ) * CONC( C,R ) & + DTDX1S & * ( RK11( C+1,R,L ) & * ( CONC( C+1,R ) - CONC( C,R ) ) & - RK11( C,R,L ) & * ( CONC( C,R ) - CONC( C-1,R ) ) ) & + DTDX2S & * ( RK22( C,R+1,L ) & * ( CONC( C,R+1 ) - CONC( C,R ) ) & - RK22( C,R,L ) & * ( CONC( C,R ) - CONC( C,R-1 ) ) ) END DO END DO 344 CONTINUE 355 CONTINUE 366 CONTINUE RETURN 1001 FORMAT( 5X, 'Negative concentrations reset to', 1PE11.3 ) 1003 FORMAT( 1X, 4I7, 9X, 1PE11.3) 1005 FORMAT( / 5X, 'H-eddy DT & integration steps: ', 1PE15.7, I8 ) END