! WRF:MODEL_LAYER:PHYSICS MODULE module_diffusion_em USE module_configure USE module_bc USE module_state_description USE module_big_step_utilities_em USE module_model_constants USE module_wrf_error CONTAINS !======================================================================= !======================================================================= SUBROUTINE cal_deform_and_div( config_flags, u, v, w, div, & defor11, defor22, defor33, & defor12, defor13, defor23, & u_base, v_base, msfux, msfuy, & msfvx, msfvy, msftx, msfty, & rdx, rdy, dn, dnw, rdz, rdzw, & fnm, fnp, cf1, cf2, cf3, zx, zy, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) ! History: Sep 2003 Changes by Jason Knievel and George Bryan, NCAR ! Oct 2001 Converted to mass core by Bill Skamarock, NCAR ! ... ... ! Purpose: This routine calculates deformation and 3-d divergence. ! References: Klemp and Wilhelmson (JAS 1978) ! Chen and Dudhia (NCAR WRF physics report 2000) !----------------------------------------------------------------------- ! Comments 10-MAR-05 ! Equations 13a-f, Chen and Dudhia 2000, Appendix A: ! Eqn 13a: D11=defor11= 2m^2 * (partial du^/dX + partial dpsi/dx * partial du^/dpsi) ! Eqn 13b: D22=defor22= 2m^2 * (partial dv^/dY + partial dpsi/dy * partial dv^/dpsi) ! Eqn 13c: D33=defor33= 2 * partial dw/dz [SIMPLER FORM] ! Eqn 13d: D12=defor12= m^2 * (partial dv^/dX + partial du^/dY + ! partial dpsi/dx * partial dv^/dpsi + ! partial dpsi/dy * partial du^/dpsi) ! Eqn 13e: D13=defor13= m^2 * (partial dw^/dX + partial dpsi/dx * partial dw^/dpsi) ! + partial du/dz [SIMPLER FORM] ! Eqn 13f: D23=defor23= m^2 * (partial dw^/dY + partial dpsi/dy * partial dw^/dpsi) ! + partial dv/dz [SIMPLER FORM] !----------------------------------------------------------------------- ! Begin declarations. IMPLICIT NONE 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, INTENT( IN ) & :: rdx, rdy, cf1, cf2, cf3 REAL, DIMENSION( kms:kme ), INTENT( IN ) & :: fnm, fnp, dn, dnw, u_base, v_base REAL, DIMENSION( ims:ime , jms:jme ), INTENT( IN ) & :: msfux, msfuy, msfvx, msfvy, msftx, msfty REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( IN ) & :: u, v, w, zx, zy, rdz, rdzw REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( INOUT ) & :: defor11, defor22, defor33, defor12, defor13, defor23, div ! Local variables. INTEGER & :: i, j, k, ktf, ktes1, ktes2, i_start, i_end, j_start, j_end REAL & :: tmp, tmpzx, tmpzy, tmpzeta_z, cft1, cft2 REAL, DIMENSION( its:ite, jts:jte ) & :: mm, zzavg, zeta_zd12 REAL, DIMENSION( its-2:ite+2, kts:kte, jts-2:jte+2 ) & :: tmp1, hat, hatavg ! End declarations. !----------------------------------------------------------------------- ! Comments 10-MAR-2005 ! Treat all differentials as 'div-style' [or 'curl-style'], ! i.e., du/dX becomes (in map coordinate space) mx*my * d(u/my)/dx, ! NB - all equations referred to here are from Chen and Dudhia 2002, from the ! WRF physics documents web pages: ! http://www.mmm.ucar.edu/wrf/users/docs/wrf-doc-physics.pdf !======================================================================= ! In the following section, calculate 3-d divergence and the first three ! (defor11, defor22, defor33) of six deformation terms. ktes1 = kte-1 ktes2 = kte-2 cft2 = - 0.5 * dnw(ktes1) / dn(ktes1) cft1 = 1.0 - cft2 ktf = MIN( kte, kde-1 ) i_start = its i_end = MIN( ite, ide-1 ) j_start = jts j_end = MIN( jte, jde-1 ) ! Square the map scale factor. DO j = j_start, j_end DO i = i_start, i_end mm(i,j) = msftx(i,j) * msfty(i,j) END DO END DO !----------------------------------------------------------------------- ! Calculate du/dx. ! Apply a coordinate transformation to zonal velocity, u. DO j = j_start, j_end DO k = kts, ktf DO i = i_start, i_end+1 hat(i,k,j) = u(i,k,j) / msfuy(i,j) END DO END DO END DO ! Average in x and z. DO j=j_start,j_end DO k=kts+1,ktf DO i=i_start,i_end hatavg(i,k,j) = 0.5 * & ( fnm(k) * ( hat(i,k ,j) + hat(i+1, k,j) ) + & fnp(k) * ( hat(i,k-1,j) + hat(i+1,k-1,j) ) ) END DO END DO END DO ! Extrapolate to top and bottom of domain (to w levels). DO j = j_start, j_end DO i = i_start, i_end hatavg(i,1,j) = 0.5 * ( & cf1 * hat(i ,1,j) + & cf2 * hat(i ,2,j) + & cf3 * hat(i ,3,j) + & cf1 * hat(i+1,1,j) + & cf2 * hat(i+1,2,j) + & cf3 * hat(i+1,3,j) ) hatavg(i,kte,j) = 0.5 * ( & cft1 * ( hat(i,ktes1,j) + hat(i+1,ktes1,j) ) + & cft2 * ( hat(i,ktes2,j) + hat(i+1,ktes2,j) ) ) END DO END DO ! Comments 10-MAR-05 ! Eqn 13a: D11=defor11= 2m^2 * (partial du^/dX + partial dpsi/dx * partial du^/dpsi) ! Below, D11 is set = 2*tmp1 ! => tmp1 = m^2 * (partial du^/dX + partial dpsi/dx * partial du^/dpsi) ! tmpzx = averaged value of dpsi/dx (=zx) DO j = j_start, j_end DO k = kts, ktf DO i = i_start, i_end tmpzx = 0.25 * ( & zx(i,k ,j) + zx(i+1,k ,j) + & zx(i,k+1,j) + zx(i+1,k+1,j) ) tmp1(i,k,j) = ( hatavg(i,k+1,j) - hatavg(i,k,j) ) *tmpzx * rdzw(i,k,j) ! tmp1 to here = partial dpsi/dx * partial du^/dpsi: END DO END DO END DO DO j = j_start, j_end DO k = kts, ktf DO i = i_start, i_end tmp1(i,k,j) = mm(i,j) * ( rdx * ( hat(i+1,k,j) - hat(i,k,j) ) - & tmp1(i,k,j)) END DO END DO END DO ! End calculation of du/dx. !----------------------------------------------------------------------- !----------------------------------------------------------------------- ! Calculate defor11 (2*du/dx). ! Comments 10-MAR-05 ! Eqn 13a: D11=defor11= 2 m^2 * (partial du^/dX + partial dpsi/dx * partial du^/dpsi) ! = 2*tmp1 DO j = j_start, j_end DO k = kts, ktf DO i = i_start, i_end defor11(i,k,j) = 2.0 * tmp1(i,k,j) END DO END DO END DO ! End calculation of defor11. !----------------------------------------------------------------------- !----------------------------------------------------------------------- ! Calculate zonal divergence (du/dx) and add it to the divergence array. DO j = j_start, j_end DO k = kts, ktf DO i = i_start, i_end div(i,k,j) = tmp1(i,k,j) END DO END DO END DO ! End calculation of zonal divergence. !----------------------------------------------------------------------- !----------------------------------------------------------------------- ! Calculate dv/dy. ! Apply a coordinate transformation to meridional velocity, v. DO j = j_start, j_end+1 DO k = kts, ktf DO i = i_start, i_end ! Because msfvx at the poles will be undefined (1./0.), we will have ! trouble. But we are OK since v at the poles is 0., and that takes ! precedence in this case. IF ((config_flags%polar) .AND. ((j == jds) .OR. (j == jde))) THEN hat(i,k,j) = 0. ELSE ! normal code hat(i,k,j) = v(i,k,j) / msfvx(i,j) ENDIF END DO END DO END DO ! Account for the slope in y of eta surfaces. DO j=j_start,j_end DO k=kts+1,ktf DO i=i_start,i_end hatavg(i,k,j) = 0.5 * ( & fnm(k) * ( hat(i,k ,j) + hat(i,k ,j+1) ) + & fnp(k) * ( hat(i,k-1,j) + hat(i,k-1,j+1) ) ) END DO END DO END DO ! Extrapolate to top and bottom of domain (to w levels). DO j = j_start, j_end DO i = i_start, i_end hatavg(i,1,j) = 0.5 * ( & cf1 * hat(i,1,j ) + & cf2 * hat(i,2,j ) + & cf3 * hat(i,3,j ) + & cf1 * hat(i,1,j+1) + & cf2 * hat(i,2,j+1) + & cf3 * hat(i,3,j+1) ) hatavg(i,kte,j) = 0.5 * ( & cft1 * ( hat(i,ktes1,j) + hat(i,ktes1,j+1) ) + & cft2 * ( hat(i,ktes2,j) + hat(i,ktes2,j+1) ) ) END DO END DO ! Comments 10-MAR-05 ! Eqn 13b: D22=defor22= 2m^2 * (partial dv^/dY + partial dpsi/dy * partial dv^/dpsi) ! Below, D22 is set = 2*tmp1 ! => tmp1 = m^2 * (partial dv^/dY + partial dpsi/dy * partial dv^/dpsi) ! tmpzy = averaged value of dpsi/dy (=zy) DO j = j_start, j_end DO k = kts, ktf DO i = i_start, i_end tmpzy = 0.25 * ( & zy(i,k ,j) + zy(i,k ,j+1) + & zy(i,k+1,j) + zy(i,k+1,j+1) ) tmp1(i,k,j) = ( hatavg(i,k+1,j) - hatavg(i,k,j) ) * tmpzy * rdzw(i,k,j) ! tmp1 to here = partial dpsi/dy * partial dv^/dpsi: END DO END DO END DO DO j = j_start, j_end DO k = kts, ktf DO i = i_start, i_end tmp1(i,k,j) = mm(i,j) * ( & rdy * ( hat(i,k,j+1) - hat(i,k,j) ) - tmp1(i,k,j) ) END DO END DO END DO ! End calculation of dv/dy. !----------------------------------------------------------------------- !----------------------------------------------------------------------- ! Calculate defor22 (2*dv/dy). ! Comments 10-MAR-05 ! Eqn 13b: D22=defor22= 2 m^2 * (partial dv^/dY + partial dpsi/dy * partial dv^/dpsi) ! = 2*tmp1 DO j = j_start, j_end DO k = kts, ktf DO i = i_start, i_end defor22(i,k,j) = 2.0 * tmp1(i,k,j) END DO END DO END DO ! End calculation of defor22. !----------------------------------------------------------------------- !----------------------------------------------------------------------- ! Calculate meridional divergence (dv/dy) and add it to the divergence ! array. DO j = j_start, j_end DO k = kts, ktf DO i = i_start, i_end div(i,k,j) = div(i,k,j) + tmp1(i,k,j) END DO END DO END DO ! End calculation of meridional divergence. !----------------------------------------------------------------------- !----------------------------------------------------------------------- ! Comments 10-MAR-05 ! Eqn 13c: D33=defor33= 2 * partial dw/dz ! Below, D33 is set = 2*tmp1 ! => tmp1 = partial dw/dz ! Calculate dw/dz. DO j = j_start, j_end DO k = kts, ktf DO i = i_start, i_end tmp1(i,k,j) = ( w(i,k+1,j) - w(i,k,j) ) * rdzw(i,k,j) END DO END DO END DO ! End calculation of dw/dz. !----------------------------------------------------------------------- !----------------------------------------------------------------------- ! Calculate defor33 (2*dw/dz). DO j = j_start, j_end DO k = kts, ktf DO i = i_start, i_end defor33(i,k,j) = 2.0 * tmp1(i,k,j) END DO END DO END DO ! End calculation of defor33. !----------------------------------------------------------------------- !----------------------------------------------------------------------- ! Calculate vertical divergence (dw/dz) and add it to the divergence ! array. DO j = j_start, j_end DO k = kts, ktf DO i = i_start, i_end div(i,k,j) = div(i,k,j) + tmp1(i,k,j) END DO END DO END DO ! End calculation of vertical divergence. !----------------------------------------------------------------------- ! Three-dimensional divergence is now finished and values are in array ! "div." Also, the first three (defor11, defor22, defor33) of six ! deformation terms are now calculated at pressure points. !======================================================================= ! Comments 10-MAR-2005 ! Treat all differentials as 'div-style' [or 'curl-style'], ! i.e., du/dY becomes (in map coordinate space) mx*my * d(u/mx)/dy, ! dv/dX becomes (in map coordinate space) mx*my * d(v/my)/dx, ! (see e.g. Haltiner and Williams p. 441) !======================================================================= ! Calculate the final three deformations (defor12, defor13, defor23) at ! vorticity points. i_start = its i_end = ite j_start = jts j_end = jte IF ( config_flags%open_xs .OR. config_flags%specified .OR. & config_flags%nested) i_start = MAX( ids+1, its ) IF ( config_flags%open_xe .OR. config_flags%specified .OR. & config_flags%nested) i_end = MIN( ide-1, ite ) IF ( config_flags%open_ys .OR. config_flags%specified .OR. & config_flags%nested) j_start = MAX( jds+1, jts ) IF ( config_flags%open_ye .OR. config_flags%specified .OR. & config_flags%nested) j_end = MIN( jde-1, jte ) IF ( config_flags%periodic_x ) i_start = its IF ( config_flags%periodic_x ) i_end = ite !----------------------------------------------------------------------- ! Calculate du/dy. ! First, calculate an average mapscale factor. ! Comments 10-MAR-05 ! du/dy => need u map scale factor in x (which is defined at u points) ! averaged over j and j-1 ! dv/dx => need v map scale factor in y (which is defined at v points) ! averaged over i and i-1 DO j = j_start, j_end DO i = i_start, i_end mm(i,j) = 0.25 * ( msfux(i,j-1) + msfux(i,j) ) * ( msfvy(i-1,j) + msfvy(i,j) ) END DO END DO ! Apply a coordinate transformation to zonal velocity, u. DO j =j_start-1, j_end DO k =kts, ktf DO i =i_start, i_end ! Fixes to set_physical_bc2/3d for polar boundary conditions ! remove issues with loop over j hat(i,k,j) = u(i,k,j) / msfux(i,j) END DO END DO END DO ! Average in y and z. DO j=j_start,j_end DO k=kts+1,ktf DO i=i_start,i_end hatavg(i,k,j) = 0.5 * ( & fnm(k) * ( hat(i,k ,j-1) + hat(i,k ,j) ) + & fnp(k) * ( hat(i,k-1,j-1) + hat(i,k-1,j) ) ) END DO END DO END DO ! Extrapolate to top and bottom of domain (to w levels). DO j = j_start, j_end DO i = i_start, i_end hatavg(i,1,j) = 0.5 * ( & cf1 * hat(i,1,j-1) + & cf2 * hat(i,2,j-1) + & cf3 * hat(i,3,j-1) + & cf1 * hat(i,1,j ) + & cf2 * hat(i,2,j ) + & cf3 * hat(i,3,j ) ) hatavg(i,kte,j) = 0.5 * ( & cft1 * ( hat(i,ktes1,j-1) + hat(i,ktes1,j) ) + & cft2 * ( hat(i,ktes2,j-1) + hat(i,ktes2,j) ) ) END DO END DO ! tmpzy = averaged value of dpsi/dy (=zy) on vorticity grid ! tmp1 = partial dpsi/dy * partial du^/dpsi DO j = j_start, j_end DO k = kts, ktf DO i = i_start, i_end tmpzy = 0.25 * ( & zy(i-1,k ,j) + zy(i,k ,j) + & zy(i-1,k+1,j) + zy(i,k+1,j) ) tmp1(i,k,j) = ( hatavg(i,k+1,j) - hatavg(i,k,j) ) * & 0.25 * tmpzy * ( rdzw(i,k,j) + rdzw(i-1,k,j) + & rdzw(i-1,k,j-1) + rdzw(i,k,j-1) ) END DO END DO END DO ! End calculation of du/dy. !---------------------------------------------------------------------- !----------------------------------------------------------------------- ! Add the first term to defor12 (du/dy+dv/dx) at vorticity points. ! Comments 10-MAR-05 ! Eqn 13d: D12=defor12= m^2 * (partial dv^/dX + partial du^/dY + ! partial dpsi/dx * partial dv^/dpsi + ! partial dpsi/dy * partial du^/dpsi) ! Here deal with m^2 * (partial du^/dY + partial dpsi/dy * partial du^/dpsi) ! Still need to add v^ terms: ! m^2 * (partial dv^/dX + partial dpsi/dx * partial dv^/dpsi) DO j = j_start, j_end DO k = kts, ktf DO i = i_start, i_end defor12(i,k,j) = mm(i,j) * ( & rdy * ( hat(i,k,j) - hat(i,k,j-1) ) - tmp1(i,k,j) ) END DO END DO END DO ! End addition of the first term to defor12. !----------------------------------------------------------------------- !----------------------------------------------------------------------- ! Calculate dv/dx. ! Apply a coordinate transformation to meridional velocity, v. DO j = j_start, j_end DO k = kts, ktf DO i = i_start-1, i_end hat(i,k,j) = v(i,k,j) / msfvy(i,j) END DO END DO END DO ! Account for the slope in x of eta surfaces. DO j = j_start, j_end DO k = kts+1, ktf DO i = i_start, i_end hatavg(i,k,j) = 0.5 * ( & fnm(k) * ( hat(i-1,k ,j) + hat(i,k ,j) ) + & fnp(k) * ( hat(i-1,k-1,j) + hat(i,k-1,j) ) ) END DO END DO END DO ! Extrapolate to top and bottom of domain (to w levels). DO j = j_start, j_end DO i = i_start, i_end hatavg(i,1,j) = 0.5 * ( & cf1 * hat(i-1,1,j) + & cf2 * hat(i-1,2,j) + & cf3 * hat(i-1,3,j) + & cf1 * hat(i ,1,j) + & cf2 * hat(i ,2,j) + & cf3 * hat(i ,3,j) ) hatavg(i,kte,j) = 0.5 * ( & cft1 * ( hat(i,ktes1,j) + hat(i-1,ktes1,j) ) + & cft2 * ( hat(i,ktes2,j) + hat(i-1,ktes2,j) ) ) END DO END DO ! Fixes to set_physical_bc2/3d have made any check for polar B.C.'s ! unnecessary in this place. zx, rdzw, and hatavg are all defined ! in places they need to be and the values at the poles are replications ! of the values one grid point in, so the averaging over j and j-1 works ! to act as just using the value at j or j-1 (with out extra code). ! ! tmpzx = averaged value of dpsi/dx (=zx) on vorticity grid ! tmp1 = partial dpsi/dx * partial dv^/dpsi DO j = j_start, j_end DO k = kts, ktf DO i = i_start, i_end tmpzx = 0.25 * ( & zx(i,k ,j-1) + zx(i,k ,j) + & zx(i,k+1,j-1) + zx(i,k+1,j) ) tmp1(i,k,j) = ( hatavg(i,k+1,j) - hatavg(i,k,j) ) * & 0.25 * tmpzx * ( rdzw(i,k,j) + rdzw(i,k,j-1) + & rdzw(i-1,k,j-1) + rdzw(i-1,k,j) ) END DO END DO END DO ! End calculation of dv/dx. !----------------------------------------------------------------------- !----------------------------------------------------------------------- ! Add the second term to defor12 (du/dy+dv/dx) at vorticity points. ! Comments 10-MAR-05 ! Eqn 13d: D12=defor12= m^2 * (partial dv^/dX + partial du^/dY + ! partial dpsi/dx * partial dv^/dpsi + ! partial dpsi/dy * partial du^/dpsi) ! Here adding v^ terms: ! m^2 * (partial dv^/dX + partial dpsi/dx * partial dv^/dpsi) DO j = j_start, j_end DO k = kts, ktf DO i = i_start, i_end defor12(i,k,j) = defor12(i,k,j) + & mm(i,j) * ( & rdx * ( hat(i,k,j) - hat(i-1,k,j) ) - tmp1(i,k,j) ) END DO END DO END DO ! End addition of the second term to defor12. !----------------------------------------------------------------------- !----------------------------------------------------------------------- ! Update the boundary for defor12 (might need to change later). IF ( .NOT. config_flags%periodic_x .AND. i_start .EQ. ids+1 ) THEN DO j = jts, jte DO k = kts, kte defor12(ids,k,j) = defor12(ids+1,k,j) END DO END DO END IF IF ( .NOT. config_flags%periodic_y .AND. j_start .EQ. jds+1) THEN DO k = kts, kte DO i = its, ite defor12(i,k,jds) = defor12(i,k,jds+1) END DO END DO END IF IF ( .NOT. config_flags%periodic_x .AND. i_end .EQ. ide-1) THEN DO j = jts, jte DO k = kts, kte defor12(ide,k,j) = defor12(ide-1,k,j) END DO END DO END IF IF ( .NOT. config_flags%periodic_y .AND. j_end .EQ. jde-1) THEN DO k = kts, kte DO i = its, ite defor12(i,k,jde) = defor12(i,k,jde-1) END DO END DO END IF ! End update of boundary for defor12. !----------------------------------------------------------------------- ! Comments 10-MAR-05 ! Further deformation terms not needed for 2-dimensional Smagorinsky diffusion, ! so those terms have not been dealt with yet. ! A "y" has simply been added to all map scale factors to allow the model to ! compile without errors. !----------------------------------------------------------------------- ! Calculate dw/dx. i_start = its i_end = MIN( ite, ide-1 ) j_start = jts j_end = MIN( jte, jde-1 ) IF ( config_flags%open_xs .OR. config_flags%specified .OR. & config_flags%nested) i_start = MAX( ids+1, its ) IF ( config_flags%open_ys .OR. config_flags%specified .OR. & config_flags%nested) j_start = MAX( jds+1, jts ) IF ( config_flags%periodic_x ) i_start = its IF ( config_flags%periodic_x ) i_end = MIN( ite, ide ) IF ( config_flags%periodic_y ) j_end = MIN( jte, jde ) ! Square the mapscale factor. DO j = jts, jte DO i = its, ite mm(i,j) = msfux(i,j) * msfuy(i,j) END DO END DO ! Apply a coordinate transformation to vertical velocity, w. This is for both ! defor13 and defor23. DO j = j_start, j_end DO k = kts, kte DO i = i_start, i_end hat(i,k,j) = w(i,k,j) / msfty(i,j) END DO END DO END DO i = i_start-1 DO j = j_start, MIN( jte, jde-1 ) DO k = kts, kte hat(i,k,j) = w(i,k,j) / msfty(i,j) END DO END DO j = j_start-1 DO k = kts, kte DO i = i_start, MIN( ite, ide-1 ) hat(i,k,j) = w(i,k,j) / msfty(i,j) END DO END DO ! QUESTION: What is this for? DO j = j_start, j_end DO k = kts, ktf DO i = i_start, i_end hatavg(i,k,j) = 0.25 * ( & hat(i ,k ,j) + & hat(i ,k+1,j) + & hat(i-1,k ,j) + & hat(i-1,k+1,j) ) END DO END DO END DO ! Calculate dw/dx. DO j = j_start, j_end DO k = kts+1, ktf DO i = i_start, i_end tmp1(i,k,j) = ( hatavg(i,k,j) - hatavg(i,k-1,j) ) * zx(i,k,j) * & 0.5 * ( rdz(i,k,j) + rdz(i-1,k,j) ) END DO END DO END DO ! End calculation of dw/dx. !----------------------------------------------------------------------- !----------------------------------------------------------------------- ! Add the first term (dw/dx) to defor13 (dw/dx+du/dz) at vorticity ! points. DO j = j_start, j_end DO k = kts+1, ktf DO i = i_start, i_end defor13(i,k,j) = mm(i,j) * ( & rdx * ( hat(i,k,j) - hat(i-1,k,j) ) - tmp1(i,k,j) ) END DO END DO END DO DO j = j_start, j_end DO i = i_start, i_end defor13(i,kts,j ) = 0.0 defor13(i,ktf+1,j) = 0.0 END DO END DO ! End addition of the first term to defor13. !----------------------------------------------------------------------- !----------------------------------------------------------------------- ! Calculate du/dz. IF ( config_flags%mix_full_fields ) THEN DO j = j_start, j_end DO k = kts+1, ktf DO i = i_start, i_end tmp1(i,k,j) = ( u(i,k,j) - u(i,k-1,j) ) * & 0.5 * ( rdz(i,k,j) + rdz(i-1,k,j) ) END DO END DO END DO ELSE DO j = j_start, j_end DO k = kts+1, ktf DO i = i_start, i_end tmp1(i,k,j) = ( u(i,k,j) - u_base(k) - u(i,k-1,j) + u_base(k-1) ) * & 0.5 * ( rdz(i,k,j) + rdz(i-1,k,j) ) END DO END DO END DO END IF !----------------------------------------------------------------------- ! Add the second term (du/dz) to defor13 (dw/dx+du/dz) at vorticity ! points. DO j = j_start, j_end DO k = kts+1, ktf DO i = i_start, i_end defor13(i,k,j) = defor13(i,k,j) + tmp1(i,k,j) END DO END DO END DO ! End addition of the second term to defor13. !----------------------------------------------------------------------- !----------------------------------------------------------------------- ! Calculate dw/dy. i_start = its i_end = MIN( ite, ide-1 ) j_start = jts j_end = MIN( jte, jde-1 ) IF ( config_flags%open_xs .OR. config_flags%specified .OR. & config_flags%nested) i_start = MAX( ids+1, its ) IF ( config_flags%open_ys .OR. config_flags%specified .OR. & config_flags%nested) j_start = MAX( jds+1, jts ) IF ( config_flags%periodic_y ) j_end = MIN( jte, jde ) IF ( config_flags%periodic_x ) i_start = its IF ( config_flags%periodic_x ) i_end = MIN( ite, ide-1 ) ! Square mapscale factor. DO j = jts, jte DO i = its, ite mm(i,j) = msfvx(i,j) * msfvy(i,j) END DO END DO ! Apply a coordinate transformation to vertical velocity, w. Added by CW 7/19/07 DO j = j_start, j_end DO k = kts, kte DO i = i_start, i_end hat(i,k,j) = w(i,k,j) / msftx(i,j) END DO END DO END DO i = i_start-1 DO j = j_start, MIN( jte, jde-1 ) DO k = kts, kte hat(i,k,j) = w(i,k,j) / msftx(i,j) END DO END DO j = j_start-1 DO k = kts, kte DO i = i_start, MIN( ite, ide-1 ) hat(i,k,j) = w(i,k,j) / msftx(i,j) END DO END DO ! QUESTION: What is this for? DO j = j_start, j_end DO k = kts, ktf DO i = i_start, i_end hatavg(i,k,j) = 0.25 * ( & hat(i,k ,j ) + & hat(i,k+1,j ) + & hat(i,k ,j-1) + & hat(i,k+1,j-1) ) END DO END DO END DO ! Calculate dw/dy and store in tmp1. DO j = j_start, j_end DO k = kts+1, ktf DO i = i_start, i_end tmp1(i,k,j) = ( hatavg(i,k,j) - hatavg(i,k-1,j) ) * zy(i,k,j) * & 0.5 * ( rdz(i,k,j) + rdz(i,k,j-1) ) END DO END DO END DO ! End calculation of dw/dy. !----------------------------------------------------------------------- !----------------------------------------------------------------------- ! Add the first term (dw/dy) to defor23 (dw/dy+dv/dz) at vorticity ! points. DO j = j_start, j_end DO k = kts+1, ktf DO i = i_start, i_end defor23(i,k,j) = mm(i,j) * ( & rdy * ( hat(i,k,j) - hat(i,k,j-1) ) - tmp1(i,k,j) ) END DO END DO END DO DO j = j_start, j_end DO i = i_start, i_end defor23(i,kts,j ) = 0.0 defor23(i,ktf+1,j) = 0.0 END DO END DO ! End addition of the first term to defor23. !----------------------------------------------------------------------- !----------------------------------------------------------------------- ! Calculate dv/dz. IF ( config_flags%mix_full_fields ) THEN DO j = j_start, j_end DO k = kts+1, ktf DO i = i_start, i_end tmp1(i,k,j) = ( v(i,k,j) - v(i,k-1,j) ) * & 0.5 * ( rdz(i,k,j) + rdz(i,k,j-1) ) END DO END DO END DO ELSE DO j = j_start, j_end DO k = kts+1, ktf DO i = i_start, i_end tmp1(i,k,j) = ( v(i,k,j) - v_base(k) - v(i,k-1,j) + v_base(k-1) ) * & 0.5 * ( rdz(i,k,j) + rdz(i,k,j-1) ) END DO END DO END DO END IF ! End calculation of dv/dz. !----------------------------------------------------------------------- !----------------------------------------------------------------------- ! Add the second term (dv/dz) to defor23 (dw/dy+dv/dz) at vorticity ! points. ! Add tmp1 to defor23. DO j = j_start, j_end DO k = kts+1, ktf DO i = i_start, i_end defor23(i,k,j) = defor23(i,k,j) + tmp1(i,k,j) END DO END DO END DO ! End addition of the second term to defor23. !----------------------------------------------------------------------- !----------------------------------------------------------------------- ! Update the boundary for defor13 and defor23 (might need to change ! later). IF ( .NOT. config_flags%periodic_x .AND. i_start .EQ. ids+1) THEN DO j = jts, jte DO k = kts, kte defor13(ids,k,j) = defor13(ids+1,k,j) defor23(ids,k,j) = defor23(ids+1,k,j) END DO END DO END IF IF ( .NOT. config_flags%periodic_y .AND. j_start .EQ. jds+1) THEN DO k = kts, kte DO i = its, ite defor13(i,k,jds) = defor13(i,k,jds+1) defor23(i,k,jds) = defor23(i,k,jds+1) END DO END DO END IF IF ( .NOT. config_flags%periodic_x .AND. i_end .EQ. ide-1) THEN DO j = jts, jte DO k = kts, kte defor13(ide,k,j) = defor13(ide-1,k,j) defor23(ide,k,j) = defor23(ide-1,k,j) END DO END DO END IF IF ( .NOT. config_flags%periodic_y .AND. j_end .EQ. jde-1) THEN DO k = kts, kte DO i = its, ite defor13(i,k,jde) = defor13(i,k,jde-1) defor23(i,k,jde) = defor23(i,k,jde-1) END DO END DO END IF ! End update of boundary for defor13 and defor23. !----------------------------------------------------------------------- ! The second three (defor12, defor13, defor23) of six deformation terms ! are now calculated at vorticity points. !======================================================================= END SUBROUTINE cal_deform_and_div !======================================================================= !======================================================================= SUBROUTINE calculate_km_kh( config_flags, dt, & dampcoef, zdamp, damp_opt, & xkmh, xkmv, xkhh, xkhv, & BN2, khdif, kvdif, div, & defor11, defor22, defor33, & defor12, defor13, defor23, & tke, p8w, t8w, theta, t, p, moist, & dn, dnw, dx, dy, rdz, rdzw, isotropic, & n_moist, cf1, cf2, cf3, warm_rain, & mix_upper_bound, & msftx, msfty, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) ! History: Sep 2003 Changes by George Bryan and Jason Knievel, NCAR ! Oct 2001 Converted to mass core by Bill Skamarock, NCAR ! ... ... ! Purpose: This routine calculates exchange coefficients for the TKE ! scheme. ! References: Klemp and Wilhelmson (JAS 1978) ! Deardorff (B-L Meteor 1980) ! Chen and Dudhia (NCAR WRF physics report 2000) !----------------------------------------------------------------------- ! Begin declarations. IMPLICIT NONE TYPE( grid_config_rec_type ), INTENT( IN ) & :: config_flags INTEGER, INTENT( IN ) & :: n_moist, damp_opt, isotropic, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte LOGICAL, INTENT( IN ) & :: warm_rain REAL, INTENT( IN ) & :: dx, dy, zdamp, dt, dampcoef, cf1, cf2, cf3, khdif, kvdif REAL, DIMENSION( kms:kme ), INTENT( IN ) & :: dnw, dn REAL, DIMENSION( ims:ime, kms:kme, jms:jme, n_moist ), INTENT( INOUT ) & :: moist REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( INOUT ) & :: xkmv, xkmh, xkhv, xkhh, BN2 REAL, DIMENSION( ims:ime , kms:kme, jms:jme ), INTENT( IN ) & :: defor11, defor22, defor33, defor12, defor13, defor23, & div, rdz, rdzw, p8w, t8w, theta, t, p REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( INOUT ) & :: tke REAL, INTENT( IN ) & :: mix_upper_bound REAL, DIMENSION( ims:ime, jms:jme ), INTENT( IN ) & :: msftx, msfty ! Local variables. INTEGER & :: i_start, i_end, j_start, j_end, ktf, i, j, k ! End declarations. !----------------------------------------------------------------------- ktf = MIN( kte, kde-1 ) i_start = its i_end = MIN( ite, ide-1 ) j_start = jts j_end = MIN( jte, jde-1 ) CALL calculate_N2( config_flags, BN2, moist, & theta, t, p, p8w, t8w, & dnw, dn, rdz, rdzw, & n_moist, cf1, cf2, cf3, warm_rain, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) ! Select a scheme for calculating diffusion coefficients. km_coef: SELECT CASE( config_flags%km_opt ) CASE (1) CALL isotropic_km( config_flags, xkmh, xkmv, & xkhh, xkhv, khdif, kvdif, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) CASE (2) CALL tke_km( config_flags, xkmh, xkmv, & xkhh, xkhv, BN2, tke, p8w, t8w, theta, & rdz, rdzw, dx, dy, dt, isotropic, & mix_upper_bound, msftx, msfty, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) CASE (3) CALL smag_km( config_flags, xkmh, xkmv, & xkhh, xkhv, BN2, div, & defor11, defor22, defor33, & defor12, defor13, defor23, & rdzw, dx, dy, dt, isotropic, & mix_upper_bound, msftx, msfty, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) CASE (4) CALL smag2d_km( config_flags, xkmh, xkmv, & xkhh, xkhv, defor11, defor22, defor12, & rdzw, dx, dy, msftx, msfty, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) CASE DEFAULT CALL wrf_error_fatal( 'Please choose diffusion coefficient scheme' ) END SELECT km_coef IF ( damp_opt .eq. 1 ) THEN CALL cal_dampkm( config_flags, xkmh, xkhh, xkmv, xkhv, & dx, dy, dt, dampcoef, rdz, rdzw, zdamp, & msftx, msfty, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) END IF END SUBROUTINE calculate_km_kh !======================================================================= SUBROUTINE cal_dampkm( config_flags,xkmh,xkhh,xkmv,xkhv, & dx,dy,dt,dampcoef, & rdz, rdzw ,zdamp, & msftx, msfty, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte ) !----------------------------------------------------------------------- ! Begin declarations. IMPLICIT NONE 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 , INTENT(IN ) :: zdamp,dx,dy,dt,dampcoef REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: xkmh , & xkhh , & xkmv , & xkhv REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(IN ) :: rdz, & rdzw REAL , DIMENSION( ims:ime, jms:jme), INTENT(IN ) :: msftx, & msfty ! LOCAL VARS INTEGER :: i_start, i_end, j_start, j_end, ktf, ktfm1, i, j, k REAL :: kmmax,kmmvmax,degrad90,dz,tmp REAL :: ds REAL , DIMENSION( its:ite ) :: deltaz REAL , DIMENSION( its:ite, kts:kte, jts:jte) :: dampk,dampkv ! End declarations. !----------------------------------------------------------------------- ktf = min(kte,kde-1) ktfm1 = ktf-1 i_start = its i_end = MIN(ite,ide-1) j_start = jts j_end = MIN(jte,jde-1) ! keep upper damping diffusion away from relaxation zones at boundaries if used IF(config_flags%specified .OR. config_flags%nested)THEN i_start = MAX(i_start,ids+config_flags%spec_bdy_width-1) i_end = MIN(i_end,ide-config_flags%spec_bdy_width) j_start = MAX(j_start,jds+config_flags%spec_bdy_width-1) j_end = MIN(j_end,jde-config_flags%spec_bdy_width) ENDIF kmmax=dx*dx/dt degrad90=DEGRAD*90. DO j = j_start, j_end k=ktf DO i = i_start, i_end ! Unmodified dx used above may produce very large diffusivities ! when msftx is very large. And the above formula ignores the fact ! that dy may now be different from dx as well. Let's fix that by ! defining a "ds" as the minimum of the "real-space" (physical ! distance) values of dx and dy, and then using that smallest value ! to calculate a point-by-point kmmax ds = MIN(dx/msftx(i,j),dy/msfty(i,j)) kmmax=ds*ds/dt ! deltaz(i)=0.5*dnw(k)/zeta_z(i,j) ! dz=dnw(k)/zeta_z(i,j) dz = 1./rdzw(i,k,j) deltaz(i) = 0.5*dz kmmvmax=dz*dz/dt tmp=min(deltaz(i)/zdamp,1.) dampk(i,k,j)=cos(degrad90*tmp)*cos(degrad90*tmp)*kmmax*dampcoef dampkv(i,k,j)=cos(degrad90*tmp)*cos(degrad90*tmp)*kmmvmax*dampcoef ! set upper limit on vertical K (based on horizontal K) dampkv(i,k,j)=min(dampkv(i,k,j),dampk(i,k,j)) ENDDO DO k = ktfm1,kts,-1 DO i = i_start, i_end ! Unmodified dx used above may produce very large diffusivities ! when msftx is very large. And the above formula ignores the fact ! that dy may now be different from dx as well. Let's fix that by ! defining a "ds" as the minimum of the "real-space" (physical ! distance) values of dx and dy, and then using that smallest value ! to calculate a point-by-point kmmax ds = MIN(dx/msftx(i,j),dy/msfty(i,j)) kmmax=ds*ds/dt ! deltaz(i)=deltaz(i)+dn(k)/zeta_z(i,j) ! dz=dnw(k)/zeta_z(i,j) dz = 1./rdz(i,k,j) deltaz(i) = deltaz(i) + dz dz = 1./rdzw(i,k,j) kmmvmax=dz*dz/dt tmp=min(deltaz(i)/zdamp,1.) dampk(i,k,j)=cos(degrad90*tmp)*cos(degrad90*tmp)*kmmax*dampcoef dampkv(i,k,j)=cos(degrad90*tmp)*cos(degrad90*tmp)*kmmvmax*dampcoef ! set upper limit on vertical K (based on horizontal K) dampkv(i,k,j)=min(dampkv(i,k,j),dampk(i,k,j)) ENDDO ENDDO ENDDO DO j = j_start, j_end DO k = kts,ktf DO i = i_start, i_end xkmh(i,k,j)=max(xkmh(i,k,j),dampk(i,k,j)) xkhh(i,k,j)=max(xkhh(i,k,j),dampk(i,k,j)) xkmv(i,k,j)=max(xkmv(i,k,j),dampkv(i,k,j)) xkhv(i,k,j)=max(xkhv(i,k,j),dampkv(i,k,j)) ENDDO ENDDO ENDDO END SUBROUTINE cal_dampkm !======================================================================= !======================================================================= SUBROUTINE calculate_N2( config_flags, BN2, moist, & theta, t, p, p8w, t8w, & dnw, dn, rdz, rdzw, & n_moist, cf1, cf2, cf3, warm_rain, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) !----------------------------------------------------------------------- ! Begin declarations. IMPLICIT NONE TYPE( grid_config_rec_type ), INTENT( IN ) & :: config_flags INTEGER, INTENT( IN ) & :: n_moist, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte LOGICAL, INTENT( IN ) & :: warm_rain REAL, INTENT( IN ) & :: cf1, cf2, cf3 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( INOUT ) & :: BN2 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( IN ) & :: rdz, rdzw, theta, t, p, p8w, t8w REAL, DIMENSION( kms:kme ), INTENT( IN ) & :: dnw, dn REAL, DIMENSION( ims:ime, kms:kme, jms:jme, n_moist), INTENT( INOUT ) & :: moist ! Local variables. INTEGER & :: i, j, k, ktf, ispe, ktes1, ktes2, & i_start, i_end, j_start, j_end REAL & :: coefa, thetaep1, thetaem1, qc_cr, es, tc, qlpqi, qsw, qsi, & tmpdz, xlvqv, thetaesfc, thetasfc, qvtop, qvsfc, thetatop, thetaetop REAL, DIMENSION( its:ite, jts:jte ) & :: tmp1sfc, tmp1top REAL, DIMENSION( its:ite, kts:kte, jts:jte ) & :: tmp1, qvs, qctmp ! End declarations. !----------------------------------------------------------------------- qc_cr = 0.00001 ! in Kg/Kg ktf = MIN( kte, kde-1 ) ktes1 = kte-1 ktes2 = kte-2 i_start = its i_end = MIN( ite, ide-1 ) j_start = jts j_end = MIN( jte, jde-1 ) IF ( config_flags%open_xs .OR. config_flags%specified .OR. & config_flags%nested) i_start = MAX( ids+1, its ) IF ( config_flags%open_xe .OR. config_flags%specified .OR. & config_flags%nested) i_end = MIN( ide-2, ite ) IF ( config_flags%open_ys .OR. config_flags%specified .OR. & config_flags%nested) j_start = MAX( jds+1, jts ) IF ( config_flags%open_ye .OR. config_flags%specified .OR. & config_flags%nested) j_end = MIN( jde-2 ,jte ) IF ( config_flags%periodic_x ) i_start = its IF ( config_flags%periodic_x ) i_end = MIN( ite, ide-1 ) IF ( P_QC .GT. PARAM_FIRST_SCALAR) THEN DO j = j_start, j_end DO k = kts, ktf DO i = i_start, i_end qctmp(i,k,j) = moist(i,k,j,P_QC) END DO END DO END DO ELSE DO j = j_start, j_end DO k = kts, ktf DO i = i_start, i_end qctmp(i,k,j) = 0.0 END DO END DO END DO END IF DO j = jts, jte DO k = kts, kte DO i = its, ite tmp1(i,k,j) = 0.0 END DO END DO END DO DO j = jts,jte DO i = its,ite tmp1sfc(i,j) = 0.0 tmp1top(i,j) = 0.0 END DO END DO DO ispe = PARAM_FIRST_SCALAR, n_moist IF ( ispe .EQ. P_QV .OR. ispe .EQ. P_QC .OR. ispe .EQ. P_QI) THEN DO j = j_start, j_end DO k = kts, ktf DO i = i_start, i_end tmp1(i,k,j) = tmp1(i,k,j) + moist(i,k,j,ispe) END DO END DO END DO DO j = j_start, j_end DO i = i_start, i_end tmp1sfc(i,j) = tmp1sfc(i,j) + & cf1 * moist(i,1,j,ispe) + & cf2 * moist(i,2,j,ispe) + & cf3 * moist(i,3,j,ispe) tmp1top(i,j) = tmp1top(i,j) + & moist(i,ktes1,j,ispe) + & ( moist(i,ktes1,j,ispe) - moist(i,ktes2,j,ispe) ) * & 0.5 * dnw(ktes1) / dn(ktes1) END DO END DO END IF END DO ! Calculate saturation mixing ratio. DO j = j_start, j_end DO k = kts, ktf DO i = i_start, i_end tc = t(i,k,j) - SVPT0 es = 1000.0 * SVP1 * EXP( SVP2 * tc / ( t(i,k,j) - SVP3 ) ) qvs(i,k,j) = EP_2 * es / ( p(i,k,j) - es ) END DO END DO END DO DO j = j_start, j_end DO k = kts+1, ktf-1 DO i = i_start, i_end tmpdz = 1.0 / rdz(i,k,j) + 1.0 / rdz(i,k+1,j) IF ( moist(i,k,j,P_QV) .GE. qvs(i,k,j) .OR. qctmp(i,k,j) .GE. qc_cr) THEN xlvqv = XLV * moist(i,k,j,P_QV) coefa = ( 1.0 + xlvqv / R_d / t(i,k,j) ) / & ( 1.0 + XLV * xlvqv / Cp / R_v / t(i,k,j) / t(i,k,j) ) / & theta(i,k,j) thetaep1 = theta(i,k+1,j) * & ( 1.0 + XLV * qvs(i,k+1,j) / Cp / t(i,k+1,j) ) thetaem1 = theta(i,k-1,j) * & ( 1.0 + XLV * qvs(i,k-1,j) / Cp / t(i,k-1,j) ) BN2(i,k,j) = g * ( coefa * ( thetaep1 - thetaem1 ) / tmpdz - & ( tmp1(i,k+1,j) - tmp1(i,k-1,j) ) / tmpdz ) ELSE BN2(i,k,j) = g * ( (theta(i,k+1,j) - theta(i,k-1,j) ) / & theta(i,k,j) / tmpdz + & 1.61 * ( moist(i,k+1,j,P_QV) - moist(i,k-1,j,P_QV) ) / & tmpdz - & ( tmp1(i,k+1,j) - tmp1(i,k-1,j) ) / tmpdz ) ENDIF END DO END DO END DO k = kts DO j = j_start, j_end DO i = i_start, i_end tmpdz = 1.0 / rdz(i,k+1,j) + 0.5 / rdzw(i,k,j) thetasfc = T8w(i,kts,j) / ( p8w(i,k,j) / p1000mb )**( R_d / Cp ) IF ( moist(i,k,j,P_QV) .GE. qvs(i,k,j) .OR. qctmp(i,k,j) .GE. qc_cr) THEN qvsfc = cf1 * qvs(i,1,j) + & cf2 * qvs(i,2,j) + & cf3 * qvs(i,3,j) xlvqv = XLV * moist(i,k,j,P_QV) coefa = ( 1.0 + xlvqv / R_d / t(i,k,j) ) / & ( 1.0 + XLV * xlvqv / Cp / R_v / t(i,k,j) / t(i,k,j) ) / & theta(i,k,j) thetaep1 = theta(i,k+1,j) * & ( 1.0 + XLV * qvs(i,k+1,j) / Cp / t(i,k+1,j) ) thetaesfc = thetasfc * & ( 1.0 + XLV * qvsfc / Cp / t8w(i,kts,j) ) BN2(i,k,j) = g * ( coefa * ( thetaep1 - thetaesfc ) / tmpdz - & ( tmp1(i,k+1,j) - tmp1sfc(i,j) ) / tmpdz ) ELSE qvsfc = cf1 * moist(i,1,j,P_QV) + & cf2 * moist(i,2,j,P_QV) + & cf3 * moist(i,3,j,P_QV) ! BN2(i,k,j) = g * ( ( theta(i,k+1,j) - thetasfc ) / & ! theta(i,k,j) / tmpdz + & ! 1.61 * ( moist(i,k+1,j,P_QV) - qvsfc ) / & ! tmpdz - & ! ( tmp1(i,k+1,j) - tmp1sfc(i,j) ) / tmpdz ) !...... MARTA: change in computation of BN2 at the surface, WCS 040331 tmpdz= 1./rdzw(i,k,j) ! controlare come calcola rdzw BN2(i,k,j) = g * ( ( theta(i,k+1,j) - theta(i,k,j)) / & theta(i,k,j) / tmpdz + & 1.61 * ( moist(i,k+1,j,P_QV) - qvsfc ) / & tmpdz - & ( tmp1(i,k+1,j) - tmp1sfc(i,j) ) / tmpdz ) ! end of MARTA/WCS change ENDIF END DO END DO !...... MARTA: change in computation of BN2 at the top, WCS 040331 DO j = j_start, j_end DO i = i_start, i_end BN2(i,ktf,j)=BN2(i,ktf-1,j) END DO END DO ! end of MARTA/WCS change END SUBROUTINE calculate_N2 !======================================================================= !======================================================================= SUBROUTINE isotropic_km( config_flags, & xkmh,xkmv,xkhh,xkhv,khdif,kvdif, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte ) !----------------------------------------------------------------------- ! Begin declarations. IMPLICIT NONE 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 , INTENT(IN ) :: khdif,kvdif REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: xkmh, & xkmv, & xkhh, & xkhv ! LOCAL VARS INTEGER :: i_start, i_end, j_start, j_end, ktf, i, j, k REAL :: khdif3,kvdif3 ! End declarations. !----------------------------------------------------------------------- ktf = kte i_start = its i_end = MIN(ite,ide-1) j_start = jts j_end = MIN(jte,jde-1) ! khdif3=khdif*3. ! kvdif3=kvdif*3. khdif3=khdif/prandtl kvdif3=kvdif/prandtl DO j = j_start, j_end DO k = kts, ktf DO i = i_start, i_end xkmh(i,k,j)=khdif xkmv(i,k,j)=kvdif xkhh(i,k,j)=khdif3 xkhv(i,k,j)=kvdif3 ENDDO ENDDO ENDDO END SUBROUTINE isotropic_km !======================================================================= !======================================================================= SUBROUTINE smag_km( config_flags,xkmh,xkmv,xkhh,xkhv,BN2, & div,defor11,defor22,defor33,defor12, & defor13,defor23, & rdzw,dx,dy,dt,isotropic, & mix_upper_bound, msftx, msfty, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte ) !----------------------------------------------------------------------- ! Begin declarations. IMPLICIT NONE 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 INTEGER , INTENT(IN ) :: isotropic REAL , INTENT(IN ) :: dx, dy, dt, mix_upper_bound REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: BN2, & rdzw REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: xkmh, & xkmv, & xkhh, & xkhv REAL , DIMENSION( ims:ime , kms:kme, jms:jme ), INTENT(IN ) :: & defor11, & defor22, & defor33, & defor12, & defor13, & defor23, & div REAL , DIMENSION( ims:ime, jms:jme), INTENT(IN ) :: msftx, & msfty ! LOCAL VARS INTEGER :: i_start, i_end, j_start, j_end, ktf, i, j, k REAL :: deltas, tmp, pr, mlen_h, mlen_v, c_s REAL, DIMENSION( its:ite , kts:kte , jts:jte ) :: def2 ! End declarations. !----------------------------------------------------------------------- ktf = min(kte,kde-1) i_start = its i_end = MIN(ite,ide-1) j_start = jts j_end = MIN(jte,jde-1) IF ( config_flags%open_xs .or. config_flags%specified .or. & config_flags%nested) i_start = MAX(ids+1,its) IF ( config_flags%open_xe .or. config_flags%specified .or. & config_flags%nested) i_end = MIN(ide-2,ite) IF ( config_flags%open_ys .or. config_flags%specified .or. & config_flags%nested) j_start = MAX(jds+1,jts) IF ( config_flags%open_ye .or. config_flags%specified .or. & config_flags%nested) j_end = MIN(jde-2,jte) IF ( config_flags%periodic_x ) i_start = its IF ( config_flags%periodic_x ) i_end = MIN( ite, ide-1 ) pr = prandtl c_s = config_flags%c_s do j=j_start,j_end do k=kts,ktf do i=i_start,i_end def2(i,k,j)=0.5*(defor11(i,k,j)*defor11(i,k,j) + & defor22(i,k,j)*defor22(i,k,j) + & defor33(i,k,j)*defor33(i,k,j)) enddo enddo enddo do j=j_start,j_end do k=kts,ktf do i=i_start,i_end tmp=0.25*(defor12(i ,k,j)+defor12(i ,k,j+1)+ & defor12(i+1,k,j)+defor12(i+1,k,j+1)) def2(i,k,j)=def2(i,k,j)+tmp*tmp enddo enddo enddo do j=j_start,j_end do k=kts,ktf do i=i_start,i_end tmp=0.25*(defor13(i ,k+1,j)+defor13(i ,k,j)+ & defor13(i+1,k+1,j)+defor13(i+1,k,j)) def2(i,k,j)=def2(i,k,j)+tmp*tmp enddo enddo enddo do j=j_start,j_end do k=kts,ktf do i=i_start,i_end tmp=0.25*(defor23(i,k+1,j )+defor23(i,k,j )+ & defor23(i,k+1,j+1)+defor23(i,k,j+1)) def2(i,k,j)=def2(i,k,j)+tmp*tmp enddo enddo enddo ! IF (isotropic .EQ. 0) THEN DO j = j_start, j_end DO k = kts, ktf DO i = i_start, i_end mlen_h=sqrt(dx/msftx(i,j) * dy/msfty(i,j)) mlen_v= 1./rdzw(i,k,j) tmp=max(0.,def2(i,k,j)-BN2(i,k,j)/pr) tmp=tmp**0.5 xkmh(i,k,j)=max(c_s*c_s*mlen_h*mlen_h*tmp, 1.0E-6*mlen_h*mlen_h ) xkmh(i,k,j)=min(xkmh(i,k,j), mix_upper_bound * mlen_h * mlen_h / dt ) xkmv(i,k,j)=max(c_s*c_s*mlen_v*mlen_v*tmp, 1.0E-6*mlen_v*mlen_v ) xkmv(i,k,j)=min(xkmv(i,k,j), mix_upper_bound * mlen_v * mlen_v / dt ) xkhh(i,k,j)=xkmh(i,k,j)/pr xkhh(i,k,j)=min(xkhh(i,k,j), mix_upper_bound * mlen_h * mlen_h / dt ) xkhv(i,k,j)=xkmv(i,k,j)/pr xkhv(i,k,j)=min(xkhv(i,k,j), mix_upper_bound * mlen_v * mlen_v / dt ) ENDDO ENDDO ENDDO ELSE DO j = j_start, j_end DO k = kts, ktf DO i = i_start, i_end deltas=(dx/msftx(i,j) * dy/msfty(i,j)/rdzw(i,k,j))**0.33333333 tmp=max(0.,def2(i,k,j)-BN2(i,k,j)/pr) tmp=tmp**0.5 xkmh(i,k,j)=max(c_s*c_s*deltas*deltas*tmp, 1.0E-6*deltas*deltas ) xkmh(i,k,j)=min(xkmh(i,k,j), mix_upper_bound * dx/msftx(i,j) * dy/msfty(i,j) / dt ) xkmv(i,k,j)=xkmh(i,k,j) xkmv(i,k,j)=min(xkmv(i,k,j), mix_upper_bound / rdzw(i,k,j) / rdzw(i,k,j) / dt ) xkhh(i,k,j)=xkmh(i,k,j)/pr xkhh(i,k,j)=min(xkhh(i,k,j), mix_upper_bound * dx/msftx(i,j) * dy/msfty(i,j) / dt ) xkhv(i,k,j)=xkmv(i,k,j)/pr xkhv(i,k,j)=min(xkhv(i,k,j), mix_upper_bound / rdzw(i,k,j) / rdzw(i,k,j) / dt ) ENDDO ENDDO ENDDO ENDIF END SUBROUTINE smag_km !======================================================================= !======================================================================= SUBROUTINE smag2d_km( config_flags,xkmh,xkmv,xkhh,xkhv, & defor11,defor22,defor12, & rdzw,dx,dy,msftx, msfty, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte ) !----------------------------------------------------------------------- ! Begin declarations. IMPLICIT NONE 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 , INTENT(IN ) :: dx, dy REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: rdzw REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: xkmh, & xkmv, & xkhh, & xkhv REAL , DIMENSION( ims:ime , kms:kme, jms:jme ), INTENT(IN ) :: & defor11, & defor22, & defor12 REAL , DIMENSION( ims:ime, jms:jme), INTENT(IN ) :: msftx, & msfty ! LOCAL VARS INTEGER :: i_start, i_end, j_start, j_end, ktf, i, j, k REAL :: deltas, tmp, pr, mlen_h, c_s REAL, DIMENSION( its:ite , kts:kte , jts:jte ) :: def2 ! End declarations. !----------------------------------------------------------------------- ktf = min(kte,kde-1) i_start = its i_end = MIN(ite,ide-1) j_start = jts j_end = MIN(jte,jde-1) IF ( config_flags%open_xs .or. config_flags%specified .or. & config_flags%nested) i_start = MAX(ids+1,its) IF ( config_flags%open_xe .or. config_flags%specified .or. & config_flags%nested) i_end = MIN(ide-2,ite) IF ( config_flags%open_ys .or. config_flags%specified .or. & config_flags%nested) j_start = MAX(jds+1,jts) IF ( config_flags%open_ye .or. config_flags%specified .or. & config_flags%nested) j_end = MIN(jde-2,jte) IF ( config_flags%periodic_x ) i_start = its IF ( config_flags%periodic_x ) i_end = MIN( ite, ide-1 ) pr=prandtl c_s = config_flags%c_s do j=j_start,j_end do k=kts,ktf do i=i_start,i_end def2(i,k,j)=0.25*((defor11(i,k,j)-defor22(i,k,j))*(defor11(i,k,j)-defor22(i,k,j))) tmp=0.25*(defor12(i ,k,j)+defor12(i ,k,j+1)+ & defor12(i+1,k,j)+defor12(i+1,k,j+1)) def2(i,k,j)=def2(i,k,j)+tmp*tmp enddo enddo enddo ! DO j = j_start, j_end DO k = kts, ktf DO i = i_start, i_end mlen_h=sqrt(dx/msftx(i,j) * dy/msfty(i,j)) tmp=sqrt(def2(i,k,j)) ! xkmh(i,k,j)=max(c_s*c_s*mlen_h*mlen_h*tmp, 1.0E-6*mlen_h*mlen_h ) xkmh(i,k,j)=c_s*c_s*mlen_h*mlen_h*tmp xkmh(i,k,j)=min(xkmh(i,k,j), 10.*mlen_h ) xkmv(i,k,j)=0. xkhh(i,k,j)=xkmh(i,k,j)/pr xkhv(i,k,j)=0. ENDDO ENDDO ENDDO END SUBROUTINE smag2d_km !======================================================================= !======================================================================= SUBROUTINE tke_km( config_flags, xkmh, xkmv, xkhh, xkhv, & bn2, tke, p8w, t8w, theta, & rdz, rdzw, dx,dy, dt, isotropic, & mix_upper_bound, msftx, msfty, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) ! History: Sep 2003 Changes by Jason Knievel and George Bryan, NCAR ! Oct 2001 Converted to mass core by Bill Skamarock, NCAR ! ... ... ! Purpose: This routine calculates the exchange coefficients for the ! TKE turbulence parameterization. ! References: Klemp and Wilhelmson (JAS 1978) ! Chen and Dudhia (NCAR WRF physics report 2000) !----------------------------------------------------------------------- ! Begin declarations. IMPLICIT NONE 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 INTEGER, INTENT( IN ) :: isotropic REAL, INTENT( IN ) & :: dx, dy, dt REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( IN ) & :: tke, p8w, t8w, theta, rdz, rdzw, bn2 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( INOUT ) & :: xkmh, xkmv, xkhh, xkhv REAL, INTENT( IN ) & :: mix_upper_bound REAL , DIMENSION( ims:ime, jms:jme), INTENT(IN ) :: msftx, & msfty ! Local variables. REAL, DIMENSION( its:ite, kts:kte, jts:jte ) & :: l_scale REAL, DIMENSION( its:ite, kts:kte, jts:jte ) & :: dthrdn REAL & :: deltas, tmp, mlen_s, mlen_h, mlen_v, tmpdz, & thetasfc, thetatop, minkx, pr_inv, pr_inv_h, pr_inv_v, c_k INTEGER & :: i_start, i_end, j_start, j_end, ktf, i, j, k REAL, PARAMETER :: tke_seed_value = 1.e-06 REAL :: tke_seed REAL, PARAMETER :: epsilon = 1.e-10 ! End declarations. !----------------------------------------------------------------------- ktf = MIN( kte, kde-1 ) i_start = its i_end = MIN( ite, ide-1 ) j_start = jts j_end = MIN( jte, jde-1 ) IF ( config_flags%open_xs .OR. config_flags%specified .OR. & config_flags%nested) i_start = MAX( ids+1, its ) IF ( config_flags%open_xe .OR. config_flags%specified .OR. & config_flags%nested) i_end = MIN( ide-2, ite ) IF ( config_flags%open_ys .OR. config_flags%specified .OR. & config_flags%nested) j_start = MAX( jds+1, jts ) IF ( config_flags%open_ye .OR. config_flags%specified .OR. & config_flags%nested) j_end = MIN( jde-2, jte) IF ( config_flags%periodic_x ) i_start = its IF ( config_flags%periodic_x ) i_end = MIN( ite, ide-1 ) ! in the absence of surface drag or a surface heat flux, there ! is no way to generate tke without pre-existing tke. Use ! tke_seed if the drag and flux are off. c_k = config_flags%c_k tke_seed = tke_seed_value if( (config_flags%tke_drag_coefficient .gt. epsilon) .or. & (config_flags%tke_heat_flux .gt. epsilon) ) tke_seed = 0. DO j = j_start, j_end DO k = kts+1, ktf-1 DO i = i_start, i_end tmpdz = 1.0 / ( rdz(i,k+1,j) + rdz(i,k,j) ) dthrdn(i,k,j) = ( theta(i,k+1,j) - theta(i,k-1,j) ) / tmpdz END DO END DO END DO k = kts DO j = j_start, j_end DO i = i_start, i_end tmpdz = 1.0 / ( rdzw(i,k+1,j) + rdzw(i,k,j) ) thetasfc = T8w(i,kts,j) / ( p8w(i,k,j) / p1000mb )**( R_d / Cp ) dthrdn(i,k,j) = ( theta(i,k+1,j) - thetasfc ) / tmpdz END DO END DO k = ktf DO j = j_start, j_end DO i = i_start, i_end tmpdz = 1.0 / rdz(i,k,j) + 0.5 / rdzw(i,k,j) thetatop = T8w(i,kde,j) / ( p8w(i,kde,j) / p1000mb )**( R_d / Cp ) dthrdn(i,k,j) = ( thetatop - theta(i,k-1,j) ) / tmpdz END DO END DO IF ( isotropic .EQ. 0 ) THEN DO j = j_start, j_end DO k = kts, ktf DO i = i_start, i_end mlen_h = SQRT( dx/msftx(i,j) * dy/msfty(i,j) ) tmp = SQRT( MAX( tke(i,k,j), tke_seed ) ) deltas = 1.0 / rdzw(i,k,j) mlen_v = deltas IF ( dthrdn(i,k,j) .GT. 0.) THEN mlen_s = 0.76 * tmp / ( ABS( g / theta(i,k,j) * dthrdn(i,k,j) ) )**0.5 mlen_v = MIN( mlen_v, mlen_s ) END IF xkmh(i,k,j) = MAX( c_k * tmp * mlen_h, 1.0E-6 * mlen_h * mlen_h ) xkmh(i,k,j) = MIN( xkmh(i,k,j), mix_upper_bound * mlen_h *mlen_h / dt ) xkmv(i,k,j) = MAX( c_k * tmp * mlen_v, 1.0E-6 * deltas * deltas ) xkmv(i,k,j) = MIN( xkmv(i,k,j), mix_upper_bound * deltas *deltas / dt ) pr_inv_h = 1./prandtl pr_inv_v = 1.0 + 2.0 * mlen_v / deltas xkhh(i,k,j) = xkmh(i,k,j) * pr_inv_h xkhv(i,k,j) = xkmv(i,k,j) * pr_inv_v END DO END DO END DO ELSE CALL calc_l_scale( config_flags, tke, BN2, l_scale, & i_start, i_end, ktf, j_start, j_end, & dx, dy, rdzw, msftx, msfty, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) DO j = j_start, j_end DO k = kts, ktf DO i = i_start, i_end tmp = SQRT( MAX( tke(i,k,j), tke_seed ) ) deltas = ( dx/msftx(i,j) * dy/msfty(i,j) / rdzw(i,k,j) )**0.33333333 xkmh(i,k,j) = c_k * tmp * l_scale(i,k,j) xkmh(i,k,j) = MIN( mix_upper_bound * dx/msftx(i,j) * dy/msfty(i,j) / dt, xkmh(i,k,j) ) xkmv(i,k,j) = c_k * tmp * l_scale(i,k,j) xkmv(i,k,j) = MIN( mix_upper_bound / rdzw(i,k,j) / rdzw(i,k,j) / dt , xkmv(i,k,j) ) pr_inv = 1.0 + 2.0 * l_scale(i,k,j) / deltas xkhh(i,k,j) = MIN( mix_upper_bound * dx/msftx(i,j) * dy/msfty(i,j) / dt, xkmh(i,k,j) * pr_inv ) xkhv(i,k,j) = MIN( mix_upper_bound / rdzw(i,k,j) / rdzw(i,k,j) / dt, xkmv(i,k,j) * pr_inv ) END DO END DO END DO END IF END SUBROUTINE tke_km !======================================================================= !======================================================================= SUBROUTINE calc_l_scale( config_flags, tke, BN2, l_scale, & i_start, i_end, ktf, j_start, j_end, & dx, dy, rdzw, msftx, msfty, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) ! History: Sep 2003 Written by Bryan and Knievel, NCAR ! Purpose: This routine calculates the length scale, based on stability, ! for TKE parameterization of subgrid-scale turbulence. !----------------------------------------------------------------------- ! Begin declarations. IMPLICIT NONE TYPE( grid_config_rec_type ), INTENT( IN ) & :: config_flags INTEGER, INTENT( IN ) & :: i_start, i_end, ktf, j_start, j_end, & 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 ) & :: BN2, tke, rdzw REAL, INTENT( IN ) & :: dx, dy REAL, DIMENSION( its:ite, kts:kte, jts:jte ), INTENT( OUT ) & :: l_scale REAL , DIMENSION( ims:ime, jms:jme), INTENT(IN ) :: msftx, & msfty ! Local variables. INTEGER & :: i, j, k REAL & :: deltas, tmp ! End declarations. !----------------------------------------------------------------------- DO j = j_start, j_end DO k = kts, ktf DO i = i_start, i_end deltas = ( dx/msftx(i,j) * dy/msfty(i,j) / rdzw(i,k,j) )**0.33333333 l_scale(i,k,j) = deltas IF ( BN2(i,k,j) .gt. 1.0e-6 ) THEN tmp = SQRT( MAX( tke(i,k,j), 1.0e-6 ) ) l_scale(i,k,j) = 0.76 * tmp / SQRT( BN2(i,k,j) ) l_scale(i,k,j) = MIN( l_scale(i,k,j), deltas) l_scale(i,k,j) = MAX( l_scale(i,k,j), 0.001 * deltas ) END IF END DO END DO END DO END SUBROUTINE calc_l_scale !======================================================================= !======================================================================= SUBROUTINE horizontal_diffusion_2 ( rt_tendf, ru_tendf, rv_tendf, rw_tendf, & tke_tendf, & moist_tendf, n_moist, & chem_tendf, n_chem, & scalar_tendf, n_scalar, & thp, theta, mu, tke, config_flags, & defor11, defor22, defor12, & defor13, defor23, div, & moist, chem, scalar, & msfux, msfuy, msfvx, msfvy, & msftx, msfty, xkmh, xkhh,km_opt, & rdx, rdy, rdz, rdzw, fnm, fnp, & cf1, cf2, cf3, zx, zy, dn, dnw, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) !----------------------------------------------------------------------- ! Begin declarations. IMPLICIT NONE 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 INTEGER , INTENT(IN ) :: n_moist, n_chem, n_scalar, km_opt REAL , INTENT(IN ) :: cf1, cf2, cf3 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fnm REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fnp REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: dnw REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: dn REAL , DIMENSION( ims:ime, jms:jme) , INTENT(IN ) :: msfux, & msfuy, & msfvx, & msfvy, & msftx, & msfty, & mu REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(INOUT) ::rt_tendf,& ru_tendf,& rv_tendf,& rw_tendf,& tke_tendf REAL , DIMENSION( ims:ime, kms:kme, jms:jme, n_moist), & INTENT(INOUT) :: moist_tendf REAL , DIMENSION( ims:ime, kms:kme, jms:jme, n_chem), & INTENT(INOUT) :: chem_tendf REAL , DIMENSION( ims:ime, kms:kme, jms:jme, n_scalar), & INTENT(INOUT) :: scalar_tendf REAL , DIMENSION( ims:ime, kms:kme, jms:jme, n_moist), & INTENT(IN ) :: moist REAL , DIMENSION( ims:ime, kms:kme, jms:jme, n_chem), & INTENT(IN ) :: chem REAL , DIMENSION( ims:ime, kms:kme, jms:jme, n_scalar) , & INTENT(IN ) :: scalar REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(IN ) ::defor11, & defor22, & defor12, & defor13, & defor23, & div, & xkmh, & xkhh, & zx, & zy, & theta, & thp, & tke, & rdz, & rdzw REAL , INTENT(IN ) :: rdx, & rdy ! LOCAL VARS INTEGER :: im, ic, is ! REAL , DIMENSION(its-1:ite+1, kts:kte, jts-1:jte+1) :: xkhh ! End declarations. !----------------------------------------------------------------------- !----------------------------------------------------------------------- ! Call diffusion subroutines. CALL horizontal_diffusion_u_2( ru_tendf, mu, config_flags, & defor11, defor12, div, & tke(ims,kms,jms), & msfux, msfuy, xkmh, rdx, rdy, fnm, fnp, & zx, zy, rdzw, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) CALL horizontal_diffusion_v_2( rv_tendf, mu, config_flags, & defor12, defor22, div, & tke(ims,kms,jms), & msfvx, msfvy, xkmh, rdx, rdy, fnm, fnp, & zx, zy, rdzw, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) CALL horizontal_diffusion_w_2( rw_tendf, mu, config_flags, & defor13, defor23, div, & tke(ims,kms,jms), & msftx, msfty, xkmh, rdx, rdy, fnm, fnp, & zx, zy, rdz, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) CALL horizontal_diffusion_s ( rt_tendf, mu, config_flags, thp, & msftx, msfty, msfux, msfuy, & msfvx, msfvy, xkhh, rdx, rdy, & fnm, fnp, cf1, cf2, cf3, & zx, zy, rdz, rdzw, dnw, dn, & .false., & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) IF (km_opt .eq. 2) & CALL horizontal_diffusion_s ( tke_tendf(ims,kms,jms), & mu, config_flags, & tke(ims,kms,jms), & msftx, msfty, msfux, msfuy, & msfvx, msfvy, xkhh, rdx, rdy, & fnm, fnp, cf1, cf2, cf3, & zx, zy, rdz, rdzw, dnw, dn, & .true., & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) IF (n_moist .ge. PARAM_FIRST_SCALAR) THEN moist_loop: do im = PARAM_FIRST_SCALAR, n_moist CALL horizontal_diffusion_s( moist_tendf(ims,kms,jms,im), & mu, config_flags, & moist(ims,kms,jms,im), & msftx, msfty, msfux, msfuy, & msfvx, msfvy, xkhh, rdx, rdy, & fnm, fnp, cf1, cf2, cf3, & zx, zy, rdz, rdzw, dnw, dn, & .false., & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) ENDDO moist_loop ENDIF IF (n_chem .ge. PARAM_FIRST_SCALAR) THEN chem_loop: do ic = PARAM_FIRST_SCALAR, n_chem CALL horizontal_diffusion_s( chem_tendf(ims,kms,jms,ic), & mu, config_flags, & chem(ims,kms,jms,ic), & msftx, msfty, msfux, msfuy, & msfvx, msfvy, xkhh, rdx, rdy, & fnm, fnp, cf1, cf2, cf3, & zx, zy, rdz, rdzw, dnw, dn, & .false., & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) ENDDO chem_loop ENDIF IF (n_scalar .ge. PARAM_FIRST_SCALAR) THEN scalar_loop: do is = PARAM_FIRST_SCALAR, n_scalar CALL horizontal_diffusion_s( scalar_tendf(ims,kms,jms,is), & mu, config_flags, & scalar(ims,kms,jms,is), & msftx, msfty, msfux, msfuy, & msfvx, msfvy, xkhh, rdx, rdy, & fnm, fnp, cf1, cf2, cf3, & zx, zy, rdz, rdzw, dnw, dn, & .false., & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) ENDDO scalar_loop ENDIF END SUBROUTINE horizontal_diffusion_2 !======================================================================= !======================================================================= SUBROUTINE horizontal_diffusion_u_2( tendency, mu, config_flags, & defor11, defor12, div, tke, & msfux, msfuy, & xkmh, rdx, rdy, fnm, fnp, & zx, zy, rdzw, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) !----------------------------------------------------------------------- ! Begin declarations. IMPLICIT NONE 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( kms:kme ) , INTENT(IN ) :: fnm REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fnp REAL , DIMENSION( ims:ime, jms:jme) , INTENT(IN ) :: msfux, & msfuy, & mu REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(INOUT) ::tendency REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(IN ) :: rdzw REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(IN ) ::defor11, & defor12, & div, & tke, & xkmh, & zx, & zy REAL , INTENT(IN ) :: rdx, & rdy ! Local data INTEGER :: i, j, k, ktf INTEGER :: i_start, i_end, j_start, j_end INTEGER :: is_ext,ie_ext,js_ext,je_ext REAL , DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1) :: titau1avg, & titau2avg, & titau1, & titau2, & xkxavg, & rravg ! new ! zxavg, & ! zyavg REAL :: mrdx, mrdy, rcoup REAL :: tmpzy, tmpzeta_z REAL :: term1, term2, term3 ! End declarations. !----------------------------------------------------------------------- ktf=MIN(kte,kde-1) !----------------------------------------------------------------------- ! u : p (.), u(|), w(-) ! ! p u p u u u ! ! p | . | . | . | k+1 | . | . | . | k+1 ! ! w - 13 - - k+1 13 k+1 ! ! p | 11 O 11 | . | k | 12 O 12 | . | k ! ! w - 13 - - k 13 k ! ! p | . | . | . | k-1 | . | . | . | k-1 ! ! i-1 i i i+1 j-1 j j j+1 j+1 ! i_start = its i_end = ite j_start = jts j_end = MIN(jte,jde-1) IF ( config_flags%open_xs .or. config_flags%specified .or. & config_flags%nested) i_start = MAX(ids+1,its) IF ( config_flags%open_xe .or. config_flags%specified .or. & config_flags%nested) i_end = MIN(ide-1,ite) IF ( config_flags%open_ys .or. config_flags%specified .or. & config_flags%nested) j_start = MAX(jds+1,jts) IF ( config_flags%open_ye .or. config_flags%specified .or. & config_flags%nested) j_end = MIN(jde-2,jte) IF ( config_flags%periodic_x ) i_start = its IF ( config_flags%periodic_x ) i_end = ite ! titau1 = titau11 is_ext=1 ie_ext=0 js_ext=0 je_ext=0 CALL cal_titau_11_22_33( config_flags, titau1, & mu, tke, xkmh, defor11, & is_ext, ie_ext, js_ext, je_ext, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) ! titau2 = titau12 is_ext=0 ie_ext=0 js_ext=0 je_ext=1 CALL cal_titau_12_21( config_flags, titau2, & mu, xkmh, defor12, & is_ext, ie_ext, js_ext, je_ext, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) ! titau1avg = titau11avg ! titau2avg = titau12avg DO j = j_start, j_end DO k = kts+1,ktf DO i = i_start, i_end titau1avg(i,k,j)=0.5*(fnm(k)*(titau1(i-1,k ,j)+titau1(i,k ,j))+ & fnp(k)*(titau1(i-1,k-1,j)+titau1(i,k-1,j))) titau2avg(i,k,j)=0.5*(fnm(k)*(titau2(i,k ,j+1)+titau2(i,k ,j))+ & fnp(k)*(titau2(i,k-1,j+1)+titau2(i,k-1,j))) tmpzy = 0.25*( zy(i-1,k,j )+zy(i,k,j )+ & zy(i-1,k,j+1)+zy(i,k,j+1) ) ! tmpzeta_z = 0.5*(zeta_z(i,j)+zeta_z(i-1,j)) ! titau1avg(i,k,j)=titau1avg(i,k,j)*zx(i,k,j)*tmpzeta_z ! titau2avg(i,k,j)=titau2avg(i,k,j)*tmpzy *tmpzeta_z titau1avg(i,k,j)=titau1avg(i,k,j)*zx(i,k,j) titau2avg(i,k,j)=titau2avg(i,k,j)*tmpzy ENDDO ENDDO ENDDO ! DO j = j_start, j_end DO i = i_start, i_end titau1avg(i,kts,j)=0. titau1avg(i,ktf+1,j)=0. titau2avg(i,kts,j)=0. titau2avg(i,ktf+1,j)=0. ENDDO ENDDO ! DO j = j_start, j_end DO k = kts,ktf DO i = i_start, i_end mrdx=msfux(i,j)*rdx mrdy=msfuy(i,j)*rdy tendency(i,k,j)=tendency(i,k,j)- & (mrdx*(titau1(i,k,j )-titau1(i-1,k,j))+ & mrdy*(titau2(i,k,j+1)-titau2(i,k,j ))- & msfuy(i,j)*rdzw(i,k,j)*((titau1avg(i,k+1,j)-titau1avg(i,k,j))+ & (titau2avg(i,k+1,j)-titau2avg(i,k,j)) & ) ) ENDDO ENDDO ENDDO END SUBROUTINE horizontal_diffusion_u_2 !======================================================================= !======================================================================= SUBROUTINE horizontal_diffusion_v_2( tendency, mu, config_flags, & defor12, defor22, div, tke, & msfvx, msfvy, & xkmh, rdx, rdy, fnm, fnp, & zx, zy, rdzw, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) !----------------------------------------------------------------------- ! Begin declarations. IMPLICIT NONE 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( kms:kme ) , INTENT(IN ) :: fnm REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fnp REAL , DIMENSION( ims:ime, jms:jme) , INTENT(IN ) :: msfvx, & msfvy, & mu REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: tendency REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(IN ) ::defor12, & defor22, & div, & tke, & xkmh, & zx, & zy, & rdzw REAL , INTENT(IN ) :: rdx, & rdy ! Local data INTEGER :: i, j, k, ktf INTEGER :: i_start, i_end, j_start, j_end INTEGER :: is_ext,ie_ext,js_ext,je_ext REAL , DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1) :: titau1avg, & titau2avg, & titau1, & titau2, & xkxavg, & rravg ! new ! zxavg, & ! zyavg REAL :: mrdx, mrdy, rcoup REAL :: tmpzx, tmpzeta_z ! End declarations. !----------------------------------------------------------------------- ktf=MIN(kte,kde-1) !----------------------------------------------------------------------- ! v : p (.), v(+), w(-) ! ! p v p v v v ! ! p + . + . + . + k+1 + . + . + . + k+1 ! ! w - 23 - - k+1 23 k+1 ! ! p + 22 O 22 + . + k + 21 O 21 + . + k ! ! w - 23 - - k 23 k ! ! p + . + . + . + k-1 + . + . + . + k-1 ! ! j-1 j j j+1 i-1 i i i+1 i+1 ! i_start = its i_end = MIN(ite,ide-1) j_start = jts j_end = jte IF ( config_flags%open_xs .or. config_flags%specified .or. & config_flags%nested) i_start = MAX(ids+1,its) IF ( config_flags%open_xe .or. config_flags%specified .or. & config_flags%nested) i_end = MIN(ide-2,ite) IF ( config_flags%open_ys .or. config_flags%specified .or. & config_flags%nested) j_start = MAX(jds+1,jts) IF ( config_flags%open_ye .or. config_flags%specified .or. & config_flags%nested) j_end = MIN(jde-1,jte) IF ( config_flags%periodic_x ) i_start = its IF ( config_flags%periodic_x ) i_end = MIN(ite,ide-1) ! titau1 = titau21 is_ext=0 ie_ext=1 js_ext=0 je_ext=0 CALL cal_titau_12_21( config_flags, titau1, & mu, xkmh, defor12, & is_ext,ie_ext,js_ext,je_ext, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) ! titau2 = titau22 is_ext=0 ie_ext=0 js_ext=1 je_ext=0 CALL cal_titau_11_22_33( config_flags, titau2, & mu, tke, xkmh, defor22, & is_ext, ie_ext, js_ext, je_ext, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) DO j = j_start, j_end DO k = kts+1,ktf DO i = i_start, i_end titau1avg(i,k,j)=0.5*(fnm(k)*(titau1(i+1,k ,j)+titau1(i,k ,j))+ & fnp(k)*(titau1(i+1,k-1,j)+titau1(i,k-1,j))) titau2avg(i,k,j)=0.5*(fnm(k)*(titau2(i,k ,j-1)+titau2(i,k ,j))+ & fnp(k)*(titau2(i,k-1,j-1)+titau2(i,k-1,j))) tmpzx = 0.25*( zx(i,k,j )+zx(i+1,k,j )+ & zx(i,k,j-1)+zx(i+1,k,j-1) ) titau1avg(i,k,j)=titau1avg(i,k,j)*tmpzx titau2avg(i,k,j)=titau2avg(i,k,j)*zy(i,k,j) ENDDO ENDDO ENDDO DO j = j_start, j_end DO i = i_start, i_end titau1avg(i,kts,j)=0. titau1avg(i,ktf+1,j)=0. titau2avg(i,kts,j)=0. titau2avg(i,ktf+1,j)=0. ENDDO ENDDO ! DO j = j_start, j_end DO k = kts,ktf DO i = i_start, i_end mrdx=msfvx(i,j)*rdx mrdy=msfvy(i,j)*rdy tendency(i,k,j)=tendency(i,k,j)- & (mrdy*(titau2(i ,k,j)-titau2(i,k,j-1))+ & mrdx*(titau1(i+1,k,j)-titau1(i,k,j ))- & msfvy(i,j)*rdzw(i,k,j)*((titau1avg(i,k+1,j)-titau1avg(i,k,j))+ & (titau2avg(i,k+1,j)-titau2avg(i,k,j)) & ) & ) ENDDO ENDDO ENDDO END SUBROUTINE horizontal_diffusion_v_2 !======================================================================= !======================================================================= SUBROUTINE horizontal_diffusion_w_2( tendency, mu, config_flags, & defor13, defor23, div, tke, & msftx, msfty, & xkmh, rdx, rdy, fnm, fnp, & zx, zy, rdz, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) !----------------------------------------------------------------------- ! Begin declarations. IMPLICIT NONE 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( kms:kme ) , INTENT(IN ) :: fnm REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fnp REAL , DIMENSION( ims:ime, jms:jme) , INTENT(IN ) :: msftx, & msfty, & mu REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: tendency REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(IN ) ::defor13, & defor23, & div, & tke, & xkmh, & zx, & zy, & rdz REAL , INTENT(IN ) :: rdx, & rdy ! Local data INTEGER :: i, j, k, ktf INTEGER :: i_start, i_end, j_start, j_end INTEGER :: is_ext,ie_ext,js_ext,je_ext REAL , DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1) :: titau1avg, & titau2avg, & titau1, & titau2, & xkxavg, & rravg ! new ! zxavg, & ! zyavg REAL :: mrdx, mrdy, rcoup REAL :: tmpzx, tmpzy, tmpzeta_z ! End declarations. !----------------------------------------------------------------------- ktf=MIN(kte,kde-1) !----------------------------------------------------------------------- ! w : p (.), u(|), v(+), w(-) ! ! p u p u p v p v ! ! w - - - k+1 w - - - k+1 ! ! p . | 33 | . k p . + 33 + . k ! ! w - 31 O 31 - k w - 32 O 32 - k ! ! p . | 33 | . k-1 p . | 33 | . k-1 ! ! w - - - k-1 w - - - k-1 ! ! i-1 i i i+1 j-1 j j j+1 ! i_start = its i_end = MIN(ite,ide-1) j_start = jts j_end = MIN(jte,jde-1) IF ( config_flags%open_xs .or. config_flags%specified .or. & config_flags%nested) i_start = MAX(ids+1,its) IF ( config_flags%open_xe .or. config_flags%specified .or. & config_flags%nested) i_end = MIN(ide-2,ite) IF ( config_flags%open_ys .or. config_flags%specified .or. & config_flags%nested) j_start = MAX(jds+1,jts) IF ( config_flags%open_ye .or. config_flags%specified .or. & config_flags%nested) j_end = MIN(jde-2,jte) IF ( config_flags%periodic_x ) i_start = its IF ( config_flags%periodic_x ) i_end = MIN(ite,ide-1) ! titau1 = titau31 is_ext=0 ie_ext=1 js_ext=0 je_ext=0 CALL cal_titau_13_31( config_flags, titau1, defor13, & mu, xkmh, fnm, fnp, & is_ext, ie_ext, js_ext, je_ext, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) ! titau2 = titau32 is_ext=0 ie_ext=0 js_ext=0 je_ext=1 CALL cal_titau_23_32( config_flags, titau2, defor23, & mu, xkmh, fnm, fnp, & is_ext, ie_ext, js_ext, je_ext, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) ! titau1avg = titau31avg * zx * zeta_z = titau13avg * zx * zeta_z ! titau2avg = titau32avg * zy * zeta_z = titau23avg * zy * zeta_z DO j = j_start, j_end DO k = kts,ktf DO i = i_start, i_end titau1avg(i,k,j)=0.25*(titau1(i+1,k+1,j)+titau1(i,k+1,j)+ & titau1(i+1,k ,j)+titau1(i,k ,j)) titau2avg(i,k,j)=0.25*(titau2(i,k+1,j+1)+titau2(i,k+1,j)+ & titau2(i,k ,j+1)+titau2(i,k ,j)) ! new tmpzx =0.25*( zx(i,k ,j)+zx(i+1,k ,j)+ & zx(i,k+1,j)+zx(i+1,k+1,j) ) tmpzy =0.25*( zy(i,k ,j)+zy(i,k ,j+1)+ & zy(i,k+1,j)+zy(i,k+1,j+1) ) titau1avg(i,k,j)=titau1avg(i,k,j)*tmpzx titau2avg(i,k,j)=titau2avg(i,k,j)*tmpzy ! titau1avg(i,k,j)=titau1avg(i,k,j)*tmpzx*zeta_z(i,j) ! titau2avg(i,k,j)=titau2avg(i,k,j)*tmpzy*zeta_z(i,j) ENDDO ENDDO ENDDO DO j = j_start, j_end DO i = i_start, i_end titau1avg(i,ktf+1,j)=0. titau2avg(i,ktf+1,j)=0. ENDDO ENDDO DO j = j_start, j_end DO k = kts+1,ktf DO i = i_start, i_end mrdx=msftx(i,j)*rdx mrdy=msfty(i,j)*rdy tendency(i,k,j)=tendency(i,k,j)- & (mrdx*(titau1(i+1,k,j)-titau1(i,k,j))+ & mrdy*(titau2(i,k,j+1)-titau2(i,k,j))- & msfty(i,j)*rdz(i,k,j)*(titau1avg(i,k,j)-titau1avg(i,k-1,j)+ & titau2avg(i,k,j)-titau2avg(i,k-1,j) & ) & ) ! msft(i,j)/dn(k)*(titau1avg(i,k,j)-titau1avg(i,k-1,j)+ & ! titau2avg(i,k,j)-titau2avg(i,k-1,j) & ! ) & ! ) ENDDO ENDDO ENDDO END SUBROUTINE horizontal_diffusion_w_2 !======================================================================= !======================================================================= SUBROUTINE horizontal_diffusion_s (tendency, mu, config_flags, var, & msftx, msfty, msfux, msfuy, & msfvx, msfvy, xkhh, rdx, rdy, & fnm, fnp, cf1, cf2, cf3, & zx, zy, rdz, rdzw, dnw, dn, & doing_tke, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) !----------------------------------------------------------------------- ! Begin declarations. IMPLICIT NONE 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 LOGICAL, INTENT(IN ) :: doing_tke REAL , INTENT(IN ) :: cf1, cf2, cf3 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fnm REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fnp REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: dn REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: dnw REAL , DIMENSION( ims:ime, jms:jme) , INTENT(IN ) :: msfux REAL , DIMENSION( ims:ime, jms:jme) , INTENT(IN ) :: msfuy REAL , DIMENSION( ims:ime, jms:jme) , INTENT(IN ) :: msfvx REAL , DIMENSION( ims:ime, jms:jme) , INTENT(IN ) :: msfvy REAL , DIMENSION( ims:ime, jms:jme) , INTENT(IN ) :: msftx REAL , DIMENSION( ims:ime, jms:jme) , INTENT(IN ) :: msfty REAL , DIMENSION( ims:ime, jms:jme) , INTENT(IN ) :: mu ! REAL , DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1), & ! INTENT(IN ) :: xkhh REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: tendency REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(IN ) :: & xkhh, & rdz, & rdzw REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(IN ) :: var, & zx, & zy REAL , INTENT(IN ) :: rdx, & rdy ! Local data INTEGER :: i, j, k, ktf INTEGER :: i_start, i_end, j_start, j_end REAL , DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1) :: H1avg, & H2avg, & H1, & H2, & xkxavg ! new ! zxavg, & ! zyavg REAL , DIMENSION( its:ite, kts:kte, jts:jte) :: tmptendf REAL :: mrdx, mrdy, rcoup REAL :: tmpzx, tmpzy, tmpzeta_z, rdzu, rdzv INTEGER :: ktes1,ktes2 ! End declarations. !----------------------------------------------------------------------- ktf=MIN(kte,kde-1) !----------------------------------------------------------------------- ! scalars: t (.), u(|), v(+), w(-) ! ! t u t u t v t v ! ! w - 3 - k+1 w - 3 - k+1 ! ! t . 1 O 1 . k t . 2 O 2 . k ! ! w - 3 - k w - 3 - k ! ! t . | . | . k-1 t . + . + . k-1 ! ! w - - - k-1 w - - - k-1 ! ! t i-1 i i i+1 j-1 j j j+1 ! ktes1=kte-1 ktes2=kte-2 i_start = its i_end = MIN(ite,ide-1) j_start = jts j_end = MIN(jte,jde-1) IF ( config_flags%open_xs .or. config_flags%specified .or. & config_flags%nested) i_start = MAX(ids+1,its) IF ( config_flags%open_xe .or. config_flags%specified .or. & config_flags%nested) i_end = MIN(ide-2,ite) IF ( config_flags%open_ys .or. config_flags%specified .or. & config_flags%nested) j_start = MAX(jds+1,jts) IF ( config_flags%open_ye .or. config_flags%specified .or. & config_flags%nested) j_end = MIN(jde-2,jte) IF ( config_flags%periodic_x ) i_start = its IF ( config_flags%periodic_x ) i_end = MIN(ite,ide-1) ! diffusion of the TKE needs mutiple 2 IF ( doing_tke ) THEN DO j = j_start, j_end DO k = kts,ktf DO i = i_start, i_end tmptendf(i,k,j)=tendency(i,k,j) ENDDO ENDDO ENDDO ENDIF ! H1 = partial var over partial x DO j = j_start, j_end DO k = kts, ktf DO i = i_start, i_end + 1 ! new ! zxavg(i,k,j) =0.5*( zx(i-1,k,j)+ zx(i,k,j)) xkxavg(i,k,j)=0.5*(xkhh(i-1,k,j)+xkhh(i,k,j)) ENDDO ENDDO ENDDO DO j = j_start, j_end DO k = kts+1, ktf DO i = i_start, i_end + 1 H1avg(i,k,j)=0.5*(fnm(k)*(var(i-1,k ,j)+var(i,k ,j))+ & fnp(k)*(var(i-1,k-1,j)+var(i,k-1,j))) ENDDO ENDDO ENDDO DO j = j_start, j_end DO i = i_start, i_end + 1 H1avg(i,kts ,j)=0.5*(cf1*var(i ,1,j)+cf2*var(i ,2,j)+ & cf3*var(i ,3,j)+cf1*var(i-1,1,j)+ & cf2*var(i-1,2,j)+cf3*var(i-1,3,j)) H1avg(i,ktf+1,j)=0.5*(var(i,ktes1,j)+(var(i,ktes1,j)- & var(i,ktes2,j))*0.5*dnw(ktes1)/dn(ktes1)+ & var(i-1,ktes1,j)+(var(i-1,ktes1,j)- & var(i-1,ktes2,j))*0.5*dnw(ktes1)/dn(ktes1)) ENDDO ENDDO DO j = j_start, j_end DO k = kts, ktf DO i = i_start, i_end + 1 ! new tmpzx = 0.5*( zx(i,k,j)+ zx(i,k+1,j)) rdzu = 2./(1./rdzw(i,k,j) + 1./rdzw(i-1,k,j)) H1(i,k,j)=-msfuy(i,j)*xkxavg(i,k,j)*( & rdx*(var(i,k,j)-var(i-1,k,j)) - tmpzx* & (H1avg(i,k+1,j)-H1avg(i,k,j))*rdzu ) ! tmpzeta_z = 0.5*(zeta_z(i,j)+zeta_z(i-1,j)) ! H1(i,k,j)=-msfuy(i,j)*xkxavg(i,k,j)*( & ! rdx*(var(i,k,j)-var(i-1,k,j)) - tmpzx*tmpzeta_z* & ! (H1avg(i,k+1,j)-H1avg(i,k,j))/dnw(k)) ENDDO ENDDO ENDDO ! H2 = partial var over partial y DO j = j_start, j_end + 1 DO k = kts, ktf DO i = i_start, i_end ! new ! zyavg(i,k,j) =0.5*( zy(i,k,j-1)+ zy(i,k,j)) xkxavg(i,k,j)=0.5*(xkhh(i,k,j-1)+xkhh(i,k,j)) ENDDO ENDDO ENDDO DO j = j_start, j_end + 1 DO k = kts+1, ktf DO i = i_start, i_end ! new H2avg(i,k,j)=0.5*(fnm(k)*(var(i,k ,j-1)+var(i,k ,j))+ & fnp(k)*(var(i,k-1,j-1)+var(i,k-1,j))) ENDDO ENDDO ENDDO DO j = j_start, j_end + 1 DO i = i_start, i_end H2avg(i,kts ,j)=0.5*(cf1*var(i,1,j )+cf2*var(i ,2,j)+ & cf3*var(i,3,j )+cf1*var(i,1,j-1)+ & cf2*var(i,2,j-1)+cf3*var(i,3,j-1)) H2avg(i,ktf+1,j)=0.5*(var(i,ktes1,j)+(var(i,ktes1,j)- & var(i,ktes2,j))*0.5*dnw(ktes1)/dn(ktes1)+ & var(i,ktes1,j-1)+(var(i,ktes1,j-1)- & var(i,ktes2,j-1))*0.5*dnw(ktes1)/dn(ktes1)) ENDDO ENDDO DO j = j_start, j_end + 1 DO k = kts, ktf DO i = i_start, i_end ! new tmpzy = 0.5*( zy(i,k,j)+ zy(i,k+1,j)) rdzv = 2./(1./rdzw(i,k,j) + 1./rdzw(i,k,j-1)) H2(i,k,j)=-msfvy(i,j)*xkxavg(i,k,j)*( & rdy*(var(i,k,j)-var(i,k,j-1)) - tmpzy* & (H2avg(i ,k+1,j)-H2avg(i,k,j))*rdzv) ! tmpzeta_z = 0.5*(zeta_z(i,j)+zeta_z(i,j-1)) ! H2(i,k,j)=-msfvy(i,j)*xkxavg(i,k,j)*( & ! rdy*(var(i,k,j)-var(i,k,j-1)) - tmpzy*tmpzeta_z* & ! (H2avg(i ,k+1,j)-H2avg(i,k,j))/dnw(k)) ENDDO ENDDO ENDDO DO j = j_start, j_end DO k = kts+1, ktf DO i = i_start, i_end H1avg(i,k,j)=0.5*(fnm(k)*(H1(i+1,k ,j)+H1(i,k ,j))+ & fnp(k)*(H1(i+1,k-1,j)+H1(i,k-1,j))) H2avg(i,k,j)=0.5*(fnm(k)*(H2(i,k ,j+1)+H2(i,k ,j))+ & fnp(k)*(H2(i,k-1,j+1)+H2(i,k-1,j))) ! new ! zxavg(i,k,j)=fnm(k)*zx(i,k,j)+fnp(k)*zx(i,k-1,j) ! zyavg(i,k,j)=fnm(k)*zy(i,k,j)+fnp(k)*zy(i,k-1,j) ! H1avg(i,k,j)=zx*H1avg*zeta_z ! H2avg(i,k,j)=zy*H2avg*zeta_z tmpzx = 0.5*( zx(i,k,j)+ zx(i+1,k,j )) tmpzy = 0.5*( zy(i,k,j)+ zy(i ,k,j+1)) H1avg(i,k,j)=H1avg(i,k,j)*tmpzx H2avg(i,k,j)=H2avg(i,k,j)*tmpzy ! H1avg(i,k,j)=H1avg(i,k,j)*tmpzx*zeta_z(i,j) ! H2avg(i,k,j)=H2avg(i,k,j)*tmpzy*zeta_z(i,j) ENDDO ENDDO ENDDO DO j = j_start, j_end DO i = i_start, i_end H1avg(i,kts ,j)=0. H1avg(i,ktf+1,j)=0. H2avg(i,kts ,j)=0. H2avg(i,ktf+1,j)=0. ENDDO ENDDO DO j = j_start, j_end DO k = kts,ktf DO i = i_start, i_end mrdx=msftx(i,j)*rdx mrdy=msfty(i,j)*rdy tendency(i,k,j)=tendency(i,k,j)- & (mrdx*0.5*((mu(i+1,j)+mu(i,j))*H1(i+1,k,j)- & (mu(i-1,j)+mu(i,j))*H1(i ,k,j))+ & mrdy*0.5*((mu(i,j+1)+mu(i,j))*H2(i,k,j+1)- & (mu(i,j-1)+mu(i,j))*H2(i,k,j ))- & msfty(i,j)*mu(i,j)*(H1avg(i,k+1,j)-H1avg(i,k,j)+ & H2avg(i,k+1,j)-H2avg(i,k,j) & )*rdzw(i,k,j) & ) ENDDO ENDDO ENDDO IF ( doing_tke ) THEN DO j = j_start, j_end DO k = kts,ktf DO i = i_start, i_end tendency(i,k,j)=tmptendf(i,k,j)+2.* & (tendency(i,k,j)-tmptendf(i,k,j)) ENDDO ENDDO ENDDO ENDIF END SUBROUTINE horizontal_diffusion_s !======================================================================= !======================================================================= SUBROUTINE vertical_diffusion_2 ( ru_tendf, rv_tendf, rw_tendf, rt_tendf, & tke_tendf, moist_tendf, n_moist, & chem_tendf, n_chem, & scalar_tendf, n_scalar, & u_2, v_2, & thp,u_base,v_base,t_base,qv_base,mu,tke, & config_flags,defor13,defor23,defor33,div, & moist, chem, scalar, xkmv, xkhv,km_opt, & fnm, fnp, dn, dnw, rdz, rdzw, & hfx, qfx, ust, rho, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) !----------------------------------------------------------------------- ! Begin declarations. IMPLICIT NONE 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 INTEGER , INTENT(IN ) :: n_moist, n_chem, n_scalar, km_opt REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fnm REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fnp REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: dnw REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: dn REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mu REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: qv_base REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: u_base REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: v_base REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: t_base REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(INOUT) ::ru_tendf,& rv_tendf,& rw_tendf,& tke_tendf,& rt_tendf REAL , DIMENSION( ims:ime, kms:kme, jms:jme, n_moist), & INTENT(INOUT) :: moist_tendf REAL , DIMENSION( ims:ime, kms:kme, jms:jme, n_chem), & INTENT(INOUT) :: chem_tendf REAL , DIMENSION( ims:ime, kms:kme, jms:jme, n_scalar) , & INTENT(INOUT) :: scalar_tendf REAL , DIMENSION( ims:ime, kms:kme, jms:jme, n_moist), & INTENT(INOUT) :: moist REAL , DIMENSION( ims:ime, kms:kme, jms:jme, n_chem), & INTENT(INOUT) :: chem REAL , DIMENSION( ims:ime, kms:kme, jms:jme, n_scalar) , & INTENT(IN ) :: scalar REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(IN ) ::defor13, & defor23, & defor33, & div, & xkmv, & xkhv, & tke, & rdz, & u_2, & v_2, & rdzw REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(IN ) :: rho REAL , DIMENSION( ims:ime, jms:jme), INTENT(IN ) :: hfx, & qfx REAL , DIMENSION( ims:ime, jms:jme), INTENT(IN ) :: ust REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: thp ! LOCAL VAR REAL , DIMENSION( ims:ime, kms:kme, jms:jme) :: var_mix INTEGER :: im, i,j,k INTEGER :: i_start, i_end, j_start, j_end ! REAL , DIMENSION( its:ite, kts:kte, jts:jte) :: xkhv !*************************************************************************** !*************************************************************************** !MODIFICA VARIABILI PER I FLUSSI ! REAL :: V0_u,V0_v,tao_xz,tao_yz,ustar,cd0 REAL :: xsfc,psi1,vk2,zrough,lnz REAL :: heat_flux, moist_flux, heat_flux0 REAL :: cpm ! !FINE MODIFICA VARIABILI PER I FLUSSI !*************************************************************************** ! ! End declarations. !----------------------------------------------------------------------- i_start = its i_end = MIN(ite,ide-1) j_start = jts j_end = MIN(jte,jde-1) ! !----------------------------------------------------------------------- CALL vertical_diffusion_u_2( ru_tendf, config_flags, mu, & defor13, xkmv, & dnw, rdzw, fnm, fnp, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) CALL vertical_diffusion_v_2( rv_tendf, config_flags, mu, & defor23, xkmv, & dnw, rdzw, fnm, fnp, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) CALL vertical_diffusion_w_2( rw_tendf, config_flags, mu, & defor33, tke(ims,kms,jms), & div, xkmv, & dn, rdz, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) !***************************************** !***************************************** ! MODIFICA al flusso di momento alla parete ! vflux: SELECT CASE( config_flags%isfflx ) CASE (0) ! Assume cd a constant, specified in namelist cd0 = config_flags%tke_drag_coefficient ! constant drag coefficient ! set in namelist.input ! !calcolo del modulo della velocita DO j = j_start, j_end DO i = i_start, ite V0_u=0. tao_xz=0. V0_u= sqrt((u_2(i,kts,j)**2) + & (((v_2(i ,kts,j )+ & v_2(i ,kts,j+1)+ & v_2(i-1,kts,j )+ & v_2(i-1,kts,j+1))/4)**2))+epsilon tao_xz=cd0*V0_u*u_2(i,kts,j) ru_tendf(i,kts,j)=ru_tendf(i,kts,j) & -0.25*(mu(i,j)+mu(i-1,j))*tao_xz*(rdzw(i,kts,j)+rdzw(i-1,kts,j)) ENDDO ENDDO ! DO j = j_start, jte DO i = i_start, i_end V0_v=0. tao_yz=0. V0_v= sqrt((v_2(i,kts,j)**2) + & (((u_2(i ,kts,j )+ & u_2(i ,kts,j-1)+ & u_2(i+1,kts,j )+ & u_2(i+1,kts,j-1))/4)**2))+epsilon tao_yz=cd0*V0_v*v_2(i,kts,j) rv_tendf(i,kts,j)=rv_tendf(i,kts,j) & -0.25*(mu(i,j)+mu(i,j-1))*tao_yz*(rdzw(i,kts,j)+rdzw(i,kts,j-1)) ENDDO ENDDO CASE (1,2) ! ustar computed from surface routine DO j = j_start, j_end DO i = i_start, ite V0_u=0. tao_xz=0. V0_u= sqrt((u_2(i,kts,j)**2) + & (((v_2(i ,kts,j )+ & v_2(i ,kts,j+1)+ & v_2(i-1,kts,j )+ & v_2(i-1,kts,j+1))/4)**2))+epsilon ustar=0.5*(ust(i,j)+ust(i-1,j)) tao_xz=ustar*ustar*u_2(i,kts,j)/V0_u ru_tendf(i,kts,j)=ru_tendf(i,kts,j) & -0.25*(mu(i,j)+mu(i-1,j))*tao_xz*(rdzw(i,kts,j)+rdzw(i-1,kts,j)) ENDDO ENDDO DO j = j_start, jte DO i = i_start, i_end V0_v=0. tao_yz=0. V0_v= sqrt((v_2(i,kts,j)**2) + & (((u_2(i ,kts,j )+ & u_2(i ,kts,j-1)+ & u_2(i+1,kts,j )+ & u_2(i+1,kts,j-1))/4)**2))+epsilon ustar=0.5*(ust(i,j)+ust(i,j-1)) tao_yz=ustar*ustar*v_2(i,kts,j)/V0_v rv_tendf(i,kts,j)=rv_tendf(i,kts,j) & -0.25*(mu(i,j)+mu(i,j-1))*tao_yz*(rdzw(i,kts,j)+rdzw(i,kts,j-1)) ENDDO ENDDO CASE DEFAULT CALL wrf_error_fatal( 'isfflx value invalid for diff_opt=2' ) END SELECT vflux ! ! FINE MODIFICA al flusso di momento alla parete !***************************************** !***************************************** IF ( config_flags%mix_full_fields ) THEN DO j=jts,min(jte,jde-1) DO k=kts,kte-1 DO i=its,min(ite,ide-1) var_mix(i,k,j) = thp(i,k,j) ENDDO ENDDO ENDDO ELSE DO j=jts,min(jte,jde-1) DO k=kts,kte-1 DO i=its,min(ite,ide-1) var_mix(i,k,j) = thp(i,k,j) - t_base(k) ENDDO ENDDO ENDDO END IF CALL vertical_diffusion_s( rt_tendf, config_flags, var_mix, mu, xkhv, & dn, dnw, rdz, rdzw, fnm, fnp, & .false., & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) !***************************************** !***************************************** !MODIFICA al flusso di calore ! ! hflux: SELECT CASE( config_flags%isfflx ) CASE (0,2) ! with fixed surface heat flux given in the namelist heat_flux = config_flags%tke_heat_flux ! constant heat flux value ! set in namelist.input DO j = j_start, j_end DO i = i_start, i_end rt_tendf(i,kts,j)=rt_tendf(i,kts,j) & +mu(i,j)*heat_flux*rdzw(i,kts,j) ENDDO ENDDO CASE (1) ! use surface heat flux computed from surface routine DO j = j_start, j_end DO i = i_start, i_end cpm = cp * (1. + 0.8 * moist(i,kts,j,P_QV)) heat_flux = hfx(i,j)/cpm/rho(i,1,j) rt_tendf(i,kts,j)=rt_tendf(i,kts,j) & +mu(i,j)*heat_flux*rdzw(i,kts,j) ENDDO ENDDO CASE DEFAULT CALL wrf_error_fatal( 'isfflx value invalid for diff_opt=2' ) END SELECT hflux ! ! FINE MODIFICA al flusso di calore !***************************************** !***************************************** If (km_opt .eq. 2) then CALL vertical_diffusion_s( tke_tendf(ims,kms,jms), & config_flags, tke(ims,kms,jms), & mu, xkhv, & dn, dnw, rdz, rdzw, fnm, fnp, & .true., & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) endif IF (n_moist .ge. PARAM_FIRST_SCALAR) THEN moist_loop: do im = PARAM_FIRST_SCALAR, n_moist IF ( (.not. config_flags%mix_full_fields) .and. (im == P_QV) ) THEN DO j=jts,min(jte,jde-1) DO k=kts,kte-1 DO i=its,min(ite,ide-1) var_mix(i,k,j) = moist(i,k,j,im) - qv_base(k) ENDDO ENDDO ENDDO ELSE DO j=jts,min(jte,jde-1) DO k=kts,kte-1 DO i=its,min(ite,ide-1) var_mix(i,k,j) = moist(i,k,j,im) ENDDO ENDDO ENDDO END IF CALL vertical_diffusion_s( moist_tendf(ims,kms,jms,im), & config_flags, var_mix, & mu, xkhv, & dn, dnw, rdz, rdzw, fnm, fnp, & .false., & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) !***************************************** !***************************************** !MODIFICATIONS for water vapor flux ! ! qflux: SELECT CASE( config_flags%isfflx ) CASE (0) ! do nothing CASE (1,2) ! with surface moisture flux IF ( im == P_QV ) THEN DO j = j_start, j_end DO i = i_start, i_end moist_flux = qfx(i,j)/rho(i,1,j) moist_tendf(i,kts,j,im)=moist_tendf(i,kts,j,im) & +mu(i,j)*moist_flux*rdzw(i,kts,j) ENDDO ENDDO ENDIF CASE DEFAULT CALL wrf_error_fatal( 'isfflx value invalid for diff_opt=2' ) END SELECT qflux ! ! END MODIFICATIONS for water vapor flux !***************************************** !***************************************** ENDDO moist_loop ENDIF IF (n_chem .ge. PARAM_FIRST_SCALAR) THEN chem_loop: do im = PARAM_FIRST_SCALAR, n_chem CALL vertical_diffusion_s( chem_tendf(ims,kms,jms,im), & config_flags, chem(ims,kms,jms,im), & mu, xkhv, & dn, dnw, rdz, rdzw, fnm, fnp, & .false., & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) ENDDO chem_loop ENDIF IF (n_scalar .ge. PARAM_FIRST_SCALAR) THEN scalar_loop: do im = PARAM_FIRST_SCALAR, n_scalar CALL vertical_diffusion_s( scalar_tendf(ims,kms,jms,im), & config_flags, scalar(ims,kms,jms,im), & mu, xkhv, & dn, dnw, rdz, rdzw, fnm, fnp, & .false., & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) ENDDO scalar_loop ENDIF END SUBROUTINE vertical_diffusion_2 !======================================================================= !======================================================================= SUBROUTINE vertical_diffusion_u_2( tendency, config_flags, mu, & defor13, xkmv, & dnw, rdzw, fnm, fnp, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) !----------------------------------------------------------------------- ! Begin declarations. IMPLICIT NONE 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( kms:kme ) , INTENT(IN ) :: fnm REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fnp REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: dnw ! REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: zeta_z REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(INOUT) ::tendency REAL , DIMENSION( ims:ime , kms:kme, jms:jme ) , & INTENT(IN ) ::defor13, & xkmv, & rdzw REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mu ! LOCAL VARS INTEGER :: i, j, k, ktf INTEGER :: i_start, i_end, j_start, j_end INTEGER :: is_ext,ie_ext,js_ext,je_ext REAL , DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1) :: titau3 REAL , DIMENSION( its:ite, jts:jte) :: zzavg REAL :: rdzu ! End declarations. !----------------------------------------------------------------------- ktf=MIN(kte,kde-1) i_start = its i_end = ite j_start = jts j_end = MIN(jte,jde-1) IF ( config_flags%open_xs .or. config_flags%specified .or. & config_flags%nested) i_start = MAX(ids+1,its) IF ( config_flags%open_xe .or. config_flags%specified .or. & config_flags%nested) i_end = MIN(ide-1,ite) IF ( config_flags%open_ys .or. config_flags%specified .or. & config_flags%nested) j_start = MAX(jds+1,jts) IF ( config_flags%open_ye .or. config_flags%specified .or. & config_flags%nested) j_end = MIN(jde-2,jte) IF ( config_flags%periodic_x ) i_start = its IF ( config_flags%periodic_x ) i_end = ite ! titau3 = titau13 is_ext=0 ie_ext=0 js_ext=0 je_ext=0 CALL cal_titau_13_31( config_flags, titau3, defor13, & mu, xkmv, fnm, fnp, & is_ext, ie_ext, js_ext, je_ext, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) ! DO j = j_start, j_end DO k=kts+1,ktf DO i = i_start, i_end rdzu = 2./(1./rdzw(i,k,j) + 1./rdzw(i-1,k,j)) tendency(i,k,j)=tendency(i,k,j)-rdzu*(titau3(i,k+1,j)-titau3(i,k,j)) ENDDO ENDDO ENDDO ! ******** MODIF... ! we will pick up the surface drag (titau3(i,kts,j)) later ! DO j = j_start, j_end k=kts DO i = i_start, i_end rdzu = 2./(1./rdzw(i,k,j) + 1./rdzw(i-1,k,j)) tendency(i,k,j)=tendency(i,k,j)-rdzu*(titau3(i,k+1,j)) ENDDO ENDDO ! ******** MODIF... END SUBROUTINE vertical_diffusion_u_2 !======================================================================= !======================================================================= SUBROUTINE vertical_diffusion_v_2( tendency, config_flags, mu, & defor23, xkmv, & dnw, rdzw, fnm, fnp, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) !----------------------------------------------------------------------- ! Begin declarations. IMPLICIT NONE 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( kms:kme ) , INTENT(IN ) :: fnm REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fnp REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: dnw ! REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: zeta_z REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(INOUT) ::tendency REAL , DIMENSION( ims:ime , kms:kme, jms:jme ) , & INTENT(IN ) ::defor23, & xkmv, & rdzw REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mu ! LOCAL VARS INTEGER :: i, j, k, ktf INTEGER :: i_start, i_end, j_start, j_end INTEGER :: is_ext,ie_ext,js_ext,je_ext REAL , DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1) :: titau3 REAL , DIMENSION( its:ite, jts:jte) :: zzavg REAL :: rdzv ! End declarations. !----------------------------------------------------------------------- ktf=MIN(kte,kde-1) i_start = its i_end = MIN(ite,ide-1) j_start = jts j_end = jte IF ( config_flags%open_xs .or. config_flags%specified .or. & config_flags%nested) i_start = MAX(ids+1,its) IF ( config_flags%open_xe .or. config_flags%specified .or. & config_flags%nested) i_end = MIN(ide-2,ite) IF ( config_flags%open_ys .or. config_flags%specified .or. & config_flags%nested) j_start = MAX(jds+1,jts) IF ( config_flags%open_ye .or. config_flags%specified .or. & config_flags%nested) j_end = MIN(jde-1,jte) IF ( config_flags%periodic_x ) i_start = its IF ( config_flags%periodic_x ) i_end = MIN(ite,ide-1) ! titau3 = titau23 is_ext=0 ie_ext=0 js_ext=0 je_ext=0 CALL cal_titau_23_32( config_flags, titau3, defor23, & mu, xkmv, fnm, fnp, & is_ext, ie_ext, js_ext, je_ext, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) DO j = j_start, j_end DO k = kts+1,ktf DO i = i_start, i_end rdzv = 2./(1./rdzw(i,k,j) + 1./rdzw(i,k,j-1)) tendency(i,k,j)=tendency(i,k,j)-rdzv*(titau3(i,k+1,j)-titau3(i,k,j)) ENDDO ENDDO ENDDO ! ******** MODIF... ! we will pick up the surface drag (titau3(i,kts,j)) later ! DO j = j_start, j_end k=kts DO i = i_start, i_end rdzv = 2./(1./rdzw(i,k,j) + 1./rdzw(i,k,j-1)) tendency(i,k,j)=tendency(i,k,j)-rdzv*(titau3(i,k+1,j)) ENDDO ENDDO ! ******** MODIF... END SUBROUTINE vertical_diffusion_v_2 !======================================================================= !======================================================================= SUBROUTINE vertical_diffusion_w_2(tendency, config_flags, mu, & defor33, tke, div, xkmv, & dn, rdz, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) !----------------------------------------------------------------------- ! Begin declarations. IMPLICIT NONE 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( kms:kme ) , INTENT(IN ) :: dn REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(INOUT) ::tendency REAL , DIMENSION( ims:ime , kms:kme, jms:jme ) , & INTENT(IN ) ::defor33, & tke, & div, & xkmv, & rdz REAL , DIMENSION( ims:ime, jms:jme), INTENT(IN ) :: mu ! LOCAL VARS INTEGER :: i, j, k, ktf INTEGER :: i_start, i_end, j_start, j_end INTEGER :: is_ext,ie_ext,js_ext,je_ext REAL , DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1) :: titau3 ! End declarations. !----------------------------------------------------------------------- ktf=MIN(kte,kde-1) i_start = its i_end = MIN(ite,ide-1) j_start = jts j_end = MIN(jte,jde-1) IF ( config_flags%open_xs .or. config_flags%specified .or. & config_flags%nested) i_start = MAX(ids+1,its) IF ( config_flags%open_xe .or. config_flags%specified .or. & config_flags%nested) i_end = MIN(ide-2,ite) IF ( config_flags%open_ys .or. config_flags%specified .or. & config_flags%nested) j_start = MAX(jds+1,jts) IF ( config_flags%open_ye .or. config_flags%specified .or. & config_flags%nested) j_end = MIN(jde-2,jte) IF ( config_flags%periodic_x ) i_start = its IF ( config_flags%periodic_x ) i_end = MIN(ite,ide-1) ! titau3 = titau33 is_ext=0 ie_ext=0 js_ext=0 je_ext=0 CALL cal_titau_11_22_33( config_flags, titau3, & mu, tke, xkmv, defor33, & is_ext, ie_ext, js_ext, je_ext, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) ! DO j = j_start, j_end ! DO k = kts+1, ktf ! DO i = i_start, i_end ! titau3(i,k,j)=titau3(i,k,j)*zeta_z(i,j) ! ENDDO ! ENDDO ! ENDDO DO j = j_start, j_end DO k = kts+1, ktf DO i = i_start, i_end tendency(i,k,j)=tendency(i,k,j)-rdz(i,k,j)*(titau3(i,k,j)-titau3(i,k-1,j)) ENDDO ENDDO ENDDO END SUBROUTINE vertical_diffusion_w_2 !======================================================================= !======================================================================= SUBROUTINE vertical_diffusion_s( tendency, config_flags, var, mu, xkhv, & dn, dnw, rdz, rdzw, fnm, fnp, & doing_tke, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) !----------------------------------------------------------------------- ! Begin declarations. IMPLICIT NONE 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 LOGICAL, INTENT(IN ) :: doing_tke REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fnm REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fnp REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: dn REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: dnw REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(INOUT) ::tendency REAL , DIMENSION( ims:ime , kms:kme, jms:jme ) , INTENT(IN) :: xkhv REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN) :: mu REAL , DIMENSION( ims:ime , kms:kme, jms:jme ) , & INTENT(IN ) :: var, & rdz, & rdzw ! LOCAL VARS INTEGER :: i, j, k, ktf INTEGER :: i_start, i_end, j_start, j_end REAL , DIMENSION( its:ite, kts:kte, jts:jte) :: H3, & xkxavg, & rravg REAL , DIMENSION( its:ite, kts:kte, jts:jte) :: tmptendf ! End declarations. !----------------------------------------------------------------------- ktf=MIN(kte,kde-1) i_start = its i_end = MIN(ite,ide-1) j_start = jts j_end = MIN(jte,jde-1) IF ( config_flags%open_xs .or. config_flags%specified .or. & config_flags%nested) i_start = MAX(ids+1,its) IF ( config_flags%open_xe .or. config_flags%specified .or. & config_flags%nested) i_end = MIN(ide-2,ite) IF ( config_flags%open_ys .or. config_flags%specified .or. & config_flags%nested) j_start = MAX(jds+1,jts) IF ( config_flags%open_ye .or. config_flags%specified .or. & config_flags%nested) j_end = MIN(jde-2,jte) IF ( config_flags%periodic_x ) i_start = its IF ( config_flags%periodic_x ) i_end = MIN(ite,ide-1) IF (doing_tke) THEN DO j = j_start, j_end DO k = kts,ktf DO i = i_start, i_end tmptendf(i,k,j)=tendency(i,k,j) ENDDO ENDDO ENDDO ENDIF ! H3 xkxavg = 0. DO j = j_start, j_end DO k = kts+1,ktf DO i = i_start, i_end xkxavg(i,k,j)=fnm(k)*xkhv(i,k,j)+fnp(k)*xkhv(i,k-1,j) H3(i,k,j)=-xkxavg(i,k,j)*(var(i,k,j)-var(i,k-1,j))*rdz(i,k,j) ! H3(i,k,j)=-xkxavg(i,k,j)*zeta_z(i,j)* & ! (var(i,k,j)-var(i,k-1,j))/dn(k) ENDDO ENDDO ENDDO DO j = j_start, j_end DO i = i_start, i_end H3(i,kts,j)=0. H3(i,ktf+1,j)=0. ! H3(i,kts,j)=H3(i,kts+1,j) ! H3(i,ktf+1,j)=H3(i,ktf,j) ENDDO ENDDO DO j = j_start, j_end DO k = kts,ktf DO i = i_start, i_end tendency(i,k,j)=tendency(i,k,j) & -mu(i,j)*(H3(i,k+1,j)-H3(i,k,j))*rdzw(i,k,j) ENDDO ENDDO ENDDO IF (doing_tke) THEN DO j = j_start, j_end DO k = kts,ktf DO i = i_start, i_end tendency(i,k,j)=tmptendf(i,k,j)+2.* & (tendency(i,k,j)-tmptendf(i,k,j)) ENDDO ENDDO ENDDO ENDIF END SUBROUTINE vertical_diffusion_s !======================================================================= !======================================================================= SUBROUTINE cal_titau_11_22_33( config_flags, titau, & mu, tke, xkx, defor, & is_ext, ie_ext, js_ext, je_ext, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) ! History: Sep 2003 Changes by George Bryan and Jason Knievel, NCAR ! Oct 2001 Converted to mass core by Bill Skamarock, NCAR ! Aug 2000 Original code by Shu-Hua Chen, UC-Davis ! Purpose: This routine calculates stress terms (taus) for use in ! the calculation of production of TKE by sheared wind ! References: Klemp and Wilhelmson (JAS 1978) ! Deardorff (B-L Meteor 1980) ! Chen and Dudhia (NCAR WRF physics report 2000) ! Key: !----------------------------------------------------------------------- ! Begin declarations. IMPLICIT NONE 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 INTEGER, INTENT( IN ) & :: is_ext, ie_ext, js_ext, je_ext REAL, DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1 ), INTENT( INOUT ) & :: titau REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( IN ) & :: defor, xkx, tke REAL, DIMENSION( ims:ime, jms:jme ), INTENT( IN ) & :: mu ! Local variables. INTEGER & :: i, j, k, ktf, i_start, i_end, j_start, j_end ! End declarations. !----------------------------------------------------------------------- ktf = MIN( kte, kde-1 ) i_start = its i_end = ite j_start = jts j_end = jte IF ( config_flags%open_xs .OR. config_flags%specified .OR. & config_flags%nested) i_start = MAX( ids+1, its ) IF ( config_flags%open_xe .OR. config_flags%specified .OR. & config_flags%nested) i_end = MIN( ide-1, ite ) IF ( config_flags%open_ys .OR. config_flags%specified .OR. & config_flags%nested) j_start = MAX( jds+1, jts ) IF ( config_flags%open_ye .OR. config_flags%specified .OR. & config_flags%nested) j_end = MIN( jde-1, jte ) IF ( config_flags%periodic_x ) i_start = its IF ( config_flags%periodic_x ) i_end = ite i_start = i_start - is_ext i_end = i_end + ie_ext j_start = j_start - js_ext j_end = j_end + je_ext IF ( config_flags%km_opt .EQ. 2) THEN DO j = j_start,j_end DO k = kts,ktf DO i = i_start,i_end titau(i,k,j) = mu(i,j) * ( - xkx(i,k,j) * ( defor(i,k,j) ) ) END DO END DO END DO ELSE DO j = j_start, j_end DO k = kts, ktf DO i = i_start, i_end titau(i,k,j) = - mu(i,j) * xkx(i,k,j) * defor(i,k,j) END DO END DO END DO END IF END SUBROUTINE cal_titau_11_22_33 !======================================================================= !======================================================================= SUBROUTINE cal_titau_12_21( config_flags, titau, & mu, xkx, defor, & is_ext, ie_ext, js_ext, je_ext, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) ! History: Sep 2003 Modifications by George Bryan and Jason Knievel, NCAR ! Oct 2001 Converted to mass core by Bill Skamarock, NCAR ! Aug 2000 Original code by Shu-Hua Chen, UC-Davis ! Pusrpose This routine calculates the stress terms (taus) for use in ! the calculation of production of TKE by sheared wind ! References: Klemp and Wilhelmson (JAS 1978) ! Deardorff (B-L Meteor 1980) ! Chen and Dudhia (NCAR WRF physics report 2000) ! Key: !----------------------------------------------------------------------- ! Begin declarations. IMPLICIT NONE 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 INTEGER, INTENT( IN ) & :: is_ext, ie_ext, js_ext, je_ext REAL, DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1 ), INTENT( INOUT ) & :: titau REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( IN ) & :: defor, xkx REAL, DIMENSION( ims:ime, jms:jme ), INTENT( IN ) & :: mu ! Local variables. INTEGER & :: i, j, k, ktf, i_start, i_end, j_start, j_end REAL, DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1 ) & :: xkxavg REAL, DIMENSION( its-1:ite+1, jts-1:jte+1 ) & :: muavg ! End declarations. !----------------------------------------------------------------------- ktf = MIN( kte, kde-1 ) ! Needs one more point in the x and y directions. i_start = its i_end = ite j_start = jts j_end = jte IF ( config_flags%open_xs .OR. config_flags%specified .OR. & config_flags%nested ) i_start = MAX( ids+1, its ) IF ( config_flags%open_xe .OR. config_flags%specified .OR. & config_flags%nested ) i_end = MIN( ide-1, ite ) IF ( config_flags%open_ys .OR. config_flags%specified .OR. & config_flags%nested ) j_start = MAX( jds+1, jts ) IF ( config_flags%open_ye .OR. config_flags%specified .OR. & config_flags%nested ) j_end = MIN( jde-1, jte ) IF ( config_flags%periodic_x ) i_start = its IF ( config_flags%periodic_x ) i_end = ite i_start = i_start - is_ext i_end = i_end + ie_ext j_start = j_start - js_ext j_end = j_end + je_ext DO j = j_start, j_end DO k = kts, ktf DO i = i_start, i_end xkxavg(i,k,j) = 0.25 * ( xkx(i-1,k,j ) + xkx(i,k,j ) + & xkx(i-1,k,j-1) + xkx(i,k,j-1) ) END DO END DO END DO DO j = j_start, j_end DO i = i_start, i_end muavg(i,j) = 0.25 * ( mu(i-1,j ) + mu(i,j ) + & mu(i-1,j-1) + mu(i,j-1) ) END DO END DO ! titau12 or titau21 DO j = j_start, j_end DO k = kts, ktf DO i = i_start, i_end titau(i,k,j) = - muavg(i,j) * xkxavg(i,k,j) * defor(i,k,j) END DO END DO END DO END SUBROUTINE cal_titau_12_21 !======================================================================= SUBROUTINE cal_titau_13_31( config_flags, titau, & defor, mu, xkx, fnm, fnp, & is_ext, ie_ext, js_ext, je_ext, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) ! History: Sep 2003 Modifications by George Bryan and Jason Knievel, NCAR ! Oct 2001 Converted to mass core by Bill Skamarock, NCAR ! Aug 2000 Original code by Shu-Hua Chen, UC-Davis ! Purpose: This routine calculates the stress terms (taus) for use in ! the calculation of production of TKE by sheared wind ! References: Klemp and Wilhelmson (JAS 1978) ! Deardorff (B-L Meteor 1980) ! Chen and Dudhia (NCAR WRF physics report 2000) ! Key: !----------------------------------------------------------------------- ! Begin declarations. IMPLICIT NONE 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 INTEGER, INTENT( IN ) & :: is_ext, ie_ext, js_ext, je_ext REAL, DIMENSION( kms:kme ), INTENT( IN ) & :: fnm, fnp REAL, DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1 ), INTENT( INOUT ) & :: titau REAL, DIMENSION( ims:ime, kms:kme, jms:jme), INTENT( IN ) & :: defor, xkx REAL, DIMENSION( ims:ime, jms:jme), INTENT( IN ) & :: mu ! Local variables. INTEGER & :: i, j, k, ktf, i_start, i_end, j_start, j_end REAL, DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1 ) & :: xkxavg REAL, DIMENSION( its-1:ite+1, jts-1:jte+1 ) & :: muavg ! End declarations. !----------------------------------------------------------------------- ktf = MIN( kte, kde-1 ) ! Find ide-1 and jde-1 for averaging to p point. i_start = its i_end = ite j_start = jts j_end = MIN( jte, jde-1 ) IF ( config_flags%open_xs .OR. config_flags%specified .OR. & config_flags%nested) i_start = MAX( ids+1, its ) IF ( config_flags%open_xe .OR. config_flags%specified .OR. & config_flags%nested) i_end = MIN( ide-1, ite ) IF ( config_flags%open_ys .OR. config_flags%specified .OR. & config_flags%nested) j_start = MAX( jds+1, jts ) IF ( config_flags%open_ye .OR. config_flags%specified .OR. & config_flags%nested) j_end = MIN( jde-2, jte ) IF ( config_flags%periodic_x ) i_start = its IF ( config_flags%periodic_x ) i_end = ite i_start = i_start - is_ext i_end = i_end + ie_ext j_start = j_start - js_ext j_end = j_end + je_ext DO j = j_start, j_end DO k = kts+1, ktf DO i = i_start, i_end xkxavg(i,k,j) = 0.5 * ( fnm(k) * ( xkx(i,k ,j) + xkx(i-1,k ,j) ) + & fnp(k) * ( xkx(i,k-1,j) + xkx(i-1,k-1,j) ) ) END DO END DO END DO DO j = j_start, j_end DO i = i_start, i_end muavg(i,j) = 0.5 * ( mu(i,j) + mu(i-1,j) ) END DO END DO DO j = j_start, j_end DO k = kts+1, ktf DO i = i_start, i_end titau(i,k,j) = - muavg(i,j) * xkxavg(i,k,j) * defor(i,k,j) ENDDO ENDDO ENDDO DO j = j_start, j_end DO i = i_start, i_end titau(i,kts ,j) = 0.0 titau(i,ktf+1,j) = 0.0 ENDDO ENDDO END SUBROUTINE cal_titau_13_31 !======================================================================= !======================================================================= SUBROUTINE cal_titau_23_32( config_flags, titau, defor, & mu, xkx, fnm, fnp, & is_ext, ie_ext, js_ext, je_ext, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) ! History: Sep 2003 Changes by George Bryan and Jason Knievel, NCAR ! Oct 2001 Converted to mass core by Bill Skamarock, NCAR ! Aug 2000 Original code by Shu-Hua Chen, UC-Davis ! Purpose: This routine calculates stress terms (taus) for use in ! the calculation of production of TKE by sheared wind ! References: Klemp and Wilhelmson (JAS 1978) ! Deardorff (B-L Meteor 1980) ! Chen and Dudhia (NCAR WRF physics report 2000) ! Key: !----------------------------------------------------------------------- ! Begin declarations. IMPLICIT NONE 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 INTEGER, INTENT( IN ) & :: is_ext,ie_ext,js_ext,je_ext REAL, DIMENSION( kms:kme ), INTENT( IN ) & :: fnm, fnp REAL, DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1 ), INTENT( INOUT ) & :: titau REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( IN ) & :: defor, xkx REAL, DIMENSION( ims:ime, jms:jme ), INTENT( IN ) & :: mu ! Local variables. INTEGER & :: i, j, k, ktf, i_start, i_end, j_start, j_end REAL, DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1 ) & :: xkxavg REAL, DIMENSION( its-1:ite+1, jts-1:jte+1 ) & :: muavg ! End declarations. !----------------------------------------------------------------------- ktf = MIN( kte, kde-1 ) ! Find ide-1 and jde-1 for averaging to p point. i_start = its i_end = MIN( ite, ide-1 ) j_start = jts j_end = jte IF ( config_flags%open_xs .OR. config_flags%specified .OR. & config_flags%nested) i_start = MAX( ids+1, its ) IF ( config_flags%open_xe .OR. config_flags%specified .OR. & config_flags%nested) i_end = MIN( ide-2, ite ) IF ( config_flags%open_ys .OR. config_flags%specified .OR. & config_flags%nested) j_start = MAX( jds+1, jts ) IF ( config_flags%open_ye .OR. config_flags%specified .OR. & config_flags%nested) j_end = MIN( jde-1, jte ) IF ( config_flags%periodic_x ) i_start = its IF ( config_flags%periodic_x ) i_end = MIN( ite, ide-1 ) i_start = i_start - is_ext i_end = i_end + ie_ext j_start = j_start - js_ext j_end = j_end + je_ext DO j = j_start, j_end DO k = kts+1, ktf DO i = i_start, i_end xkxavg(i,k,j) = 0.5 * ( fnm(k) * ( xkx(i,k ,j) + xkx(i,k ,j-1) ) + & fnp(k) * ( xkx(i,k-1,j) + xkx(i,k-1,j-1) ) ) END DO END DO END DO DO j = j_start, j_end DO i = i_start, i_end muavg(i,j) = 0.5 * ( mu(i,j) + mu(i,j-1) ) END DO END DO DO j = j_start, j_end DO k = kts+1, ktf DO i = i_start, i_end titau(i,k,j) = - muavg(i,j) * xkxavg(i,k,j) * defor(i,k,j) END DO END DO END DO DO j = j_start, j_end DO i = i_start, i_end titau(i,kts ,j) = 0.0 titau(i,ktf+1,j) = 0.0 END DO END DO END SUBROUTINE cal_titau_23_32 !======================================================================= !======================================================================= SUBROUTINE phy_bc ( config_flags,div,defor11,defor22,defor33, & defor12,defor13,defor23,xkmh,xkmv,xkhh,xkhv,tke, & RUBLTEN, RVBLTEN, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & ips, ipe, jps, jpe, kps, kpe, & its, ite, jts, jte, kts, kte ) !------------------------------------------------------------------------------ ! Begin declarations. IMPLICIT NONE TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & ips, ipe, jps, jpe, kps, kpe, & its, ite, jts, jte, kts, kte REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(INOUT) ::RUBLTEN, & RVBLTEN, & defor11, & defor22, & defor33, & defor12, & defor13, & defor23, & xkmh, & xkmv, & xkhh, & xkhv, & tke, & div ! End declarations. !----------------------------------------------------------------------- IF( (config_flags%bl_pbl_physics .GT. 0) & .OR. (config_flags%modif_wrf) ) THEN !****MARS: Always do that ?? > yes, for periodic simu -- grep in share CALL set_physical_bc3d( RUBLTEN , 't', config_flags, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & ips, ipe, jps, jpe, kps, kpe, & its, ite, jts, jte, kts, kte ) CALL set_physical_bc3d( RVBLTEN , 't', config_flags, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & ips, ipe, jps, jpe, kps, kpe, & its, ite, jts, jte, kts, kte ) ENDIF ! move out of the conditional, below; horiz coeffs needed for ! all diff_opt cases. JM CALL set_physical_bc3d( xkmh , 't', config_flags, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & ips, ipe, jps, jpe, kps, kpe, & its, ite, jts, jte, kts, kte ) CALL set_physical_bc3d( xkhh , 't', config_flags, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & ips, ipe, jps, jpe, kps, kpe, & its, ite, jts, jte, kts, kte ) IF(config_flags%diff_opt .eq. 2) THEN CALL set_physical_bc3d( xkmv , 't', config_flags, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & ips, ipe, jps, jpe, kps, kpe, & its, ite, jts, jte, kts, kte ) CALL set_physical_bc3d( xkhv , 't', config_flags, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & ips, ipe, jps, jpe, kps, kpe, & its, ite, jts, jte, kts, kte ) CALL set_physical_bc3d( tke , 't', config_flags, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & ips, ipe, jps, jpe, kps, kpe, & its, ite, jts, jte, kts, kte ) CALL set_physical_bc3d( div , 't', config_flags, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & ips, ipe, jps, jpe, kps, kpe, & its, ite, jts, jte, kts, kte ) CALL set_physical_bc3d( defor11 , 't', config_flags, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & ips, ipe, jps, jpe, kps, kpe, & its, ite, jts, jte, kts, kte ) CALL set_physical_bc3d( defor22 , 't', config_flags, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & ips, ipe, jps, jpe, kps, kpe, & its, ite, jts, jte, kts, kte ) CALL set_physical_bc3d( defor33 , 't', config_flags, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & ips, ipe, jps, jpe, kps, kpe, & its, ite, jts, jte, kts, kte ) CALL set_physical_bc3d( defor12 , 'd', config_flags, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & ips, ipe, jps, jpe, kps, kpe, & its, ite, jts, jte, kts, kte ) CALL set_physical_bc3d( defor13 , 'e', config_flags, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & ips, ipe, jps, jpe, kps, kpe, & its, ite, jts, jte, kts, kte ) CALL set_physical_bc3d( defor23 , 'f', config_flags, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & ips, ipe, jps, jpe, kps, kpe, & its, ite, jts, jte, kts, kte ) ENDIF END SUBROUTINE phy_bc !======================================================================= !======================================================================= SUBROUTINE tke_rhs( tendency, BN2, config_flags, & defor11, defor22, defor33, & defor12, defor13, defor23, & u, v, w, div, tke, mu, & theta, p, p8w, t8w, z, fnm, fnp, & cf1, cf2, cf3, msftx, msfty, & xkmh, xkmv, xkhv, & rdx, rdy, dx, dy, dt, zx, zy, & rdz, rdzw, dn, dnw, isotropic, & hfx, qfx, qv, ust, rho, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) !----------------------------------------------------------------------- ! Begin declarations. IMPLICIT NONE 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 INTEGER, INTENT( IN ) :: isotropic REAL, INTENT( IN ) & :: cf1, cf2, cf3, dt, rdx, rdy, dx, dy REAL, DIMENSION( kms:kme ), INTENT( IN ) & :: fnm, fnp, dnw, dn REAL, DIMENSION( ims:ime, jms:jme ), INTENT( IN ) & :: msftx, msfty REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( INOUT ) & :: tendency REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( IN ) & :: defor11, defor22, defor33, defor12, defor13, defor23, & div, BN2, tke, xkmh, xkmv, xkhv, zx, zy, u, v, w, theta, & p, p8w, t8w, z, rdz, rdzw REAL, DIMENSION( ims:ime, jms:jme ), INTENT( IN ) & :: mu REAL, DIMENSION ( ims:ime, jms:jme ), INTENT( IN ) & :: hfx, ust, qfx REAL, DIMENSION ( ims:ime, kms:kme, jms:jme ), INTENT ( IN ) & :: qv, rho ! Local variables. INTEGER & :: i, j, k, ktf, i_start, i_end, j_start, j_end ! End declarations. !----------------------------------------------------------------------- CALL tke_shear( tendency, config_flags, & defor11, defor22, defor33, & defor12, defor13, defor23, & u, v, w, tke, ust, mu, fnm, fnp, & cf1, cf2, cf3, msftx, msfty, & xkmh, xkmv, & rdx, rdy, zx, zy, rdz, rdzw, dnw, dn, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) CALL tke_buoyancy( tendency, config_flags, mu, & tke, xkhv, BN2, theta, dt, & hfx, qfx, qv, rho, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) CALL tke_dissip( tendency, config_flags, & mu, tke, bn2, theta, p8w, t8w, z, & dx, dy,rdz, rdzw, isotropic, & msftx, msfty, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) ! Set a lower limit on TKE. ktf = MIN( kte, kde-1 ) i_start = its i_end = MIN( ite, ide-1 ) j_start = jts j_end = MIN( jte, jde-1 ) IF ( config_flags%open_xs .or. config_flags%specified .or. & config_flags%nested) i_start = MAX(ids+1,its) IF ( config_flags%open_xe .or. config_flags%specified .or. & config_flags%nested) i_end = MIN(ide-2,ite) IF ( config_flags%open_ys .or. config_flags%specified .or. & config_flags%nested) j_start = MAX(jds+1,jts) IF ( config_flags%open_ye .or. config_flags%specified .or. & config_flags%nested) j_end = MIN(jde-2,jte) IF ( config_flags%periodic_x ) i_start = its IF ( config_flags%periodic_x ) i_end = MIN( ite, ide-1 ) DO j = j_start, j_end DO k = kts, ktf DO i = i_start, i_end tendency(i,k,j) = max( tendency(i,k,j), -mu(i,j) * max( 0.0 , tke(i,k,j) ) / dt ) END DO END DO END DO END SUBROUTINE tke_rhs !======================================================================= !======================================================================= SUBROUTINE tke_buoyancy( tendency, config_flags, mu, & tke, xkhv, BN2, theta, dt, & hfx, qfx, qv, rho, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) !----------------------------------------------------------------------- ! Begin declarations. IMPLICIT NONE 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, INTENT( IN ) & :: dt REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( INOUT ) & :: tendency REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( IN ) & :: xkhv, tke, BN2, theta REAL, DIMENSION( ims:ime, jms:jme ), INTENT( IN ) & :: mu REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT ( IN ) & :: qv, rho REAL, DIMENSION(ims:ime, jms:jme ), INTENT ( IN ) :: hfx, qfx ! Local variables. INTEGER & :: i, j, k, ktf INTEGER & :: i_start, i_end, j_start, j_end REAL :: heat_flux, heat_flux0 REAL :: cpm ! End declarations. !----------------------------------------------------------------------- !----------------------------------------------------------------------- ! Add to the TKE tendency the term that accounts for production of TKE ! due to buoyant motions. ktf = MIN( kte, kde-1 ) i_start = its i_end = MIN( ite, ide-1 ) j_start = jts j_end = MIN( jte, jde-1 ) IF ( config_flags%open_xs .OR. config_flags%specified .OR. & config_flags%nested ) i_start = MAX( ids+1, its ) IF ( config_flags%open_xe .OR. config_flags%specified .OR. & config_flags%nested ) i_end = MIN( ide-2, ite ) IF ( config_flags%open_ys .OR. config_flags%specified .OR. & config_flags%nested ) j_start = MAX( jds+1, jts ) IF ( config_flags%open_ye .OR. config_flags%specified .OR. & config_flags%nested ) j_end = MIN( jde-2, jte ) IF ( config_flags%periodic_x ) i_start = its IF ( config_flags%periodic_x ) i_end = MIN( ite, ide-1 ) DO j = j_start, j_end DO k = kts+1, ktf DO i = i_start, i_end tendency(i,k,j) = tendency(i,k,j) - mu(i,j) * xkhv(i,k,j) * BN2(i,k,j) END DO END DO END DO ! MARTA: change in the computation of the tke's tendency at the surface. ! the buoyancy flux is the average of the surface heat flux (0.06) and the ! flux at the first w level ! ! WCS 040331 hflux: SELECT CASE( config_flags%isfflx ) CASE (0,2) ! with fixed surface heat flux given in the namelist heat_flux0 = config_flags%tke_heat_flux ! constant heat flux value ! set in namelist.input ! LES mods K=KTS PRINT *, 'USE constant H/rho/cp : e.g. value of Hs ', heat_flux0*cp*rho(i_start,k,j_start) DO j = j_start, j_end DO i = i_start, i_end heat_flux = heat_flux0 tendency(i,k,j)= tendency(i,k,j) - & mu(i,j)*((xkhv(i,k,j)*BN2(i,k,j))- (g/theta(i,k,j))*heat_flux)/2. ENDDO ENDDO CASE (1) ! use surface heat flux computed from surface routine PRINT *, 'USE surface heat flux from PHYSICS', hfx(i_start,j_start) K=KTS DO j = j_start, j_end DO i = i_start, i_end cpm = cp * (1. + 0.8*qv(i,k,j)) heat_flux = (hfx(i,j)/cpm)/rho(i,k,j) tendency(i,k,j)= tendency(i,k,j) - & mu(i,j)*((xkhv(i,k,j)*BN2(i,k,j))- (g/theta(i,k,j))*heat_flux)/2. ENDDO ENDDO CASE DEFAULT CALL wrf_error_fatal( 'isfflx value invalid for diff_opt=2' ) END SELECT hflux ! end of MARTA/WCS change ! The tendency array now includes production of TKE from buoyant ! motions. !----------------------------------------------------------------------- END SUBROUTINE tke_buoyancy !======================================================================= !======================================================================= SUBROUTINE tke_dissip( tendency, config_flags, & mu, tke, bn2, theta, p8w, t8w, z, & dx, dy, rdz, rdzw, isotropic, & msftx, msfty, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) ! History: Sep 2003 Changes by George Bryan and Jason Knievel, NCAR ! Oct 2001 Converted to mass core by Bill Skamarock, NCAR ! Aug 2000 Original code by Shu-Hua Chen, UC-Davis ! Purpose: This routine calculates dissipation of turbulent kinetic ! energy. ! References: Klemp and Wilhelmson (JAS 1978) ! Deardorff (B-L Meteor 1980) ! Chen and Dudhia (NCAR WRF physics report 2000) !----------------------------------------------------------------------- ! Begin declarations. IMPLICIT NONE 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 INTEGER, INTENT( IN ) :: isotropic REAL, INTENT( IN ) & :: dx, dy REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( INOUT ) & :: tendency REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( IN ) & :: tke, bn2, theta, p8w, t8w, z, rdz, rdzw REAL, DIMENSION( ims:ime, jms:jme ), INTENT( IN ) & :: mu REAL, DIMENSION( ims:ime, jms:jme ), INTENT( IN ) & :: msftx, msfty ! Local variables. REAL, DIMENSION( its:ite, kts:kte, jts:jte ) & :: dthrdn REAL, DIMENSION( its:ite, kts:kte, jts:jte ) & :: l_scale REAL, DIMENSION( its:ite ) & :: sumtke, sumtkez INTEGER & :: i, j, k, ktf, i_start, i_end, j_start, j_end REAL & :: disp_len, deltas, coefc, tmpdz, len_s, thetasfc, & thetatop, len_0, tketmp, tmp, ce1, ce2, c_k ! End declarations. !----------------------------------------------------------------------- c_k = config_flags%c_k ce1 = ( c_k / 0.10 ) * 0.19 ce2 = max( 0.0 , 0.93 - ce1 ) ktf = MIN( kte, kde-1 ) i_start = its i_end = MIN(ite,ide-1) j_start = jts j_end = MIN(jte,jde-1) IF ( config_flags%open_xs .OR. config_flags%specified .OR. & config_flags%nested) i_start = MAX( ids+1, its ) IF ( config_flags%open_xe .OR. config_flags%specified .OR. & config_flags%nested) i_end = MIN( ide-2, ite ) IF ( config_flags%open_ys .OR. config_flags%specified .OR. & config_flags%nested) j_start = MAX( jds+1, jts ) IF ( config_flags%open_ye .OR. config_flags%specified .OR. & config_flags%nested) j_end = MIN( jde-2, jte ) IF ( config_flags%periodic_x ) i_start = its IF ( config_flags%periodic_x ) i_end = MIN( ite, ide-1 ) CALL calc_l_scale( config_flags, tke, BN2, l_scale, & i_start, i_end, ktf, j_start, j_end, & dx, dy, rdzw, msftx, msfty, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) DO j = j_start, j_end DO k = kts, ktf DO i = i_start, i_end deltas = ( dx/msftx(i,j) * dy/msfty(i,j) / rdzw(i,k,j) )**0.33333333 tketmp = MAX( tke(i,k,j), 1.0e-6 ) ! Apply Deardorff's (1980) "wall effect" at the bottom of the domain. ! For LES with fine grid, no need for this wall effect! IF ( k .eq. kts .or. k .eq. ktf ) then coefc = 3.9 ELSE coefc = ce1 + ce2 * l_scale(i,k,j) / deltas END IF tendency(i,k,j) = tendency(i,k,j) - & mu(i,j) * coefc * tketmp**1.5 / l_scale(i,k,j) END DO END DO END DO END SUBROUTINE tke_dissip !======================================================================= !======================================================================= SUBROUTINE tke_shear( tendency, config_flags, & defor11, defor22, defor33, & defor12, defor13, defor23, & u, v, w, tke, ust, mu, fnm, fnp, & cf1, cf2, cf3, msftx, msfty, & xkmh, xkmv, & rdx, rdy, zx, zy, rdz, rdzw, dn, dnw, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) ! History: Sep 2003 Rewritten by George Bryan and Jason Knievel, ! NCAR ! Oct 2001 Converted to mass core by Bill Skamarock, NCAR ! Aug 2000 Original code by Shu-Hua Chen, UC-Davis ! Purpose: This routine calculates the production of turbulent ! kinetic energy by stresses due to sheared wind. ! References: Klemp and Wilhelmson (JAS 1978) ! Deardorff (B-L Meteor 1980) ! Chen and Dudhia (NCAR WRF physics report 2000) ! Key: ! avg temporary working array ! cf1 ! cf2 ! cf3 ! defor11 deformation term ( du/dx + du/dx ) ! defor12 deformation term ( dv/dx + du/dy ); same as defor21 ! defor13 deformation term ( dw/dx + du/dz ); same as defor31 ! defor22 deformation term ( dv/dy + dv/dy ) ! defor23 deformation term ( dw/dy + dv/dz ); same as defor32 ! defor33 deformation term ( dw/dz + dw/dz ) ! div 3-d divergence ! dn ! dnw ! fnm ! fnp ! msftx ! msfty ! rdx ! rdy ! tendency ! titau tau (stress tensor) with a tilde, indicating division by ! a map-scale factor and the fraction of the total modeled ! atmosphere beneath a given altitude (titau = tau/m/zeta) ! tke turbulent kinetic energy !----------------------------------------------------------------------- ! Begin declarations. IMPLICIT NONE 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, INTENT( IN ) & :: cf1, cf2, cf3, rdx, rdy REAL, DIMENSION( kms:kme ), INTENT( IN ) & :: fnm, fnp, dn, dnw REAL, DIMENSION( ims:ime, jms:jme ), INTENT( IN ) & :: msftx, msfty REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( INOUT ) & :: tendency REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( IN ) & :: defor11, defor22, defor33, defor12, defor13, defor23, & tke, xkmh, xkmv, zx, zy, u, v, w, rdz, rdzw REAL, DIMENSION( ims:ime, jms:jme ), INTENT( IN ) & :: mu REAL, DIMENSION( ims:ime, jms:jme ), INTENT( IN ) & :: ust ! Local variables. INTEGER & :: i, j, k, ktf, ktes1, ktes2, & i_start, i_end, j_start, j_end, & is_ext, ie_ext, js_ext, je_ext REAL & :: mtau REAL, DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1 ) & :: avg, titau, tmp2 REAL, DIMENSION( its:ite, kts:kte, jts:jte ) & :: titau12, tmp1, zxavg, zyavg REAL :: absU, cd0, Cd ! End declarations. !----------------------------------------------------------------------- ktf = MIN( kte, kde-1 ) ktes1 = kte-1 ktes2 = kte-2 i_start = its i_end = MIN( ite, ide-1 ) j_start = jts j_end = MIN( jte, jde-1 ) IF ( config_flags%open_xs .OR. config_flags%specified .OR. & config_flags%nested ) i_start = MAX( ids+1, its ) IF ( config_flags%open_xe .OR. config_flags%specified .OR. & config_flags%nested ) i_end = MIN( ide-2, ite ) IF ( config_flags%open_ys .OR. config_flags%specified .OR. & config_flags%nested ) j_start = MAX( jds+1, jts ) IF ( config_flags%open_ye .OR. config_flags%specified .OR. & config_flags%nested ) j_end = MIN( jde-2, jte ) IF ( config_flags%periodic_x ) i_start = its IF ( config_flags%periodic_x ) i_end = MIN( ite, ide-1 ) DO j = j_start, j_end DO k = kts, ktf DO i = i_start, i_end zxavg(i,k,j) = 0.25 * ( zx(i,k ,j) + zx(i+1,k ,j) + & zx(i,k+1,j) + zx(i+1,k+1,j) ) zyavg(i,k,j) = 0.25 * ( zy(i,k ,j) + zy(i,k ,j+1) + & zy(i,k+1,j) + zy(i,k+1,j+1) ) END DO END DO END DO ! Begin calculating production of turbulence due to shear. The approach ! is to add together contributions from six terms, each of which is the ! square of a deformation that is then multiplied by an exchange ! coefficiant. The same exchange coefficient is assumed for horizontal ! and vertical coefficients for some of the terms (the vertical value is ! the one used). ! For defor11. DO j = j_start, j_end DO k = kts, ktf DO i = i_start, i_end tendency(i,k,j) = tendency(i,k,j) + 0.5 * & mu(i,j) * xkmh(i,k,j) * ( ( defor11(i,k,j) )**2 ) END DO END DO END DO ! For defor22. DO j = j_start, j_end DO k = kts, ktf DO i = i_start, i_end tendency(i,k,j) = tendency(i,k,j) + 0.5 * & mu(i,j) * xkmh(i,k,j) * ( ( defor22(i,k,j) )**2 ) END DO END DO END DO ! For defor33. DO j = j_start, j_end DO k = kts, ktf DO i = i_start, i_end tendency(i,k,j) = tendency(i,k,j) + 0.5 * & mu(i,j) * xkmv(i,k,j) * ( ( defor33(i,k,j) )**2 ) END DO END DO END DO ! For defor12. DO j = j_start, j_end DO k = kts, ktf DO i = i_start, i_end avg(i,k,j) = 0.25 * & ( ( defor12(i ,k,j)**2 ) + ( defor12(i ,k,j+1)**2 ) + & ( defor12(i+1,k,j)**2 ) + ( defor12(i+1,k,j+1)**2 ) ) END DO END DO END DO DO j = j_start, j_end DO k = kts, ktf DO i = i_start, i_end tendency(i,k,j) = tendency(i,k,j) + mu(i,j) * xkmh(i,k,j) * avg(i,k,j) END DO END DO END DO ! For defor13. DO j = j_start, j_end DO k = kts+1, ktf DO i = i_start, i_end+1 tmp2(i,k,j) = defor13(i,k,j) END DO END DO END DO DO j = j_start, j_end DO i = i_start, i_end+1 tmp2(i,kts ,j) = 0.0 tmp2(i,ktf+1,j) = 0.0 END DO END DO DO j = j_start, j_end DO k = kts, ktf DO i = i_start, i_end avg(i,k,j) = 0.25 * & ( ( tmp2(i ,k+1,j)**2 ) + ( tmp2(i ,k,j)**2 ) + & ( tmp2(i+1,k+1,j)**2 ) + ( tmp2(i+1,k,j)**2 ) ) END DO END DO END DO DO j = j_start, j_end DO k = kts, ktf DO i = i_start, i_end tendency(i,k,j) = tendency(i,k,j) + mu(i,j) * xkmv(i,k,j) * avg(i,k,j) END DO END DO END DO !MARTA: add the drag at the surface; WCS 040331 K=KTS uflux: SELECT CASE( config_flags%isfflx ) CASE (0) ! Assume cd a constant, specified in namelist cd0 = config_flags%tke_drag_coefficient ! drag coefficient set ! in namelist.input DO j = j_start, j_end DO i = i_start, i_end absU=0.5*sqrt((u(i,k,j)+u(i+1,k,j))**2+(v(i,k,j)+v(i,k,j+1))**2) Cd = cd0 tendency(i,k,j) = tendency(i,k,j) + & mu(i,j)*( (u(i,k,j)+u(i+1,k,j))*0.5* & Cd*absU*(defor13(i,kts+1,j)+defor13(i+1,kts+1,j))*0.5 ) END DO END DO CASE (1,2) ! ustar computed from surface routine PRINT *, 'USE u star (surf drag) from PHYSICS', ust(i_start,j_start) DO j = j_start, j_end DO i = i_start, i_end absU=0.5*sqrt((u(i,k,j)+u(i+1,k,j))**2+(v(i,k,j)+v(i,k,j+1))**2)+epsilon Cd = (ust(i,j)**2)/(absU**2) tendency(i,k,j) = tendency(i,k,j) + & mu(i,j)*( (u(i,k,j)+u(i+1,k,j))*0.5* & Cd*absU*(defor13(i,kts+1,j)+defor13(i+1,kts+1,j))*0.5 ) END DO END DO CASE DEFAULT CALL wrf_error_fatal( 'isfflx value invalid for diff_opt=2' ) END SELECT uflux ! end of MARTA/WCS change ! For defor23. DO j = j_start, j_end+1 DO k = kts+1, ktf DO i = i_start, i_end tmp2(i,k,j) = defor23(i,k,j) END DO END DO END DO DO j = j_start, j_end+1 DO i = i_start, i_end tmp2(i,kts, j) = 0.0 tmp2(i,ktf+1,j) = 0.0 END DO END DO DO j = j_start, j_end DO k = kts, ktf DO i = i_start, i_end avg(i,k,j) = 0.25 * & ( ( tmp2(i,k+1,j )**2 ) + ( tmp2(i,k,j )**2) + & ( tmp2(i,k+1,j+1)**2 ) + ( tmp2(i,k,j+1)**2) ) END DO END DO END DO DO j = j_start, j_end DO k = kts, ktf DO i = i_start, i_end tendency(i,k,j) = tendency(i,k,j) + mu(i,j) * xkmv(i,k,j) * avg(i,k,j) END DO END DO END DO !MARTA: add the drag at the surface; WCS 040331 K=KTS vflux: SELECT CASE( config_flags%isfflx ) CASE (0) ! Assume cd a constant, specified in namelist cd0 = config_flags%tke_drag_coefficient ! constant drag coefficient ! set in namelist.input DO j = j_start, j_end DO i = i_start, i_end absU=0.5*sqrt((u(i,k,j)+u(i+1,k,j))**2+(v(i,k,j)+v(i,k,j+1))**2) Cd = cd0 tendency(i,k,j) = tendency(i,k,j) + & mu(i,j)*( (v(i,k,j)+v(i,k,j+1))*0.5* & Cd*absU*(defor23(i,kts+1,j)+defor23(i,kts+1,j+1))*0.5 ) END DO END DO CASE (1,2) ! ustar computed from surface routine DO j = j_start, j_end DO i = i_start, i_end absU=0.5*sqrt((u(i,k,j)+u(i+1,k,j))**2+(v(i,k,j)+v(i,k,j+1))**2)+epsilon Cd = (ust(i,j)**2)/(absU**2) tendency(i,k,j) = tendency(i,k,j) + & mu(i,j)*( (v(i,k,j)+v(i,k,j+1))*0.5* & Cd*absU*(defor23(i,kts+1,j)+defor23(i,kts+1,j+1))*0.5 ) END DO END DO CASE DEFAULT CALL wrf_error_fatal( 'isfflx value invalid for diff_opt=2' ) END SELECT vflux ! end of MARTA/WCS change END SUBROUTINE tke_shear !======================================================================= !======================================================================= SUBROUTINE compute_diff_metrics( config_flags, ph, phb, z, rdz, rdzw, & zx, zy, rdx, rdy, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) !----------------------------------------------------------------------- ! Begin declarations. IMPLICIT NONE 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 ) & :: ph, phb REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( OUT ) & :: rdz, rdzw, zx, zy, z REAL, INTENT( IN ) & :: rdx, rdy ! Local variables. INTEGER & :: i, j, k, i_start, i_end, j_start, j_end, ktf ! End declarations. !----------------------------------------------------------------------- ktf = MIN( kte, kde-1 ) ! Bug fix, WCS, 22 april 2002. ! We need rdzw in halo for average to u and v points. j_start = jts-1 j_end = jte ! Begin with dz computations. DO j = j_start, j_end IF ( ( j_start >= jts ) .AND. ( j_end <= MIN( jte, jde-1 ) ) ) THEN i_start = its-1 i_end = ite ELSE i_start = its i_end = MIN( ite, ide-1 ) END IF ! Compute z at w points for rdz and rdzw computations. We'll switch z ! to z at p points before returning DO k = 1, kte ! Bug fix, WCS, 22 april 2002 DO i = i_start, i_end z(i,k,j) = ( ph(i,k,j) + phb(i,k,j) ) / g END DO END DO DO k = 1, ktf DO i = i_start, i_end rdzw(i,k,j) = 1.0 / ( z(i,k+1,j) - z(i,k,j) ) END DO END DO DO k = 2, ktf DO i = i_start, i_end rdz(i,k,j) = 2.0 / ( z(i,k+1,j) - z(i,k-1,j) ) END DO END DO ! Bug fix, WCS, 22 april 2002; added the following code DO i = i_start, i_end rdz(i,1,j) = 2./(z(i,2,j)-z(i,1,j)) END DO END DO ! End bug fix. ! Now compute zx and zy; we'll assume that the halo for ph and phb is ! properly filled. 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 = 1, kte DO i = MAX( ids+1, its ), i_end zx(i,k,j) = rdx * ( phb(i,k,j) - phb(i-1,k,j) ) / g END DO END DO END DO DO j = j_start, j_end DO k = 1, kte DO i = MAX( ids+1, its ), i_end zx(i,k,j) = zx(i,k,j) + rdx * ( ph(i,k,j) - ph(i-1,k,j) ) / g END DO END DO END DO DO j = MAX( jds+1, jts ), j_end DO k = 1, kte DO i = i_start, i_end zy(i,k,j) = rdy * ( phb(i,k,j) - phb(i,k,j-1) ) / g END DO END DO END DO DO j = MAX( jds+1, jts ), j_end DO k = 1, kte DO i = i_start, i_end zy(i,k,j) = zy(i,k,j) + rdy * ( ph(i,k,j) - ph(i,k,j-1) ) / g END DO END DO END DO ! Some b.c. on zx and zy. IF ( .NOT. config_flags%periodic_x ) THEN IF ( ite == ide ) THEN DO j = j_start, j_end DO k = 1, ktf zx(ide,k,j) = 0.0 END DO END DO END IF IF ( its == ids ) THEN DO j = j_start, j_end DO k = 1, ktf zx(ids,k,j) = 0.0 END DO END DO END IF ELSE IF ( ite == ide ) THEN DO j=j_start,j_end DO k=1,ktf zx(ide,k,j) = rdx * ( phb(ide,k,j) - phb(ide-1,k,j) ) / g END DO END DO DO j = j_start, j_end DO k = 1, ktf zx(ide,k,j) = zx(ide,k,j) + rdx * ( ph(ide,k,j) - ph(ide-1,k,j) ) / g END DO END DO END IF IF ( its == ids ) THEN DO j = j_start, j_end DO k = 1, ktf zx(ids,k,j) = rdx * ( phb(ids,k,j) - phb(ids-1,k,j) ) / g END DO END DO DO j =j_start,j_end DO k =1,ktf zx(ids,k,j) = zx(ids,k,j) + rdx * ( ph(ids,k,j) - ph(ids-1,k,j) ) / g END DO END DO END IF END IF IF ( .NOT. config_flags%periodic_y ) THEN IF ( jte == jde ) THEN DO k =1, ktf DO i =i_start, i_end zy(i,k,jde) = 0.0 END DO END DO END IF IF ( jts == jds ) THEN DO k =1, ktf DO i =i_start, i_end zy(i,k,jds) = 0.0 END DO END DO END IF ELSE IF ( jte == jde ) THEN DO j=j_start, j_end DO k=1, ktf zy(i,k,jde) = rdy * ( phb(i,k,jde) - phb(i,k,jde-1) ) / g END DO END DO DO j = j_start, j_end DO k = 1, ktf zy(i,k,jde) = zy(i,k,jde) + rdy * ( ph(i,k,jde) - ph(i,k,jde-1) ) / g END DO END DO END IF IF ( jts == jds ) THEN DO j = j_start, j_end DO k = 1, ktf zy(i,k,jds) = rdy * ( phb(i,k,jds) - phb(i,k,jds-1) ) / g END DO END DO DO j = j_start, j_end DO k = 1, ktf zy(i,k,jds) = zy(i,k,jds) + rdy * ( ph(i,k,jds) - ph(i,k,jds-1) ) / g END DO END DO END IF END IF ! Calculate z at p points. DO j = j_start, j_end DO k = 1, ktf DO i = i_start, i_end z(i,k,j) = 0.5 * & ( ph(i,k,j) + phb(i,k,j) + ph(i,k+1,j) + phb(i,k+1,j) ) / g END DO END DO END DO END SUBROUTINE compute_diff_metrics !======================================================================= !======================================================================= END MODULE module_diffusion_em !======================================================================= !=======================================================================