!WRF:MODEL_LAYER:DYNAMICS ! MODULE module_advect_em USE module_bc USE module_model_constants USE module_wrf_error CONTAINS SUBROUTINE mass_flux_divergence ( field, field_old, tendency, & ru, rv, rom, & mut, config_flags, & msfu, msfv, msft, & fzm, fzp, & rdx, rdy, rdzw, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) IMPLICIT NONE ! Input data TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: field, & field_old, & ru, & rv, & rom REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mut REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: msfu, & msfv, & msft REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fzm, & fzp, & rdzw REAL , INTENT(IN ) :: rdx, & rdy ! Local data INTEGER :: i, j, k, itf, jtf, ktf INTEGER :: i_start, i_end, j_start, j_end INTEGER :: imin, imax, jmin, jmax REAL :: mrdx, mrdy, ub, vb, uw, vw REAL , DIMENSION(its:ite,kts:kte) :: vflux LOGICAL :: specified !--------------- horizontal flux specified = .false. if(config_flags%specified .or. config_flags%nested) specified = .true. ktf=MIN(kte,kde-1) i_start = its i_end = MIN(ite,ide-1) j_start = jts j_end = MIN(jte,jde-1) DO j = j_start, j_end DO k = kts, ktf DO i = i_start, i_end mrdx=msft(i,j)*rdx tendency(i,k,j)=tendency(i,k,j)-mrdx*0.5 & *(ru(i+1,k,j)*(field(i+1,k,j)+field(i ,k,j)) & -ru(i ,k,j)*(field(i ,k,j)+field(i-1,k,j))) ENDDO ENDDO ENDDO DO j = j_start, j_end DO k = kts, ktf DO i = i_start, i_end mrdy=msft(i,j)*rdy tendency(i,k,j)=tendency(i,k,j) -mrdy*0.5 & *(rv(i,k,j+1)*(field(i,k,j+1)+field(i,k,j )) & -rv(i,k,j )*(field(i,k,j )+field(i,k,j-1))) ENDDO ENDDO ENDDO !---------------- vertical flux divergence DO i = i_start, i_end vflux(i,kts)=0. vflux(i,kte)=0. ENDDO DO j = j_start, j_end DO k = kts+1, ktf DO i = i_start, i_end vflux(i,k)=rom(i,k,j)*(fzm(k)*field(i,k,j)+fzp(k)*field(i,k-1,j)) ENDDO ENDDO DO k = kts, ktf DO i = i_start, i_end tendency(i,k,j)=tendency(i,k,j)-rdzw(k)*(vflux(i,k+1)-vflux(i,k)) ENDDO ENDDO ENDDO END SUBROUTINE mass_flux_divergence !------------------------------------------------------------------------------- SUBROUTINE advect_u ( u, u_old, tendency, & ru, rv, rom, & mut, config_flags, & msfu, msfv, msft, & fzm, fzp, & rdx, rdy, rdzw, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) IMPLICIT NONE ! Input data TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: u, & u_old, & ru, & rv, & rom REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mut REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: msfu, & msfv, & msft REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fzm, & fzp, & rdzw REAL , INTENT(IN ) :: rdx, & rdy ! Local data INTEGER :: i, j, k, itf, jtf, ktf INTEGER :: i_start, i_end, j_start, j_end INTEGER :: i_start_f, i_end_f, j_start_f, j_end_f INTEGER :: jmin, jmax, jp, jm, imin, imax, im, ip INTEGER :: jp1, jp0, jtmp INTEGER :: horz_order, vert_order REAL :: mrdx, mrdy, ub, vb, uw, vw, dvm, dvp REAL , DIMENSION(its:ite, kts:kte) :: vflux REAL, DIMENSION( its-1:ite+1, kts:kte ) :: fqx REAL, DIMENSION( its:ite, kts:kte, 2) :: fqy LOGICAL :: degrade_xs, degrade_ys LOGICAL :: degrade_xe, degrade_ye ! definition of flux operators, 3rd, 4th, 5th or 6th order REAL :: flux3, flux4, flux5, flux6 REAL :: q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua, vel flux4(q_im2, q_im1, q_i, q_ip1, ua) = & ( 7.*(q_i + q_im1) - (q_ip1 + q_im2) )/12.0 flux3(q_im2, q_im1, q_i, q_ip1, ua) = & flux4(q_im2, q_im1, q_i, q_ip1, ua) + & sign(1.,ua)*((q_ip1 - q_im2)-3.*(q_i-q_im1))/12.0 flux6(q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua) = & ( 37.*(q_i+q_im1) - 8.*(q_ip1+q_im2) & +(q_ip2+q_im3) )/60.0 flux5(q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua) = & flux6(q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua) & -sign(1.,ua)*( & (q_ip2-q_im3)-5.*(q_ip1-q_im2)+10.*(q_i-q_im1) )/60.0 LOGICAL :: specified specified = .false. if(config_flags%specified .or. config_flags%nested) specified = .true. ! set order for vertical and horzontal flux operators horz_order = config_flags%h_mom_adv_order vert_order = config_flags%v_mom_adv_order ktf=MIN(kte,kde-1) ! begin with horizontal flux divergence horizontal_order_test : IF( horz_order == 6 ) THEN ! determine boundary mods for flux operators ! We degrade the flux operators from 3rd/4th order ! to second order one gridpoint in from the boundaries for ! all boundary conditions except periodic and symmetry - these ! conditions have boundary zone data fill for correct application ! of the higher order flux stencils degrade_xs = .true. degrade_xe = .true. degrade_ys = .true. degrade_ye = .true. IF( config_flags%periodic_x .or. & config_flags%symmetric_xs .or. & (its > ids+2) ) degrade_xs = .false. IF( config_flags%periodic_x .or. & config_flags%symmetric_xe .or. & (ite < ide-2) ) degrade_xe = .false. IF( config_flags%periodic_y .or. & config_flags%symmetric_ys .or. & (jts > jds+2) ) degrade_ys = .false. IF( config_flags%periodic_y .or. & config_flags%symmetric_ye .or. & (jte < jde-3) ) degrade_ye = .false. !--------------- y - advection first i_start = its i_end = ite IF ( config_flags%open_xs .or. specified ) i_start = MAX(ids+1,its) IF ( config_flags%open_xe .or. specified ) i_end = MIN(ide-1,ite) IF ( config_flags%periodic_x ) i_start = its IF ( config_flags%periodic_x ) i_end = ite j_start = jts j_end = MIN(jte,jde-1) ! higher order flux has a 5 or 7 point stencil, so compute ! bounds so we can switch to second order flux close to the boundary j_start_f = j_start j_end_f = j_end+1 IF(degrade_ys) then j_start = MAX(jts,jds+1) j_start_f = jds+3 ENDIF IF(degrade_ye) then j_end = MIN(jte,jde-2) j_end_f = jde-3 ENDIF ! compute fluxes, 5th or 6th order jp1 = 2 jp0 = 1 j_loop_y_flux_6 : DO j = j_start, j_end+1 IF( (j >= j_start_f ) .and. (j <= j_end_f) ) THEN ! use full stencil DO k=kts,ktf DO i = i_start, i_end vel = 0.5*(rv(i,k,j)+rv(i-1,k,j)) fqy( i, k, jp1 ) = vel*flux6( & u(i,k,j-3), u(i,k,j-2), u(i,k,j-1), & u(i,k,j ), u(i,k,j+1), u(i,k,j+2), vel ) ENDDO ENDDO ! we must be close to some boundary where we need to reduce the order of the stencil ELSE IF ( j == jds+1 ) THEN ! 2nd order flux next to south boundary DO k=kts,ktf DO i = i_start, i_end fqy(i, k, jp1) = 0.25*(rv(i,k,j)+rv(i-1,k,j)) & *(u(i,k,j)+u(i,k,j-1)) ENDDO ENDDO ELSE IF ( j == jds+2 ) THEN ! third of 4th order flux 2 in from south boundary DO k=kts,ktf DO i = i_start, i_end vel = 0.5*(rv(i,k,j)+rv(i-1,k,j)) fqy( i, k, jp1 ) = vel*flux4( & u(i,k,j-2),u(i,k,j-1), u(i,k,j),u(i,k,j+1),vel ) ENDDO ENDDO ELSE IF ( j == jde-1 ) THEN ! 2nd order flux next to north boundary DO k=kts,ktf DO i = i_start, i_end fqy(i, k, jp1) = 0.25*(rv(i,k,j)+rv(i-1,k,j)) & *(u(i,k,j)+u(i,k,j-1)) ENDDO ENDDO ELSE IF ( j == jde-2 ) THEN ! 3rd or 4th order flux 2 in from north boundary DO k=kts,ktf DO i = i_start, i_end vel = 0.5*(rv(i,k,j)+rv(i-1,k,j)) fqy( i, k, jp1 ) = vel*flux4( & u(i,k,j-2),u(i,k,j-1), & u(i,k,j),u(i,k,j+1),vel ) ENDDO ENDDO END IF !stopped ! y flux-divergence into tendency IF(j > j_start) THEN DO k=kts,ktf DO i = i_start, i_end mrdy=msfu(i,j-1)*rdy tendency(i,k,j-1) = tendency(i,k,j-1) - mrdy*(fqy(i,k,jp1)-fqy(i,k,jp0)) ENDDO ENDDO ENDIF jtmp = jp1 jp1 = jp0 jp0 = jtmp ENDDO j_loop_y_flux_6 ! next, x - flux divergence i_start = its i_end = ite j_start = jts j_end = MIN(jte,jde-1) ! higher order flux has a 5 or 7 point stencil, so compute ! bounds so we can switch to second order flux close to the boundary i_start_f = i_start i_end_f = i_end+1 IF(degrade_xs) then i_start = MAX(ids+1,its) i_start_f = ids+3 ENDIF IF(degrade_xe) then i_end = MIN(ide-1,ite) i_end_f = ide-2 ENDIF ! compute fluxes DO j = j_start, j_end ! 5th or 6th order flux DO k=kts,ktf DO i = i_start_f, i_end_f vel = 0.5*(ru(i,k,j)+ru(i-1,k,j)) fqx( i,k ) = vel*flux6( u(i-3,k,j), u(i-2,k,j), & u(i-1,k,j), u(i ,k,j), & u(i+1,k,j), u(i+2,k,j), & vel ) ENDDO ENDDO ! lower order fluxes close to boundaries (if not periodic or symmetric) ! specified uses upstream normal wind at boundaries IF( degrade_xs ) THEN IF( i_start == ids+1 ) THEN ! second order flux next to the boundary i = ids+1 DO k=kts,ktf ub = u(i-1,k,j) IF (specified .AND. u(i,k,j) .LT. 0.)ub = u(i,k,j) fqx(i, k) = 0.25*(ru(i,k,j)+ru(i-1,k,j)) & *(u(i,k,j)+ub) ENDDO END IF i = ids+2 DO k=kts,ktf vel = 0.5*(ru(i,k,j)+ru(i-1,k,j)) fqx( i, k ) = vel*flux4( u(i-2,k,j), u(i-1,k,j), & u(i ,k,j), u(i+1,k,j), & vel ) ENDDO ENDIF IF( degrade_xe ) THEN IF( i_end == ide-1 ) THEN ! second order flux next to the boundary i = ide DO k=kts,ktf ub = u(i,k,j) IF (specified .AND. u(i-1,k,j) .GT. 0.)ub = u(i-1,k,j) fqx(i, k) = 0.25*(ru(i,k,j)+ru(i-1,k,j)) & *(u(i-1,k,j)+ub) ENDDO ENDIF DO k=kts,ktf i = ide-1 vel = 0.5*(ru(i,k,j)+ru(i-1,k,j)) fqx( i,k ) = vel*flux4( u(i-2,k,j), u(i-1,k,j), & u(i ,k,j), u(i+1,k,j), & vel ) ENDDO ENDIF ! x flux-divergence into tendency DO k=kts,ktf DO i = i_start, i_end mrdx=msfu(i,j)*rdx tendency(i,k,j) = tendency(i,k,j) - mrdx*(fqx(i+1,k)-fqx(i,k)) ENDDO ENDDO ENDDO ELSE IF( horz_order == 5 ) THEN ! 5th order horizontal flux calculation ! This code is EXACTLY the same as the 6th order code ! EXCEPT the 5th order and 3rd operators are used in ! place of the 6th and 4th order operators ! determine boundary mods for flux operators ! We degrade the flux operators from 3rd/4th order ! to second order one gridpoint in from the boundaries for ! all boundary conditions except periodic and symmetry - these ! conditions have boundary zone data fill for correct application ! of the higher order flux stencils degrade_xs = .true. degrade_xe = .true. degrade_ys = .true. degrade_ye = .true. IF( config_flags%periodic_x .or. & config_flags%symmetric_xs .or. & (its > ids+2) ) degrade_xs = .false. IF( config_flags%periodic_x .or. & config_flags%symmetric_xe .or. & (ite < ide-2) ) degrade_xe = .false. IF( config_flags%periodic_y .or. & config_flags%symmetric_ys .or. & (jts > jds+2) ) degrade_ys = .false. IF( config_flags%periodic_y .or. & config_flags%symmetric_ye .or. & (jte < jde-3) ) degrade_ye = .false. !--------------- y - advection first i_start = its i_end = ite IF ( config_flags%open_xs .or. specified ) i_start = MAX(ids+1,its) IF ( config_flags%open_xe .or. specified ) i_end = MIN(ide-1,ite) IF ( config_flags%periodic_x ) i_start = its IF ( config_flags%periodic_x ) i_end = ite j_start = jts j_end = MIN(jte,jde-1) ! higher order flux has a 5 or 7 point stencil, so compute ! bounds so we can switch to second order flux close to the boundary j_start_f = j_start j_end_f = j_end+1 IF(degrade_ys) then j_start = MAX(jts,jds+1) j_start_f = jds+3 ENDIF IF(degrade_ye) then j_end = MIN(jte,jde-2) j_end_f = jde-3 ENDIF ! compute fluxes, 5th or 6th order jp1 = 2 jp0 = 1 j_loop_y_flux_5 : DO j = j_start, j_end+1 IF( (j >= j_start_f ) .and. (j <= j_end_f) ) THEN ! use full stencil DO k=kts,ktf DO i = i_start, i_end vel = 0.5*(rv(i,k,j)+rv(i-1,k,j)) fqy( i, k, jp1 ) = vel*flux5( & u(i,k,j-3), u(i,k,j-2), u(i,k,j-1), & u(i,k,j ), u(i,k,j+1), u(i,k,j+2), vel ) ENDDO ENDDO ! we must be close to some boundary where we need to reduce the order of the stencil ELSE IF ( j == jds+1 ) THEN ! 2nd order flux next to south boundary DO k=kts,ktf DO i = i_start, i_end fqy(i, k, jp1) = 0.25*(rv(i,k,j)+rv(i-1,k,j)) & *(u(i,k,j)+u(i,k,j-1)) ENDDO ENDDO ELSE IF ( j == jds+2 ) THEN ! third of 4th order flux 2 in from south boundary DO k=kts,ktf DO i = i_start, i_end vel = 0.5*(rv(i,k,j)+rv(i-1,k,j)) fqy( i, k, jp1 ) = vel*flux3( & u(i,k,j-2),u(i,k,j-1), u(i,k,j),u(i,k,j+1),vel ) ENDDO ENDDO ELSE IF ( j == jde-1 ) THEN ! 2nd order flux next to north boundary DO k=kts,ktf DO i = i_start, i_end fqy(i, k, jp1) = 0.25*(rv(i,k,j)+rv(i-1,k,j)) & *(u(i,k,j)+u(i,k,j-1)) ENDDO ENDDO ELSE IF ( j == jde-2 ) THEN ! 3rd or 4th order flux 2 in from north boundary DO k=kts,ktf DO i = i_start, i_end vel = 0.5*(rv(i,k,j)+rv(i-1,k,j)) fqy( i, k, jp1 ) = vel*flux3( & u(i,k,j-2),u(i,k,j-1), & u(i,k,j),u(i,k,j+1),vel ) ENDDO ENDDO END IF ! y flux-divergence into tendency IF(j > j_start) THEN DO k=kts,ktf DO i = i_start, i_end mrdy=msfu(i,j-1)*rdy tendency(i,k,j-1) = tendency(i,k,j-1) - mrdy*(fqy(i,k,jp1)-fqy(i,k,jp0)) ENDDO ENDDO ENDIF jtmp = jp1 jp1 = jp0 jp0 = jtmp ENDDO j_loop_y_flux_5 ! next, x - flux divergence i_start = its i_end = ite j_start = jts j_end = MIN(jte,jde-1) ! higher order flux has a 5 or 7 point stencil, so compute ! bounds so we can switch to second order flux close to the boundary i_start_f = i_start i_end_f = i_end+1 IF(degrade_xs) then i_start = MAX(ids+1,its) i_start_f = ids+3 ENDIF IF(degrade_xe) then i_end = MIN(ide-1,ite) i_end_f = ide-2 ENDIF ! compute fluxes DO j = j_start, j_end ! 5th or 6th order flux DO k=kts,ktf DO i = i_start_f, i_end_f vel = 0.5*(ru(i,k,j)+ru(i-1,k,j)) fqx( i,k ) = vel*flux5( u(i-3,k,j), u(i-2,k,j), & u(i-1,k,j), u(i ,k,j), & u(i+1,k,j), u(i+2,k,j), & vel ) ENDDO ENDDO ! lower order fluxes close to boundaries (if not periodic or symmetric) ! specified uses upstream normal wind at boundaries IF( degrade_xs ) THEN IF( i_start == ids+1 ) THEN ! second order flux next to the boundary i = ids+1 DO k=kts,ktf ub = u(i-1,k,j) IF (specified .AND. u(i,k,j) .LT. 0.)ub = u(i,k,j) fqx(i, k) = 0.25*(ru(i,k,j)+ru(i-1,k,j)) & *(u(i,k,j)+ub) ENDDO END IF i = ids+2 DO k=kts,ktf vel = 0.5*(ru(i,k,j)+ru(i-1,k,j)) fqx( i, k ) = vel*flux3( u(i-2,k,j), u(i-1,k,j), & u(i ,k,j), u(i+1,k,j), & vel ) ENDDO ENDIF IF( degrade_xe ) THEN IF( i_end == ide-1 ) THEN ! second order flux next to the boundary i = ide DO k=kts,ktf ub = u(i,k,j) IF (specified .AND. u(i-1,k,j) .GT. 0.)ub = u(i-1,k,j) fqx(i, k) = 0.25*(ru(i,k,j)+ru(i-1,k,j)) & *(u(i-1,k,j)+ub) ENDDO ENDIF DO k=kts,ktf i = ide-1 vel = 0.5*(ru(i,k,j)+ru(i-1,k,j)) fqx( i,k ) = vel*flux3( u(i-2,k,j), u(i-1,k,j), & u(i ,k,j), u(i+1,k,j), & vel ) ENDDO ENDIF ! x flux-divergence into tendency DO k=kts,ktf DO i = i_start, i_end mrdx=msfu(i,j)*rdx tendency(i,k,j) = tendency(i,k,j) - mrdx*(fqx(i+1,k)-fqx(i,k)) ENDDO ENDDO ENDDO ELSE IF( horz_order == 4 ) THEN ! determine boundary mods for flux operators ! We degrade the flux operators from 3rd/4th order ! to second order one gridpoint in from the boundaries for ! all boundary conditions except periodic and symmetry - these ! conditions have boundary zone data fill for correct application ! of the higher order flux stencils degrade_xs = .true. degrade_xe = .true. degrade_ys = .true. degrade_ye = .true. IF( config_flags%periodic_x .or. & config_flags%symmetric_xs .or. & (its > ids+1) ) degrade_xs = .false. IF( config_flags%periodic_x .or. & config_flags%symmetric_xe .or. & (ite < ide-1) ) degrade_xe = .false. IF( config_flags%periodic_y .or. & config_flags%symmetric_ys .or. & (jts > jds+1) ) degrade_ys = .false. IF( config_flags%periodic_y .or. & config_flags%symmetric_ye .or. & (jte < jde-2) ) degrade_ye = .false. !--------------- x - advection first i_start = its i_end = ite j_start = jts j_end = MIN(jte,jde-1) ! 3rd or 4th order flux has a 5 point stencil, so compute ! bounds so we can switch to second order flux close to the boundary i_start_f = i_start i_end_f = i_end+1 IF(degrade_xs) then i_start = ids+1 i_start_f = i_start+1 ENDIF IF(degrade_xe) then i_end = ide-1 i_end_f = ide-1 ENDIF ! compute fluxes DO j = j_start, j_end DO k=kts,ktf DO i = i_start_f, i_end_f vel = 0.5*(ru(i,k,j)+ru(i-1,k,j)) fqx( i, k ) = vel*flux4( u(i-2,k,j), u(i-1,k,j), & u(i ,k,j), u(i+1,k,j), vel ) ENDDO ENDDO ! second order flux close to boundaries (if not periodic or symmetric) ! specified uses upstream normal wind at boundaries IF( degrade_xs ) THEN i = i_start DO k=kts,ktf ub = u(i-1,k,j) IF (specified .AND. u(i,k,j) .LT. 0.)ub = u(i,k,j) fqx(i, k) = 0.25*(ru(i,k,j)+ru(i-1,k,j)) & *(u(i,k,j)+ub) ENDDO ENDIF IF( degrade_xe ) THEN i = i_end+1 DO k=kts,ktf ub = u(i,k,j) IF (specified .AND. u(i-1,k,j) .GT. 0.)ub = u(i-1,k,j) fqx(i, k) = 0.25*(ru(i,k,j)+ru(i-1,k,j)) & *(u(i-1,k,j)+ub) ENDDO ENDIF ! x flux-divergence into tendency DO k=kts,ktf DO i = i_start, i_end mrdx=msfu(i,j)*rdx tendency(i,k,j) = tendency(i,k,j) - mrdx*(fqx(i+1,k)-fqx(i,k)) ENDDO ENDDO ENDDO ! y flux divergence i_start = its i_end = ite IF ( config_flags%open_xs .or. specified ) i_start = MAX(ids+1,its) IF ( config_flags%open_xe .or. specified ) i_end = MIN(ide-1,ite) IF ( config_flags%periodic_x ) i_start = its IF ( config_flags%periodic_x ) i_end = ite j_start = jts j_end = MIN(jte,jde-1) ! 3rd or 4th order flux has a 5 point stencil, so compute ! bounds so we can switch to second order flux close to the boundary j_start_f = j_start j_end_f = j_end+1 !CJM these may not work with tiling because they define j_start and end in terms of domain dim IF(degrade_ys) then j_start = jds+1 j_start_f = j_start+1 ENDIF IF(degrade_ye) then j_end = jde-2 j_end_f = jde-2 ENDIF ! j flux loop for v flux of u momentum jp1 = 2 jp0 = 1 DO j = j_start, j_end+1 IF ( (j < j_start_f) .and. degrade_ys) THEN DO k = kts, ktf DO i = i_start, i_end fqy(i, k, jp1) = 0.25*(rv(i,k,j_start)+rv(i-1,k,j_start)) & *(u(i,k,j_start)+u(i,k,j_start-1)) ENDDO ENDDO ELSE IF ((j > j_end_f) .and. degrade_ye) THEN DO k = kts, ktf DO i = i_start, i_end fqy(i, k, jp1) = 0.25*(rv(i,k,j_end+1)+rv(i-1,k,j_end+1)) & *(u(i,k,j_end+1)+u(i,k,j_end)) ENDDO ENDDO ELSE ! 3rd or 4th order flux DO k = kts, ktf DO i = i_start, i_end vel = 0.5*(rv(i,k,j)+rv(i-1,k,j)) fqy( i, k, jp1 ) = vel*flux4( u(i,k,j-2), u(i,k,j-1), & u(i,k,j ), u(i,k,j+1), & vel ) ENDDO ENDDO END IF IF (j > j_start) THEN ! y flux-divergence into tendency DO k=kts,ktf DO i = i_start, i_end mrdy=msfu(i,j-1)*rdy tendency(i,k,j-1) = tendency(i,k,j-1) - mrdy*(fqy(i,k,jp1)-fqy(i,k,jp0)) ENDDO ENDDO END IF jtmp = jp1 jp1 = jp0 jp0 = jtmp ENDDO ELSE IF ( horz_order == 3 ) THEN ! As with the 5th and 6th order flux chioces, the 3rd and 4th order ! code is EXACTLY the same EXCEPT for the flux operator. ! determine boundary mods for flux operators ! We degrade the flux operators from 3rd/4th order ! to second order one gridpoint in from the boundaries for ! all boundary conditions except periodic and symmetry - these ! conditions have boundary zone data fill for correct application ! of the higher order flux stencils degrade_xs = .true. degrade_xe = .true. degrade_ys = .true. degrade_ye = .true. IF( config_flags%periodic_x .or. & config_flags%symmetric_xs .or. & (its > ids+1) ) degrade_xs = .false. IF( config_flags%periodic_x .or. & config_flags%symmetric_xe .or. & (ite < ide-1) ) degrade_xe = .false. IF( config_flags%periodic_y .or. & config_flags%symmetric_ys .or. & (jts > jds+1) ) degrade_ys = .false. IF( config_flags%periodic_y .or. & config_flags%symmetric_ye .or. & (jte < jde-2) ) degrade_ye = .false. !--------------- x - advection first i_start = its i_end = ite j_start = jts j_end = MIN(jte,jde-1) ! 3rd or 4th order flux has a 5 point stencil, so compute ! bounds so we can switch to second order flux close to the boundary i_start_f = i_start i_end_f = i_end+1 IF(degrade_xs) then i_start = ids+1 i_start_f = i_start+1 ENDIF IF(degrade_xe) then i_end = ide-1 i_end_f = ide-1 ENDIF ! compute fluxes DO j = j_start, j_end DO k=kts,ktf DO i = i_start_f, i_end_f vel = 0.5*(ru(i,k,j)+ru(i-1,k,j)) fqx( i, k ) = vel*flux3( u(i-2,k,j), u(i-1,k,j), & u(i ,k,j), u(i+1,k,j), vel ) ENDDO ENDDO ! second order flux close to boundaries (if not periodic or symmetric) ! specified uses upstream normal wind at boundaries IF( degrade_xs ) THEN i = i_start DO k=kts,ktf ub = u(i-1,k,j) IF (specified .AND. u(i,k,j) .LT. 0.)ub = u(i,k,j) fqx(i, k) = 0.25*(ru(i,k,j)+ru(i-1,k,j)) & *(u(i,k,j)+ub) ENDDO ENDIF IF( degrade_xe ) THEN i = i_end+1 DO k=kts,ktf ub = u(i,k,j) IF (specified .AND. u(i-1,k,j) .GT. 0.)ub = u(i-1,k,j) fqx(i, k) = 0.25*(ru(i,k,j)+ru(i-1,k,j)) & *(u(i-1,k,j)+ub) ENDDO ENDIF ! x flux-divergence into tendency DO k=kts,ktf DO i = i_start, i_end mrdx=msfu(i,j)*rdx tendency(i,k,j) = tendency(i,k,j) - mrdx*(fqx(i+1,k)-fqx(i,k)) ENDDO ENDDO ENDDO ! y flux divergence i_start = its i_end = ite IF ( config_flags%open_xs .or. specified ) i_start = MAX(ids+1,its) IF ( config_flags%open_xe .or. specified ) i_end = MIN(ide-1,ite) IF ( config_flags%periodic_x ) i_start = its IF ( config_flags%periodic_x ) i_end = ite j_start = jts j_end = MIN(jte,jde-1) ! 3rd or 4th order flux has a 5 point stencil, so compute ! bounds so we can switch to second order flux close to the boundary j_start_f = j_start j_end_f = j_end+1 !CJM these may not work with tiling because they define j_start and end in terms of domain dim IF(degrade_ys) then j_start = jds+1 j_start_f = j_start+1 ENDIF IF(degrade_ye) then j_end = jde-2 j_end_f = jde-2 ENDIF ! j flux loop for v flux of u momentum jp1 = 2 jp0 = 1 DO j = j_start, j_end+1 IF ( (j < j_start_f) .and. degrade_ys) THEN DO k = kts, ktf DO i = i_start, i_end fqy(i, k, jp1) = 0.25*(rv(i,k,j_start)+rv(i-1,k,j_start)) & *(u(i,k,j_start)+u(i,k,j_start-1)) ENDDO ENDDO ELSE IF ((j > j_end_f) .and. degrade_ye) THEN DO k = kts, ktf DO i = i_start, i_end fqy(i, k, jp1) = 0.25*(rv(i,k,j_end+1)+rv(i-1,k,j_end+1)) & *(u(i,k,j_end+1)+u(i,k,j_end)) ENDDO ENDDO ELSE ! 3rd or 4th order flux DO k = kts, ktf DO i = i_start, i_end vel = 0.5*(rv(i,k,j)+rv(i-1,k,j)) fqy( i, k, jp1 ) = vel*flux3( u(i,k,j-2), u(i,k,j-1), & u(i,k,j ), u(i,k,j+1), & vel ) ENDDO ENDDO END IF IF (j > j_start) THEN ! y flux-divergence into tendency DO k=kts,ktf DO i = i_start, i_end mrdy=msfu(i,j-1)*rdy tendency(i,k,j-1) = tendency(i,k,j-1) - mrdy*(fqy(i,k,jp1)-fqy(i,k,jp0)) ENDDO ENDDO END IF jtmp = jp1 jp1 = jp0 jp0 = jtmp ENDDO ELSE IF ( horz_order == 2 ) THEN i_start = its i_end = ite j_start = jts j_end = MIN(jte,jde-1) IF ( config_flags%open_xs ) i_start = MAX(ids+1,its) IF ( config_flags%open_xe ) i_end = MIN(ide-1,ite) IF ( specified ) i_start = MAX(ids+2,its) IF ( specified ) i_end = MIN(ide-2,ite) IF ( config_flags%periodic_x ) i_start = its IF ( config_flags%periodic_x ) i_end = ite DO j = j_start, j_end DO k=kts,ktf DO i = i_start, i_end mrdx=msfu(i,j)*rdx tendency(i,k,j)=tendency(i,k,j)-mrdx*0.25 & *((ru(i+1,k,j)+ru(i,k,j))*(u(i+1,k,j)+u(i,k,j)) & -(ru(i,k,j)+ru(i-1,k,j))*(u(i,k,j)+u(i-1,k,j))) ENDDO ENDDO ENDDO IF ( specified .AND. its .LE. ids+1 .AND. .NOT. config_flags%periodic_x ) THEN DO j = j_start, j_end DO k=kts,ktf i = ids+1 mrdx=msfu(i,j)*rdx ub = u(i-1,k,j) IF (u(i,k,j) .LT. 0.) ub = u(i,k,j) tendency(i,k,j)=tendency(i,k,j)-mrdx*0.25 & *((ru(i+1,k,j)+ru(i,k,j))*(u(i+1,k,j)+u(i,k,j)) & -(ru(i,k,j)+ru(i-1,k,j))*(u(i,k,j)+ub)) ENDDO ENDDO ENDIF IF ( specified .AND. ite .GE. ide-1 .AND. .NOT. config_flags%periodic_x ) THEN DO j = j_start, j_end DO k=kts,ktf i = ide-1 mrdx=msfu(i,j)*rdx ub = u(i+1,k,j) IF (u(i,k,j) .GT. 0.) ub = u(i,k,j) tendency(i,k,j)=tendency(i,k,j)-mrdx*0.25 & *((ru(i+1,k,j)+ru(i,k,j))*(ub+u(i,k,j)) & -(ru(i,k,j)+ru(i-1,k,j))*(u(i,k,j)+u(i-1,k,j))) ENDDO ENDDO ENDIF IF ( config_flags%open_ys .or. specified ) j_start = MAX(jds+1,jts) IF ( config_flags%open_ye .or. specified ) j_end = MIN(jde-2,jte) DO j = j_start, j_end DO k=kts,ktf DO i = i_start, i_end mrdy=msfu(i,j)*rdy tendency(i,k,j)=tendency(i,k,j)-mrdy*0.25 & *((rv(i,k,j+1)+rv(i-1,k,j+1))*(u(i,k,j+1)+u(i,k,j)) & -(rv(i,k,j)+rv(i-1,k,j))*(u(i,k,j)+u(i,k,j-1))) ENDDO ENDDO ENDDO ELSE WRITE ( wrf_err_message , * ) 'module_advect: advect_u_6a: h_order not known ',horz_order CALL wrf_error_fatal ( TRIM( wrf_err_message ) ) ENDIF horizontal_order_test ! radiative lateral boundary condition in x for normal velocity (u) IF ( (config_flags%open_xs) .and. its == ids ) THEN j_start = jts j_end = MIN(jte,jde-1) DO j = j_start, j_end DO k = kts, ktf ub = MIN(ru(its,k,j)-cb*mut(its,j), 0.) tendency(its,k,j) = tendency(its,k,j) & - rdx*ub*(u_old(its+1,k,j) - u_old(its,k,j)) ENDDO ENDDO ENDIF IF ( (config_flags%open_xe) .and. ite == ide ) THEN j_start = jts j_end = MIN(jte,jde-1) DO j = j_start, j_end DO k = kts, ktf ub = MAX(ru(ite,k,j)+cb*mut(ite-1,j), 0.) tendency(ite,k,j) = tendency(ite,k,j) & - rdx*ub*(u_old(ite,k,j) - u_old(ite-1,k,j)) ENDDO ENDDO ENDIF ! pick up the rest of the horizontal radiation boundary conditions. ! (these are the computations that don't require 'cb') ! first, set to index ranges i_start = its i_end = MIN(ite,ide) imin = ids imax = ide-1 IF (config_flags%open_xs) THEN i_start = MAX(ids+1, its) imin = ids ENDIF IF (config_flags%open_xe) THEN i_end = MIN(ite,ide-1) imax = ide-1 ENDIF IF( (config_flags%open_ys) .and. (jts == jds)) THEN DO i = i_start, i_end mrdy=msfu(i,jts)*rdy ip = MIN( imax, i ) im = MAX( imin, i-1 ) DO k=kts,ktf vw = 0.5*(rv(ip,k,jts)+rv(im,k,jts)) vb = MIN( vw, 0. ) dvm = rv(ip,k,jts+1)-rv(ip,k,jts) dvp = rv(im,k,jts+1)-rv(im,k,jts) tendency(i,k,jts)=tendency(i,k,jts)-mrdy*( & vb*(u_old(i,k,jts+1)-u_old(i,k,jts)) & +0.5*u(i,k,jts)*(dvm+dvp)) ENDDO ENDDO ENDIF IF( (config_flags%open_ye) .and. (jte == jde)) THEN DO i = i_start, i_end mrdy=msfu(i,jte-1)*rdy ip = MIN( imax, i ) im = MAX( imin, i-1 ) DO k=kts,ktf vw = 0.5*(rv(ip,k,jte)+rv(im,k,jte)) vb = MAX( vw, 0. ) dvm = rv(ip,k,jte)-rv(ip,k,jte-1) dvp = rv(im,k,jte)-rv(im,k,jte-1) tendency(i,k,jte-1)=tendency(i,k,jte-1)-mrdy*( & vb*(u_old(i,k,jte-1)-u_old(i,k,jte-2)) & +0.5*u(i,k,jte-1)*(dvm+dvp)) ENDDO ENDDO ENDIF !-------------------- vertical advection i_start = its i_end = ite j_start = jts j_end = min(jte,jde-1) ! IF ( config_flags%open_xs ) i_start = MAX(ids+1,its) ! IF ( config_flags%open_xe ) i_end = MIN(ide-1,ite) IF ( config_flags%open_ys .or. specified ) i_start = MAX(ids+1,its) IF ( config_flags%open_ye .or. specified ) i_end = MIN(ide-1,ite) IF ( config_flags%periodic_x ) i_start = its IF ( config_flags%periodic_x ) i_end = ite DO i = i_start, i_end vflux(i,kts)=0. vflux(i,kte)=0. ENDDO vert_order_test : IF (vert_order == 6) THEN DO j = j_start, j_end DO k=kts+3,ktf-2 DO i = i_start, i_end vel=0.5*(rom(i-1,k,j)+rom(i,k,j)) vflux(i,k) = vel*flux6( & u(i,k-3,j), u(i,k-2,j), u(i,k-1,j), & u(i,k ,j), u(i,k+1,j), u(i,k+2,j), -vel ) ENDDO ENDDO DO i = i_start, i_end k=kts+1 vflux(i,k)=0.5*(rom(i,k,j)+rom(i-1,k,j)) & *(fzm(k)*u(i,k,j)+fzp(k)*u(i,k-1,j)) k = kts+2 vel=0.5*(rom(i,k,j)+rom(i-1,k,j)) vflux(i,k) = vel*flux4( & u(i,k-2,j), u(i,k-1,j), & u(i,k ,j), u(i,k+1,j), -vel ) k = ktf-1 vel=0.5*(rom(i,k,j)+rom(i-1,k,j)) vflux(i,k) = vel*flux4( & u(i,k-2,j), u(i,k-1,j), & u(i,k ,j), u(i,k+1,j), -vel ) k=ktf vflux(i,k)=0.5*(rom(i,k,j)+rom(i-1,k,j)) & *(fzm(k)*u(i,k,j)+fzp(k)*u(i,k-1,j)) ENDDO DO k=kts,ktf DO i = i_start, i_end tendency(i,k,j)=tendency(i,k,j)-rdzw(k)*(vflux(i,k+1)-vflux(i,k)) ENDDO ENDDO ENDDO ELSE IF (vert_order == 5) THEN DO j = j_start, j_end DO k=kts+3,ktf-2 DO i = i_start, i_end vel=0.5*(rom(i-1,k,j)+rom(i,k,j)) vflux(i,k) = vel*flux5( & u(i,k-3,j), u(i,k-2,j), u(i,k-1,j), & u(i,k ,j), u(i,k+1,j), u(i,k+2,j), -vel ) ENDDO ENDDO DO i = i_start, i_end k=kts+1 vflux(i,k)=0.5*(rom(i,k,j)+rom(i-1,k,j)) & *(fzm(k)*u(i,k,j)+fzp(k)*u(i,k-1,j)) k = kts+2 vel=0.5*(rom(i,k,j)+rom(i-1,k,j)) vflux(i,k) = vel*flux3( & u(i,k-2,j), u(i,k-1,j), & u(i,k ,j), u(i,k+1,j), -vel ) k = ktf-1 vel=0.5*(rom(i,k,j)+rom(i-1,k,j)) vflux(i,k) = vel*flux3( & u(i,k-2,j), u(i,k-1,j), & u(i,k ,j), u(i,k+1,j), -vel ) k=ktf vflux(i,k)=0.5*(rom(i,k,j)+rom(i-1,k,j)) & *(fzm(k)*u(i,k,j)+fzp(k)*u(i,k-1,j)) ENDDO DO k=kts,ktf DO i = i_start, i_end tendency(i,k,j)=tendency(i,k,j)-rdzw(k)*(vflux(i,k+1)-vflux(i,k)) ENDDO ENDDO ENDDO ELSE IF (vert_order == 4) THEN DO j = j_start, j_end DO k=kts+2,ktf-1 DO i = i_start, i_end vel=0.5*(rom(i-1,k,j)+rom(i,k,j)) vflux(i,k) = vel*flux4( & u(i,k-2,j), u(i,k-1,j), & u(i,k ,j), u(i,k+1,j), -vel ) ENDDO ENDDO DO i = i_start, i_end k=kts+1 vflux(i,k)=0.5*(rom(i,k,j)+rom(i-1,k,j)) & *(fzm(k)*u(i,k,j)+fzp(k)*u(i,k-1,j)) k=ktf vflux(i,k)=0.5*(rom(i,k,j)+rom(i-1,k,j)) & *(fzm(k)*u(i,k,j)+fzp(k)*u(i,k-1,j)) ENDDO DO k=kts,ktf DO i = i_start, i_end tendency(i,k,j)=tendency(i,k,j)-rdzw(k)*(vflux(i,k+1)-vflux(i,k)) ENDDO ENDDO ENDDO ELSE IF (vert_order == 3) THEN DO j = j_start, j_end DO k=kts+2,ktf-1 DO i = i_start, i_end vel=0.5*(rom(i-1,k,j)+rom(i,k,j)) vflux(i,k) = vel*flux3( & u(i,k-2,j), u(i,k-1,j), & u(i,k ,j), u(i,k+1,j), -vel ) ENDDO ENDDO DO i = i_start, i_end k=kts+1 vflux(i,k)=0.5*(rom(i,k,j)+rom(i-1,k,j)) & *(fzm(k)*u(i,k,j)+fzp(k)*u(i,k-1,j)) k=ktf vflux(i,k)=0.5*(rom(i,k,j)+rom(i-1,k,j)) & *(fzm(k)*u(i,k,j)+fzp(k)*u(i,k-1,j)) ENDDO DO k=kts,ktf DO i = i_start, i_end tendency(i,k,j)=tendency(i,k,j)-rdzw(k)*(vflux(i,k+1)-vflux(i,k)) ENDDO ENDDO ENDDO ELSE IF (vert_order == 2) THEN DO j = j_start, j_end DO k=kts+1,ktf DO i = i_start, i_end vflux(i,k)=0.5*(rom(i,k,j)+rom(i-1,k,j)) & *(fzm(k)*u(i,k,j)+fzp(k)*u(i,k-1,j)) ENDDO ENDDO DO k=kts,ktf DO i = i_start, i_end tendency(i,k,j)=tendency(i,k,j)-rdzw(k)*(vflux(i,k+1)-vflux(i,k)) ENDDO ENDDO ENDDO ELSE WRITE ( wrf_err_message , * ) 'module_advect: advect_u_6a: v_order not known ',vert_order CALL wrf_error_fatal ( TRIM( wrf_err_message ) ) ENDIF vert_order_test END SUBROUTINE advect_u !------------------------------------------------------------------------------- SUBROUTINE advect_v ( v, v_old, tendency, & ru, rv, rom, & mut, config_flags, & msfu, msfv, msft, & fzm, fzp, & rdx, rdy, rdzw, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) IMPLICIT NONE ! Input data TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: v, & v_old, & ru, & rv, & rom REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mut REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: msfu, & msfv, & msft REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fzm, & fzp, & rdzw REAL , INTENT(IN ) :: rdx, & rdy ! Local data INTEGER :: i, j, k, itf, jtf, ktf INTEGER :: i_start, i_end, j_start, j_end INTEGER :: i_start_f, i_end_f, j_start_f, j_end_f INTEGER :: jmin, jmax, jp, jm, imin, imax REAL :: mrdx, mrdy, ub, vb, uw, vw, dup, dum REAL , DIMENSION(its:ite, kts:kte) :: vflux REAL, DIMENSION( its:ite+1, kts:kte ) :: fqx REAL, DIMENSION( its:ite, kts:kte, 2 ) :: fqy INTEGER :: horz_order INTEGER :: vert_order LOGICAL :: degrade_xs, degrade_ys LOGICAL :: degrade_xe, degrade_ye INTEGER :: jp1, jp0, jtmp ! definition of flux operators, 3rd, 4th, 5th or 6th order REAL :: flux3, flux4, flux5, flux6 REAL :: q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua, vel flux4(q_im2, q_im1, q_i, q_ip1, ua) = & ( 7.*(q_i + q_im1) - (q_ip1 + q_im2) )/12.0 flux3(q_im2, q_im1, q_i, q_ip1, ua) = & flux4(q_im2, q_im1, q_i, q_ip1, ua) + & sign(1.,ua)*((q_ip1 - q_im2)-3.*(q_i-q_im1))/12.0 flux6(q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua) = & ( 37.*(q_i+q_im1) - 8.*(q_ip1+q_im2) & +(q_ip2+q_im3) )/60.0 flux5(q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua) = & flux6(q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua) & -sign(1.,ua)*( & (q_ip2-q_im3)-5.*(q_ip1-q_im2)+10.*(q_i-q_im1) )/60.0 LOGICAL :: specified specified = .false. if(config_flags%specified .or. config_flags%nested) specified = .true. ! set order for the advection schemes ktf=MIN(kte,kde-1) horz_order = config_flags%h_mom_adv_order vert_order = config_flags%v_mom_adv_order ! here is the choice of flux operators horizontal_order_test : IF( horz_order == 6 ) THEN ! determine boundary mods for flux operators ! We degrade the flux operators from 3rd/4th order ! to second order one gridpoint in from the boundaries for ! all boundary conditions except periodic and symmetry - these ! conditions have boundary zone data fill for correct application ! of the higher order flux stencils degrade_xs = .true. degrade_xe = .true. degrade_ys = .true. degrade_ye = .true. IF( config_flags%periodic_x .or. & config_flags%symmetric_xs .or. & (its > ids+2) ) degrade_xs = .false. IF( config_flags%periodic_x .or. & config_flags%symmetric_xe .or. & (ite < ide-3) ) degrade_xe = .false. IF( config_flags%periodic_y .or. & config_flags%symmetric_ys .or. & (jts > jds+2) ) degrade_ys = .false. IF( config_flags%periodic_y .or. & config_flags%symmetric_ye .or. & (jte < jde-2) ) degrade_ye = .false. !--------------- y - advection first ktf=MIN(kte,kde-1) i_start = its i_end = MIN(ite,ide-1) j_start = jts j_end = jte ! higher order flux has a 5 or 7 point stencil, so compute ! bounds so we can switch to second order flux close to the boundary j_start_f = j_start j_end_f = j_end+1 IF(degrade_ys) then j_start = MAX(jts,jds+1) j_start_f = jds+3 ENDIF IF(degrade_ye) then j_end = MIN(jte,jde-1) j_end_f = jde-2 ENDIF ! compute fluxes, 5th or 6th order jp1 = 2 jp0 = 1 j_loop_y_flux_6 : DO j = j_start, j_end+1 IF( (j >= j_start_f ) .and. (j <= j_end_f) ) THEN DO k=kts,ktf DO i = i_start, i_end vel = 0.5*(rv(i,k,j)+rv(i,k,j-1)) fqy( i, k, jp1 ) = vel*flux6( & v(i,k,j-3), v(i,k,j-2), v(i,k,j-1), & v(i,k,j ), v(i,k,j+1), v(i,k,j+2), vel ) ENDDO ENDDO ! we must be close to some boundary where we need to reduce the order of the stencil ! specified uses upstream normal wind at boundaries ELSE IF ( j == jds+1 ) THEN ! 2nd order flux next to south boundary DO k=kts,ktf DO i = i_start, i_end vb = v(i,k,j-1) IF (specified .AND. v(i,k,j) .LT. 0.)vb = v(i,k,j) fqy(i, k, jp1) = 0.25*(rv(i,k,j)+rv(i,k,j-1)) & *(v(i,k,j)+vb) ENDDO ENDDO ELSE IF ( j == jds+2 ) THEN ! third of 4th order flux 2 in from south boundary DO k=kts,ktf DO i = i_start, i_end vel = 0.5*(rv(i,k,j)+rv(i,k,j-1)) fqy( i, k, jp1 ) = vel*flux4( & v(i,k,j-2),v(i,k,j-1),v(i,k,j),v(i,k,j+1),vel ) ENDDO ENDDO ELSE IF ( j == jde ) THEN ! 2nd order flux next to north boundary DO k=kts,ktf DO i = i_start, i_end vb = v(i,k,j) IF (specified .AND. v(i,k,j-1) .GT. 0.)vb = v(i,k,j-1) fqy(i, k, jp1) = 0.25*(rv(i,k,j)+rv(i,k,j-1)) & *(vb+v(i,k,j-1)) ENDDO ENDDO ELSE IF ( j == jde-1 ) THEN ! 3rd or 4th order flux 2 in from north boundary DO k=kts,ktf DO i = i_start, i_end vel = 0.5*(rv(i,k,j)+rv(i,k,j-1)) fqy( i, k, jp1 ) = vel*flux4( & v(i,k,j-2),v(i,k,j-1),v(i,k,j),v(i,k,j+1),vel ) ENDDO ENDDO END IF ! y flux-divergence into tendency IF(j > j_start) THEN DO k=kts,ktf DO i = i_start, i_end mrdy=msfv(i,j-1)*rdy tendency(i,k,j-1) = tendency(i,k,j-1) - mrdy*(fqy(i,k,jp1)-fqy(i,k,jp0)) ENDDO ENDDO ENDIF jtmp = jp1 jp1 = jp0 jp0 = jtmp ENDDO j_loop_y_flux_6 ! next, x - flux divergence i_start = its i_end = MIN(ite,ide-1) j_start = jts j_end = jte IF ( config_flags%open_ys .or. specified ) j_start = MAX(jds+1,jts) IF ( config_flags%open_ye .or. specified ) j_end = MIN(jde-1,jte) ! higher order flux has a 5 or 7 point stencil, so compute ! bounds so we can switch to second order flux close to the boundary i_start_f = i_start i_end_f = i_end+1 IF(degrade_xs) then i_start = MAX(ids+1,its) i_start_f = i_start+2 ENDIF IF(degrade_xe) then i_end = MIN(ide-2,ite) i_end_f = ide-3 ENDIF ! compute fluxes DO j = j_start, j_end ! 5th or 6th order flux DO k=kts,ktf DO i = i_start_f, i_end_f vel = 0.5*(ru(i,k,j)+ru(i,k,j-1)) fqx( i, k ) = vel*flux6( v(i-3,k,j), v(i-2,k,j), & v(i-1,k,j), v(i ,k,j), & v(i+1,k,j), v(i+2,k,j), & vel ) ENDDO ENDDO ! lower order fluxes close to boundaries (if not periodic or symmetric) IF( degrade_xs ) THEN IF( i_start == ids+1 ) THEN ! second order flux next to the boundary i = ids+1 DO k=kts,ktf fqx(i,k) = 0.25*(ru(i,k,j)+ru(i,k,j-1)) & *(v(i,k,j)+v(i-1,k,j)) ENDDO ENDIF i = ids+2 DO k=kts,ktf vel = 0.5*(ru(i,k,j)+ru(i,k,j-1)) fqx( i,k ) = vel*flux4( v(i-2,k,j), v(i-1,k,j), & v(i ,k,j), v(i+1,k,j), & vel ) ENDDO ENDIF IF( degrade_xe ) THEN IF( i_end == ide-2 ) THEN ! second order flux next to the boundary i = ide-1 DO k=kts,ktf fqx(i,k) = 0.25*(ru(i_end+1,k,j)+ru(i_end+1,k,j-1)) & *(v(i_end+1,k,j)+v(i_end,k,j)) ENDDO ENDIF i = ide-2 DO k=kts,ktf vel = 0.5*(ru(i,k,j)+ru(i,k,j-1)) fqx( i,k ) = vel*flux4( v(i-2,k,j), v(i-1,k,j), & v(i ,k,j), v(i+1,k,j), & vel ) ENDDO ENDIF ! x flux-divergence into tendency DO k=kts,ktf DO i = i_start, i_end mrdx=msfv(i,j)*rdx tendency(i,k,j) = tendency(i,k,j) - mrdx*(fqx(i+1,k)-fqx(i,k)) ENDDO ENDDO ENDDO ELSE IF( horz_order == 5 ) THEN ! 5th order horizontal flux calculation ! This code is EXACTLY the same as the 6th order code ! EXCEPT the 5th order and 3rd operators are used in ! place of the 6th and 4th order operators ! determine boundary mods for flux operators ! We degrade the flux operators from 3rd/4th order ! to second order one gridpoint in from the boundaries for ! all boundary conditions except periodic and symmetry - these ! conditions have boundary zone data fill for correct application ! of the higher order flux stencils degrade_xs = .true. degrade_xe = .true. degrade_ys = .true. degrade_ye = .true. IF( config_flags%periodic_x .or. & config_flags%symmetric_xs .or. & (its > ids+2) ) degrade_xs = .false. IF( config_flags%periodic_x .or. & config_flags%symmetric_xe .or. & (ite < ide-3) ) degrade_xe = .false. IF( config_flags%periodic_y .or. & config_flags%symmetric_ys .or. & (jts > jds+2) ) degrade_ys = .false. IF( config_flags%periodic_y .or. & config_flags%symmetric_ye .or. & (jte < jde-2) ) degrade_ye = .false. !--------------- y - advection first i_start = its i_end = MIN(ite,ide-1) j_start = jts j_end = jte ! higher order flux has a 5 or 7 point stencil, so compute ! bounds so we can switch to second order flux close to the boundary j_start_f = j_start j_end_f = j_end+1 IF(degrade_ys) then j_start = MAX(jts,jds+1) j_start_f = jds+3 ENDIF IF(degrade_ye) then j_end = MIN(jte,jde-1) j_end_f = jde-2 ENDIF ! compute fluxes, 5th or 6th order jp1 = 2 jp0 = 1 j_loop_y_flux_5 : DO j = j_start, j_end+1 IF( (j >= j_start_f ) .and. (j <= j_end_f) ) THEN DO k=kts,ktf DO i = i_start, i_end vel = 0.5*(rv(i,k,j)+rv(i,k,j-1)) fqy( i, k, jp1 ) = vel*flux5( & v(i,k,j-3), v(i,k,j-2), v(i,k,j-1), & v(i,k,j ), v(i,k,j+1), v(i,k,j+2), vel ) ENDDO ENDDO ! we must be close to some boundary where we need to reduce the order of the stencil ! specified uses upstream normal wind at boundaries ELSE IF ( j == jds+1 ) THEN ! 2nd order flux next to south boundary DO k=kts,ktf DO i = i_start, i_end vb = v(i,k,j-1) IF (specified .AND. v(i,k,j) .LT. 0.)vb = v(i,k,j) fqy(i, k, jp1) = 0.25*(rv(i,k,j)+rv(i,k,j-1)) & *(v(i,k,j)+vb) ENDDO ENDDO ELSE IF ( j == jds+2 ) THEN ! third of 4th order flux 2 in from south boundary DO k=kts,ktf DO i = i_start, i_end vel = 0.5*(rv(i,k,j)+rv(i,k,j-1)) fqy( i, k, jp1 ) = vel*flux3( & v(i,k,j-2),v(i,k,j-1),v(i,k,j),v(i,k,j+1),vel ) ENDDO ENDDO ELSE IF ( j == jde ) THEN ! 2nd order flux next to north boundary DO k=kts,ktf DO i = i_start, i_end vb = v(i,k,j) IF (specified .AND. v(i,k,j-1) .GT. 0.)vb = v(i,k,j-1) fqy(i, k, jp1) = 0.25*(rv(i,k,j)+rv(i,k,j-1)) & *(vb+v(i,k,j-1)) ENDDO ENDDO ELSE IF ( j == jde-1 ) THEN ! 3rd or 4th order flux 2 in from north boundary DO k=kts,ktf DO i = i_start, i_end vel = 0.5*(rv(i,k,j)+rv(i,k,j-1)) fqy( i, k, jp1 ) = vel*flux3( & v(i,k,j-2),v(i,k,j-1),v(i,k,j),v(i,k,j+1),vel ) ENDDO ENDDO END IF ! y flux-divergence into tendency IF(j > j_start) THEN DO k=kts,ktf DO i = i_start, i_end mrdy=msfv(i,j-1)*rdy tendency(i,k,j-1) = tendency(i,k,j-1) - mrdy*(fqy(i,k,jp1)-fqy(i,k,jp0)) ENDDO ENDDO ENDIF jtmp = jp1 jp1 = jp0 jp0 = jtmp ENDDO j_loop_y_flux_5 ! next, x - flux divergence i_start = its i_end = MIN(ite,ide-1) j_start = jts j_end = jte IF ( config_flags%open_ys .or. specified ) j_start = MAX(jds+1,jts) IF ( config_flags%open_ye .or. specified ) j_end = MIN(jde-1,jte) ! higher order flux has a 5 or 7 point stencil, so compute ! bounds so we can switch to second order flux close to the boundary i_start_f = i_start i_end_f = i_end+1 IF(degrade_xs) then i_start = MAX(ids+1,its) i_start_f = i_start+2 ENDIF IF(degrade_xe) then i_end = MIN(ide-2,ite) i_end_f = ide-3 ENDIF ! compute fluxes DO j = j_start, j_end ! 5th or 6th order flux DO k=kts,ktf DO i = i_start_f, i_end_f vel = 0.5*(ru(i,k,j)+ru(i,k,j-1)) fqx( i, k ) = vel*flux5( v(i-3,k,j), v(i-2,k,j), & v(i-1,k,j), v(i ,k,j), & v(i+1,k,j), v(i+2,k,j), & vel ) ENDDO ENDDO ! lower order fluxes close to boundaries (if not periodic or symmetric) IF( degrade_xs ) THEN IF( i_start == ids+1 ) THEN ! second order flux next to the boundary i = ids+1 DO k=kts,ktf fqx(i,k) = 0.25*(ru(i,k,j)+ru(i,k,j-1)) & *(v(i,k,j)+v(i-1,k,j)) ENDDO ENDIF i = ids+2 DO k=kts,ktf vel = 0.5*(ru(i,k,j)+ru(i,k,j-1)) fqx( i,k ) = vel*flux3( v(i-2,k,j), v(i-1,k,j), & v(i ,k,j), v(i+1,k,j), & vel ) ENDDO ENDIF IF( degrade_xe ) THEN IF( i_end == ide-2 ) THEN ! second order flux next to the boundary i = ide-1 DO k=kts,ktf fqx(i,k) = 0.25*(ru(i_end+1,k,j)+ru(i_end+1,k,j-1)) & *(v(i_end+1,k,j)+v(i_end,k,j)) ENDDO ENDIF i = ide-2 DO k=kts,ktf vel = 0.5*(ru(i,k,j)+ru(i,k,j-1)) fqx( i,k ) = vel*flux3( v(i-2,k,j), v(i-1,k,j), & v(i ,k,j), v(i+1,k,j), & vel ) ENDDO ENDIF ! x flux-divergence into tendency DO k=kts,ktf DO i = i_start, i_end mrdx=msfv(i,j)*rdx tendency(i,k,j) = tendency(i,k,j) - mrdx*(fqx(i+1,k)-fqx(i,k)) ENDDO ENDDO ENDDO ELSE IF( horz_order == 4 ) THEN ! determine boundary mods for flux operators ! We degrade the flux operators from 3rd/4th order ! to second order one gridpoint in from the boundaries for ! all boundary conditions except periodic and symmetry - these ! conditions have boundary zone data fill for correct application ! of the higher order flux stencils degrade_xs = .true. degrade_xe = .true. degrade_ys = .true. degrade_ye = .true. IF( config_flags%periodic_x .or. & config_flags%symmetric_xs .or. & (its > ids+1) ) degrade_xs = .false. IF( config_flags%periodic_x .or. & config_flags%symmetric_xe .or. & (ite < ide-2) ) degrade_xe = .false. IF( config_flags%periodic_y .or. & config_flags%symmetric_ys .or. & (jts > jds+1) ) degrade_ys = .false. IF( config_flags%periodic_y .or. & config_flags%symmetric_ye .or. & (jte < jde-1) ) degrade_ye = .false. !--------------- y - advection first ktf=MIN(kte,kde-1) i_start = its i_end = MIN(ite,ide-1) j_start = jts j_end = jte ! 3rd or 4th order flux has a 5 point stencil, so compute ! bounds so we can switch to second order flux close to the boundary j_start_f = j_start j_end_f = j_end+1 !CJM May not work with tiling because defined in terms of domain dims IF(degrade_ys) then j_start = jds+1 j_start_f = j_start+1 ENDIF IF(degrade_ye) then j_end = jde-1 j_end_f = jde-1 ENDIF ! compute fluxes ! specified uses upstream normal wind at boundaries jp0 = 1 jp1 = 2 DO j = j_start, j_end+1 IF ((j == j_start) .and. degrade_ys) THEN DO k = kts,ktf DO i = i_start, i_end vb = v(i,k,j-1) IF (specified .AND. v(i,k,j) .LT. 0.)vb = v(i,k,j) fqy(i, k, jp1) = 0.25*(rv(i,k,j)+rv(i,k,j-1)) & *(v(i,k,j)+vb) ENDDO ENDDO ELSE IF ((j == j_end+1) .and. degrade_ye) THEN DO k = kts, ktf DO i = i_start, i_end vb = v(i,k,j) IF (specified .AND. v(i,k,j-1) .GT. 0.)vb = v(i,k,j-1) fqy(i, k, jp1) = 0.25*(rv(i,k,j)+rv(i,k,j-1)) & *(vb+v(i,k,j-1)) ENDDO ENDDO ELSE DO k = kts, ktf DO i = i_start, i_end vel = 0.5*(rv(i,k,j)+rv(i,k,j-1)) fqy( i,k,jp1 ) = vel*flux4( v(i,k,j-2), v(i,k,j-1), & v(i,k,j ), v(i,k,j+1), & vel ) ENDDO ENDDO END IF IF( j > j_start) THEN DO k = kts, ktf DO i = i_start, i_end mrdy=msfv(i,j-1)*rdy tendency(i,k,j-1) = tendency(i,k,j-1) - mrdy*(fqy(i,k,jp1)-fqy(i,k,jp0)) ENDDO ENDDO END IF jtmp = jp1 jp1 = jp0 jp0 = jtmp ENDDO ! next, x - flux divergence i_start = its i_end = MIN(ite,ide-1) j_start = jts j_end = jte IF ( config_flags%open_ys .or. specified ) j_start = MAX(jds+1,jts) IF ( config_flags%open_ye .or. specified ) j_end = MIN(jde-1,jte) ! 3rd or 4th order flux has a 5 point stencil, so compute ! bounds so we can switch to second order flux close to the boundary i_start_f = i_start i_end_f = i_end+1 IF(degrade_xs) then i_start = ids+1 i_start_f = i_start+1 ENDIF IF(degrade_xe) then i_end = ide-2 i_end_f = ide-2 ENDIF ! compute fluxes DO j = j_start, j_end ! 3rd or 4th order flux DO k=kts,ktf DO i = i_start_f, i_end_f vel = 0.5*(ru(i,k,j)+ru(i,k,j-1)) fqx( i, k ) = vel*flux4( v(i-2,k,j), v(i-1,k,j), & v(i ,k,j), v(i+1,k,j), & vel ) ENDDO ENDDO ! second order flux close to boundaries (if not periodic or symmetric) IF( degrade_xs ) THEN DO k=kts,ktf fqx(i_start,k) = 0.25*(ru(i_start,k,j)+ru(i_start,k,j-1)) & *(v(i_start,k,j)+v(i_start-1,k,j)) ENDDO ENDIF IF( degrade_xe ) THEN DO k=kts,ktf fqx(i_end+1,k) = 0.25*(ru(i_end+1,k,j)+ru(i_end+1,k,j-1)) & *(v(i_end+1,k,j)+v(i_end,k,j)) ENDDO ENDIF ! x flux-divergence into tendency DO k=kts,ktf DO i = i_start, i_end mrdx=msfv(i,j)*rdx tendency(i,k,j) = tendency(i,k,j) - mrdx*(fqx(i+1,k)-fqx(i,k)) ENDDO ENDDO ENDDO ELSE IF( horz_order == 3 ) THEN ! determine boundary mods for flux operators ! We degrade the flux operators from 3rd/4th order ! to second order one gridpoint in from the boundaries for ! all boundary conditions except periodic and symmetry - these ! conditions have boundary zone data fill for correct application ! of the higher order flux stencils degrade_xs = .true. degrade_xe = .true. degrade_ys = .true. degrade_ye = .true. IF( config_flags%periodic_x .or. & config_flags%symmetric_xs .or. & (its > ids+1) ) degrade_xs = .false. IF( config_flags%periodic_x .or. & config_flags%symmetric_xe .or. & (ite < ide-2) ) degrade_xe = .false. IF( config_flags%periodic_y .or. & config_flags%symmetric_ys .or. & (jts > jds+1) ) degrade_ys = .false. IF( config_flags%periodic_y .or. & config_flags%symmetric_ye .or. & (jte < jde-1) ) degrade_ye = .false. !--------------- y - advection first ktf=MIN(kte,kde-1) i_start = its i_end = MIN(ite,ide-1) j_start = jts j_end = jte ! 3rd or 4th order flux has a 5 point stencil, so compute ! bounds so we can switch to second order flux close to the boundary j_start_f = j_start j_end_f = j_end+1 !CJM May not work with tiling because defined in terms of domain dims IF(degrade_ys) then j_start = jds+1 j_start_f = j_start+1 ENDIF IF(degrade_ye) then j_end = jde-1 j_end_f = jde-1 ENDIF ! compute fluxes ! specified uses upstream normal wind at boundaries jp0 = 1 jp1 = 2 DO j = j_start, j_end+1 IF ((j == j_start) .and. degrade_ys) THEN DO k = kts,ktf DO i = i_start, i_end vb = v(i,k,j-1) IF (specified .AND. v(i,k,j) .LT. 0.)vb = v(i,k,j) fqy(i, k, jp1) = 0.25*(rv(i,k,j)+rv(i,k,j-1)) & *(v(i,k,j)+vb) ENDDO ENDDO ELSE IF ((j == j_end+1) .and. degrade_ye) THEN DO k = kts, ktf DO i = i_start, i_end vb = v(i,k,j) IF (specified .AND. v(i,k,j-1) .GT. 0.)vb = v(i,k,j-1) fqy(i, k, jp1) = 0.25*(rv(i,k,j)+rv(i,k,j-1)) & *(vb+v(i,k,j-1)) ENDDO ENDDO ELSE DO k = kts, ktf DO i = i_start, i_end vel = 0.5*(rv(i,k,j)+rv(i,k,j-1)) fqy( i,k,jp1 ) = vel*flux3( v(i,k,j-2), v(i,k,j-1), & v(i,k,j ), v(i,k,j+1), & vel ) ENDDO ENDDO END IF IF( j > j_start) THEN DO k = kts, ktf DO i = i_start, i_end mrdy=msfv(i,j-1)*rdy tendency(i,k,j-1) = tendency(i,k,j-1) - mrdy*(fqy(i,k,jp1)-fqy(i,k,jp0)) ENDDO ENDDO END IF jtmp = jp1 jp1 = jp0 jp0 = jtmp ENDDO ! next, x - flux divergence i_start = its i_end = MIN(ite,ide-1) j_start = jts j_end = jte IF ( config_flags%open_ys .or. specified ) j_start = MAX(jds+1,jts) IF ( config_flags%open_ye .or. specified ) j_end = MIN(jde-1,jte) ! 3rd or 4th order flux has a 5 point stencil, so compute ! bounds so we can switch to second order flux close to the boundary i_start_f = i_start i_end_f = i_end+1 IF(degrade_xs) then i_start = ids+1 i_start_f = i_start+1 ENDIF IF(degrade_xe) then i_end = ide-2 i_end_f = ide-2 ENDIF ! compute fluxes DO j = j_start, j_end ! 3rd or 4th order flux DO k=kts,ktf DO i = i_start_f, i_end_f vel = 0.5*(ru(i,k,j)+ru(i,k,j-1)) fqx( i, k ) = vel*flux3( v(i-2,k,j), v(i-1,k,j), & v(i ,k,j), v(i+1,k,j), & vel ) ENDDO ENDDO ! second order flux close to boundaries (if not periodic or symmetric) IF( degrade_xs ) THEN DO k=kts,ktf fqx(i_start,k) = 0.25*(ru(i_start,k,j)+ru(i_start,k,j-1)) & *(v(i_start,k,j)+v(i_start-1,k,j)) ENDDO ENDIF IF( degrade_xe ) THEN DO k=kts,ktf fqx(i_end+1,k) = 0.25*(ru(i_end+1,k,j)+ru(i_end+1,k,j-1)) & *(v(i_end+1,k,j)+v(i_end,k,j)) ENDDO ENDIF ! x flux-divergence into tendency DO k=kts,ktf DO i = i_start, i_end mrdx=msfv(i,j)*rdx tendency(i,k,j) = tendency(i,k,j) - mrdx*(fqx(i+1,k)-fqx(i,k)) ENDDO ENDDO ENDDO ELSE IF( horz_order == 2 ) THEN i_start = its i_end = MIN(ite,ide-1) j_start = jts j_end = jte IF ( config_flags%open_ys ) j_start = MAX(jds+1,jts) IF ( config_flags%open_ye ) j_end = MIN(jde-1,jte) IF ( specified ) j_start = MAX(jds+2,jts) IF ( specified ) j_end = MIN(jde-2,jte) DO j = j_start, j_end DO k=kts,ktf DO i = i_start, i_end mrdy=msfv(i,j)*rdy tendency(i,k,j)=tendency(i,k,j) -mrdy*0.25 & *((rv(i,k,j+1)+rv(i,k,j ))*(v(i,k,j+1)+v(i,k,j )) & -(rv(i,k,j )+rv(i,k,j-1))*(v(i,k,j )+v(i,k,j-1))) ENDDO ENDDO ENDDO ! specified uses upstream normal wind at boundaries IF ( specified .AND. jts .LE. jds+1 ) THEN j = jds+1 DO k=kts,ktf DO i = i_start, i_end mrdy=msfv(i,j)*rdy vb = v(i,k,j-1) IF (v(i,k,j) .LT. 0.) vb = v(i,k,j) tendency(i,k,j)=tendency(i,k,j) -mrdy*0.25 & *((rv(i,k,j+1)+rv(i,k,j ))*(v(i,k,j+1)+v(i,k,j )) & -(rv(i,k,j )+rv(i,k,j-1))*(v(i,k,j )+vb)) ENDDO ENDDO ENDIF IF ( specified .AND. jte .GE. jde-1 ) THEN j = jde-1 DO k=kts,ktf DO i = i_start, i_end mrdy=msfv(i,j)*rdy vb = v(i,k,j+1) IF (v(i,k,j) .GT. 0.) vb = v(i,k,j) tendency(i,k,j)=tendency(i,k,j) -mrdy*0.25 & *((rv(i,k,j+1)+rv(i,k,j ))*(vb+v(i,k,j )) & -(rv(i,k,j )+rv(i,k,j-1))*(v(i,k,j )+v(i,k,j-1))) ENDDO ENDDO ENDIF IF ( .NOT. config_flags%periodic_x ) THEN IF ( config_flags%open_xs .or. specified ) i_start = MAX(ids+1,its) IF ( config_flags%open_xe .or. specified ) i_end = MIN(ide-2,ite) ENDIF DO j = j_start, j_end DO k=kts,ktf DO i = i_start, i_end mrdx=msfv(i,j)*rdx tendency(i,k,j)=tendency(i,k,j)-mrdx*0.25 & *((ru(i+1,k,j)+ru(i+1,k,j-1))*(v(i+1,k,j)+v(i ,k,j)) & -(ru(i ,k,j)+ru(i ,k,j-1))*(v(i ,k,j)+v(i-1,k,j))) ENDDO ENDDO ENDDO ELSE WRITE ( wrf_err_message , * ) 'module_advect: advect_v_6a: h_order not known ',horz_order CALL wrf_error_fatal ( TRIM( wrf_err_message ) ) ENDIF horizontal_order_test ! radiative lateral boundary condition in y for normal velocity (v) IF ( (config_flags%open_ys) .and. jts == jds ) THEN i_start = its i_end = MIN(ite,ide-1) DO i = i_start, i_end DO k = kts, ktf vb = MIN(rv(i,k,jts)-cb*mut(i,jts), 0.) tendency(i,k,jts) = tendency(i,k,jts) & - rdy*vb*(v_old(i,k,jts+1) - v_old(i,k,jts)) ENDDO ENDDO ENDIF IF ( (config_flags%open_ye) .and. jte == jde ) THEN i_start = its i_end = MIN(ite,ide-1) DO i = i_start, i_end DO k = kts, ktf vb = MAX(rv(i,k,jte)+cb*mut(i,jte-1), 0.) tendency(i,k,jte) = tendency(i,k,jte) & - rdy*vb*(v_old(i,k,jte) - v_old(i,k,jte-1)) ENDDO ENDDO ENDIF ! pick up the rest of the horizontal radiation boundary conditions. ! (these are the computations that don't require 'cb'. ! first, set to index ranges j_start = jts j_end = MIN(jte,jde) jmin = jds jmax = jde-1 IF (config_flags%open_ys) THEN j_start = MAX(jds+1, jts) jmin = jds ENDIF IF (config_flags%open_ye) THEN j_end = MIN(jte,jde-1) jmax = jde-1 ENDIF ! compute x (u) conditions for v, w, or scalar IF( (config_flags%open_xs) .and. (its == ids)) THEN DO j = j_start, j_end mrdx=msfv(its,j)*rdx jp = MIN( jmax, j ) jm = MAX( jmin, j-1 ) DO k=kts,ktf uw = 0.5*(ru(its,k,jp)+ru(its,k,jm)) ub = MIN( uw, 0. ) dup = ru(its+1,k,jp)-ru(its,k,jp) dum = ru(its+1,k,jm)-ru(its,k,jm) tendency(its,k,j)=tendency(its,k,j)-mrdx*( & ub*(v_old(its+1,k,j)-v_old(its,k,j)) & +0.5*v(its,k,j)*(dup+dum)) ENDDO ENDDO ENDIF IF( (config_flags%open_xe) .and. (ite == ide) ) THEN DO j = j_start, j_end mrdx=msfv(ite-1,j)*rdx jp = MIN( jmax, j ) jm = MAX( jmin, j-1 ) DO k=kts,ktf uw = 0.5*(ru(ite,k,jp)+ru(ite,k,jm)) ub = MAX( uw, 0. ) dup = ru(ite,k,jp)-ru(ite-1,k,jp) dum = ru(ite,k,jm)-ru(ite-1,k,jm) ! tendency(ite-1,k,j)=tendency(ite-1,k,j)-mrdx*( & ! ub*(v_old(ite-1,k,j)-v_old(ite-2,k,j)) & ! +0.5*v(ite-1,k,j)* & ! ( ru(ite,k,jp)-ru(ite-1,k,jp) & ! +ru(ite,k,jm)-ru(ite-1,k,jm)) ) tendency(ite-1,k,j)=tendency(ite-1,k,j)-mrdx*( & ub*(v_old(ite-1,k,j)-v_old(ite-2,k,j)) & +0.5*v(ite-1,k,j)*(dup+dum)) ENDDO ENDDO ENDIF !-------------------- vertical advection i_start = its i_end = MIN(ite,ide-1) j_start = jts j_end = jte DO i = i_start, i_end vflux(i,kts)=0. vflux(i,kte)=0. ENDDO IF ( config_flags%open_ys .or. specified ) j_start = MAX(jds+1,jts) IF ( config_flags%open_ye .or. specified ) j_end = MIN(jde-1,jte) vert_order_test : IF (vert_order == 6) THEN DO j = j_start, j_end DO k=kts+3,ktf-2 DO i = i_start, i_end vel=0.5*(rom(i,k,j)+rom(i,k,j-1)) vflux(i,k) = vel*flux6( & v(i,k-3,j), v(i,k-2,j), v(i,k-1,j), & v(i,k ,j), v(i,k+1,j), v(i,k+2,j), -vel ) ENDDO ENDDO DO i = i_start, i_end k=kts+1 vflux(i,k)=0.5*(rom(i,k,j)+rom(i,k,j-1)) & *(fzm(k)*v(i,k,j)+fzp(k)*v(i,k-1,j)) k = kts+2 vel=0.5*(rom(i,k,j)+rom(i,k,j-1)) vflux(i,k) = vel*flux4( & v(i,k-2,j), v(i,k-1,j), & v(i,k ,j), v(i,k+1,j), -vel ) k = ktf-1 vel=0.5*(rom(i,k,j)+rom(i,k,j-1)) vflux(i,k) = vel*flux4( & v(i,k-2,j), v(i,k-1,j), & v(i,k ,j), v(i,k+1,j), -vel ) k=ktf vflux(i,k)=0.5*(rom(i,k,j)+rom(i,k,j-1)) & *(fzm(k)*v(i,k,j)+fzp(k)*v(i,k-1,j)) ENDDO DO k=kts,ktf DO i = i_start, i_end tendency(i,k,j)=tendency(i,k,j)-rdzw(k)*(vflux(i,k+1)-vflux(i,k)) ENDDO ENDDO ENDDO ELSE IF (vert_order == 5) THEN DO j = j_start, j_end DO k=kts+3,ktf-2 DO i = i_start, i_end vel=0.5*(rom(i,k,j)+rom(i,k,j-1)) vflux(i,k) = vel*flux5( & v(i,k-3,j), v(i,k-2,j), v(i,k-1,j), & v(i,k ,j), v(i,k+1,j), v(i,k+2,j), -vel ) ENDDO ENDDO DO i = i_start, i_end k=kts+1 vflux(i,k)=0.5*(rom(i,k,j)+rom(i,k,j-1)) & *(fzm(k)*v(i,k,j)+fzp(k)*v(i,k-1,j)) k = kts+2 vel=0.5*(rom(i,k,j)+rom(i,k,j-1)) vflux(i,k) = vel*flux3( & v(i,k-2,j), v(i,k-1,j), & v(i,k ,j), v(i,k+1,j), -vel ) k = ktf-1 vel=0.5*(rom(i,k,j)+rom(i,k,j-1)) vflux(i,k) = vel*flux3( & v(i,k-2,j), v(i,k-1,j), & v(i,k ,j), v(i,k+1,j), -vel ) k=ktf vflux(i,k)=0.5*(rom(i,k,j)+rom(i,k,j-1)) & *(fzm(k)*v(i,k,j)+fzp(k)*v(i,k-1,j)) ENDDO DO k=kts,ktf DO i = i_start, i_end tendency(i,k,j)=tendency(i,k,j)-rdzw(k)*(vflux(i,k+1)-vflux(i,k)) ENDDO ENDDO ENDDO ELSE IF (vert_order == 4) THEN DO j = j_start, j_end DO k=kts+2,ktf-1 DO i = i_start, i_end vel=0.5*(rom(i,k,j)+rom(i,k,j-1)) vflux(i,k) = vel*flux4( & v(i,k-2,j), v(i,k-1,j), & v(i,k ,j), v(i,k+1,j), -vel ) ENDDO ENDDO DO i = i_start, i_end k=kts+1 vflux(i,k)=0.5*(rom(i,k,j)+rom(i,k,j-1)) & *(fzm(k)*v(i,k,j)+fzp(k)*v(i,k-1,j)) k=ktf vflux(i,k)=0.5*(rom(i,k,j)+rom(i,k,j-1)) & *(fzm(k)*v(i,k,j)+fzp(k)*v(i,k-1,j)) ENDDO DO k=kts,ktf DO i = i_start, i_end tendency(i,k,j)=tendency(i,k,j)-rdzw(k)*(vflux(i,k+1)-vflux(i,k)) ENDDO ENDDO ENDDO ELSE IF (vert_order == 3) THEN DO j = j_start, j_end DO k=kts+2,ktf-1 DO i = i_start, i_end vel=0.5*(rom(i,k,j)+rom(i,k,j-1)) vflux(i,k) = vel*flux3( & v(i,k-2,j), v(i,k-1,j), & v(i,k ,j), v(i,k+1,j), -vel ) ENDDO ENDDO DO i = i_start, i_end k=kts+1 vflux(i,k)=0.5*(rom(i,k,j)+rom(i,k,j-1)) & *(fzm(k)*v(i,k,j)+fzp(k)*v(i,k-1,j)) k=ktf vflux(i,k)=0.5*(rom(i,k,j)+rom(i,k,j-1)) & *(fzm(k)*v(i,k,j)+fzp(k)*v(i,k-1,j)) ENDDO DO k=kts,ktf DO i = i_start, i_end tendency(i,k,j)=tendency(i,k,j)-rdzw(k)*(vflux(i,k+1)-vflux(i,k)) ENDDO ENDDO ENDDO ELSE IF (vert_order == 2) THEN DO j = j_start, j_end DO k=kts+1,ktf DO i = i_start, i_end vflux(i,k)=0.5*(rom(i,k,j)+rom(i,k,j-1)) & *(fzm(k)*v(i,k,j)+fzp(k)*v(i,k-1,j)) ENDDO ENDDO DO k=kts,ktf DO i = i_start, i_end tendency(i,k,j)=tendency(i,k,j)-rdzw(k)*(vflux(i,k+1)-vflux(i,k)) ENDDO ENDDO ENDDO ELSE WRITE ( wrf_err_message , * ) 'module_advect: advect_v_6a: v_order not known ',vert_order CALL wrf_error_fatal ( TRIM( wrf_err_message ) ) ENDIF vert_order_test END SUBROUTINE advect_v !------------------------------------------------------------------- SUBROUTINE advect_scalar ( field, field_old, tendency, & ru, rv, rom, & mut, config_flags, & msfu, msfv, msft, & fzm, fzp, & rdx, rdy, rdzw, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) IMPLICIT NONE ! Input data TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: field, & field_old, & ru, & rv, & rom REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mut REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: msfu, & msfv, & msft REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fzm, & fzp, & rdzw REAL , INTENT(IN ) :: rdx, & rdy ! Local data INTEGER :: i, j, k, itf, jtf, ktf INTEGER :: i_start, i_end, j_start, j_end INTEGER :: i_start_f, i_end_f, j_start_f, j_end_f INTEGER :: jmin, jmax, jp, jm, imin, imax REAL :: mrdx, mrdy, ub, vb, uw, vw REAL , DIMENSION(its:ite, kts:kte) :: vflux REAL, DIMENSION( its:ite+1, kts:kte ) :: fqx REAL, DIMENSION( its:ite, kts:kte, 2 ) :: fqy INTEGER :: horz_order, vert_order LOGICAL :: degrade_xs, degrade_ys LOGICAL :: degrade_xe, degrade_ye INTEGER :: jp1, jp0, jtmp ! definition of flux operators, 3rd, 4th, 5th or 6th order REAL :: flux3, flux4, flux5, flux6 REAL :: q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua, vel flux4(q_im2, q_im1, q_i, q_ip1, ua) = & ( 7.*(q_i + q_im1) - (q_ip1 + q_im2) )/12.0 flux3(q_im2, q_im1, q_i, q_ip1, ua) = & flux4(q_im2, q_im1, q_i, q_ip1, ua) + & sign(1.,ua)*((q_ip1 - q_im2)-3.*(q_i-q_im1))/12.0 flux6(q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua) = & ( 37.*(q_i+q_im1) - 8.*(q_ip1+q_im2) & +(q_ip2+q_im3) )/60.0 flux5(q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua) = & flux6(q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua) & -sign(1.,ua)*( & (q_ip2-q_im3)-5.*(q_ip1-q_im2)+10.*(q_i-q_im1) )/60.0 LOGICAL :: specified specified = .false. if(config_flags%specified .or. config_flags%nested) specified = .true. ! set order for the advection schemes ktf=MIN(kte,kde-1) horz_order = config_flags%h_sca_adv_order vert_order = config_flags%v_sca_adv_order ! begin with horizontal flux divergence ! here is the choice of flux operators horizontal_order_test : IF( horz_order == 6 ) THEN ! determine boundary mods for flux operators ! We degrade the flux operators from 3rd/4th order ! to second order one gridpoint in from the boundaries for ! all boundary conditions except periodic and symmetry - these ! conditions have boundary zone data fill for correct application ! of the higher order flux stencils degrade_xs = .true. degrade_xe = .true. degrade_ys = .true. degrade_ye = .true. IF( config_flags%periodic_x .or. & config_flags%symmetric_xs .or. & (its > ids+2) ) degrade_xs = .false. IF( config_flags%periodic_x .or. & config_flags%symmetric_xe .or. & (ite < ide-3) ) degrade_xe = .false. IF( config_flags%periodic_y .or. & config_flags%symmetric_ys .or. & (jts > jds+2) ) degrade_ys = .false. IF( config_flags%periodic_y .or. & config_flags%symmetric_ye .or. & (jte < jde-3) ) degrade_ye = .false. !--------------- y - advection first ktf=MIN(kte,kde-1) i_start = its i_end = MIN(ite,ide-1) j_start = jts j_end = MIN(jte,jde-1) ! higher order flux has a 5 or 7 point stencil, so compute ! bounds so we can switch to second order flux close to the boundary j_start_f = j_start j_end_f = j_end+1 IF(degrade_ys) then j_start = MAX(jts,jds+1) j_start_f = jds+3 ENDIF IF(degrade_ye) then j_end = MIN(jte,jde-2) j_end_f = jde-3 ENDIF ! compute fluxes, 5th or 6th order jp1 = 2 jp0 = 1 j_loop_y_flux_6 : DO j = j_start, j_end+1 IF( (j >= j_start_f ) .and. (j <= j_end_f) ) THEN ! use full stencil DO k=kts,ktf DO i = i_start, i_end vel = rv(i,k,j) fqy( i, k, jp1 ) = vel*flux6( & field(i,k,j-3), field(i,k,j-2), field(i,k,j-1), & field(i,k,j ), field(i,k,j+1), field(i,k,j+2), vel ) ENDDO ENDDO ELSE IF ( j == jds+1 ) THEN ! 2nd order flux next to south boundary DO k=kts,ktf DO i = i_start, i_end fqy(i,k, jp1) = 0.5*rv(i,k,j)* & (field(i,k,j)+field(i,k,j-1)) ENDDO ENDDO ELSE IF ( j == jds+2 ) THEN ! third of 4th order flux 2 in from south boundary DO k=kts,ktf DO i = i_start, i_end vel = rv(i,k,j) fqy( i, k, jp1 ) = vel*flux4( & field(i,k,j-2),field(i,k,j-1),field(i,k,j),field(i,k,j+1),vel ) ENDDO ENDDO ELSE IF ( j == jde-1 ) THEN ! 2nd order flux next to north boundary DO k=kts,ktf DO i = i_start, i_end fqy(i, k, jp1) = 0.5*rv(i,k,j)* & (field(i,k,j)+field(i,k,j-1)) ENDDO ENDDO ELSE IF ( j == jde-2 ) THEN ! 3rd or 4th order flux 2 in from north boundary DO k=kts,ktf DO i = i_start, i_end vel = rv(i,k,j) fqy( i, k, jp1) = vel*flux4( & field(i,k,j-2),field(i,k,j-1), & field(i,k,j),field(i,k,j+1),vel ) ENDDO ENDDO ENDIF ! y flux-divergence into tendency IF(j > j_start) THEN DO k=kts,ktf DO i = i_start, i_end mrdy=msft(i,j-1)*rdy tendency(i,k,j-1) = tendency(i,k,j-1) - mrdy*(fqy(i,k,jp1)-fqy(i,k,jp0)) ENDDO ENDDO ENDIF jtmp = jp1 jp1 = jp0 jp0 = jtmp ENDDO j_loop_y_flux_6 ! next, x - flux divergence i_start = its i_end = MIN(ite,ide-1) j_start = jts j_end = MIN(jte,jde-1) ! higher order flux has a 5 or 7 point stencil, so compute ! bounds so we can switch to second order flux close to the boundary i_start_f = i_start i_end_f = i_end+1 IF(degrade_xs) then i_start = MAX(ids+1,its) i_start_f = i_start+2 ENDIF IF(degrade_xe) then i_end = MIN(ide-2,ite) i_end_f = ide-3 ENDIF ! compute fluxes DO j = j_start, j_end ! 5th or 6th order flux DO k=kts,ktf DO i = i_start_f, i_end_f vel = ru(i,k,j) fqx( i,k ) = vel*flux6( field(i-3,k,j), field(i-2,k,j), & field(i-1,k,j), field(i ,k,j), & field(i+1,k,j), field(i+2,k,j), & vel ) ENDDO ENDDO ! lower order fluxes close to boundaries (if not periodic or symmetric) IF( degrade_xs ) THEN IF( i_start == ids+1 ) THEN ! second order flux next to the boundary i = ids+1 DO k=kts,ktf fqx(i,k) = 0.5*(ru(i,k,j)) & *(field(i,k,j)+field(i-1,k,j)) ENDDO ENDIF i = ids+2 DO k=kts,ktf vel = ru(i,k,j) fqx( i,k ) = vel*flux4( field(i-2,k,j), field(i-1,k,j), & field(i ,k,j), field(i+1,k,j), & vel ) ENDDO ENDIF IF( degrade_xe ) THEN IF( i_end == ide-2 ) THEN ! second order flux next to the boundary i = ide-1 DO k=kts,ktf fqx(i,k) = 0.5*(ru(i,k,j)) & *(field(i,k,j)+field(i-1,k,j)) ENDDO ENDIF i = ide-2 DO k=kts,ktf vel = ru(i,k,j) fqx( i,k ) = vel*flux4( field(i-2,k,j), field(i-1,k,j), & field(i ,k,j), field(i+1,k,j), & vel ) ENDDO ENDIF ! x flux-divergence into tendency DO k=kts,ktf DO i = i_start, i_end mrdx=msft(i,j)*rdx tendency(i,k,j) = tendency(i,k,j) - mrdx*(fqx(i+1,k)-fqx(i,k)) ENDDO ENDDO ENDDO ELSE IF( horz_order == 5 ) THEN ! determine boundary mods for flux operators ! We degrade the flux operators from 3rd/4th order ! to second order one gridpoint in from the boundaries for ! all boundary conditions except periodic and symmetry - these ! conditions have boundary zone data fill for correct application ! of the higher order flux stencils degrade_xs = .true. degrade_xe = .true. degrade_ys = .true. degrade_ye = .true. IF( config_flags%periodic_x .or. & config_flags%symmetric_xs .or. & (its > ids+2) ) degrade_xs = .false. IF( config_flags%periodic_x .or. & config_flags%symmetric_xe .or. & (ite < ide-3) ) degrade_xe = .false. IF( config_flags%periodic_y .or. & config_flags%symmetric_ys .or. & (jts > jds+2) ) degrade_ys = .false. IF( config_flags%periodic_y .or. & config_flags%symmetric_ye .or. & (jte < jde-3) ) degrade_ye = .false. !--------------- y - advection first ktf=MIN(kte,kde-1) i_start = its i_end = MIN(ite,ide-1) j_start = jts j_end = MIN(jte,jde-1) ! higher order flux has a 5 or 7 point stencil, so compute ! bounds so we can switch to second order flux close to the boundary j_start_f = j_start j_end_f = j_end+1 IF(degrade_ys) then j_start = MAX(jts,jds+1) j_start_f = jds+3 ENDIF IF(degrade_ye) then j_end = MIN(jte,jde-2) j_end_f = jde-3 ENDIF ! compute fluxes, 5th or 6th order jp1 = 2 jp0 = 1 j_loop_y_flux_5 : DO j = j_start, j_end+1 IF( (j >= j_start_f ) .and. (j <= j_end_f) ) THEN ! use full stencil DO k=kts,ktf DO i = i_start, i_end vel = rv(i,k,j) fqy( i, k, jp1 ) = vel*flux5( & field(i,k,j-3), field(i,k,j-2), field(i,k,j-1), & field(i,k,j ), field(i,k,j+1), field(i,k,j+2), vel ) ENDDO ENDDO ELSE IF ( j == jds+1 ) THEN ! 2nd order flux next to south boundary DO k=kts,ktf DO i = i_start, i_end fqy(i,k, jp1) = 0.5*rv(i,k,j)* & (field(i,k,j)+field(i,k,j-1)) ENDDO ENDDO ELSE IF ( j == jds+2 ) THEN ! third of 4th order flux 2 in from south boundary DO k=kts,ktf DO i = i_start, i_end vel = rv(i,k,j) fqy( i, k, jp1 ) = vel*flux3( & field(i,k,j-2),field(i,k,j-1),field(i,k,j),field(i,k,j+1),vel ) ENDDO ENDDO ELSE IF ( j == jde-1 ) THEN ! 2nd order flux next to north boundary DO k=kts,ktf DO i = i_start, i_end fqy(i, k, jp1) = 0.5*rv(i,k,j)* & (field(i,k,j)+field(i,k,j-1)) ENDDO ENDDO ELSE IF ( j == jde-2 ) THEN ! 3rd or 4th order flux 2 in from north boundary DO k=kts,ktf DO i = i_start, i_end vel = rv(i,k,j) fqy( i, k, jp1) = vel*flux3( & field(i,k,j-2),field(i,k,j-1), & field(i,k,j),field(i,k,j+1),vel ) ENDDO ENDDO ENDIF ! y flux-divergence into tendency IF(j > j_start) THEN DO k=kts,ktf DO i = i_start, i_end mrdy=msft(i,j-1)*rdy tendency(i,k,j-1) = tendency(i,k,j-1) - mrdy*(fqy(i,k,jp1)-fqy(i,k,jp0)) ENDDO ENDDO ENDIF jtmp = jp1 jp1 = jp0 jp0 = jtmp ENDDO j_loop_y_flux_5 ! next, x - flux divergence i_start = its i_end = MIN(ite,ide-1) j_start = jts j_end = MIN(jte,jde-1) ! higher order flux has a 5 or 7 point stencil, so compute ! bounds so we can switch to second order flux close to the boundary i_start_f = i_start i_end_f = i_end+1 IF(degrade_xs) then i_start = MAX(ids+1,its) i_start_f = i_start+2 ENDIF IF(degrade_xe) then i_end = MIN(ide-2,ite) i_end_f = ide-3 ENDIF ! compute fluxes DO j = j_start, j_end ! 5th or 6th order flux DO k=kts,ktf DO i = i_start_f, i_end_f vel = ru(i,k,j) fqx( i,k ) = vel*flux5( field(i-3,k,j), field(i-2,k,j), & field(i-1,k,j), field(i ,k,j), & field(i+1,k,j), field(i+2,k,j), & vel ) ENDDO ENDDO ! lower order fluxes close to boundaries (if not periodic or symmetric) IF( degrade_xs ) THEN IF( i_start == ids+1 ) THEN ! second order flux next to the boundary i = ids+1 DO k=kts,ktf fqx(i,k) = 0.5*(ru(i,k,j)) & *(field(i,k,j)+field(i-1,k,j)) ENDDO ENDIF i = ids+2 DO k=kts,ktf vel = ru(i,k,j) fqx( i,k ) = vel*flux3( field(i-2,k,j), field(i-1,k,j), & field(i ,k,j), field(i+1,k,j), & vel ) ENDDO ENDIF IF( degrade_xe ) THEN IF( i_end == ide-2 ) THEN ! second order flux next to the boundary i = ide-1 DO k=kts,ktf fqx(i,k) = 0.5*(ru(i,k,j)) & *(field(i,k,j)+field(i-1,k,j)) ENDDO ENDIF i = ide-2 DO k=kts,ktf vel = ru(i,k,j) fqx( i,k ) = vel*flux3( field(i-2,k,j), field(i-1,k,j), & field(i ,k,j), field(i+1,k,j), & vel ) ENDDO ENDIF ! x flux-divergence into tendency DO k=kts,ktf DO i = i_start, i_end mrdx=msft(i,j)*rdx tendency(i,k,j) = tendency(i,k,j) - mrdx*(fqx(i+1,k)-fqx(i,k)) ENDDO ENDDO ENDDO ELSE IF( horz_order == 4 ) THEN degrade_xs = .true. degrade_xe = .true. degrade_ys = .true. degrade_ye = .true. IF( config_flags%periodic_x .or. & config_flags%symmetric_xs .or. & (its > ids+1) ) degrade_xs = .false. IF( config_flags%periodic_x .or. & config_flags%symmetric_xe .or. & (ite < ide-2) ) degrade_xe = .false. IF( config_flags%periodic_y .or. & config_flags%symmetric_ys .or. & (jts > jds+1) ) degrade_ys = .false. IF( config_flags%periodic_y .or. & config_flags%symmetric_ye .or. & (jte < jde-2) ) degrade_ye = .false. ! begin flux computations ! start with x flux divergence ktf=MIN(kte,kde-1) i_start = its i_end = MIN(ite,ide-1) j_start = jts j_end = MIN(jte,jde-1) ! 3rd or 4th order flux has a 5 point stencil, so compute ! bounds so we can switch to second order flux close to the boundary i_start_f = i_start i_end_f = i_end+1 IF(degrade_xs) then i_start = ids+1 i_start_f = i_start+1 ENDIF IF(degrade_xe) then i_end = ide-2 i_end_f = ide-2 ENDIF ! compute fluxes DO j = j_start, j_end ! 3rd or 4th order flux DO k=kts,ktf DO i = i_start_f, i_end_f fqx( i, k) = ru(i,k,j)*flux4( field(i-2,k,j), field(i-1,k,j), & field(i ,k,j), field(i+1,k,j), & ru(i,k,j) ) ENDDO ENDDO ! second order flux close to boundaries (if not periodic or symmetric) IF( degrade_xs ) THEN DO k=kts,ktf fqx(i_start, k) = 0.5*ru(i_start,k,j) & *(field(i_start,k,j)+field(i_start-1,k,j)) ENDDO ENDIF IF( degrade_xe ) THEN DO k=kts,ktf fqx(i_end+1,k ) = 0.5*ru(i_end+1,k,j) & *(field(i_end+1,k,j)+field(i_end,k,j)) ENDDO ENDIF ! x flux-divergence into tendency DO k=kts,ktf DO i = i_start, i_end mrdx=msft(i,j)*rdx tendency(i,k,j) = tendency(i,k,j) - mrdx*(fqx(i+1,k)-fqx(i,k)) ENDDO ENDDO ENDDO ! next -> y flux divergence calculation i_start = its i_end = MIN(ite,ide-1) j_start = jts j_end = MIN(jte,jde-1) ! 3rd or 4th order flux has a 5 point stencil, so compute ! bounds so we can switch to second order flux close to the boundary j_start_f = j_start j_end_f = j_end+1 IF(degrade_ys) then j_start = jds+1 j_start_f = j_start+1 ENDIF IF(degrade_ye) then j_end = jde-2 j_end_f = jde-2 ENDIF jp1 = 2 jp0 = 1 DO j = j_start, j_end+1 IF ((j < j_start_f) .and. degrade_ys) THEN DO k = kts, ktf DO i = i_start, i_end fqy(i,k,jp1) = 0.5*rv(i,k,j_start) & *(field(i,k,j_start)+field(i,k,j_start-1)) ENDDO ENDDO ELSE IF ((j > j_end_f) .and. degrade_ye) THEN DO k = kts, ktf DO i = i_start, i_end fqy(i,k,jp1) = 0.5*rv(i,k,j_end+1) & *(field(i,k,j_end+1)+field(i,k,j_end)) ENDDO ENDDO ELSE ! 3rd or 4th order flux DO k = kts, ktf DO i = i_start, i_end fqy( i, k, jp1 ) = rv(i,k,j)*flux4( field(i,k,j-2), field(i,k,j-1), & field(i,k,j ), field(i,k,j+1), & rv(i,k,j) ) ENDDO ENDDO END IF IF ( j > j_start ) THEN ! y flux-divergence into tendency DO k=kts,ktf DO i = i_start, i_end mrdy=msft(i,j-1)*rdy tendency(i,k,j-1) = tendency(i,k,j-1) - mrdy*(fqy(i,k,jp1)-fqy(i,k,jp0)) ENDDO ENDDO END IF jtmp = jp1 jp1 = jp0 jp0 = jtmp ENDDO ELSE IF( horz_order == 3 ) THEN degrade_xs = .true. degrade_xe = .true. degrade_ys = .true. degrade_ye = .true. IF( config_flags%periodic_x .or. & config_flags%symmetric_xs .or. & (its > ids+1) ) degrade_xs = .false. IF( config_flags%periodic_x .or. & config_flags%symmetric_xe .or. & (ite < ide-2) ) degrade_xe = .false. IF( config_flags%periodic_y .or. & config_flags%symmetric_ys .or. & (jts > jds+1) ) degrade_ys = .false. IF( config_flags%periodic_y .or. & config_flags%symmetric_ye .or. & (jte < jde-2) ) degrade_ye = .false. ! begin flux computations ! start with x flux divergence ktf=MIN(kte,kde-1) i_start = its i_end = MIN(ite,ide-1) j_start = jts j_end = MIN(jte,jde-1) ! 3rd or 4th order flux has a 5 point stencil, so compute ! bounds so we can switch to second order flux close to the boundary i_start_f = i_start i_end_f = i_end+1 IF(degrade_xs) then i_start = ids+1 i_start_f = i_start+1 ENDIF IF(degrade_xe) then i_end = ide-2 i_end_f = ide-2 ENDIF ! compute fluxes DO j = j_start, j_end ! 3rd or 4th order flux DO k=kts,ktf DO i = i_start_f, i_end_f fqx( i, k) = ru(i,k,j)*flux3( field(i-2,k,j), field(i-1,k,j), & field(i ,k,j), field(i+1,k,j), & ru(i,k,j) ) ENDDO ENDDO ! second order flux close to boundaries (if not periodic or symmetric) IF( degrade_xs ) THEN DO k=kts,ktf fqx(i_start, k) = 0.5*ru(i_start,k,j) & *(field(i_start,k,j)+field(i_start-1,k,j)) ENDDO ENDIF IF( degrade_xe ) THEN DO k=kts,ktf fqx(i_end+1,k ) = 0.5*ru(i_end+1,k,j) & *(field(i_end+1,k,j)+field(i_end,k,j)) ENDDO ENDIF ! x flux-divergence into tendency DO k=kts,ktf DO i = i_start, i_end mrdx=msft(i,j)*rdx tendency(i,k,j) = tendency(i,k,j) - mrdx*(fqx(i+1,k)-fqx(i,k)) ENDDO ENDDO ENDDO ! next -> y flux divergence calculation i_start = its i_end = MIN(ite,ide-1) j_start = jts j_end = MIN(jte,jde-1) ! 3rd or 4th order flux has a 5 point stencil, so compute ! bounds so we can switch to second order flux close to the boundary j_start_f = j_start j_end_f = j_end+1 IF(degrade_ys) then j_start = jds+1 j_start_f = j_start+1 ENDIF IF(degrade_ye) then j_end = jde-2 j_end_f = jde-2 ENDIF jp1 = 2 jp0 = 1 DO j = j_start, j_end+1 IF ((j < j_start_f) .and. degrade_ys) THEN DO k = kts, ktf DO i = i_start, i_end fqy(i,k,jp1) = 0.5*rv(i,k,j_start) & *(field(i,k,j_start)+field(i,k,j_start-1)) ENDDO ENDDO ELSE IF ((j > j_end_f) .and. degrade_ye) THEN DO k = kts, ktf DO i = i_start, i_end fqy(i,k,jp1) = 0.5*rv(i,k,j_end+1) & *(field(i,k,j_end+1)+field(i,k,j_end)) ENDDO ENDDO ELSE ! 3rd or 4th order flux DO k = kts, ktf DO i = i_start, i_end fqy( i, k, jp1 ) = rv(i,k,j)*flux3( field(i,k,j-2), field(i,k,j-1), & field(i,k,j ), field(i,k,j+1), & rv(i,k,j) ) ENDDO ENDDO END IF IF ( j > j_start ) THEN ! y flux-divergence into tendency DO k=kts,ktf DO i = i_start, i_end mrdy=msft(i,j-1)*rdy tendency(i,k,j-1) = tendency(i,k,j-1) - mrdy*(fqy(i,k,jp1)-fqy(i,k,jp0)) ENDDO ENDDO END IF jtmp = jp1 jp1 = jp0 jp0 = jtmp ENDDO ELSE IF( horz_order == 2 ) THEN i_start = its i_end = MIN(ite,ide-1) j_start = jts j_end = MIN(jte,jde-1) IF ( .NOT. config_flags%periodic_x ) THEN IF ( config_flags%open_xs .or. specified ) i_start = MAX(ids+1,its) IF ( config_flags%open_xe .or. specified ) i_end = MIN(ide-2,ite) ENDIF DO j = j_start, j_end DO k = kts, ktf DO i = i_start, i_end mrdx=msft(i,j)*rdx tendency(i,k,j)=tendency(i,k,j)-mrdx*0.5 & *(ru(i+1,k,j)*(field(i+1,k,j)+field(i ,k,j)) & -ru(i ,k,j)*(field(i ,k,j)+field(i-1,k,j))) ENDDO ENDDO ENDDO i_start = its i_end = MIN(ite,ide-1) IF ( config_flags%open_ys .or. specified ) j_start = MAX(jds+1,jts) IF ( config_flags%open_ye .or. specified ) j_end = MIN(jde-2,jte) DO j = j_start, j_end DO k = kts, ktf DO i = i_start, i_end mrdy=msft(i,j)*rdy tendency(i,k,j)=tendency(i,k,j) -mrdy*0.5 & *(rv(i,k,j+1)*(field(i,k,j+1)+field(i,k,j )) & -rv(i,k,j )*(field(i,k,j )+field(i,k,j-1))) ENDDO ENDDO ENDDO ELSE WRITE ( wrf_err_message , * ) 'module_advect: advect_scalar_6a, h_order not known ',horz_order CALL wrf_error_fatal ( TRIM( wrf_err_message ) ) ENDIF horizontal_order_test ! pick up the rest of the horizontal radiation boundary conditions. ! (these are the computations that don't require 'cb'. ! first, set to index ranges i_start = its i_end = MIN(ite,ide-1) j_start = jts j_end = MIN(jte,jde-1) ! compute x (u) conditions for v, w, or scalar IF( (config_flags%open_xs) .and. (its == ids) ) THEN DO j = j_start, j_end DO k = kts, ktf ub = MIN( 0.5*(ru(its,k,j)+ru(its+1,k,j)), 0. ) tendency(its,k,j) = tendency(its,k,j) & - rdx*( & ub*( field_old(its+1,k,j) & - field_old(its ,k,j) ) + & field(its,k,j)*(ru(its+1,k,j)-ru(its,k,j)) & ) ENDDO ENDDO ENDIF IF( (config_flags%open_xe) .and. (ite == ide) ) THEN DO j = j_start, j_end DO k = kts, ktf ub = MAX( 0.5*(ru(ite-1,k,j)+ru(ite,k,j)), 0. ) tendency(i_end,k,j) = tendency(i_end,k,j) & - rdx*( & ub*( field_old(i_end ,k,j) & - field_old(i_end-1,k,j) ) + & field(i_end,k,j)*(ru(ite,k,j)-ru(ite-1,k,j)) & ) ENDDO ENDDO ENDIF IF( (config_flags%open_ys) .and. (jts == jds) ) THEN DO i = i_start, i_end DO k = kts, ktf vb = MIN( 0.5*(rv(i,k,jts)+rv(i,k,jts+1)), 0. ) tendency(i,k,jts) = tendency(i,k,jts) & - rdy*( & vb*( field_old(i,k,jts+1) & - field_old(i,k,jts ) ) + & field(i,k,jts)*(rv(i,k,jts+1)-rv(i,k,jts)) & ) ENDDO ENDDO ENDIF IF( (config_flags%open_ye) .and. (jte == jde)) THEN DO i = i_start, i_end DO k = kts, ktf vb = MAX( 0.5*(rv(i,k,jte-1)+rv(i,k,jte)), 0. ) tendency(i,k,j_end) = tendency(i,k,j_end) & - rdy*( & vb*( field_old(i,k,j_end ) & - field_old(i,k,j_end-1) ) + & field(i,k,j_end)*(rv(i,k,jte)-rv(i,k,jte-1)) & ) ENDDO ENDDO ENDIF !-------------------- vertical advection i_start = its i_end = MIN(ite,ide-1) j_start = jts j_end = MIN(jte,jde-1) DO i = i_start, i_end vflux(i,kts)=0. vflux(i,kte)=0. ENDDO vert_order_test : IF (vert_order == 6) THEN DO j = j_start, j_end DO k=kts+3,ktf-2 DO i = i_start, i_end vel=rom(i,k,j) vflux(i,k) = vel*flux6( & field(i,k-3,j), field(i,k-2,j), field(i,k-1,j), & field(i,k ,j), field(i,k+1,j), field(i,k+2,j), -vel ) ENDDO ENDDO DO i = i_start, i_end k=kts+1 vflux(i,k)=rom(i,k,j)*(fzm(k)*field(i,k,j)+fzp(k)*field(i,k-1,j)) k = kts+2 vel=rom(i,k,j) vflux(i,k) = vel*flux4( & field(i,k-2,j), field(i,k-1,j), & field(i,k ,j), field(i,k+1,j), -vel ) k = ktf-1 vel=rom(i,k,j) vflux(i,k) = vel*flux4( & field(i,k-2,j), field(i,k-1,j), & field(i,k ,j), field(i,k+1,j), -vel ) k=ktf vflux(i,k)=rom(i,k,j)*(fzm(k)*field(i,k,j)+fzp(k)*field(i,k-1,j)) ENDDO DO k=kts,ktf DO i = i_start, i_end tendency(i,k,j)=tendency(i,k,j)-rdzw(k)*(vflux(i,k+1)-vflux(i,k)) ENDDO ENDDO ENDDO ELSE IF (vert_order == 5) THEN DO j = j_start, j_end DO k=kts+3,ktf-2 DO i = i_start, i_end vel=rom(i,k,j) vflux(i,k) = vel*flux5( & field(i,k-3,j), field(i,k-2,j), field(i,k-1,j), & field(i,k ,j), field(i,k+1,j), field(i,k+2,j), -vel ) ENDDO ENDDO DO i = i_start, i_end k=kts+1 vflux(i,k)=rom(i,k,j)*(fzm(k)*field(i,k,j)+fzp(k)*field(i,k-1,j)) k = kts+2 vel=rom(i,k,j) vflux(i,k) = vel*flux3( & field(i,k-2,j), field(i,k-1,j), & field(i,k ,j), field(i,k+1,j), -vel ) k = ktf-1 vel=rom(i,k,j) vflux(i,k) = vel*flux3( & field(i,k-2,j), field(i,k-1,j), & field(i,k ,j), field(i,k+1,j), -vel ) k=ktf vflux(i,k)=rom(i,k,j)*(fzm(k)*field(i,k,j)+fzp(k)*field(i,k-1,j)) ENDDO DO k=kts,ktf DO i = i_start, i_end tendency(i,k,j)=tendency(i,k,j)-rdzw(k)*(vflux(i,k+1)-vflux(i,k)) ENDDO ENDDO ENDDO ELSE IF (vert_order == 4) THEN DO j = j_start, j_end DO k=kts+2,ktf-1 DO i = i_start, i_end vel=rom(i,k,j) vflux(i,k) = vel*flux4( & field(i,k-2,j), field(i,k-1,j), & field(i,k ,j), field(i,k+1,j), -vel ) ENDDO ENDDO DO i = i_start, i_end k=kts+1 vflux(i,k)=rom(i,k,j)*(fzm(k)*field(i,k,j)+fzp(k)*field(i,k-1,j)) k=ktf vflux(i,k)=rom(i,k,j)*(fzm(k)*field(i,k,j)+fzp(k)*field(i,k-1,j)) ENDDO DO k=kts,ktf DO i = i_start, i_end tendency(i,k,j)=tendency(i,k,j)-rdzw(k)*(vflux(i,k+1)-vflux(i,k)) ENDDO ENDDO ENDDO ELSE IF (vert_order == 3) THEN DO j = j_start, j_end DO k=kts+2,ktf-1 DO i = i_start, i_end vel=rom(i,k,j) vflux(i,k) = vel*flux3( & field(i,k-2,j), field(i,k-1,j), & field(i,k ,j), field(i,k+1,j), -vel ) ENDDO ENDDO DO i = i_start, i_end k=kts+1 vflux(i,k)=rom(i,k,j)*(fzm(k)*field(i,k,j)+fzp(k)*field(i,k-1,j)) k=ktf vflux(i,k)=rom(i,k,j)*(fzm(k)*field(i,k,j)+fzp(k)*field(i,k-1,j)) ENDDO DO k=kts,ktf DO i = i_start, i_end tendency(i,k,j)=tendency(i,k,j)-rdzw(k)*(vflux(i,k+1)-vflux(i,k)) ENDDO ENDDO ENDDO ELSE IF (vert_order == 2) THEN DO j = j_start, j_end DO k = kts+1, ktf DO i = i_start, i_end vflux(i,k)=rom(i,k,j)*(fzm(k)*field(i,k,j)+fzp(k)*field(i,k-1,j)) ENDDO ENDDO DO k = kts, ktf DO i = i_start, i_end tendency(i,k,j)=tendency(i,k,j)-rdzw(k)*(vflux(i,k+1)-vflux(i,k)) ENDDO ENDDO ENDDO ELSE WRITE (wrf_err_message,*) ' advect_scalar_6a, v_order not known ',vert_order CALL wrf_error_fatal ( wrf_err_message ) ENDIF vert_order_test END SUBROUTINE advect_scalar !--------------------------------------------------------------------------------- SUBROUTINE advect_w ( w, w_old, tendency, & ru, rv, rom, & mut, config_flags, & msfu, msfv, msft, & fzm, fzp, & rdx, rdy, rdzu, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) IMPLICIT NONE ! Input data TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: w, & w_old, & ru, & rv, & rom REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mut REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: msfu, & msfv, & msft REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fzm, & fzp, & rdzu REAL , INTENT(IN ) :: rdx, & rdy ! Local data INTEGER :: i, j, k, itf, jtf, ktf INTEGER :: i_start, i_end, j_start, j_end INTEGER :: i_start_f, i_end_f, j_start_f, j_end_f INTEGER :: jmin, jmax, jp, jm, imin, imax REAL :: mrdx, mrdy, ub, vb, uw, vw REAL , DIMENSION(its:ite, kts:kte) :: vflux INTEGER :: horz_order, vert_order REAL, DIMENSION( its:ite+1, kts:kte ) :: fqx REAL, DIMENSION( its:ite, kts:kte, 2 ) :: fqy LOGICAL :: degrade_xs, degrade_ys LOGICAL :: degrade_xe, degrade_ye INTEGER :: jp1, jp0, jtmp ! definition of flux operators, 3rd, 4th, 5th or 6th order REAL :: flux3, flux4, flux5, flux6 REAL :: q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua, vel flux4(q_im2, q_im1, q_i, q_ip1, ua) = & ( 7.*(q_i + q_im1) - (q_ip1 + q_im2) )/12.0 flux3(q_im2, q_im1, q_i, q_ip1, ua) = & flux4(q_im2, q_im1, q_i, q_ip1, ua) + & sign(1.,ua)*((q_ip1 - q_im2)-3.*(q_i-q_im1))/12.0 flux6(q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua) = & ( 37.*(q_i+q_im1) - 8.*(q_ip1+q_im2) & +(q_ip2+q_im3) )/60.0 flux5(q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua) = & flux6(q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua) & -sign(1.,ua)*( & (q_ip2-q_im3)-5.*(q_ip1-q_im2)+10.*(q_i-q_im1) )/60.0 LOGICAL :: specified specified = .false. if(config_flags%specified .or. config_flags%nested) specified = .true. ! set order for the advection scheme ktf=MIN(kte,kde-1) horz_order = config_flags%h_sca_adv_order vert_order = config_flags%v_sca_adv_order ! here is the choice of flux operators ! begin with horizontal flux divergence horizontal_order_test : IF( horz_order == 6 ) THEN ! determine boundary mods for flux operators ! We degrade the flux operators from 3rd/4th order ! to second order one gridpoint in from the boundaries for ! all boundary conditions except periodic and symmetry - these ! conditions have boundary zone data fill for correct application ! of the higher order flux stencils degrade_xs = .true. degrade_xe = .true. degrade_ys = .true. degrade_ye = .true. IF( config_flags%periodic_x .or. & config_flags%symmetric_xs .or. & (its > ids+2) ) degrade_xs = .false. IF( config_flags%periodic_x .or. & config_flags%symmetric_xe .or. & (ite < ide-3) ) degrade_xe = .false. IF( config_flags%periodic_y .or. & config_flags%symmetric_ys .or. & (jts > jds+2) ) degrade_ys = .false. IF( config_flags%periodic_y .or. & config_flags%symmetric_ye .or. & (jte < jde-3) ) degrade_ye = .false. !--------------- y - advection first i_start = its i_end = MIN(ite,ide-1) j_start = jts j_end = MIN(jte,jde-1) ! higher order flux has a 5 or 7 point stencil, so compute ! bounds so we can switch to second order flux close to the boundary j_start_f = j_start j_end_f = j_end+1 IF(degrade_ys) then j_start = MAX(jts,jds+1) j_start_f = jds+3 ENDIF IF(degrade_ye) then j_end = MIN(jte,jde-2) j_end_f = jde-3 ENDIF ! compute fluxes, 5th or 6th order jp1 = 2 jp0 = 1 j_loop_y_flux_6 : DO j = j_start, j_end+1 IF( (j >= j_start_f ) .and. (j <= j_end_f) ) THEN DO k=kts+1,ktf DO i = i_start, i_end vel = fzm(k)*rv(i,k,j)+fzp(k)*rv(i,k-1,j) fqy( i, k, jp1 ) = vel*flux6( & w(i,k,j-3), w(i,k,j-2), w(i,k,j-1), & w(i,k,j ), w(i,k,j+1), w(i,k,j+2), vel ) ENDDO ENDDO k = ktf+1 DO i = i_start, i_end vel = (2.-fzm(k-1))*rv(i,k-1,j)-fzp(k-1)*rv(i,k-2,j) fqy( i, k, jp1 ) = vel*flux6( & w(i,k,j-3), w(i,k,j-2), w(i,k,j-1), & w(i,k,j ), w(i,k,j+1), w(i,k,j+2), vel ) ENDDO ELSE IF ( j == jds+1 ) THEN ! 2nd order flux next to south boundary DO k=kts+1,ktf DO i = i_start, i_end fqy(i, k, jp1) = 0.5*(fzm(k)*rv(i,k,j)+fzp(k)*rv(i,k-1,j))* & (w(i,k,j)+w(i,k,j-1)) ENDDO ENDDO k = ktf+1 DO i = i_start, i_end fqy(i, k, jp1) = 0.5*((2.-fzm(k-1))*rv(i,k-1,j)-fzp(k-1)*rv(i,k-2,j))* & (w(i,k,j)+w(i,k,j-1)) ENDDO ELSE IF ( j == jds+2 ) THEN ! third of 4th order flux 2 in from south boundary DO k=kts+1,ktf DO i = i_start, i_end vel = fzm(k)*rv(i,k,j)+fzp(k)*rv(i,k-1,j) fqy( i, k, jp1 ) = vel*flux4( & w(i,k,j-2),w(i,k,j-1),w(i,k,j),w(i,k,j+1),vel ) ENDDO ENDDO k = ktf+1 DO i = i_start, i_end vel = (2.-fzm(k-1))*rv(i,k-1,j)-fzp(k-1)*rv(i,k-2,j) fqy( i, k, jp1 ) = vel*flux4( & w(i,k,j-2),w(i,k,j-1),w(i,k,j),w(i,k,j+1),vel ) ENDDO ELSE IF ( j == jde-1 ) THEN ! 2nd order flux next to north boundary DO k=kts+1,ktf DO i = i_start, i_end fqy(i, k, jp1) = 0.5*(fzm(k)*rv(i,k,j)+fzp(k)*rv(i,k-1,j))* & (w(i,k,j)+w(i,k,j-1)) ENDDO ENDDO k = ktf+1 DO i = i_start, i_end fqy(i, k, jp1) = 0.5*((2.-fzm(k-1))*rv(i,k-1,j)-fzp(k-1)*rv(i,k-2,j))* & (w(i,k,j)+w(i,k,j-1)) ENDDO ELSE IF ( j == jde-2 ) THEN ! 3rd or 4th order flux 2 in from north boundary DO k=kts+1,ktf DO i = i_start, i_end vel = fzm(k)*rv(i,k,j)+fzp(k)*rv(i,k-1,j) fqy( i, k, jp1 ) = vel*flux4( & w(i,k,j-2),w(i,k,j-1), & w(i,k,j),w(i,k,j+1),vel ) ENDDO ENDDO k = ktf+1 DO i = i_start, i_end vel = (2.-fzm(k-1))*rv(i,k-1,j)-fzp(k-1)*rv(i,k-2,j) fqy( i, k, jp1 ) = vel*flux4( & w(i,k,j-2),w(i,k,j-1), & w(i,k,j),w(i,k,j+1),vel ) ENDDO ENDIF ! y flux-divergence into tendency IF(j > j_start) THEN DO k=kts+1,ktf+1 DO i = i_start, i_end mrdy=msft(i,j-1)*rdy tendency(i,k,j-1) = tendency(i,k,j-1) - mrdy*(fqy(i,k,jp1)-fqy(i,k,jp0)) ENDDO ENDDO ENDIF jtmp = jp1 jp1 = jp0 jp0 = jtmp ENDDO j_loop_y_flux_6 ! next, x - flux divergence i_start = its i_end = MIN(ite,ide-1) j_start = jts j_end = MIN(jte,jde-1) ! higher order flux has a 5 or 7 point stencil, so compute ! bounds so we can switch to second order flux close to the boundary i_start_f = i_start i_end_f = i_end+1 IF(degrade_xs) then i_start = MAX(ids+1,its) i_start_f = i_start+2 ENDIF IF(degrade_xe) then i_end = MIN(ide-2,ite) i_end_f = ide-3 ENDIF ! compute fluxes DO j = j_start, j_end ! 5th or 6th order flux DO k=kts+1,ktf DO i = i_start_f, i_end_f vel = fzm(k)*ru(i,k,j)+fzp(k)*ru(i,k-1,j) fqx( i,k ) = vel*flux6( w(i-3,k,j), w(i-2,k,j), & w(i-1,k,j), w(i ,k,j), & w(i+1,k,j), w(i+2,k,j), & vel ) ENDDO ENDDO k = ktf+1 DO i = i_start_f, i_end_f vel = (2.-fzm(k-1))*ru(i,k-1,j)-fzp(k-1)*ru(i,k-2,j) fqx( i,k ) = vel*flux6( w(i-3,k,j), w(i-2,k,j), & w(i-1,k,j), w(i ,k,j), & w(i+1,k,j), w(i+2,k,j), & vel ) ENDDO ! lower order fluxes close to boundaries (if not periodic or symmetric) IF( degrade_xs ) THEN IF( i_start == ids+1 ) THEN ! second order flux next to the boundary i = ids+1 DO k=kts+1,ktf fqx(i,k) = 0.5*(fzm(k)*ru(i,k,j)+fzp(k)*ru(i,k-1,j)) & *(w(i,k,j)+w(i-1,k,j)) ENDDO k = ktf+1 fqx(i,k) = 0.5*((2.-fzm(k-1))*ru(i,k-1,j)-fzp(k-1)*ru(i,k-2,j)) & *(w(i,k,j)+w(i-1,k,j)) ENDIF DO k=kts+1,ktf i = i_start+1 vel = fzm(k)*ru(i,k,j)+fzp(k)*ru(i,k-1,j) fqx( i,k ) = vel*flux4( w(i-2,k,j), w(i-1,k,j), & w(i ,k,j), w(i+1,k,j), & vel ) ENDDO k = ktf+1 vel = (2.-fzm(k-1))*ru(i,k-1,j)-fzp(k-1)*ru(i,k-2,j) fqx( i,k ) = vel*flux4( w(i-2,k,j), w(i-1,k,j), & w(i ,k,j), w(i+1,k,j), & vel ) ENDIF IF( degrade_xe ) THEN IF( i_end == ide-2 ) THEN ! second order flux next to the boundary i = ide-1 DO k=kts+1,ktf fqx(i,k) = 0.5*(fzm(k)*ru(i,k,j)+fzp(k)*ru(i,k-1,j)) & *(w(i,k,j)+w(i-1,k,j)) ENDDO k = ktf+1 fqx(i,k) = 0.5*((2.-fzm(k-1))*ru(i,k-1,j)-fzp(k-1)*ru(i,k-2,j)) & *(w(i,k,j)+w(i-1,k,j)) ENDIF i = ide-2 DO k=kts+1,ktf vel = fzm(k)*ru(i,k,j)+fzp(k)*ru(i,k-1,j) fqx( i,k ) = vel*flux4( w(i-2,k,j), w(i-1,k,j), & w(i ,k,j), w(i+1,k,j), & vel ) ENDDO k = ktf+1 vel = (2.-fzm(k-1))*ru(i,k-1,j)-fzp(k-1)*ru(i,k-2,j) fqx( i,k ) = vel*flux4( w(i-2,k,j), w(i-1,k,j), & w(i ,k,j), w(i+1,k,j), & vel ) ENDIF ! x flux-divergence into tendency DO k=kts+1,ktf+1 DO i = i_start, i_end mrdx=msft(i,j)*rdx tendency(i,k,j) = tendency(i,k,j) - mrdx*(fqx(i+1,k)-fqx(i,k)) ENDDO ENDDO ENDDO ELSE IF (horz_order == 5 ) THEN ! determine boundary mods for flux operators ! We degrade the flux operators from 3rd/4th order ! to second order one gridpoint in from the boundaries for ! all boundary conditions except periodic and symmetry - these ! conditions have boundary zone data fill for correct application ! of the higher order flux stencils degrade_xs = .true. degrade_xe = .true. degrade_ys = .true. degrade_ye = .true. IF( config_flags%periodic_x .or. & config_flags%symmetric_xs .or. & (its > ids+2) ) degrade_xs = .false. IF( config_flags%periodic_x .or. & config_flags%symmetric_xe .or. & (ite < ide-3) ) degrade_xe = .false. IF( config_flags%periodic_y .or. & config_flags%symmetric_ys .or. & (jts > jds+2) ) degrade_ys = .false. IF( config_flags%periodic_y .or. & config_flags%symmetric_ye .or. & (jte < jde-3) ) degrade_ye = .false. !--------------- y - advection first i_start = its i_end = MIN(ite,ide-1) j_start = jts j_end = MIN(jte,jde-1) ! higher order flux has a 5 or 7 point stencil, so compute ! bounds so we can switch to second order flux close to the boundary j_start_f = j_start j_end_f = j_end+1 IF(degrade_ys) then j_start = MAX(jts,jds+1) j_start_f = jds+3 ENDIF IF(degrade_ye) then j_end = MIN(jte,jde-2) j_end_f = jde-3 ENDIF ! compute fluxes, 5th or 6th order jp1 = 2 jp0 = 1 j_loop_y_flux_5 : DO j = j_start, j_end+1 IF( (j >= j_start_f ) .and. (j <= j_end_f) ) THEN DO k=kts+1,ktf DO i = i_start, i_end vel = fzm(k)*rv(i,k,j)+fzp(k)*rv(i,k-1,j) fqy( i, k, jp1 ) = vel*flux5( & w(i,k,j-3), w(i,k,j-2), w(i,k,j-1), & w(i,k,j ), w(i,k,j+1), w(i,k,j+2), vel ) ENDDO ENDDO k = ktf+1 DO i = i_start, i_end vel = (2.-fzm(k-1))*rv(i,k-1,j)-fzp(k-1)*rv(i,k-2,j) fqy( i, k, jp1 ) = vel*flux5( & w(i,k,j-3), w(i,k,j-2), w(i,k,j-1), & w(i,k,j ), w(i,k,j+1), w(i,k,j+2), vel ) ENDDO ELSE IF ( j == jds+1 ) THEN ! 2nd order flux next to south boundary DO k=kts+1,ktf DO i = i_start, i_end fqy(i, k, jp1) = 0.5*(fzm(k)*rv(i,k,j)+fzp(k)*rv(i,k-1,j))* & (w(i,k,j)+w(i,k,j-1)) ENDDO ENDDO k = ktf+1 DO i = i_start, i_end fqy(i, k, jp1) = 0.5*((2.-fzm(k-1))*rv(i,k-1,j)-fzp(k-1)*rv(i,k-2,j))* & (w(i,k,j)+w(i,k,j-1)) ENDDO ELSE IF ( j == jds+2 ) THEN ! third of 4th order flux 2 in from south boundary DO k=kts+1,ktf DO i = i_start, i_end vel = fzm(k)*rv(i,k,j)+fzp(k)*rv(i,k-1,j) fqy( i, k, jp1 ) = vel*flux3( & w(i,k,j-2),w(i,k,j-1),w(i,k,j),w(i,k,j+1),vel ) ENDDO ENDDO k = ktf+1 DO i = i_start, i_end vel = (2.-fzm(k-1))*rv(i,k-1,j)-fzp(k-1)*rv(i,k-2,j) fqy( i, k, jp1 ) = vel*flux3( & w(i,k,j-2),w(i,k,j-1),w(i,k,j),w(i,k,j+1),vel ) ENDDO ELSE IF ( j == jde-1 ) THEN ! 2nd order flux next to north boundary DO k=kts+1,ktf DO i = i_start, i_end fqy(i, k, jp1) = 0.5*(fzm(k)*rv(i,k,j)+fzp(k)*rv(i,k-1,j))* & (w(i,k,j)+w(i,k,j-1)) ENDDO ENDDO k = ktf+1 DO i = i_start, i_end fqy(i, k, jp1) = 0.5*((2.-fzm(k-1))*rv(i,k-1,j)-fzp(k-1)*rv(i,k-2,j))* & (w(i,k,j)+w(i,k,j-1)) ENDDO ELSE IF ( j == jde-2 ) THEN ! 3rd or 4th order flux 2 in from north boundary DO k=kts+1,ktf DO i = i_start, i_end vel = fzm(k)*rv(i,k,j)+fzp(k)*rv(i,k-1,j) fqy( i, k, jp1 ) = vel*flux3( & w(i,k,j-2),w(i,k,j-1), & w(i,k,j),w(i,k,j+1),vel ) ENDDO ENDDO k = ktf+1 DO i = i_start, i_end vel = (2.-fzm(k-1))*rv(i,k-1,j)-fzp(k-1)*rv(i,k-2,j) fqy( i, k, jp1 ) = vel*flux3( & w(i,k,j-2),w(i,k,j-1), & w(i,k,j),w(i,k,j+1),vel ) ENDDO ENDIF ! y flux-divergence into tendency IF(j > j_start) THEN DO k=kts+1,ktf+1 DO i = i_start, i_end mrdy=msft(i,j-1)*rdy tendency(i,k,j-1) = tendency(i,k,j-1) - mrdy*(fqy(i,k,jp1)-fqy(i,k,jp0)) ENDDO ENDDO ENDIF jtmp = jp1 jp1 = jp0 jp0 = jtmp ENDDO j_loop_y_flux_5 ! next, x - flux divergence i_start = its i_end = MIN(ite,ide-1) j_start = jts j_end = MIN(jte,jde-1) ! higher order flux has a 5 or 7 point stencil, so compute ! bounds so we can switch to second order flux close to the boundary i_start_f = i_start i_end_f = i_end+1 IF(degrade_xs) then i_start = MAX(ids+1,its) i_start_f = i_start+2 ENDIF IF(degrade_xe) then i_end = MIN(ide-2,ite) i_end_f = ide-3 ENDIF ! compute fluxes DO j = j_start, j_end ! 5th or 6th order flux DO k=kts+1,ktf DO i = i_start_f, i_end_f vel = fzm(k)*ru(i,k,j)+fzp(k)*ru(i,k-1,j) fqx( i,k ) = vel*flux5( w(i-3,k,j), w(i-2,k,j), & w(i-1,k,j), w(i ,k,j), & w(i+1,k,j), w(i+2,k,j), & vel ) ENDDO ENDDO k = ktf+1 DO i = i_start_f, i_end_f vel = (2.-fzm(k-1))*ru(i,k-1,j)-fzp(k-1)*ru(i,k-2,j) fqx( i,k ) = vel*flux5( w(i-3,k,j), w(i-2,k,j), & w(i-1,k,j), w(i ,k,j), & w(i+1,k,j), w(i+2,k,j), & vel ) ENDDO ! lower order fluxes close to boundaries (if not periodic or symmetric) IF( degrade_xs ) THEN IF( i_start == ids+1 ) THEN ! second order flux next to the boundary i = ids+1 DO k=kts+1,ktf fqx(i,k) = 0.5*(fzm(k)*ru(i,k,j)+fzp(k)*ru(i,k-1,j)) & *(w(i,k,j)+w(i-1,k,j)) ENDDO k = ktf+1 fqx(i,k) = 0.5*((2.-fzm(k-1))*ru(i,k-1,j)-fzp(k-1)*ru(i,k-2,j)) & *(w(i,k,j)+w(i-1,k,j)) ENDIF i = i_start+1 DO k=kts+1,ktf vel = fzm(k)*ru(i,k,j)+fzp(k)*ru(i,k-1,j) fqx( i,k ) = vel*flux3( w(i-2,k,j), w(i-1,k,j), & w(i ,k,j), w(i+1,k,j), & vel ) ENDDO k = ktf+1 vel = (2.-fzm(k-1))*ru(i,k-1,j)-fzp(k-1)*ru(i,k-2,j) fqx( i,k ) = vel*flux3( w(i-2,k,j), w(i-1,k,j), & w(i ,k,j), w(i+1,k,j), & vel ) ENDIF IF( degrade_xe ) THEN IF( i_end == ide-2 ) THEN ! second order flux next to the boundary i = ide-1 DO k=kts+1,ktf fqx(i,k) = 0.5*(fzm(k)*ru(i,k,j)+fzp(k)*ru(i,k-1,j)) & *(w(i,k,j)+w(i-1,k,j)) ENDDO k = ktf+1 fqx(i,k) = 0.5*((2.-fzm(k-1))*ru(i,k-1,j)-fzp(k-1)*ru(i,k-2,j)) & *(w(i,k,j)+w(i-1,k,j)) ENDIF i = ide-2 DO k=kts+1,ktf vel = fzm(k)*ru(i,k,j)+fzp(k)*ru(i,k-1,j) fqx( i,k ) = vel*flux3( w(i-2,k,j), w(i-1,k,j), & w(i ,k,j), w(i+1,k,j), & vel ) ENDDO k = ktf+1 vel = (2.-fzm(k-1))*ru(i,k-1,j)-fzp(k-1)*ru(i,k-2,j) fqx( i,k ) = vel*flux3( w(i-2,k,j), w(i-1,k,j), & w(i ,k,j), w(i+1,k,j), & vel ) ENDIF ! x flux-divergence into tendency DO k=kts+1,ktf+1 DO i = i_start, i_end mrdx=msft(i,j)*rdx tendency(i,k,j) = tendency(i,k,j) - mrdx*(fqx(i+1,k)-fqx(i,k)) ENDDO ENDDO ENDDO ELSE IF ( horz_order == 4 ) THEN degrade_xs = .true. degrade_xe = .true. degrade_ys = .true. degrade_ye = .true. IF( config_flags%periodic_x .or. & config_flags%symmetric_xs .or. & (its > ids+1) ) degrade_xs = .false. IF( config_flags%periodic_x .or. & config_flags%symmetric_xe .or. & (ite < ide-2) ) degrade_xe = .false. IF( config_flags%periodic_y .or. & config_flags%symmetric_ys .or. & (jts > jds+1) ) degrade_ys = .false. IF( config_flags%periodic_y .or. & config_flags%symmetric_ye .or. & (jte < jde-2) ) degrade_ye = .false. ! begin flux computations ! start with x flux divergence !--------------- ktf=MIN(kte,kde-1) i_start = its i_end = MIN(ite,ide-1) j_start = jts j_end = MIN(jte,jde-1) ! 3rd or 4th order flux has a 5 point stencil, so compute ! bounds so we can switch to second order flux close to the boundary i_start_f = i_start i_end_f = i_end+1 IF(degrade_xs) then i_start = ids+1 i_start_f = i_start+1 ENDIF IF(degrade_xe) then i_end = ide-2 i_end_f = ide-2 ENDIF ! compute fluxes DO j = j_start, j_end DO k=kts+1,ktf DO i = i_start_f, i_end_f vel = fzm(k)*ru(i,k,j)+fzp(k)*ru(i,k-1,j) fqx( i, k ) = vel*flux4( w(i-2,k,j), w(i-1,k,j), & w(i ,k,j), w(i+1,k,j), & vel ) ENDDO ENDDO k = ktf+1 DO i = i_start_f, i_end_f vel = (2.-fzm(k-1))*ru(i,k-1,j)-fzp(k-1)*ru(i,k-2,j) fqx( i, k ) = vel*flux4( w(i-2,k,j), w(i-1,k,j), & w(i ,k,j), w(i+1,k,j), & vel ) ENDDO ! second order flux close to boundaries (if not periodic or symmetric) IF( degrade_xs ) THEN DO k=kts+1,ktf fqx(i_start, k) = & 0.5*(fzm(k)*ru(i_start,k,j)+fzp(k)*ru(i_start,k-1,j)) & *(w(i_start,k,j)+w(i_start-1,k,j)) ENDDO k = ktf+1 fqx(i_start, k) = & 0.5*((2.-fzm(k-1))*ru(i_start,k-1,j)-fzp(k-1)*ru(i_start,k-2,j)) & *(w(i_start,k,j)+w(i_start-1,k,j)) ENDIF IF( degrade_xe ) THEN DO k=kts+1,ktf fqx(i_end+1, k) = & 0.5*(fzm(k)*ru(i_end+1,k,j)+fzp(k)*ru(i_end+1,k-1,j)) & *(w(i_end+1,k,j)+w(i_end,k,j)) ENDDO k = ktf+1 fqx(i_end+1, k) = & 0.5*((2.-fzm(k-1))*ru(i_end+1,k-1,j)-fzp(k-1)*ru(i_end+1,k-2,j)) & *(w(i_end+1,k,j)+w(i_end,k,j)) ENDIF ! x flux-divergence into tendency DO k=kts+1,ktf+1 DO i = i_start, i_end mrdx=msft(i,j)*rdx tendency(i,k,j) = tendency(i,k,j) - mrdx*(fqx(i+1,k)-fqx(i,k)) ENDDO ENDDO ENDDO ! next -> y flux divergence calculation i_start = its i_end = MIN(ite,ide-1) j_start = jts j_end = MIN(jte,jde-1) ! 3rd or 4th order flux has a 5 point stencil, so compute ! bounds so we can switch to second order flux close to the boundary j_start_f = j_start j_end_f = j_end+1 IF(degrade_ys) then j_start = jds+1 j_start_f = j_start+1 ENDIF IF(degrade_ye) then j_end = jde-2 j_end_f = jde-2 ENDIF jp1 = 2 jp0 = 1 DO j = j_start, j_end+1 IF ((j < j_start_f) .and. degrade_ys) THEN DO k = kts+1, ktf DO i = i_start, i_end fqy(i, k, jp1) = & 0.5*(fzm(k)*rv(i,k,j_start)+fzp(k)*rv(i,k-1,j_start)) & *(w(i,k,j_start)+w(i,k,j_start-1)) ENDDO ENDDO k = ktf+1 DO i = i_start, i_end fqy(i, k, jp1) = & 0.5*((2.-fzm(k-1))*rv(i,k-1,j_start)-fzp(k-1)*rv(i,k-2,j_start)) & *(w(i,k,j_start)+w(i,k,j_start-1)) ENDDO ELSE IF ((j > j_end_f) .and. degrade_ye) THEN DO k = kts+1, ktf DO i = i_start, i_end fqy(i, k, jp1) = & 0.5*(fzm(k)*rv(i,k,j_end+1)+fzp(k)*rv(i,k-1,j_end+1)) & *(w(i,k,j_end+1)+w(i,k,j_end)) ENDDO ENDDO k = ktf+1 DO i = i_start, i_end fqy(i, k, jp1) = & 0.5*((2.-fzm(k-1))*rv(i,k-1,j_end+1)-fzp(k-1)*rv(i,k-2,j_end+1)) & *(w(i,k,j_end+1)+w(i,k,j_end)) ENDDO ELSE ! 3rd or 4th order flux DO k = kts+1, ktf DO i = i_start, i_end vel = fzm(k)*rv(i,k,j)+fzp(k)*rv(i,k-1,j) fqy( i, k, jp1 ) = vel*flux4( w(i,k,j-2), w(i,k,j-1), & w(i,k,j ), w(i,k,j+1), & vel ) ENDDO ENDDO k = ktf+1 DO i = i_start, i_end vel = (2.-fzm(k-1))*rv(i,k-1,j)-fzp(k-1)*rv(i,k-2,j) fqy( i, k, jp1 ) = vel*flux4( w(i,k,j-2), w(i,k,j-1), & w(i,k,j ), w(i,k,j+1), & vel ) ENDDO END IF IF( j > j_start ) THEN ! y flux-divergence into tendency DO k = kts+1, ktf+1 DO i = i_start, i_end mrdy=msft(i,j-1)*rdy tendency(i,k,j-1) = tendency(i,k,j-1) - mrdy*(fqy(i,k,jp1)-fqy(i,k,jp0)) ENDDO ENDDO END IF jtmp = jp1 jp1 = jp0 jp0 = jtmp ENDDO ELSE IF ( horz_order == 3 ) THEN degrade_xs = .true. degrade_xe = .true. degrade_ys = .true. degrade_ye = .true. IF( config_flags%periodic_x .or. & config_flags%symmetric_xs .or. & (its > ids+1) ) degrade_xs = .false. IF( config_flags%periodic_x .or. & config_flags%symmetric_xe .or. & (ite < ide-2) ) degrade_xe = .false. IF( config_flags%periodic_y .or. & config_flags%symmetric_ys .or. & (jts > jds+1) ) degrade_ys = .false. IF( config_flags%periodic_y .or. & config_flags%symmetric_ye .or. & (jte < jde-2) ) degrade_ye = .false. ! begin flux computations ! start with x flux divergence !--------------- ktf=MIN(kte,kde-1) i_start = its i_end = MIN(ite,ide-1) j_start = jts j_end = MIN(jte,jde-1) ! 3rd or 4th order flux has a 5 point stencil, so compute ! bounds so we can switch to second order flux close to the boundary i_start_f = i_start i_end_f = i_end+1 IF(degrade_xs) then i_start = ids+1 i_start_f = i_start+1 ENDIF IF(degrade_xe) then i_end = ide-2 i_end_f = ide-2 ENDIF ! compute fluxes DO j = j_start, j_end DO k=kts+1,ktf DO i = i_start_f, i_end_f vel = fzm(k)*ru(i,k,j)+fzp(k)*ru(i,k-1,j) fqx( i, k ) = vel*flux3( w(i-2,k,j), w(i-1,k,j), & w(i ,k,j), w(i+1,k,j), & vel ) ENDDO ENDDO k = ktf+1 DO i = i_start_f, i_end_f vel = (2.-fzm(k-1))*ru(i,k-1,j)-fzp(k-1)*ru(i,k-2,j) fqx( i, k ) = vel*flux3( w(i-2,k,j), w(i-1,k,j), & w(i ,k,j), w(i+1,k,j), & vel ) ENDDO ! second order flux close to boundaries (if not periodic or symmetric) IF( degrade_xs ) THEN DO k=kts+1,ktf fqx(i_start, k) = & 0.5*(fzm(k)*ru(i_start,k,j)+fzp(k)*ru(i_start,k-1,j)) & *(w(i_start,k,j)+w(i_start-1,k,j)) ENDDO k = ktf+1 fqx(i_start, k) = & 0.5*((2.-fzm(k-1))*ru(i_start,k-1,j)-fzp(k-1)*ru(i_start,k-2,j)) & *(w(i_start,k,j)+w(i_start-1,k,j)) ENDIF IF( degrade_xe ) THEN DO k=kts+1,ktf fqx(i_end+1, k) = & 0.5*(fzm(k)*ru(i_end+1,k,j)+fzp(k)*ru(i_end+1,k-1,j)) & *(w(i_end+1,k,j)+w(i_end,k,j)) ENDDO k = ktf+1 fqx(i_end+1, k) = & 0.5*((2.-fzm(k-1))*ru(i_end+1,k-1,j)-fzp(k-1)*ru(i_end+1,k-2,j)) & *(w(i_end+1,k,j)+w(i_end,k,j)) ENDIF ! x flux-divergence into tendency DO k=kts+1,ktf+1 DO i = i_start, i_end mrdx=msft(i,j)*rdx tendency(i,k,j) = tendency(i,k,j) - mrdx*(fqx(i+1,k)-fqx(i,k)) ENDDO ENDDO ENDDO ! next -> y flux divergence calculation i_start = its i_end = MIN(ite,ide-1) j_start = jts j_end = MIN(jte,jde-1) ! 3rd or 4th order flux has a 5 point stencil, so compute ! bounds so we can switch to second order flux close to the boundary j_start_f = j_start j_end_f = j_end+1 IF(degrade_ys) then j_start = jds+1 j_start_f = j_start+1 ENDIF IF(degrade_ye) then j_end = jde-2 j_end_f = jde-2 ENDIF jp1 = 2 jp0 = 1 DO j = j_start, j_end+1 IF ((j < j_start_f) .and. degrade_ys) THEN DO k = kts+1, ktf DO i = i_start, i_end fqy(i, k, jp1) = & 0.5*(fzm(k)*rv(i,k,j_start)+fzp(k)*rv(i,k-1,j_start)) & *(w(i,k,j_start)+w(i,k,j_start-1)) ENDDO ENDDO k = ktf+1 DO i = i_start, i_end fqy(i, k, jp1) = & 0.5*((2.-fzm(k-1))*rv(i,k-1,j_start)-fzp(k-1)*rv(i,k-2,j_start)) & *(w(i,k,j_start)+w(i,k,j_start-1)) ENDDO ELSE IF ((j > j_end_f) .and. degrade_ye) THEN DO k = kts+1, ktf DO i = i_start, i_end fqy(i, k, jp1) = & 0.5*(fzm(k)*rv(i,k,j_end+1)+fzp(k)*rv(i,k-1,j_end+1)) & *(w(i,k,j_end+1)+w(i,k,j_end)) ENDDO ENDDO k = ktf+1 DO i = i_start, i_end fqy(i, k, jp1) = & 0.5*((2.-fzm(k-1))*rv(i,k-1,j_end+1)-fzp(k-1)*rv(i,k-2,j_end+1)) & *(w(i,k,j_end+1)+w(i,k,j_end)) ENDDO ELSE ! 3rd or 4th order flux DO k = kts+1, ktf DO i = i_start, i_end vel = fzm(k)*rv(i,k,j)+fzp(k)*rv(i,k-1,j) fqy( i, k, jp1 ) = vel*flux3( w(i,k,j-2), w(i,k,j-1), & w(i,k,j ), w(i,k,j+1), & vel ) ENDDO ENDDO k = ktf+1 DO i = i_start, i_end vel = (2.-fzm(k-1))*rv(i,k-1,j)-fzp(k-1)*rv(i,k-2,j) fqy( i, k, jp1 ) = vel*flux3( w(i,k,j-2), w(i,k,j-1), & w(i,k,j ), w(i,k,j+1), & vel ) ENDDO END IF IF( j > j_start ) THEN ! y flux-divergence into tendency DO k = kts+1, ktf+1 DO i = i_start, i_end mrdy=msft(i,j-1)*rdy tendency(i,k,j-1) = tendency(i,k,j-1) - mrdy*(fqy(i,k,jp1)-fqy(i,k,jp0)) ENDDO ENDDO END IF jtmp = jp1 jp1 = jp0 jp0 = jtmp ENDDO ELSE IF (horz_order == 2 ) THEN i_start = its i_end = MIN(ite,ide-1) j_start = jts j_end = MIN(jte,jde-1) IF ( .NOT. config_flags%periodic_x ) THEN IF ( config_flags%open_xs .or. specified ) i_start = MAX(ids+1,its) IF ( config_flags%open_xe .or. specified ) i_end = MIN(ide-2,ite) ENDIF DO j = j_start, j_end DO k=kts+1,ktf DO i = i_start, i_end mrdx=msft(i,j)*rdx tendency(i,k,j)=tendency(i,k,j)-mrdx*0.5 & *((fzm(k)*ru(i+1,k,j)+fzp(k)*ru(i+1,k-1,j)) & *(w(i+1,k,j)+w(i,k,j)) & -(fzm(k)*ru(i,k,j)+fzp(k)*ru(i,k-1,j)) & *(w(i,k,j)+w(i-1,k,j))) ENDDO ENDDO k = ktf+1 DO i = i_start, i_end mrdx=msft(i,j)*rdx tendency(i,k,j)=tendency(i,k,j)-mrdx*0.5 & *(((2.-fzm(k-1))*ru(i+1,k-1,j)-fzp(k-1)*ru(i+1,k-2,j)) & *(w(i+1,k,j)+w(i,k,j)) & -((2.-fzm(k-1))*ru(i,k-1,j)-fzp(k-1)*ru(i,k-2,j)) & *(w(i,k,j)+w(i-1,k,j))) ENDDO ENDDO i_start = its i_end = MIN(ite,ide-1) IF ( config_flags%open_ys .or. specified ) j_start = MAX(jds+1,jts) IF ( config_flags%open_ye .or. specified ) j_end = MIN(jde-2,jte) DO j = j_start, j_end DO k=kts+1,ktf DO i = i_start, i_end mrdy=msft(i,j)*rdy tendency(i,k,j)=tendency(i,k,j) -mrdy*0.5 & *((fzm(k)*rv(i,k,j+1)+fzp(k)*rv(i,k-1,j+1))* & (w(i,k,j+1)+w(i,k,j)) & -(fzm(k)*rv(i,k,j)+fzp(k)*rv(i,k-1,j)) & *(w(i,k,j)+w(i,k,j-1))) ENDDO ENDDO k = ktf+1 DO i = i_start, i_end mrdy=msft(i,j)*rdy tendency(i,k,j)=tendency(i,k,j) -mrdy*0.5 & *(((2.-fzm(k-1))*rv(i,k-1,j+1)-fzp(k-1)*rv(i,k-2,j+1))* & (w(i,k,j+1)+w(i,k,j)) & -((2.-fzm(k-1))*rv(i,k-1,j)-fzp(k-1)*rv(i,k-2,j)) & *(w(i,k,j)+w(i,k,j-1))) ENDDO ENDDO ELSE WRITE ( wrf_err_message ,*) ' advect_w_6a, h_order not known ',horz_order CALL wrf_error_fatal ( wrf_err_message ) ENDIF horizontal_order_test ! pick up the the horizontal radiation boundary conditions. ! (these are the computations that don't require 'cb'. ! first, set to index ranges i_start = its i_end = MIN(ite,ide-1) j_start = jts j_end = MIN(jte,jde-1) IF( (config_flags%open_xs) .and. (its == ids)) THEN DO j = j_start, j_end DO k = kts+1, ktf uw = 0.5*(fzm(k)*(ru(its,k ,j)+ru(its+1,k ,j)) + & fzp(k)*(ru(its,k-1,j)+ru(its+1,k-1,j)) ) ub = MIN( uw, 0. ) tendency(its,k,j) = tendency(its,k,j) & - rdx*( & ub*(w_old(its+1,k,j) - w_old(its,k,j)) + & w(its,k,j)*( & fzm(k)*(ru(its+1,k ,j)-ru(its,k ,j))+ & fzp(k)*(ru(its+1,k-1,j)-ru(its,k-1,j))) & ) ENDDO ENDDO k = ktf+1 DO j = j_start, j_end uw = 0.5*( (2.-fzm(k-1))*(ru(its,k-1,j)+ru(its+1,k-1,j)) & -fzp(k-1)*(ru(its,k-2,j)+ru(its+1,k-2,j)) ) ub = MIN( uw, 0. ) tendency(its,k,j) = tendency(its,k,j) & - rdx*( & ub*(w_old(its+1,k,j) - w_old(its,k,j)) + & w(its,k,j)*( & (2.-fzm(k-1))*(ru(its+1,k-1,j)-ru(its,k-1,j))- & fzp(k-1)*(ru(its+1,k-2,j)-ru(its,k-2,j))) & ) ENDDO ENDIF IF( (config_flags%open_xe) .and. (ite == ide)) THEN DO j = j_start, j_end DO k = kts+1, ktf uw = 0.5*(fzm(k)*(ru(ite-1,k ,j)+ru(ite,k ,j)) + & fzp(k)*(ru(ite-1,k-1,j)+ru(ite,k-1,j)) ) ub = MAX( uw, 0. ) tendency(i_end,k,j) = tendency(i_end,k,j) & - rdx*( & ub*(w_old(i_end,k,j) - w_old(i_end-1,k,j)) + & w(i_end,k,j)*( & fzm(k)*(ru(ite,k ,j)-ru(ite-1,k ,j)) + & fzp(k)*(ru(ite,k-1,j)-ru(ite-1,k-1,j))) & ) ENDDO ENDDO k = ktf+1 DO j = j_start, j_end uw = 0.5*( (2.-fzm(k-1))*(ru(ite-1,k-1,j)+ru(ite,k-1,j)) & -fzp(k-1)*(ru(ite-1,k-2,j)+ru(ite,k-2,j)) ) ub = MAX( uw, 0. ) tendency(i_end,k,j) = tendency(i_end,k,j) & - rdx*( & ub*(w_old(i_end,k,j) - w_old(i_end-1,k,j)) + & w(i_end,k,j)*( & (2.-fzm(k-1))*(ru(ite,k-1,j)-ru(ite-1,k-1,j)) - & fzp(k-1)*(ru(ite,k-2,j)-ru(ite-1,k-2,j))) & ) ENDDO ENDIF IF( (config_flags%open_ys) .and. (jts == jds)) THEN DO i = i_start, i_end DO k = kts+1, ktf vw = 0.5*( fzm(k)*(rv(i,k ,jts)+rv(i,k ,jts+1)) + & fzp(k)*(rv(i,k-1,jts)+rv(i,k-1,jts+1)) ) vb = MIN( vw, 0. ) tendency(i,k,jts) = tendency(i,k,jts) & - rdy*( & vb*(w_old(i,k,jts+1) - w_old(i,k,jts)) + & w(i,k,jts)*( & fzm(k)*(rv(i,k ,jts+1)-rv(i,k ,jts))+ & fzp(k)*(rv(i,k-1,jts+1)-rv(i,k-1,jts))) & ) ENDDO ENDDO k = ktf+1 DO i = i_start, i_end vw = 0.5*( (2.-fzm(k-1))*(rv(i,k-1,jts)+rv(i,k-1,jts+1)) & -fzp(k-1)*(rv(i,k-2,jts)+rv(i,k-2,jts+1)) ) vb = MIN( vw, 0. ) tendency(i,k,jts) = tendency(i,k,jts) & - rdy*( & vb*(w_old(i,k,jts+1) - w_old(i,k,jts)) + & w(i,k,jts)*( & (2.-fzm(k-1))*(rv(i,k-1,jts+1)-rv(i,k-1,jts))- & fzp(k-1)*(rv(i,k-2,jts+1)-rv(i,k-2,jts))) & ) ENDDO ENDIF IF( (config_flags%open_ye) .and. (jte == jde) ) THEN DO i = i_start, i_end DO k = kts+1, ktf vw = 0.5*( fzm(k)*(rv(i,k ,jte-1)+rv(i,k ,jte)) + & fzp(k)*(rv(i,k-1,jte-1)+rv(i,k-1,jte)) ) vb = MAX( vw, 0. ) tendency(i,k,j_end) = tendency(i,k,j_end) & - rdy*( & vb*(w_old(i,k,j_end) - w_old(i,k,j_end-1)) + & w(i,k,j_end)*( & fzm(k)*(rv(i,k ,jte)-rv(i,k ,jte-1))+ & fzp(k)*(rv(i,k-1,jte)-rv(i,k-1,jte-1))) & ) ENDDO ENDDO k = ktf+1 DO i = i_start, i_end vw = 0.5*( (2.-fzm(k-1))*(rv(i,k-1,jte-1)+rv(i,k-1,jte)) & -fzp(k-1)*(rv(i,k-2,jte-1)+rv(i,k-2,jte)) ) vb = MAX( vw, 0. ) tendency(i,k,j_end) = tendency(i,k,j_end) & - rdy*( & vb*(w_old(i,k,j_end) - w_old(i,k,j_end-1)) + & w(i,k,j_end)*( & (2.-fzm(k-1))*(rv(i,k-1,jte)-rv(i,k-1,jte-1))- & fzp(k-1)*(rv(i,k-2,jte)-rv(i,k-2,jte-1))) & ) ENDDO ENDIF !-------------------- vertical advection i_start = its i_end = MIN(ite,ide-1) j_start = jts j_end = MIN(jte,jde-1) DO i = i_start, i_end vflux(i,kts)=0. vflux(i,kte)=0. ENDDO vert_order_test : IF (vert_order == 6) THEN DO j = j_start, j_end DO k=kts+3,ktf-1 DO i = i_start, i_end vel=0.5*(rom(i,k,j)+rom(i,k-1,j)) vflux(i,k) = vel*flux6( & w(i,k-3,j), w(i,k-2,j), w(i,k-1,j), & w(i,k ,j), w(i,k+1,j), w(i,k+2,j), -vel ) ENDDO ENDDO DO i = i_start, i_end k=kts+1 vflux(i,k)=0.25*(rom(i,k,j)+rom(i,k-1,j))*(w(i,k,j)+w(i,k-1,j)) k = kts+2 vel=0.5*(rom(i,k,j)+rom(i,k-1,j)) vflux(i,k) = vel*flux4( & w(i,k-2,j), w(i,k-1,j), & w(i,k ,j), w(i,k+1,j), -vel ) k = ktf vel=0.5*(rom(i,k,j)+rom(i,k-1,j)) vflux(i,k) = vel*flux4( & w(i,k-2,j), w(i,k-1,j), & w(i,k ,j), w(i,k+1,j), -vel ) k=ktf+1 vflux(i,k)=0.25*(rom(i,k,j)+rom(i,k-1,j))*(w(i,k,j)+w(i,k-1,j)) ENDDO DO k=kts+1,ktf DO i = i_start, i_end tendency(i,k,j)=tendency(i,k,j)-rdzu(k)*(vflux(i,k+1)-vflux(i,k)) ENDDO ENDDO ! pick up flux contribution for w at the lid. wcs, 13 march 2004 k = ktf+1 DO i = i_start, i_end tendency(i,k,j)=tendency(i,k,j)+2.*rdzu(k-1)*(vflux(i,k)) ENDDO ENDDO ELSE IF (vert_order == 5) THEN DO j = j_start, j_end DO k=kts+3,ktf-1 DO i = i_start, i_end vel=0.5*(rom(i,k,j)+rom(i,k-1,j)) vflux(i,k) = vel*flux5( & w(i,k-3,j), w(i,k-2,j), w(i,k-1,j), & w(i,k ,j), w(i,k+1,j), w(i,k+2,j), -vel ) ENDDO ENDDO DO i = i_start, i_end k=kts+1 vflux(i,k)=0.25*(rom(i,k,j)+rom(i,k-1,j))*(w(i,k,j)+w(i,k-1,j)) k = kts+2 vel=0.5*(rom(i,k,j)+rom(i,k-1,j)) vflux(i,k) = vel*flux3( & w(i,k-2,j), w(i,k-1,j), & w(i,k ,j), w(i,k+1,j), -vel ) k = ktf vel=0.5*(rom(i,k,j)+rom(i,k-1,j)) vflux(i,k) = vel*flux3( & w(i,k-2,j), w(i,k-1,j), & w(i,k ,j), w(i,k+1,j), -vel ) k=ktf+1 vflux(i,k)=0.25*(rom(i,k,j)+rom(i,k-1,j))*(w(i,k,j)+w(i,k-1,j)) ENDDO DO k=kts+1,ktf DO i = i_start, i_end tendency(i,k,j)=tendency(i,k,j)-rdzu(k)*(vflux(i,k+1)-vflux(i,k)) ENDDO ENDDO ! pick up flux contribution for w at the lid, wcs. 13 march 2004 k = ktf+1 DO i = i_start, i_end tendency(i,k,j)=tendency(i,k,j)+2.*rdzu(k-1)*(vflux(i,k)) ENDDO ENDDO ELSE IF (vert_order == 4) THEN DO j = j_start, j_end DO k=kts+2,ktf DO i = i_start, i_end vel=0.5*(rom(i,k,j)+rom(i,k-1,j)) vflux(i,k) = vel*flux4( & w(i,k-2,j), w(i,k-1,j), & w(i,k ,j), w(i,k+1,j), -vel ) ENDDO ENDDO DO i = i_start, i_end k=kts+1 vflux(i,k)=0.25*(rom(i,k,j)+rom(i,k-1,j))*(w(i,k,j)+w(i,k-1,j)) k=ktf+1 vflux(i,k)=0.25*(rom(i,k,j)+rom(i,k-1,j))*(w(i,k,j)+w(i,k-1,j)) ENDDO DO k=kts+1,ktf DO i = i_start, i_end tendency(i,k,j)=tendency(i,k,j)-rdzu(k)*(vflux(i,k+1)-vflux(i,k)) ENDDO ENDDO ! pick up flux contribution for w at the lid, wcs. 13 march 2004 k = ktf+1 DO i = i_start, i_end tendency(i,k,j)=tendency(i,k,j)+2.*rdzu(k-1)*(vflux(i,k)) ENDDO ENDDO ELSE IF (vert_order == 3) THEN DO j = j_start, j_end DO k=kts+2,ktf DO i = i_start, i_end vel=0.5*(rom(i,k,j)+rom(i,k-1,j)) vflux(i,k) = vel*flux3( & w(i,k-2,j), w(i,k-1,j), & w(i,k ,j), w(i,k+1,j), -vel ) ENDDO ENDDO DO i = i_start, i_end k=kts+1 vflux(i,k)=0.25*(rom(i,k,j)+rom(i,k-1,j))*(w(i,k,j)+w(i,k-1,j)) k=ktf+1 vflux(i,k)=0.25*(rom(i,k,j)+rom(i,k-1,j))*(w(i,k,j)+w(i,k-1,j)) ENDDO DO k=kts+1,ktf DO i = i_start, i_end tendency(i,k,j)=tendency(i,k,j)-rdzu(k)*(vflux(i,k+1)-vflux(i,k)) ENDDO ENDDO ! pick up flux contribution for w at the lid, wcs. 13 march 2004 k = ktf+1 DO i = i_start, i_end tendency(i,k,j)=tendency(i,k,j)+2.*rdzu(k-1)*(vflux(i,k)) ENDDO ENDDO ELSE IF (vert_order == 2) THEN DO j = j_start, j_end DO k=kts+1,ktf+1 DO i = i_start, i_end vflux(i,k)=0.25*(rom(i,k,j)+rom(i,k-1,j))*(w(i,k,j)+w(i,k-1,j)) ENDDO ENDDO DO k=kts+1,ktf DO i = i_start, i_end tendency(i,k,j)=tendency(i,k,j)-rdzu(k)*(vflux(i,k+1)-vflux(i,k)) ENDDO ENDDO ! pick up flux contribution for w at the lid, wcs. 13 march 2004 k = ktf+1 DO i = i_start, i_end tendency(i,k,j)=tendency(i,k,j)+2.*rdzu(k-1)*(vflux(i,k)) ENDDO ENDDO ELSE WRITE (wrf_err_message ,*) ' advect_w, v_order not known ',vert_order CALL wrf_error_fatal ( wrf_err_message ) ENDIF vert_order_test END SUBROUTINE advect_w !---------------------------------------------------------------- SUBROUTINE advect_scalar_pd ( field, field_old, tendency, & ru, rv, rom, & mut, mub, mu_old, & config_flags, & msfu, msfv, msft, & fzm, fzp, & rdx, rdy, rdzw, dt, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) ! this is a first cut at a positive definite advection option ! for scalars in WRF. This version is memory intensive -> ! we save 3d arrays of x, y and z both high and low order fluxes ! (six in all). Alternatively, we could sweep in a direction ! and lower the cost considerably. ! uses the Smolarkiewicz MWR 1989 approach, with addition of first-order ! fluxes initially ! WCS, 3 December 2002, 24 February 2003 IMPLICIT NONE ! Input data TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: field, & field_old, & ru, & rv, & rom REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mut, mub, mu_old REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: msfu, & msfv, & msft REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fzm, & fzp, & rdzw REAL , INTENT(IN ) :: rdx, & rdy, & dt ! Local data INTEGER :: i, j, k, itf, jtf, ktf INTEGER :: i_start, i_end, j_start, j_end INTEGER :: i_start_f, i_end_f, j_start_f, j_end_f INTEGER :: jmin, jmax, jp, jm, imin, imax REAL :: mrdx, mrdy, ub, vb, uw, vw, mu ! storage for high and low order fluxes REAL, DIMENSION( its-1:ite+2, kts:kte, jts-1:jte+2 ) :: fqx, fqy, fqz REAL, DIMENSION( its-1:ite+2, kts:kte, jts-1:jte+2 ) :: fqxl, fqyl, fqzl INTEGER :: horz_order, vert_order LOGICAL :: degrade_xs, degrade_ys LOGICAL :: degrade_xe, degrade_ye INTEGER :: jp1, jp0, jtmp REAL :: flux_out, ph_low, scale REAL, PARAMETER :: eps=1.e-20 ! definition of flux operators, 3rd, 4th, 5th or 6th order REAL :: flux3, flux4, flux5, flux6, flux_upwind REAL :: q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua, vel, cr flux4(q_im2, q_im1, q_i, q_ip1, ua) = & (7./12.)*(q_i + q_im1) - (1./12.)*(q_ip1 + q_im2) flux3(q_im2, q_im1, q_i, q_ip1, ua) = & flux4(q_im2, q_im1, q_i, q_ip1, ua) + & sign(1.,ua)*(1./12.)*((q_ip1 - q_im2)-3.*(q_i-q_im1)) flux6(q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua) = & (37./60.)*(q_i+q_im1) - (2./15.)*(q_ip1+q_im2) & +(1./60.)*(q_ip2+q_im3) flux5(q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua) = & flux6(q_im3, q_im2, q_im1, q_i, q_ip1, q_ip2, ua) & -sign(1.,ua)*(1./60.)*( & (q_ip2-q_im3)-5.*(q_ip1-q_im2)+10.*(q_i-q_im1) ) flux_upwind(q_im1, q_i, cr ) = 0.5*min( 1.0,(cr+abs(cr)))*q_im1 & +0.5*max(-1.0,(cr-abs(cr)))*q_i ! flux_upwind(q_im1, q_i, cr ) = 0. REAL :: dx,dy,dz LOGICAL, PARAMETER :: pd_limit = .true. ! set order for the advection schemes ! write(6,*) ' in pd advection routine ' ktf=MIN(kte,kde-1) horz_order = config_flags%h_sca_adv_order vert_order = config_flags%v_sca_adv_order ! determine boundary mods for flux operators ! We degrade the flux operators from 3rd/4th order ! to second order one gridpoint in from the boundaries for ! all boundary conditions except periodic and symmetry - these ! conditions have boundary zone data fill for correct application ! of the higher order flux stencils degrade_xs = .true. degrade_xe = .true. degrade_ys = .true. degrade_ye = .true. ! begin with horizontal flux divergence ! here is the choice of flux operators horizontal_order_test : IF( horz_order == 6 ) THEN IF( config_flags%periodic_x .or. & config_flags%symmetric_xs .or. & (its > ids+2) ) degrade_xs = .false. IF( config_flags%periodic_x .or. & config_flags%symmetric_xe .or. & (ite < ide-3) ) degrade_xe = .false. IF( config_flags%periodic_y .or. & config_flags%symmetric_ys .or. & (jts > jds+2) ) degrade_ys = .false. IF( config_flags%periodic_y .or. & config_flags%symmetric_ye .or. & (jte < jde-3) ) degrade_ye = .false. !--------------- y - advection first !-- y flux compute; these bounds are for periodic and sym b.c. ktf=MIN(kte,kde-1) i_start = its-1 i_end = MIN(ite,ide-1)+1 j_start = jts-1 j_end = MIN(jte,jde-1)+1 j_start_f = j_start j_end_f = j_end+1 !-- modify loop bounds if open or specified IF(degrade_xs) i_start = its IF(degrade_xe) i_end = MIN(ite,ide-1) IF(degrade_ys) then j_start = MAX(jts,jds+1) j_start_f = jds+3 ENDIF IF(degrade_ye) then j_end = MIN(jte,jde-2) j_end_f = jde-3 ENDIF ! compute fluxes, 6th order j_loop_y_flux_6 : DO j = j_start, j_end+1 IF( (j >= j_start_f ) .and. (j <= j_end_f) ) THEN ! use full stencil DO k=kts,ktf DO i = i_start, i_end dy = 2./(msft(i,j)+msft(i,j-1))/rdy mu = 0.5*(mut(i,j)+mut(i,j-1)) vel = rv(i,k,j) cr = vel*dt/dy/mu fqyl(i,k,j) = mu*(dy/dt)*flux_upwind(field_old(i,k,j-1), field_old(i,k,j ), cr) fqy( i, k, j ) = vel*flux6( & field(i,k,j-3), field(i,k,j-2), field(i,k,j-1), & field(i,k,j ), field(i,k,j+1), field(i,k,j+2), vel ) fqy(i,k,j) = fqy(i,k,j) - fqyl(i,k,j) ENDDO ENDDO ELSE IF ( j == jds+1 ) THEN ! 2nd order flux next to south boundary DO k=kts,ktf DO i = i_start, i_end dy = 2./(msft(i,j)+msft(i,j-1))/rdy mu = 0.5*(mut(i,j)+mut(i,j-1)) vel = rv(i,k,j) cr = vel*dt/dy/mu fqyl(i,k,j) = mu*(dy/dt)*flux_upwind(field_old(i,k,j-1), field_old(i,k,j ), cr) fqy(i,k, j) = 0.5*rv(i,k,j)* & (field(i,k,j)+field(i,k,j-1)) fqy(i,k,j) = fqy(i,k,j) - fqyl(i,k,j) ENDDO ENDDO ELSE IF ( j == jds+2 ) THEN ! third of 4th order flux 2 in from south boundary DO k=kts,ktf DO i = i_start, i_end dy = 2./(msft(i,j)+msft(i,j-1))/rdy mu = 0.5*(mut(i,j)+mut(i,j-1)) vel = rv(i,k,j) cr = vel*dt/dy/mu fqyl(i,k,j) = mu*(dy/dt)*flux_upwind(field_old(i,k,j-1), field_old(i,k,j ), cr) fqy( i, k, j ) = vel*flux4( & field(i,k,j-2),field(i,k,j-1),field(i,k,j),field(i,k,j+1),vel ) fqy(i,k,j) = fqy(i,k,j) - fqyl(i,k,j) ENDDO ENDDO ELSE IF ( j == jde-1 ) THEN ! 2nd order flux next to north boundary DO k=kts,ktf DO i = i_start, i_end dy = 2./(msft(i,j)+msft(i,j-1))/rdy mu = 0.5*(mut(i,j)+mut(i,j-1)) vel = rv(i,k,j) cr = vel*dt/dy/mu fqyl(i,k,j) = mu*(dy/dt)*flux_upwind(field_old(i,k,j-1), field_old(i,k,j ), cr) fqy(i, k, j ) = 0.5*rv(i,k,j)* & (field(i,k,j)+field(i,k,j-1)) fqy(i,k,j) = fqy(i,k,j) - fqyl(i,k,j) ENDDO ENDDO ELSE IF ( j == jde-2 ) THEN ! 4th order flux 2 in from north boundary DO k=kts,ktf DO i = i_start, i_end dy = 2./(msft(i,j)+msft(i,j-1))/rdy mu = 0.5*(mut(i,j)+mut(i,j-1)) vel = rv(i,k,j) cr = vel*dt/dy/mu fqyl(i,k,j) = mu*(dy/dt)*flux_upwind(field_old(i,k,j-1), field_old(i,k,j ), cr) fqy( i, k, j) = vel*flux4( & field(i,k,j-2),field(i,k,j-1), & field(i,k,j),field(i,k,j+1),vel ) fqy(i,k,j) = fqy(i,k,j) - fqyl(i,k,j) ENDDO ENDDO ENDIF ENDDO j_loop_y_flux_6 ! next, x flux !-- these bounds are for periodic and sym conditions i_start = its-1 i_end = MIN(ite,ide-1)+1 i_start_f = i_start i_end_f = i_end+1 j_start = jts-1 j_end = MIN(jte,jde-1)+1 !-- modify loop bounds for open and specified b.c IF(degrade_ys) j_start = jts IF(degrade_ye) j_end = MIN(jte,jde-1) IF(degrade_xs) then i_start = MAX(ids+1,its) i_start_f = i_start+2 ENDIF IF(degrade_xe) then i_end = MIN(ide-2,ite) i_end_f = ide-3 ENDIF ! compute fluxes DO j = j_start, j_end ! 6th order flux DO k=kts,ktf DO i = i_start_f, i_end_f dx = 2./(msft(i,j)+msft(i-1,j))/rdx mu = 0.5*(mut(i,j)+mut(i-1,j)) vel = ru(i,k,j) cr = vel*dt/dx/mu fqxl(i,k,j) = mu*(dx/dt)*flux_upwind(field_old(i-1,k,j), field_old(i,k,j ), cr) fqx( i,k,j ) = vel*flux6( field(i-3,k,j), field(i-2,k,j), & field(i-1,k,j), field(i ,k,j), & field(i+1,k,j), field(i+2,k,j), & vel ) fqx(i,k,j) = fqx(i,k,j) - fqxl(i,k,j) ENDDO ENDDO ! lower order fluxes close to boundaries (if not periodic or symmetric) IF( degrade_xs ) THEN IF( i_start == ids+1 ) THEN ! second order flux next to the boundary i = ids+1 DO k=kts,ktf dx = 2./(msft(i,j)+msft(i-1,j))/rdx mu = 0.5*(mut(i,j)+mut(i-1,j)) vel = ru(i,k,j)/mu cr = vel*dt/dx fqxl(i,k,j) = mu*(dx/dt)*flux_upwind(field_old(i-1,k,j), field_old(i,k,j ), cr) fqx(i,k,j) = 0.5*(ru(i,k,j)) & *(field(i,k,j)+field(i-1,k,j)) fqx(i,k,j) = fqx(i,k,j) - fqxl(i,k,j) ENDDO ENDIF i = ids+2 DO k=kts,ktf dx = 2./(msft(i,j)+msft(i-1,j))/rdx mu = 0.5*(mut(i,j)+mut(i-1,j)) vel = ru(i,k,j) cr = vel*dt/dx/mu fqxl(i,k,j) = mu*(dx/dt)*flux_upwind(field_old(i-1,k,j), field_old(i,k,j ), cr) fqx( i,k,j ) = vel*flux4( field(i-2,k,j), field(i-1,k,j), & field(i ,k,j), field(i+1,k,j), & vel ) fqx(i,k,j) = fqx(i,k,j) - fqxl(i,k,j) ENDDO ENDIF IF( degrade_xe ) THEN IF( i_end == ide-2 ) THEN ! second order flux next to the boundary i = ide-1 DO k=kts,ktf dx = 2./(msft(i,j)+msft(i-1,j))/rdx mu = 0.5*(mut(i,j)+mut(i-1,j)) vel = ru(i,k,j) cr = vel*dt/dx/mu fqxl(i,k,j) = mu*(dx/dt)*flux_upwind(field_old(i-1,k,j), field_old(i,k,j ), cr) fqx(i,k,j) = 0.5*(ru(i,k,j)) & *(field(i,k,j)+field(i-1,k,j)) fqx(i,k,j) = fqx(i,k,j) - fqxl(i,k,j) ENDDO ENDIF i = ide-2 DO k=kts,ktf dx = 2./(msft(i,j)+msft(i-1,j))/rdx mu = 0.5*(mut(i,j)+mut(i-1,j)) vel = ru(i,k,j) cr = vel*dt/dx/mu fqxl(i,k,j) = mu*(dx/dt)*flux_upwind(field_old(i-1,k,j), field_old(i,k,j ), cr) fqx( i,k,j ) = vel*flux4( field(i-2,k,j), field(i-1,k,j), & field(i ,k,j), field(i+1,k,j), & vel ) fqx(i,k,j) = fqx(i,k,j) - fqxl(i,k,j) ENDDO ENDIF ENDDO ! enddo for outer J loop !--- end of 6th order horizontal flux calculation ELSE IF( horz_order == 5 ) THEN IF( config_flags%periodic_x .or. & config_flags%symmetric_xs .or. & (its > ids+2) ) degrade_xs = .false. IF( config_flags%periodic_x .or. & config_flags%symmetric_xe .or. & (ite < ide-3) ) degrade_xe = .false. IF( config_flags%periodic_y .or. & config_flags%symmetric_ys .or. & (jts > jds+2) ) degrade_ys = .false. IF( config_flags%periodic_y .or. & config_flags%symmetric_ye .or. & (jte < jde-3) ) degrade_ye = .false. !--------------- y - advection first !-- y flux compute; these bounds are for periodic and sym b.c. ktf=MIN(kte,kde-1) i_start = its-1 i_end = MIN(ite,ide-1)+1 j_start = jts-1 j_end = MIN(jte,jde-1)+1 j_start_f = j_start j_end_f = j_end+1 !-- modify loop bounds if open or specified IF(degrade_xs) i_start = its IF(degrade_xe) i_end = MIN(ite,ide-1) IF(degrade_ys) then j_start = MAX(jts,jds+1) j_start_f = jds+3 ENDIF IF(degrade_ye) then j_end = MIN(jte,jde-2) j_end_f = jde-3 ENDIF ! compute fluxes, 5th order j_loop_y_flux_5 : DO j = j_start, j_end+1 IF( (j >= j_start_f ) .and. (j <= j_end_f) ) THEN ! use full stencil DO k=kts,ktf DO i = i_start, i_end dy = 2./(msft(i,j)+msft(i,j-1))/rdy mu = 0.5*(mut(i,j)+mut(i,j-1)) vel = rv(i,k,j) cr = vel*dt/dy/mu fqyl(i,k,j) = mu*(dy/dt)*flux_upwind(field_old(i,k,j-1), field_old(i,k,j ), cr) fqy( i, k, j ) = vel*flux5( & field(i,k,j-3), field(i,k,j-2), field(i,k,j-1), & field(i,k,j ), field(i,k,j+1), field(i,k,j+2), vel ) fqy(i,k,j) = fqy(i,k,j) - fqyl(i,k,j) ENDDO ENDDO ELSE IF ( j == jds+1 ) THEN ! 2nd order flux next to south boundary DO k=kts,ktf DO i = i_start, i_end dy = 2./(msft(i,j)+msft(i,j-1))/rdy mu = 0.5*(mut(i,j)+mut(i,j-1)) vel = rv(i,k,j) cr = vel*dt/dy/mu fqyl(i,k,j) = mu*(dy/dt)*flux_upwind(field_old(i,k,j-1), field_old(i,k,j ), cr) fqy(i,k, j) = 0.5*rv(i,k,j)* & (field(i,k,j)+field(i,k,j-1)) fqy(i,k,j) = fqy(i,k,j) - fqyl(i,k,j) ENDDO ENDDO ELSE IF ( j == jds+2 ) THEN ! third of 4th order flux 2 in from south boundary DO k=kts,ktf DO i = i_start, i_end dy = 2./(msft(i,j)+msft(i,j-1))/rdy mu = 0.5*(mut(i,j)+mut(i,j-1)) vel = rv(i,k,j) cr = vel*dt/dy/mu fqyl(i,k,j) = mu*(dy/dt)*flux_upwind(field_old(i,k,j-1), field_old(i,k,j ), cr) fqy( i, k, j ) = vel*flux3( & field(i,k,j-2),field(i,k,j-1),field(i,k,j),field(i,k,j+1),vel ) fqy(i,k,j) = fqy(i,k,j) - fqyl(i,k,j) ENDDO ENDDO ELSE IF ( j == jde-1 ) THEN ! 2nd order flux next to north boundary DO k=kts,ktf DO i = i_start, i_end dy = 2./(msft(i,j)+msft(i,j-1))/rdy mu = 0.5*(mut(i,j)+mut(i,j-1)) vel = rv(i,k,j) cr = vel*dt/dy/mu fqyl(i,k,j) = mu*(dy/dt)*flux_upwind(field_old(i,k,j-1), field_old(i,k,j ), cr) fqy(i, k, j ) = 0.5*rv(i,k,j)* & (field(i,k,j)+field(i,k,j-1)) fqy(i,k,j) = fqy(i,k,j) - fqyl(i,k,j) ENDDO ENDDO ELSE IF ( j == jde-2 ) THEN ! 3rd or 4th order flux 2 in from north boundary DO k=kts,ktf DO i = i_start, i_end dy = 2./(msft(i,j)+msft(i,j-1))/rdy mu = 0.5*(mut(i,j)+mut(i,j-1)) vel = rv(i,k,j) cr = vel*dt/dy/mu fqyl(i,k,j) = mu*(dy/dt)*flux_upwind(field_old(i,k,j-1), field_old(i,k,j ), cr) fqy( i, k, j) = vel*flux3( & field(i,k,j-2),field(i,k,j-1), & field(i,k,j),field(i,k,j+1),vel ) fqy(i,k,j) = fqy(i,k,j) - fqyl(i,k,j) ENDDO ENDDO ENDIF ENDDO j_loop_y_flux_5 ! next, x flux !-- these bounds are for periodic and sym conditions i_start = its-1 i_end = MIN(ite,ide-1)+1 i_start_f = i_start i_end_f = i_end+1 j_start = jts-1 j_end = MIN(jte,jde-1)+1 !-- modify loop bounds for open and specified b.c IF(degrade_ys) j_start = jts IF(degrade_ye) j_end = MIN(jte,jde-1) IF(degrade_xs) then i_start = MAX(ids+1,its) i_start_f = i_start+2 ENDIF IF(degrade_xe) then i_end = MIN(ide-2,ite) i_end_f = ide-3 ENDIF ! compute fluxes DO j = j_start, j_end ! 5th order flux DO k=kts,ktf DO i = i_start_f, i_end_f dx = 2./(msft(i,j)+msft(i-1,j))/rdx mu = 0.5*(mut(i,j)+mut(i-1,j)) vel = ru(i,k,j) cr = vel*dt/dx/mu fqxl(i,k,j) = mu*(dx/dt)*flux_upwind(field_old(i-1,k,j), field_old(i,k,j ), cr) fqx( i,k,j ) = vel*flux5( field(i-3,k,j), field(i-2,k,j), & field(i-1,k,j), field(i ,k,j), & field(i+1,k,j), field(i+2,k,j), & vel ) fqx(i,k,j) = fqx(i,k,j) - fqxl(i,k,j) ENDDO ENDDO ! lower order fluxes close to boundaries (if not periodic or symmetric) IF( degrade_xs ) THEN IF( i_start == ids+1 ) THEN ! second order flux next to the boundary i = ids+1 DO k=kts,ktf dx = 2./(msft(i,j)+msft(i-1,j))/rdx mu = 0.5*(mut(i,j)+mut(i-1,j)) vel = ru(i,k,j)/mu cr = vel*dt/dx fqxl(i,k,j) = mu*(dx/dt)*flux_upwind(field_old(i-1,k,j), field_old(i,k,j ), cr) fqx(i,k,j) = 0.5*(ru(i,k,j)) & *(field(i,k,j)+field(i-1,k,j)) fqx(i,k,j) = fqx(i,k,j) - fqxl(i,k,j) ENDDO ENDIF i = ids+2 DO k=kts,ktf dx = 2./(msft(i,j)+msft(i-1,j))/rdx mu = 0.5*(mut(i,j)+mut(i-1,j)) vel = ru(i,k,j) cr = vel*dt/dx/mu fqxl(i,k,j) = mu*(dx/dt)*flux_upwind(field_old(i-1,k,j), field_old(i,k,j ), cr) fqx( i,k,j ) = vel*flux3( field(i-2,k,j), field(i-1,k,j), & field(i ,k,j), field(i+1,k,j), & vel ) fqx(i,k,j) = fqx(i,k,j) - fqxl(i,k,j) ENDDO ENDIF IF( degrade_xe ) THEN IF( i_end == ide-2 ) THEN ! second order flux next to the boundary i = ide-1 DO k=kts,ktf dx = 2./(msft(i,j)+msft(i-1,j))/rdx mu = 0.5*(mut(i,j)+mut(i-1,j)) vel = ru(i,k,j) cr = vel*dt/dx/mu fqxl(i,k,j) = mu*(dx/dt)*flux_upwind(field_old(i-1,k,j), field_old(i,k,j ), cr) fqx(i,k,j) = 0.5*(ru(i,k,j)) & *(field(i,k,j)+field(i-1,k,j)) fqx(i,k,j) = fqx(i,k,j) - fqxl(i,k,j) ENDDO ENDIF i = ide-2 DO k=kts,ktf dx = 2./(msft(i,j)+msft(i-1,j))/rdx mu = 0.5*(mut(i,j)+mut(i-1,j)) vel = ru(i,k,j) cr = vel*dt/dx/mu fqxl(i,k,j) = mu*(dx/dt)*flux_upwind(field_old(i-1,k,j), field_old(i,k,j ), cr) fqx( i,k,j ) = vel*flux3( field(i-2,k,j), field(i-1,k,j), & field(i ,k,j), field(i+1,k,j), & vel ) fqx(i,k,j) = fqx(i,k,j) - fqxl(i,k,j) ENDDO ENDIF ENDDO ! enddo for outer J loop !--- end of 5th order horizontal flux calculation ELSE IF( horz_order == 4 ) THEN IF( config_flags%periodic_x .or. & config_flags%symmetric_xs .or. & (its > ids+1) ) degrade_xs = .false. IF( config_flags%periodic_x .or. & config_flags%symmetric_xe .or. & (ite < ide-2) ) degrade_xe = .false. IF( config_flags%periodic_y .or. & config_flags%symmetric_ys .or. & (jts > jds+1) ) degrade_ys = .false. IF( config_flags%periodic_y .or. & config_flags%symmetric_ye .or. & (jte < jde-2) ) degrade_ye = .false. !--------------- y - advection first !-- y flux compute; these bounds are for periodic and sym b.c. ktf=MIN(kte,kde-1) i_start = its-1 i_end = MIN(ite,ide-1)+1 j_start = jts-1 j_end = MIN(jte,jde-1)+1 j_start_f = j_start j_end_f = j_end+1 !-- modify loop bounds if open or specified IF(degrade_xs) i_start = its IF(degrade_xe) i_end = MIN(ite,ide-1) IF(degrade_ys) then j_start = MAX(jts,jds+1) j_start_f = jds+2 ENDIF IF(degrade_ye) then j_end = MIN(jte,jde-2) j_end_f = jde-2 ENDIF ! compute fluxes, 4th order j_loop_y_flux_4 : DO j = j_start, j_end+1 IF( (j >= j_start_f ) .and. (j <= j_end_f) ) THEN ! use full stencil DO k=kts,ktf DO i = i_start, i_end dy = 2./(msft(i,j)+msft(i,j-1))/rdy mu = 0.5*(mut(i,j)+mut(i,j-1)) vel = rv(i,k,j) cr = vel*dt/dy/mu fqyl(i,k,j) = mu*(dy/dt)*flux_upwind(field_old(i,k,j-1), field_old(i,k,j ), cr) fqy( i, k, j ) = vel*flux4( field(i,k,j-2), field(i,k,j-1), & field(i,k,j ), field(i,k,j+1), vel ) fqy(i,k,j) = fqy(i,k,j) - fqyl(i,k,j) ENDDO ENDDO ELSE IF ( j == jds+1 ) THEN ! 2nd order flux next to south boundary DO k=kts,ktf DO i = i_start, i_end dy = 2./(msft(i,j)+msft(i,j-1))/rdy mu = 0.5*(mut(i,j)+mut(i,j-1)) vel = rv(i,k,j) cr = vel*dt/dy/mu fqyl(i,k,j) = mu*(dy/dt)*flux_upwind(field_old(i,k,j-1), field_old(i,k,j ), cr) fqy(i,k, j) = 0.5*rv(i,k,j)* & (field(i,k,j)+field(i,k,j-1)) fqy(i,k,j) = fqy(i,k,j) - fqyl(i,k,j) ENDDO ENDDO ELSE IF ( j == jde-1 ) THEN ! 2nd order flux next to north boundary DO k=kts,ktf DO i = i_start, i_end dy = 2./(msft(i,j)+msft(i,j-1))/rdy mu = 0.5*(mut(i,j)+mut(i,j-1)) vel = rv(i,k,j) cr = vel*dt/dy/mu fqyl(i,k,j) = mu*(dy/dt)*flux_upwind(field_old(i,k,j-1), field_old(i,k,j ), cr) fqy(i, k, j ) = 0.5*rv(i,k,j)* & (field(i,k,j)+field(i,k,j-1)) fqy(i,k,j) = fqy(i,k,j) - fqyl(i,k,j) ENDDO ENDDO ENDIF ENDDO j_loop_y_flux_4 ! next, x flux !-- these bounds are for periodic and sym conditions i_start = its-1 i_end = MIN(ite,ide-1)+1 i_start_f = i_start i_end_f = i_end+1 j_start = jts-1 j_end = MIN(jte,jde-1)+1 !-- modify loop bounds for open and specified b.c IF(degrade_ys) j_start = jts IF(degrade_ye) j_end = MIN(jte,jde-1) IF(degrade_xs) then i_start = MAX(ids+1,its) i_start_f = i_start+1 ENDIF IF(degrade_xe) then i_end = MIN(ide-2,ite) i_end_f = ide-2 ENDIF ! compute fluxes DO j = j_start, j_end ! 4th order flux DO k=kts,ktf DO i = i_start_f, i_end_f dx = 2./(msft(i,j)+msft(i-1,j))/rdx mu = 0.5*(mut(i,j)+mut(i-1,j)) vel = ru(i,k,j) cr = vel*dt/dx/mu fqxl(i,k,j) = mu*(dx/dt)*flux_upwind(field_old(i-1,k,j), field_old(i,k,j ), cr) fqx( i,k,j ) = vel*flux4( field(i-2,k,j), field(i-1,k,j), & field(i ,k,j), field(i+1,k,j), vel ) fqx(i,k,j) = fqx(i,k,j) - fqxl(i,k,j) ENDDO ENDDO ! lower order fluxes close to boundaries (if not periodic or symmetric) IF( degrade_xs ) THEN IF( i_start == ids+1 ) THEN ! second order flux next to the boundary i = ids+1 DO k=kts,ktf dx = 2./(msft(i,j)+msft(i-1,j))/rdx mu = 0.5*(mut(i,j)+mut(i-1,j)) vel = ru(i,k,j)/mu cr = vel*dt/dx fqxl(i,k,j) = mu*(dx/dt)*flux_upwind(field_old(i-1,k,j), field_old(i,k,j ), cr) fqx(i,k,j) = 0.5*(ru(i,k,j)) & *(field(i,k,j)+field(i-1,k,j)) fqx(i,k,j) = fqx(i,k,j) - fqxl(i,k,j) ENDDO ENDIF ENDIF IF( degrade_xe ) THEN IF( i_end == ide-2 ) THEN ! second order flux next to the boundary i = ide-1 DO k=kts,ktf dx = 2./(msft(i,j)+msft(i-1,j))/rdx mu = 0.5*(mut(i,j)+mut(i-1,j)) vel = ru(i,k,j) cr = vel*dt/dx/mu fqxl(i,k,j) = mu*(dx/dt)*flux_upwind(field_old(i-1,k,j), field_old(i,k,j ), cr) fqx(i,k,j) = 0.5*(ru(i,k,j)) & *(field(i,k,j)+field(i-1,k,j)) fqx(i,k,j) = fqx(i,k,j) - fqxl(i,k,j) ENDDO ENDIF ENDIF ENDDO ! enddo for outer J loop !--- end of 4th order horizontal flux calculation ELSE IF( horz_order == 3 ) THEN IF( config_flags%periodic_x .or. & config_flags%symmetric_xs .or. & (its > ids+1) ) degrade_xs = .false. IF( config_flags%periodic_x .or. & config_flags%symmetric_xe .or. & (ite < ide-2) ) degrade_xe = .false. IF( config_flags%periodic_y .or. & config_flags%symmetric_ys .or. & (jts > jds+1) ) degrade_ys = .false. IF( config_flags%periodic_y .or. & config_flags%symmetric_ye .or. & (jte < jde-2) ) degrade_ye = .false. !--------------- y - advection first !-- y flux compute; these bounds are for periodic and sym b.c. ktf=MIN(kte,kde-1) i_start = its-1 i_end = MIN(ite,ide-1)+1 j_start = jts-1 j_end = MIN(jte,jde-1)+1 j_start_f = j_start j_end_f = j_end+1 !-- modify loop bounds if open or specified IF(degrade_xs) i_start = its IF(degrade_xe) i_end = MIN(ite,ide-1) IF(degrade_ys) then j_start = MAX(jts,jds+1) j_start_f = jds+2 ENDIF IF(degrade_ye) then j_end = MIN(jte,jde-2) j_end_f = jde-2 ENDIF ! compute fluxes, 3rd order j_loop_y_flux_3 : DO j = j_start, j_end+1 IF( (j >= j_start_f ) .and. (j <= j_end_f) ) THEN ! use full stencil DO k=kts,ktf DO i = i_start, i_end dy = 2./(msft(i,j)+msft(i,j-1))/rdy mu = 0.5*(mut(i,j)+mut(i,j-1)) vel = rv(i,k,j) cr = vel*dt/dy/mu fqyl(i,k,j) = mu*(dy/dt)*flux_upwind(field_old(i,k,j-1), field_old(i,k,j ), cr) fqy( i, k, j ) = vel*flux3( field(i,k,j-2), field(i,k,j-1), & field(i,k,j ), field(i,k,j+1), vel ) fqy(i,k,j) = fqy(i,k,j) - fqyl(i,k,j) ENDDO ENDDO ELSE IF ( j == jds+1 ) THEN ! 2nd order flux next to south boundary DO k=kts,ktf DO i = i_start, i_end dy = 2./(msft(i,j)+msft(i,j-1))/rdy mu = 0.5*(mut(i,j)+mut(i,j-1)) vel = rv(i,k,j) cr = vel*dt/dy/mu fqyl(i,k,j) = mu*(dy/dt)*flux_upwind(field_old(i,k,j-1), field_old(i,k,j ), cr) fqy(i,k, j) = 0.5*rv(i,k,j)* & (field(i,k,j)+field(i,k,j-1)) fqy(i,k,j) = fqy(i,k,j) - fqyl(i,k,j) ENDDO ENDDO ELSE IF ( j == jde-1 ) THEN ! 2nd order flux next to north boundary DO k=kts,ktf DO i = i_start, i_end dy = 2./(msft(i,j)+msft(i,j-1))/rdy mu = 0.5*(mut(i,j)+mut(i,j-1)) vel = rv(i,k,j) cr = vel*dt/dy/mu fqyl(i,k,j) = mu*(dy/dt)*flux_upwind(field_old(i,k,j-1), field_old(i,k,j ), cr) fqy(i, k, j ) = 0.5*rv(i,k,j)* & (field(i,k,j)+field(i,k,j-1)) fqy(i,k,j) = fqy(i,k,j) - fqyl(i,k,j) ENDDO ENDDO ENDIF ENDDO j_loop_y_flux_3 ! next, x flux !-- these bounds are for periodic and sym conditions i_start = its-1 i_end = MIN(ite,ide-1)+1 i_start_f = i_start i_end_f = i_end+1 j_start = jts-1 j_end = MIN(jte,jde-1)+1 !-- modify loop bounds for open and specified b.c IF(degrade_ys) j_start = jts IF(degrade_ye) j_end = MIN(jte,jde-1) IF(degrade_xs) then i_start = MAX(ids+1,its) i_start_f = i_start+1 ENDIF IF(degrade_xe) then i_end = MIN(ide-2,ite) i_end_f = ide-2 ENDIF ! compute fluxes DO j = j_start, j_end ! 4th order flux DO k=kts,ktf DO i = i_start_f, i_end_f dx = 2./(msft(i,j)+msft(i-1,j))/rdx mu = 0.5*(mut(i,j)+mut(i-1,j)) vel = ru(i,k,j) cr = vel*dt/dx/mu fqxl(i,k,j) = mu*(dx/dt)*flux_upwind(field_old(i-1,k,j), field_old(i,k,j ), cr) fqx( i,k,j ) = vel*flux3( field(i-2,k,j), field(i-1,k,j), & field(i ,k,j), field(i+1,k,j), vel ) fqx(i,k,j) = fqx(i,k,j) - fqxl(i,k,j) ENDDO ENDDO ! lower order fluxes close to boundaries (if not periodic or symmetric) IF( degrade_xs ) THEN IF( i_start == ids+1 ) THEN ! second order flux next to the boundary i = ids+1 DO k=kts,ktf dx = 2./(msft(i,j)+msft(i-1,j))/rdx mu = 0.5*(mut(i,j)+mut(i-1,j)) vel = ru(i,k,j)/mu cr = vel*dt/dx fqxl(i,k,j) = mu*(dx/dt)*flux_upwind(field_old(i-1,k,j), field_old(i,k,j ), cr) fqx(i,k,j) = 0.5*(ru(i,k,j)) & *(field(i,k,j)+field(i-1,k,j)) fqx(i,k,j) = fqx(i,k,j) - fqxl(i,k,j) ENDDO ENDIF ENDIF IF( degrade_xe ) THEN IF( i_end == ide-2 ) THEN ! second order flux next to the boundary i = ide-1 DO k=kts,ktf dx = 2./(msft(i,j)+msft(i-1,j))/rdx mu = 0.5*(mut(i,j)+mut(i-1,j)) vel = ru(i,k,j) cr = vel*dt/dx/mu fqxl(i,k,j) = mu*(dx/dt)*flux_upwind(field_old(i-1,k,j), field_old(i,k,j ), cr) fqx(i,k,j) = 0.5*(ru(i,k,j)) & *(field(i,k,j)+field(i-1,k,j)) fqx(i,k,j) = fqx(i,k,j) - fqxl(i,k,j) ENDDO ENDIF ENDIF ENDDO ! enddo for outer J loop !--- end of 3rd order horizontal flux calculation ELSE IF( horz_order == 2 ) THEN IF( config_flags%periodic_x .or. & config_flags%symmetric_xs .or. & (its > ids) ) degrade_xs = .false. IF( config_flags%periodic_x .or. & config_flags%symmetric_xe .or. & (ite < ide-1) ) degrade_xe = .false. IF( config_flags%periodic_y .or. & config_flags%symmetric_ys .or. & (jts > jds) ) degrade_ys = .false. IF( config_flags%periodic_y .or. & config_flags%symmetric_ye .or. & (jte < jde-1) ) degrade_ye = .false. !-- y flux compute; these bounds are for periodic and sym b.c. ktf=MIN(kte,kde-1) i_start = its-1 i_end = MIN(ite,ide-1)+1 j_start = jts-1 j_end = MIN(jte,jde-1)+1 !-- modify loop bounds if open or specified IF(degrade_xs) i_start = its IF(degrade_xe) i_end = MIN(ite,ide-1) IF(degrade_ys) j_start = MAX(jts,jds+1) IF(degrade_ye) j_end = MIN(jte,jde-2) ! compute fluxes, 2nd order, y flux DO j = j_start, j_end+1 DO k=kts,ktf DO i = i_start, i_end dy = 2./(msft(i,j)+msft(i,j-1))/rdy mu = 0.5*(mut(i,j)+mut(i,j-1)) vel = rv(i,k,j) cr = vel*dt/dy/mu fqyl(i,k,j) = mu*(dy/dt)*flux_upwind(field_old(i,k,j-1), field_old(i,k,j ), cr) fqy(i,k, j) = 0.5*rv(i,k,j)* & (field(i,k,j)+field(i,k,j-1)) fqy(i,k,j) = fqy(i,k,j) - fqyl(i,k,j) ENDDO ENDDO ENDDO ! next, x flux DO j = j_start, j_end DO k=kts,ktf DO i = i_start, i_end+1 dx = 2./(msft(i,j)+msft(i-1,j))/rdx mu = 0.5*(mut(i,j)+mut(i-1,j)) vel = ru(i,k,j) cr = vel*dt/dx/mu fqxl(i,k,j) = mu*(dx/dt)*flux_upwind(field_old(i-1,k,j), field_old(i,k,j ), cr) fqx( i,k,j ) = vel*flux3( field(i-2,k,j), field(i-1,k,j), & field(i ,k,j), field(i+1,k,j), & vel ) fqx(i,k,j) = fqx(i,k,j) - fqxl(i,k,j) ENDDO ENDDO ENDDO !--- end of 3nd order horizontal flux calculation ELSE WRITE ( wrf_err_message , * ) 'module_advect: advect_scalar_pd, h_order not known ',horz_order CALL wrf_error_fatal ( TRIM( wrf_err_message ) ) ENDIF horizontal_order_test ! pick up the rest of the horizontal radiation boundary conditions. ! (these are the computations that don't require 'cb'. ! first, set to index ranges i_start = its i_end = MIN(ite,ide-1) j_start = jts j_end = MIN(jte,jde-1) ! compute x (u) conditions for v, w, or scalar IF( (config_flags%open_xs) .and. (its == ids) ) THEN DO j = j_start, j_end DO k = kts, ktf ub = MIN( 0.5*(ru(its,k,j)+ru(its+1,k,j)), 0. ) tendency(its,k,j) = tendency(its,k,j) & - rdx*( & ub*( field_old(its+1,k,j) & - field_old(its ,k,j) ) + & field(its,k,j)*(ru(its+1,k,j)-ru(its,k,j)) & ) ENDDO ENDDO ENDIF IF( (config_flags%open_xe) .and. (ite == ide) ) THEN DO j = j_start, j_end DO k = kts, ktf ub = MAX( 0.5*(ru(ite-1,k,j)+ru(ite,k,j)), 0. ) tendency(i_end,k,j) = tendency(i_end,k,j) & - rdx*( & ub*( field_old(i_end ,k,j) & - field_old(i_end-1,k,j) ) + & field(i_end,k,j)*(ru(ite,k,j)-ru(ite-1,k,j)) & ) ENDDO ENDDO ENDIF IF( (config_flags%open_ys) .and. (jts == jds) ) THEN DO i = i_start, i_end DO k = kts, ktf vb = MIN( 0.5*(rv(i,k,jts)+rv(i,k,jts+1)), 0. ) tendency(i,k,jts) = tendency(i,k,jts) & - rdy*( & vb*( field_old(i,k,jts+1) & - field_old(i,k,jts ) ) + & field(i,k,jts)*(rv(i,k,jts+1)-rv(i,k,jts)) & ) ENDDO ENDDO ENDIF IF( (config_flags%open_ye) .and. (jte == jde)) THEN DO i = i_start, i_end DO k = kts, ktf vb = MAX( 0.5*(rv(i,k,jte-1)+rv(i,k,jte)), 0. ) tendency(i,k,j_end) = tendency(i,k,j_end) & - rdy*( & vb*( field_old(i,k,j_end ) & - field_old(i,k,j_end-1) ) + & field(i,k,j_end)*(rv(i,k,jte)-rv(i,k,jte-1)) & ) ENDDO ENDDO ENDIF !-------------------- vertical advection !-- loop bounds for periodic or sym conditions i_start = its-1 i_end = MIN(ite,ide-1)+1 j_start = jts-1 j_end = MIN(jte,jde-1)+1 !-- loop bounds for open or specified conditions IF(degrade_xs) i_start = its IF(degrade_xe) i_end = MIN(ite,ide-1) IF(degrade_ys) j_start = jts IF(degrade_ye) j_end = MIN(jte,jde-1) vert_order_test : IF (vert_order == 6) THEN DO j = j_start, j_end DO i = i_start, i_end fqz(i,1,j) = 0. fqzl(i,1,j) = 0. fqz(i,kde,j) = 0. fqzl(i,kde,j) = 0. ENDDO DO k=kts+3,ktf-2 DO i = i_start, i_end dz = 2./(rdzw(k)+rdzw(k-1)) mu = 0.5*(mut(i,j)+mut(i,j)) vel = rom(i,k,j) cr = vel*dt/dz/mu fqzl(i,k,j) = mu*(dz/dt)*flux_upwind(field_old(i,k-1,j), field_old(i,k,j ), cr) fqz(i,k,j) = vel*flux6( field(i,k-3,j), field(i,k-2,j), field(i,k-1,j), & field(i,k ,j), field(i,k+1,j), field(i,k+2,j), -vel ) fqz(i,k,j) = fqz(i,k,j) - fqzl(i,k,j) ENDDO ENDDO DO i = i_start, i_end k=kts+1 dz = 2./(rdzw(k)+rdzw(k-1)) mu = 0.5*(mut(i,j)+mut(i,j)) vel = rom(i,k,j) cr = vel*dt/dz/mu fqzl(i,k,j) = mu*(dz/dt)*flux_upwind(field_old(i,k-1,j), field_old(i,k,j ), cr) fqz(i,k,j)=rom(i,k,j)*(fzm(k)*field(i,k,j)+fzp(k)*field(i,k-1,j)) fqz(i,k,j) = fqz(i,k,j) - fqzl(i,k,j) k=kts+2 dz = 2./(rdzw(k)+rdzw(k-1)) mu = 0.5*(mut(i,j)+mut(i,j)) vel = rom(i,k,j) cr = vel*dt/dz/mu fqzl(i,k,j) = mu*(dz/dt)*flux_upwind(field_old(i,k-1,j), field_old(i,k,j ), cr) fqz(i,k,j) = vel*flux4( & field(i,k-2,j), field(i,k-1,j), & field(i,k ,j), field(i,k+1,j), -vel ) fqz(i,k,j) = fqz(i,k,j) - fqzl(i,k,j) k=ktf-1 dz = 2./(rdzw(k)+rdzw(k-1)) mu = 0.5*(mut(i,j)+mut(i,j)) vel = rom(i,k,j) cr = vel*dt/dz/mu fqzl(i,k,j) = mu*(dz/dt)*flux_upwind(field_old(i,k-1,j), field_old(i,k,j ), cr) fqz(i,k,j) = vel*flux4( & field(i,k-2,j), field(i,k-1,j), & field(i,k ,j), field(i,k+1,j), -vel ) fqz(i,k,j) = fqz(i,k,j) - fqzl(i,k,j) k=ktf dz = 2./(rdzw(k)+rdzw(k-1)) mu = 0.5*(mut(i,j)+mut(i,j)) vel = rom(i,k,j) cr = vel*dt/dz/mu fqzl(i,k,j) = mu*(dz/dt)*flux_upwind(field_old(i,k-1,j), field_old(i,k,j ), cr) fqz(i,k,j)=rom(i,k,j)*(fzm(k)*field(i,k,j)+fzp(k)*field(i,k-1,j)) fqz(i,k,j) = fqz(i,k,j) - fqzl(i,k,j) ENDDO ENDDO ELSE IF (vert_order == 5) THEN DO j = j_start, j_end DO i = i_start, i_end fqz(i,1,j) = 0. fqzl(i,1,j) = 0. fqz(i,kde,j) = 0. fqzl(i,kde,j) = 0. ENDDO DO k=kts+3,ktf-2 DO i = i_start, i_end dz = 2./(rdzw(k)+rdzw(k-1)) mu = 0.5*(mut(i,j)+mut(i,j)) vel = rom(i,k,j) cr = vel*dt/dz/mu fqzl(i,k,j) = mu*(dz/dt)*flux_upwind(field_old(i,k-1,j), field_old(i,k,j ), cr) fqz(i,k,j) = vel*flux5( field(i,k-3,j), field(i,k-2,j), field(i,k-1,j), & field(i,k ,j), field(i,k+1,j), field(i,k+2,j), -vel ) fqz(i,k,j) = fqz(i,k,j) - fqzl(i,k,j) ENDDO ENDDO DO i = i_start, i_end k=kts+1 dz = 2./(rdzw(k)+rdzw(k-1)) mu = 0.5*(mut(i,j)+mut(i,j)) vel = rom(i,k,j) cr = vel*dt/dz/mu fqzl(i,k,j) = mu*(dz/dt)*flux_upwind(field_old(i,k-1,j), field_old(i,k,j ), cr) fqz(i,k,j)=rom(i,k,j)*(fzm(k)*field(i,k,j)+fzp(k)*field(i,k-1,j)) fqz(i,k,j) = fqz(i,k,j) - fqzl(i,k,j) k=kts+2 dz = 2./(rdzw(k)+rdzw(k-1)) mu = 0.5*(mut(i,j)+mut(i,j)) vel = rom(i,k,j) cr = vel*dt/dz/mu fqzl(i,k,j) = mu*(dz/dt)*flux_upwind(field_old(i,k-1,j), field_old(i,k,j ), cr) fqz(i,k,j) = vel*flux3( & field(i,k-2,j), field(i,k-1,j), & field(i,k ,j), field(i,k+1,j), -vel ) fqz(i,k,j) = fqz(i,k,j) - fqzl(i,k,j) k=ktf-1 dz = 2./(rdzw(k)+rdzw(k-1)) mu = 0.5*(mut(i,j)+mut(i,j)) vel = rom(i,k,j) cr = vel*dt/dz/mu fqzl(i,k,j) = mu*(dz/dt)*flux_upwind(field_old(i,k-1,j), field_old(i,k,j ), cr) fqz(i,k,j) = vel*flux3( & field(i,k-2,j), field(i,k-1,j), & field(i,k ,j), field(i,k+1,j), -vel ) fqz(i,k,j) = fqz(i,k,j) - fqzl(i,k,j) k=ktf dz = 2./(rdzw(k)+rdzw(k-1)) mu = 0.5*(mut(i,j)+mut(i,j)) vel = rom(i,k,j) cr = vel*dt/dz/mu fqzl(i,k,j) = mu*(dz/dt)*flux_upwind(field_old(i,k-1,j), field_old(i,k,j ), cr) fqz(i,k,j)=rom(i,k,j)*(fzm(k)*field(i,k,j)+fzp(k)*field(i,k-1,j)) fqz(i,k,j) = fqz(i,k,j) - fqzl(i,k,j) ENDDO ENDDO ELSE IF (vert_order == 4) THEN DO j = j_start, j_end DO i = i_start, i_end fqz(i,1,j) = 0. fqzl(i,1,j) = 0. fqz(i,kde,j) = 0. fqzl(i,kde,j) = 0. ENDDO DO k=kts+2,ktf-1 DO i = i_start, i_end dz = 2./(rdzw(k)+rdzw(k-1)) mu = 0.5*(mut(i,j)+mut(i,j)) vel = rom(i,k,j) cr = vel*dt/dz/mu fqzl(i,k,j) = mu*(dz/dt)*flux_upwind(field_old(i,k-1,j), field_old(i,k,j ), cr) fqz(i,k,j) = vel*flux4( & field(i,k-2,j), field(i,k-1,j), & field(i,k ,j), field(i,k+1,j), -vel ) fqz(i,k,j) = fqz(i,k,j) - fqzl(i,k,j) ENDDO ENDDO DO i = i_start, i_end k=kts+1 dz = 2./(rdzw(k)+rdzw(k-1)) mu = 0.5*(mut(i,j)+mut(i,j)) vel = rom(i,k,j) cr = vel*dt/dz/mu fqzl(i,k,j) = mu*(dz/dt)*flux_upwind(field_old(i,k-1,j), field_old(i,k,j ), cr) fqz(i,k,j)=rom(i,k,j)*(fzm(k)*field(i,k,j)+fzp(k)*field(i,k-1,j)) fqz(i,k,j) = fqz(i,k,j) - fqzl(i,k,j) k=ktf dz = 2./(rdzw(k)+rdzw(k-1)) mu = 0.5*(mut(i,j)+mut(i,j)) vel = rom(i,k,j) cr = vel*dt/dz/mu fqzl(i,k,j) = mu*(dz/dt)*flux_upwind(field_old(i,k-1,j), field_old(i,k,j ), cr) fqz(i,k,j)=rom(i,k,j)*(fzm(k)*field(i,k,j)+fzp(k)*field(i,k-1,j)) fqz(i,k,j) = fqz(i,k,j) - fqzl(i,k,j) ENDDO ENDDO ELSE IF (vert_order == 3) THEN DO j = j_start, j_end DO i = i_start, i_end fqz(i,1,j) = 0. fqzl(i,1,j) = 0. fqz(i,kde,j) = 0. fqzl(i,kde,j) = 0. ENDDO DO k=kts+2,ktf-1 DO i = i_start, i_end dz = 2./(rdzw(k)+rdzw(k-1)) mu = 0.5*(mut(i,j)+mut(i,j)) vel = rom(i,k,j) cr = vel*dt/dz/mu fqzl(i,k,j) = mu*(dz/dt)*flux_upwind(field_old(i,k-1,j), field_old(i,k,j ), cr) fqz(i,k,j) = vel*flux3( & field(i,k-2,j), field(i,k-1,j), & field(i,k ,j), field(i,k+1,j), -vel ) fqz(i,k,j) = fqz(i,k,j) - fqzl(i,k,j) ENDDO ENDDO DO i = i_start, i_end k=kts+1 dz = 2./(rdzw(k)+rdzw(k-1)) mu = 0.5*(mut(i,j)+mut(i,j)) vel = rom(i,k,j) cr = vel*dt/dz/mu fqzl(i,k,j) = mu*(dz/dt)*flux_upwind(field_old(i,k-1,j), field_old(i,k,j ), cr) fqz(i,k,j)=rom(i,k,j)*(fzm(k)*field(i,k,j)+fzp(k)*field(i,k-1,j)) fqz(i,k,j) = fqz(i,k,j) - fqzl(i,k,j) k=ktf dz = 2./(rdzw(k)+rdzw(k-1)) mu = 0.5*(mut(i,j)+mut(i,j)) vel = rom(i,k,j) cr = vel*dt/dz/mu fqzl(i,k,j) = mu*(dz/dt)*flux_upwind(field_old(i,k-1,j), field_old(i,k,j ), cr) fqz(i,k,j)=rom(i,k,j)*(fzm(k)*field(i,k,j)+fzp(k)*field(i,k-1,j)) fqz(i,k,j) = fqz(i,k,j) - fqzl(i,k,j) ENDDO ENDDO ELSE IF (vert_order == 2) THEN DO j = j_start, j_end DO i = i_start, i_end fqz(i,1,j) = 0. fqzl(i,1,j) = 0. fqz(i,kde,j) = 0. fqzl(i,kde,j) = 0. ENDDO DO k=kts+1,ktf DO i = i_start, i_end dz = 2./(rdzw(k)+rdzw(k-1)) mu = 0.5*(mut(i,j)+mut(i,j)) vel = rom(i,k,j) cr = vel*dt/dz/mu fqzl(i,k,j) = mu*(dz/dt)*flux_upwind(field_old(i,k-1,j), field_old(i,k,j ), cr) fqz(i,k,j)=rom(i,k,j)*(fzm(k)*field(i,k,j)+fzp(k)*field(i,k-1,j)) fqz(i,k,j) = fqz(i,k,j) - fqzl(i,k,j) ENDDO ENDDO ENDDO ELSE WRITE (wrf_err_message,*) ' advect_scalar_pd, v_order not known ',vert_order CALL wrf_error_fatal ( wrf_err_message ) ENDIF vert_order_test IF (pd_limit) THEN ! positive definite filter i_start = its-1 i_end = MIN(ite,ide-1)+1 j_start = jts-1 j_end = MIN(jte,jde-1)+1 !-- loop bounds for open or specified conditions IF(degrade_xs) i_start = its IF(degrade_xe) i_end = MIN(ite,ide-1) IF(degrade_ys) j_start = jts IF(degrade_ye) j_end = MIN(jte,jde-1) IF(config_flags%specified .or. config_flags%nested) THEN IF (degrade_xs) i_start = MAX(its,ids+1) IF (degrade_xe) i_end = MIN(ite,ide-2) IF (degrade_ys) j_start = MAX(jts,jds+1) IF (degrade_ye) j_end = MIN(jte,jde-2) END IF IF(config_flags%open_xs) THEN IF (degrade_xs) i_start = MAX(its,ids+1) END IF IF(config_flags%open_xe) THEN IF (degrade_xe) i_end = MIN(ite,ide-2) END IF IF(config_flags%open_ys) THEN IF (degrade_ys) j_start = MAX(jts,jds+1) END IF IF(config_flags%open_ye) THEN IF (degrade_ye) j_end = MIN(jte,jde-2) END IF !-- here is the limiter... DO j=j_start, j_end DO k=kts, ktf DO i=i_start, i_end ph_low = (mub(i,j)+mu_old(i,j))*field_old(i,k,j) & - dt*( msft(i,j)*( rdx*(fqxl(i+1,k,j)-fqxl(i,k,j)) & +rdy*(fqyl(i,k,j+1)-fqyl(i,k,j)) ) & +rdzw(k)*(fqzl(i,k+1,j)-fqzl(i,k,j)) ) flux_out = dt*(msft(i,j)*( rdx*( max(0.,fqx (i+1,k,j)) & -min(0.,fqx (i ,k,j)) ) & +rdy*( max(0.,fqy (i,k,j+1)) & -min(0.,fqy (i,k,j )) ) ) & +rdzw(k)*( min(0.,fqz (i,k+1,j)) & -max(0.,fqz (i,k ,j)) ) ) IF( flux_out .gt. ph_low ) THEN scale = max(0.,ph_low/(flux_out+eps)) IF( fqx (i+1,k,j) .gt. 0.) fqx(i+1,k,j) = scale*fqx(i+1,k,j) IF( fqx (i ,k,j) .lt. 0.) fqx(i ,k,j) = scale*fqx(i ,k,j) IF( fqy (i,k,j+1) .gt. 0.) fqy(i,k,j+1) = scale*fqy(i,k,j+1) IF( fqy (i,k,j ) .lt. 0.) fqy(i,k,j ) = scale*fqy(i,k,j ) ! note: z flux is opposite sign in mass coordinate because ! vertical coordinate decreases with increasing k IF( fqz (i,k+1,j) .lt. 0.) fqz(i,k+1,j) = scale*fqz(i,k+1,j) IF( fqz (i,k ,j) .gt. 0.) fqz(i,k ,j) = scale*fqz(i,k ,j) END IF ENDDO ENDDO ENDDO END IF ! add in the pd-limited flux divergence i_start = its i_end = MIN(ite,ide-1) j_start = jts j_end = MIN(jte,jde-1) DO j = j_start, j_end DO k = kts, ktf DO i = i_start, i_end tendency (i,k,j) = tendency(i,k,j) & -rdzw(k)*( fqz (i,k+1,j)-fqz (i,k,j) & +fqzl(i,k+1,j)-fqzl(i,k,j)) ENDDO ENDDO ENDDO ! x flux divergence ! IF(degrade_xs) i_start = i_start + 1 IF(degrade_xe) i_end = i_end - 1 DO j = j_start, j_end DO k = kts, ktf DO i = i_start, i_end tendency (i,k,j) = tendency(i,k,j) & - msft(i,j)*( rdx*( fqx (i+1,k,j)-fqx (i,k,j) & +fqxl(i+1,k,j)-fqxl(i,k,j)) ) ENDDO ENDDO ENDDO ! y flux divergence ! i_start = its i_end = MIN(ite,ide-1) IF(degrade_ys) j_start = j_start + 1 IF(degrade_ye) j_end = j_end - 1 DO j = j_start, j_end DO k = kts, ktf DO i = i_start, i_end tendency (i,k,j) = tendency(i,k,j) & - msft(i,j)*( rdy*( fqy (i,k,j+1)-fqy (i,k,j) & +fqyl(i,k,j+1)-fqyl(i,k,j)) ) ENDDO ENDDO ENDDO END SUBROUTINE advect_scalar_pd !---------------------------------------------------------------- END MODULE module_advect_em