/wrfv2_fire/dyn_em/module_diffusion_em.F
FORTRAN Legacy | 6015 lines | 3731 code | 1233 blank | 1051 comment | 10 complexity | 4a8d72e31f9f2c2d1120d5720275b558 MD5 | raw file
Possible License(s): AGPL-1.0
- ! WRF:MODEL_LAYER:PHYSICS
-
- MODULE module_diffusion_em
- USE module_bc, only: set_physical_bc3d
- USE module_state_description, only: p_m23, p_m13, p_m22, p_m33, p_r23, p_r13, p_r12, p_m12, p_m11
- USE module_big_step_utilities_em, only: grid_config_rec_type, param_first_scalar, p_qv, p_qi, p_qc
- USE module_model_constants
- CONTAINS
- !=======================================================================
- !=======================================================================
- SUBROUTINE cal_deform_and_div( config_flags, u, v, w, div, &
- defor11, defor22, defor33, &
- defor12, defor13, defor23, &
- nba_rij, n_nba_rij, & !JDM
- 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
- INTEGER, INTENT( IN ) :: n_nba_rij !JDM
- REAL , DIMENSION(ims:ime, kms:kme, jms:jme, n_nba_rij), INTENT(INOUT) & !JDM
- :: nba_rij
- ! 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)
- IF ( config_flags%sfs_opt .GT. 0 ) THEN ! NBA--
- !JDM____________________________________________________________________
- !
- ! s12 = du/dy + dv/dx
- ! = (du/dy - dz/dy*du/dz) + (dv/dx - dz/dx*dv/dz)
- ! ______defor12______ ___tmp1___
- !
- ! r12 = du/dy - dv/dx
- ! = (du/dy - dz/dy*du/dz) - (dv/dx - dz/dx*dv/dz)
- ! ______defor12______ ___tmp1___
- !_______________________________________________________________________
- DO j = j_start, j_end
- DO k = kts, ktf
- DO i = i_start, i_end
- nba_rij(i,k,j,P_r12) = defor12(i,k,j) - &
- mm(i,j) * ( &
- rdx * ( hat(i,k,j) - hat(i-1,k,j) ) - tmp1(i,k,j) )
- 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)
- nba_rij(ids,k,j,P_r12) = nba_rij(ids+1,k,j,P_r12)
- 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)
- nba_rij(i,k,jds,P_r12) = nba_rij(i,k,jds+1,P_r12)
- 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)
- nba_rij(ide,k,j,P_r12) = nba_rij(ide-1,k,j,P_r12)
- 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)
- nba_rij(i,k,jde,P_r12) = nba_rij(i,k,jde-1,P_r12)
- END DO
- END DO
- END IF
- ELSE ! NOT NBA--------------------------------------------------------
- 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
- ENDIF ! NBA-----------------------------------------------------------
- ! 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.
- IF ( config_flags%sfs_opt .GT. 0 ) THEN ! NBA--
- !JDM____________________________________________________________________
- !
- ! s13 = du/dz + dw/dx
- ! = du/dz + (dw/dx - dz/dx*dw/dz)
- ! = tmp1 + ______defor13______
- !
- ! r13 = du/dz - dw/dx
- ! = du/dz - (dw/dx - dz/dx*dw/dz)
- ! = tmp1 - ______defor13______
- !_______________________________________________________________________
- DO j = j_start, j_end
- DO k = kts+1, ktf
- DO i = i_start, i_end
- nba_rij(i,k,j,P_r13) = tmp1(i,k,j) - defor13(i,k,j)
- defor13(i,k,j) = defor13(i,k,j) + tmp1(i,k,j)
- END DO
- END DO
- END DO
- DO j = j_start, j_end !change for different surface B. C.
- DO i = i_start, i_end
- nba_rij(i,kts ,j,P_r13) = 0.0
- nba_rij(i,ktf+1,j,P_r13) = 0.0
- END DO
- END DO
- ELSE ! NOT NBA--------------------------------------------------------
- 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
- ENDIF ! NBA-----------------------------------------------------------
- ! 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.
- IF ( config_flags%sfs_opt .GT. 0 ) THEN ! NBA--
- !JDM___________________________________________________________________
- !
- ! s23 = dv/dz + dw/dy
- ! = dv/dz + (dw/dy - dz/dy*dw/dz)
- ! tmp1 + ______defor23______
- !
- ! r23 = dv/dz - dw/dy
- ! = dv/dz - (dw/dy - dz/dy*dw/dz)
- ! = tmp1 - ______defor23______
- ! Add tmp1 to defor23.
- DO j = j_start, j_end
- DO k = kts+1, ktf
- DO i = i_start, i_end
- nba_rij(i,k,j,P_r23) = tmp1(i,k,j) - defor23(i,k,j)
- defor23(i,k,j) = defor23(i,k,j) + tmp1(i,k,j)
- END DO
- END DO
- END DO
- DO j = j_start, j_end
- DO i = i_start, i_end
- nba_rij(i,kts ,j,P_r23) = 0.0
- nba_rij(i,ktf+1,j,P_r23) = 0.0
- 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)
- nba_rij(ids,k,j,P_r13) = nba_rij(ids+1,k,j,P_r13)
- nba_rij(ids,k,j,P_r23) = nba_rij(ids+1,k,j,P_r23)
- 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)
- nba_rij(i,k,jds,P_r13) = nba_rij(i,k,jds+1,P_r13)
- nba_rij(i,k,jds,P_r23) = nba_rij(i,k,jds+1,P_r23)
- 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)
- nba_rij(ide,k,j,P_r13) = nba_rij(ide-1,k,j,P_r13)
- nba_rij(ide,k,j,P_r23) = nba_rij(ide-1,k,j,P_r23)
- 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)
- nba_rij(i,k,jde,P_r13) = nba_rij(i,k,jde-1,P_r13)
- nba_rij(i,k,jde,P_r23) = nba_rij(i,k,jde-1,P_r23)
- END DO
- END DO
- END IF
- ELSE ! NOT NBA--------------------------------------------------------
- ! 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
- ENDIF ! NBA-----------------------------------------------------------
- ! 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, &
- tracer_tendf, n_tracer, &
- thp, theta, mu, tke, config_flags, &
- defor11, defor22, defor12, &
- defor13, defor23, &
- nba_mij, n_nba_mij, & !JDM
- div, &
- moist, chem, scalar,tracer, &
- 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,n_tracer,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_tracer), &
- INTENT(INOUT) :: tracer_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, n_tracer) , &
- INTENT(IN ) :: tracer
- 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
- INTEGER, INTENT( IN ) :: n_nba_mij !JDM
- REAL , DIMENSION(ims:ime, kms:kme, jms:jme, n_nba_mij), INTENT(INOUT) & !JDM
- :: nba_mij
- ! 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, &
- nba_mij, n_nba_mij, & !JDM
- 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, &
- nba_mij, n_nba_mij, & !JDM
- 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, &
- nba_mij, n_nba_mij, & !JDM
- 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_tracer .ge. PARAM_FIRST_SCALAR) THEN
- tracer_loop: do ic = PARAM_FIRST_SCALAR, n_tracer
- CALL horizontal_diffusion_s( tracer_tendf(ims,kms,jms,ic), &
- mu, config_flags, &
- tracer(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 tracer_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, &
- nba_mij, n_nba_mij, & !JDM
- 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
- INTEGER, INTENT( IN ) :: n_nba_mij !JDM
- REAL , DIMENSION(ims:ime, kms:kme, jms:jme, n_nba_mij), INTENT(INOUT) & !JDM
- :: nba_mij
- 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, &
- nba_mij(ims,kms,jms,P_m11), & !JDM
- 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, &
- nba_mij(ims,kms,jms,P_m12), & !JDM
- 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, &
- nba_mij, n_nba_mij, & !JDM
- 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
- INTEGER, INTENT( IN ) :: n_nba_mij !JDM
- REAL , DIMENSION(ims:ime, kms:kme, jms:jme, n_nba_mij), INTENT(INOUT) & !JDM
- :: nba_mij
- 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, &
- nba_mij(ims,kms,jms,P_m12), & !JDM
- 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, &
- nba_mij(ims,kms,jms,P_m22), & !JDM
- 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, &
- nba_mij, n_nba_mij, & !JDM
- 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
- INTEGER, INTENT( IN ) :: n_nba_mij !JDM
- REAL , DIMENSION(ims:ime, kms:kme, jms:jme, n_nba_mij), INTENT(INOUT) & !JDM
- :: nba_mij
- 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, &
- nba_mij(ims,kms,jms,P_m13), & !JDM
- 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, &
- nba_mij(ims,kms,jms,P_m23), & !JDM
- 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, &
- tracer_tendf, n_tracer, &
- u_2, v_2, &
- thp,u_base,v_base,t_base,qv_base,mu,tke, &
- config_flags,defor13,defor23,defor33, &
- nba_mij, n_nba_mij, & !JDM
- div, &
- moist,chem,scalar,tracer,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,n_tracer,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_tracer) , &
- INTENT(INOUT) :: tracer_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, n_tracer) , &
- INTENT(IN ) :: tracer
- REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(IN ) ::defor13, &
- defor23, &
- defor33, &
- div, &
- xkmv, &
- xkhv, &
- tke, &
- rdz, &
- u_2, &
- v_2, &
- rdzw
- INTEGER, INTENT( IN ) :: n_nba_mij !JDM
- REAL , DIMENSION(ims:ime, kms:kme, jms:jme, n_nba_mij), INTENT(INOUT) & !JDM
- :: nba_mij
- REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(IN ) :: rho
- REAL , DIMENSION( ims:ime, jms:jme), INTENT(INOUT) :: 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, &
- nba_mij, n_nba_mij, & !JDM
- 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, &
- nba_mij, n_nba_mij, & !JDM
- 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), &
- nba_mij, n_nba_mij, & !JDM
- 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
- cpm = cp * (1. + 0.8 * moist(i,kts,j,P_QV))
- hfx(i,j)=heat_flux*cp*rho(i,1,j) ! provided for output only
- 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)/(1.+moist(i,kts,j,P_QV))
- 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_tracer .ge. PARAM_FIRST_SCALAR) THEN
- tracer_loop: do im = PARAM_FIRST_SCALAR, n_tracer
- CALL vertical_diffusion_s( tracer_tendf(ims,kms,jms,im), &
- config_flags, tracer(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 tracer_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, &
- nba_mij, n_nba_mij, & !JDM
- 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
- INTEGER, INTENT( IN ) :: n_nba_mij !JDM
- REAL , DIMENSION(ims:ime, kms:kme, jms:jme, n_nba_mij), INTENT(INOUT) & !JDM
- :: nba_mij
- 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, &
- nba_mij(ims,kms,jms,P_m13), & !JDM
- 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, &
- nba_mij, n_nba_mij, & !JDM
- 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
- INTEGER, INTENT( IN ) :: n_nba_mij !JDM
- REAL , DIMENSION(ims:ime, kms:kme, jms:jme, n_nba_mij), INTENT(INOUT) & !JDM
- :: nba_mij
- 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, &
- nba_mij(ims,kms,jms,P_m23), & !JDM
- 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, &
- nba_mij, n_nba_mij, & !JDM
- 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
- INTEGER, INTENT( IN ) :: n_nba_mij !JDM
- REAL , DIMENSION(ims:ime, kms:kme, jms:jme, n_nba_mij), INTENT(INOUT) & !JDM
- :: nba_mij
- 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, &
- nba_mij(ims,kms,jms,P_m33), & !JDM
- 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, &
- mtau, & !JDM
- 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, kms:kme, jms:jme ), INTENT( INOUT ) & !JDM
- :: mtau
- 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%sfs_opt .GT. 0 ) THEN ! USE NBA MODEL SFS STRESSES
- DO j = j_start, j_end
- DO k = kts, ktf
- DO i = i_start, i_end
- titau(i,k,j) = mu(i,j) * mtau(i,k,j)
- END DO
- END DO
- END DO
- ELSE !NOT NBA
- IF ( config_flags%m_opt .EQ. 1 ) THEN ! ASSIGN STRESS TO MTAU FOR OUTPUT
- 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)
- mtau(i,k,j) = - xkx(i,k,j) * defor(i,k,j)
-
- END DO
- END DO
- END DO
- ELSE !NO STRESS OUTPUT
- 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
- ENDIF
- ENDIF
- END SUBROUTINE cal_titau_11_22_33
- !=======================================================================
- !=======================================================================
- SUBROUTINE cal_titau_12_21( config_flags, titau, &
- mu, xkx, defor, &
- mtau, & !JDM
- 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, kms:kme, jms:jme ), INTENT( INOUT ) & !JDM
- :: mtau
- 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
- IF ( config_flags%sfs_opt .GT. 0 ) THEN ! USE NBA MODEL SFS STRESSES
-
- DO j = j_start, j_end
- DO k = kts, ktf
- DO i = i_start, i_end
- titau(i,k,j) = muavg(i,j) * mtau(i,k,j)
- END DO
- END DO
- END DO
- ELSE ! NOT NBA
-
- IF ( config_flags%m_opt .EQ. 1 ) THEN ! ASSIGN STRESS TO MTAU FOR OUTPUT
- 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)
- mtau(i,k,j) = - xkxavg(i,k,j) * defor(i,k,j)
- END DO
- END DO
- END DO
- ELSE ! NO STRESS OUTPUT
- 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
- ENDIF
- ENDIF
- END SUBROUTINE cal_titau_12_21
- !=======================================================================
- SUBROUTINE cal_titau_13_31( config_flags, titau, &
- defor, &
- mtau, & !JDM
- 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, kms:kme, jms:jme ), INTENT( INOUT ) & !JDM
- :: mtau
- 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
- IF ( config_flags%sfs_opt .GT. 0 ) THEN ! USE NBA MODEL SFS STRESSES
-
- DO j = j_start, j_end
- DO k = kts+1, ktf
- DO i = i_start, i_end
- titau(i,k,j) = muavg(i,j) * mtau(i,k,j)
- ENDDO
- ENDDO
- ENDDO
- ELSE ! NOT NBA
-
- IF ( config_flags%m_opt .EQ. 1 ) THEN ! ASSIGN STRESS TO MTAU FOR OUTPUT
- 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)
- mtau(i,k,j) = - xkxavg(i,k,j) * defor(i,k,j)
-
- ENDDO
- ENDDO
- ENDDO
- ELSE ! NO STRESS OUTPUT
- 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
- ENDIF
- ENDIF
- 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, &
- mtau, & !JDM
- 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, kms:kme, jms:jme ), INTENT( INOUT ) & !JDM
- :: mtau
- 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
-
- IF ( config_flags%sfs_opt .EQ. 1 ) THEN ! USE NBA MODEL SFS STRESSES
- DO j = j_start, j_end
- DO k = kts+1, ktf
- DO i = i_start, i_end
- titau(i,k,j) = muavg(i,j) * mtau(i,k,j)
- END DO
- END DO
- END DO
- ELSE ! NOT NBA
- IF ( config_flags%m_opt .EQ. 1 ) THEN ! ASSIGN STRESS TO MTAU FOR OUTPUT
- 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)
- mtau(i,k,j) = - xkxavg(i,k,j) * defor(i,k,j)
- END DO
- END DO
- END DO
- ELSE ! NO STRESS OUTPUT
- 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
- ENDIF
- ENDIF
- 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, &
- RUCUTEN, RVCUTEN, &
- RUSHTEN, RVSHTEN, &
- 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, &
- RUCUTEN, &
- RVCUTEN, &
- RUSHTEN, &
- RVSHTEN, &
- defor11, &
- defor22, &
- defor33, &
- defor12, &
- defor13, &
- defor23, &
- xkmh, &
- xkmv, &
- xkhh, &
- xkhv, &
- tke, &
- div
- ! End declarations.
- !-----------------------------------------------------------------------
- IF(config_flags%bl_pbl_physics .GT. 0) THEN
- 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
- IF(config_flags%cu_physics .GT. 0) THEN
- CALL set_physical_bc3d( RUCUTEN , '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( RVCUTEN , '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
- IF(config_flags%shcu_physics .GT. 0) THEN
- CALL set_physical_bc3d( RUSHTEN , '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( RVSHTEN , '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( 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
- 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
- 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
- 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
- !=======================================================================
- !=======================================================================