pkg/atm_phys/vert_diff_mod.F90

b2ea1d2979 Jean*0001 module vert_diff_mod
                0002
                0003 !=======================================================================
                0004 !
                0005 !                         VERTICAL DIFFUSION MODULE
                0006 !
                0007 !      Routines for computing the tendencies due to vertical diffusion
                0008 !
                0009 !=======================================================================
                0010
                0011 use   constants_mod, only:  GRAV, RDGAS, RVGAS, CP_air
                0012
                0013 !use         fms_mod, only:  error_mesg, FATAL, uppercase, &
                0014 !                            write_version_number, stdlog, &
                0015 !                            mpp_pe, mpp_root_pe
                0016
                0017 !use   field_manager_mod, only: MODEL_ATMOS
                0018 !use  tracer_manager_mod, only: query_method, get_tracer_index
                0019
                0020 implicit none
                0021 private
                0022
                0023
                0024 ! public interfaces
                0025 !=======================================================================
                0026 public :: gcm_vert_diff_init,          &
                0027           gcm_vert_diff_end,           &
                0028           gcm_vert_diff,               &
                0029           gcm_vert_diff_down,          &
                0030           gcm_vert_diff_up,            &
                0031           vert_diff,                   &
                0032           surf_diff_type
                0033
                0034 !=======================================================================
                0035
                0036 ! form of interfaces
                0037 !=======================================================================
                0038
                0039 type surf_diff_type
                0040
                0041   real, pointer, dimension(:,:) :: dtmass,    &
                0042                                    dflux_t,   &
                0043                                    dflux_q,   &
                0044                                    delta_t,   &
                0045                                    delta_q
                0046
                0047 end type surf_diff_type
                0048
                0049 !real,    allocatable, dimension(:,:,:) :: e_global, f_t_global, f_q_global
                0050
                0051 logical :: do_init = .true.
                0052 logical :: do_conserve_energy = .true.
                0053 logical :: use_virtual_temp_vert_diff, do_mcm_plev
                0054 integer :: sphum
                0055
                0056 !--------------------- version number ---------------------------------
                0057
                0058 character(len=128) :: version = '$Id: vert_diff_mod.F90,v 1.1 2013/05/08 22:14:15 jmc Exp $'
                0059 character(len=128) :: tag = '$Name:  $'
                0060
                0061 real, parameter :: d608 = (RVGAS-RDGAS)/RDGAS
                0062
                0063 contains
                0064
                0065 !#######################################################################
                0066
                0067 subroutine gcm_vert_diff_init (idim, jdim, kdim,              &
                0068                                do_conserve_energy_in, myThid, &
                0069                                use_virtual_temp_vert_diff_in, &
                0070                                do_mcm_plev_in )
                0071
                0072 !type(surf_diff_type), intent(inout) :: Tri_surf
                0073  integer,              intent(in)    :: idim, jdim, kdim
                0074  logical,              intent(in)    :: do_conserve_energy_in
                0075  integer, intent(in)                 :: myThid
                0076  logical, optional,    intent(in)    :: use_virtual_temp_vert_diff_in
                0077  logical, optional,    intent(in)    :: do_mcm_plev_in
                0078
                0079 !    _BARRIER
                0080 !    _BEGIN_MASTER(myThid)
                0081      CALL BARRIER(myThid)
                0082      IF ( myThid.EQ.1 ) THEN
                0083
                0084 !   call write_version_number ( version, tag )
                0085
                0086 ! get the tracer number for specific humidity
                0087 ! jmc: --- need to do something ---
                0088     sphum = 2
                0089 !   sphum = get_tracer_index( MODEL_ATMOS, 'sphum')
                0090 !   if (mpp_pe() == mpp_root_pe()) &
                0091 !   write (stdlog(),'(a,i4)') 'Tracer number for specific humidity =',sphum
                0092
                0093     if(present(use_virtual_temp_vert_diff_in)) then
                0094       use_virtual_temp_vert_diff = use_virtual_temp_vert_diff_in
                0095     else
                0096       use_virtual_temp_vert_diff = .false.
                0097     endif
                0098     if(present(do_mcm_plev_in)) then
                0099       do_mcm_plev = do_mcm_plev_in
                0100     else
                0101       do_mcm_plev = .false.
                0102     endif
                0103
                0104  if (do_init) then
                0105
                0106 !   if (allocated(  e_global ))    deallocate (  e_global )
                0107 !   if (allocated(f_t_global ))    deallocate (f_t_global )
                0108 !   if (allocated(f_q_global ))    deallocate (f_q_global )
                0109
                0110 !   allocate(  e_global (idim, jdim, kdim-1))
                0111 !   allocate(f_t_global (idim, jdim, kdim-1))
                0112 !   allocate(f_q_global (idim, jdim, kdim-1))
                0113
                0114     do_init = .false.
                0115
                0116  endif
                0117
                0118 !call alloc_surf_diff_type ( Tri_surf, idim, jdim, myThid )
                0119
                0120  do_conserve_energy = do_conserve_energy_in
                0121
                0122      ENDIF
                0123      CALL BARRIER(myThid)
                0124
                0125 end subroutine gcm_vert_diff_init
                0126
                0127 !#######################################################################
                0128
                0129 subroutine alloc_surf_diff_type ( Tri_surf, idim, jdim, myThid )
                0130
                0131 type(surf_diff_type), intent(inout) :: Tri_surf
                0132 integer,              intent(in)    :: idim, jdim
                0133 integer,              intent(in)    :: myThid
                0134
                0135     allocate( Tri_surf%dtmass    (idim, jdim) )
                0136     allocate( Tri_surf%dflux_t   (idim, jdim) )
                0137     allocate( Tri_surf%dflux_q   (idim, jdim) )
                0138     allocate( Tri_surf%delta_t   (idim, jdim) )
                0139     allocate( Tri_surf%delta_q   (idim, jdim) )
                0140
                0141 end subroutine alloc_surf_diff_type
                0142
                0143 !#######################################################################
                0144
                0145 subroutine dealloc_surf_diff_type ( Tri_surf, myThid )
                0146
                0147 type(surf_diff_type), intent(inout) :: Tri_surf
                0148 integer,intent(in)                  :: myThid
                0149
                0150       deallocate( Tri_surf%dtmass    )
                0151       deallocate( Tri_surf%dflux_t   )
                0152       deallocate( Tri_surf%dflux_q   )
                0153       deallocate( Tri_surf%delta_t   )
                0154       deallocate( Tri_surf%delta_q   )
                0155
                0156 end subroutine dealloc_surf_diff_type
                0157
                0158 !#######################################################################
                0159
                0160 subroutine gcm_vert_diff_end (myThid)
                0161
                0162 integer, intent(in)  ::myThid
                0163
                0164   if (.not.do_init) then
                0165
                0166 !   if (allocated(   e_global ))    deallocate (   e_global)
                0167 !   if (allocated( f_t_global ))    deallocate ( f_t_global)
                0168 !   if (allocated( f_q_global ))    deallocate ( f_q_global)
                0169
                0170   endif
                0171
                0172 end subroutine gcm_vert_diff_end
                0173
                0174 !#######################################################################
                0175
                0176 subroutine gcm_vert_diff_down (is, js, delt,                 &
                0177                          u, v, t, q,                         &
                0178 !                        tr,                                 &
                0179                          diff_m, diff_t, p_half, p_full,     &
                0180                          z_full, tau_u, tau_v, dtau_datmos,  &
                0181 !                        flux_tr,                            &
                0182                          dt_u, dt_v, dt_t, dt_q,             &
                0183 !                        dt_tr,                              &
                0184                          dissipative_heat, tri_surf_dtmass,  &
                0185                          tri_surf_dflux_t, tri_surf_dflux_q, &
                0186                          tri_surf_delta_t, tri_surf_delta_q, &
                0187                          e_global, f_t_global, f_q_global,   &
                0188                          myThid, kbot )
                0189
                0190 integer, intent(in)                        :: is, js
                0191 real,    intent(in)                        :: delt
                0192 real,    intent(in)   , dimension(:,:,:)   :: u, v, t, q,     &
                0193                                               diff_m, diff_t, &
                0194                                               p_half, p_full, &
                0195                                               z_full
                0196 !real,    intent(in)   , dimension(:,:,:,:) :: tr
                0197 real,    intent(in)   , dimension(:,:)     :: dtau_datmos
                0198 real,    intent(inout), dimension(:,:)     :: tau_u, tau_v
                0199 !real,    intent(inout), dimension(:,:,:)   :: flux_tr
                0200 real,    intent(inout), dimension(:,:,:)   :: dt_u, dt_v, dt_t
                0201 real,    intent(in),    dimension(:,:,:)   :: dt_q
                0202 !real,    intent(inout), dimension(:,:,:,:) :: dt_tr
                0203 real,    intent(out)  , dimension(:,:,:)   :: dissipative_heat
                0204 !type(surf_diff_type), intent(inout)        :: Tri_surf
                0205 real,    intent(inout), dimension(:,:)     :: tri_surf_dtmass
                0206 real,    intent(inout), dimension(:,:)     :: tri_surf_dflux_t, tri_surf_dflux_q
                0207 real,    intent(inout), dimension(:,:)     :: tri_surf_delta_t, tri_surf_delta_q
                0208 real,    intent(out),   dimension(:,:,:)   :: e_global, f_t_global, f_q_global
                0209 integer, intent(in)                        :: myThid
                0210
                0211 integer, intent(in)   , dimension(:,:), optional :: kbot
                0212
                0213 real, dimension(size(u,1),size(u,2),size(u,3)) :: tt, mu, nu
                0214
                0215 real, dimension(size(u,1),size(u,2)) :: f_t_delt_n1, f_q_delt_n1, &
                0216             mu_delt_n, nu_n, e_n1, delta_t_n, delta_q_n
                0217
                0218 real    :: gcp
                0219 integer :: i, j, kb, ie, je
                0220
                0221 !-----------------------------------------------------------------------
                0222
                0223   if(do_init) call PRINT_ERROR ( &
                0224        'gcm_vert_diff_down in vert_diff_mod'//  &
                0225       'the initialization routine gcm_vert_diff_init has not been called', &
                0226        myThid )
                0227 ! if(do_init) call error_mesg ('gcm_vert_diff_down in vert_diff_mod',  &
                0228 !     'the initialization routine gcm_vert_diff_init has not been called', &
                0229 !      FATAL)
                0230
                0231  ie = is + size(t,1) -1
                0232  je = js + size(t,2) -1
                0233
                0234  gcp       = GRAV/CP_air
                0235  tt  = t + z_full*gcp   ! the vertical gradient of tt determines the
                0236                         ! diffusive flux of temperature
                0237
                0238  call compute_mu (p_half, mu, myThid)
                0239  call compute_nu (diff_m, p_half, p_full, z_full, t, q, nu, myThid)
                0240
                0241 !  diffuse u-momentum and v_momentum
                0242
                0243  call uv_vert_diff (delt, mu, nu, u, v, dtau_datmos, tau_u, tau_v,  &
                0244                     dt_u, dt_v, dt_t, dissipative_heat, myThid, kbot)
                0245
                0246 !  recompute nu for a different diffusivity
                0247  call compute_nu   (diff_t, p_half, p_full, z_full, t, q, nu, myThid)
                0248
                0249 !  diffuse tracers
                0250 !call tr_vert_diff (delt, mu, nu, tr, flux_tr, dt_tr, myThid, kbot )
                0251
                0252 !  downward sweep of tridiagonal solver for temperature and specific humidity
                0253  call vert_diff_down_2                            &
                0254          (delt, mu, nu, tt, q, dt_t, dt_q,        &
                0255          e_global             (is:ie,js:je,:),    &
                0256          f_t_global           (is:ie,js:je,:),    &
                0257          f_q_global           (is:ie,js:je,:),    &
                0258          mu_delt_n, nu_n, e_n1, f_t_delt_n1, f_q_delt_n1, &
                0259          delta_t_n, delta_q_n, myThid, kbot)
                0260
                0261 ! store information needed by flux_exchange module
                0262
                0263     tri_surf_delta_t (is:ie,js:je) = delta_t_n + mu_delt_n*nu_n*f_t_delt_n1
                0264     tri_surf_delta_q (is:ie,js:je) = delta_q_n + mu_delt_n*nu_n*f_q_delt_n1
                0265     tri_surf_dflux_t (is:ie,js:je) = -nu_n*(1.0 - e_n1)
                0266     tri_surf_dflux_q (is:ie,js:je) = -nu_n*(1.0 - e_n1)
                0267     tri_surf_dtmass  (is:ie,js:je) = mu_delt_n
                0268
                0269 !-----------------------------------------------------------------------
                0270
                0271 end subroutine gcm_vert_diff_down
                0272
                0273 !#######################################################################
                0274
                0275 subroutine gcm_vert_diff_up (is, js, delt,                   &
                0276                          tri_surf_delta_t, tri_surf_delta_q, &
                0277                          e_global, f_t_global, f_q_global,   &
                0278                          dt_t, dt_q, myThid, kbot)
                0279
                0280 integer, intent(in)                      :: is, js
                0281 real,    intent(in)                      :: delt
                0282 integer, intent(in)                      :: myThid
                0283 !type(surf_diff_type), intent(in)         :: Tri_surf
                0284 real,    intent(in),    dimension(:,:)   :: tri_surf_delta_t, tri_surf_delta_q
                0285 real,    intent(in),    dimension(:,:,:) :: e_global, f_t_global, f_q_global
                0286 real,    intent(out),   dimension(:,:,:) :: dt_t, dt_q
                0287 integer, intent(in),    dimension(:,:), optional :: kbot
                0288
                0289 integer :: ie, je
                0290
                0291  ie = is + size(dt_t,1) -1
                0292  je = js + size(dt_t,2) -1
                0293
                0294
                0295  call vert_diff_up (delt ,                              &
                0296                     e_global          (is:ie,js:je,:) , &
                0297                     f_t_global        (is:ie,js:je,:) , &
                0298                     tri_surf_delta_t  (is:ie,js:je) ,   &
                0299                     dt_t, myThid, kbot )
                0300
                0301  call vert_diff_up (delt ,                              &
                0302                     e_global          (is:ie,js:je,:) , &
                0303                     f_q_global        (is:ie,js:je,:) , &
                0304                     tri_surf_delta_q  (is:ie,js:je) ,   &
                0305                     dt_q, myThid, kbot )
                0306
                0307
                0308 end subroutine gcm_vert_diff_up
                0309
                0310 !#######################################################################
                0311
                0312 subroutine gcm_vert_diff (delt, u, v, t, q, tr,                    &
                0313                           diff_m, diff_t, p_half, p_full, z_full,  &
                0314                           dtau_datmos, dsens_datmos, devap_datmos, &
                0315                           sens, evap, tau_u, tau_v, flux_tr,       &
                0316                           dt_u, dt_v, dt_t, dt_q, dt_tr,           &
                0317                           dissipative_heat, myThid, kbot      )
                0318
                0319 !  one-step diffusion call for gcm in which there is no implicit dependence of
                0320 !    surface fluxes on surface temperature
                0321
                0322 real,    intent(in)                          :: delt
                0323 real,    intent(in)   , dimension(:,:,:)     :: u, v, t, q, p_half, p_full, &
                0324                                                 z_full, diff_m, diff_t
                0325 real,    intent(in)   , dimension(:,:,:,:)   :: tr
                0326 real,    intent(in)   , dimension(:,:)       :: dtau_datmos, dsens_datmos, &
                0327                                                 devap_datmos
                0328 real,    intent(inout), dimension(:,:,:)     :: flux_tr
                0329 real,    intent(inout), dimension(:,:)       :: tau_u, tau_v, sens, evap
                0330 real,    intent(inout), dimension(:,:,:)     :: dt_u, dt_v, dt_t, dt_q
                0331 real,    intent(inout), dimension(:,:,:,:)   :: dt_tr
                0332 real,    intent(out)  , dimension(:,:,:)     :: dissipative_heat
                0333 integer, intent(in)                          :: myThid
                0334
                0335 integer, intent(in)   , dimension(:,:), optional :: kbot
                0336
                0337 real, dimension(size(u,1),size(u,2),size(u,3)) :: mu, nu
                0338
                0339
                0340 !-----------------------------------------------------------------------
                0341
                0342  call compute_mu (p_half, mu, myThid)
                0343
                0344  call compute_nu (diff_m, p_half, p_full, z_full, t, q, nu, myThid)
                0345
                0346  call uv_vert_diff (delt, mu, nu, u, v, dtau_datmos, tau_u, tau_v, &
                0347                     dt_u, dt_v, dt_t, dissipative_heat, myThid, kbot )
                0348
                0349  call compute_nu   (diff_t, p_half, p_full, z_full, t, q, nu, myThid)
                0350
                0351  call tq_vert_diff (delt, mu, nu, t, q, z_full,  &
                0352                     dsens_datmos, devap_datmos,  &
                0353                     sens, evap, dt_t, dt_q, myThid, kbot )
                0354
                0355  call tr_vert_diff (delt, mu, nu, tr, flux_tr, dt_tr, myThid, kbot )
                0356
                0357 end subroutine gcm_vert_diff
                0358
                0359 !#######################################################################
                0360
                0361 subroutine vert_diff (delt, xi, t, q, diff, p_half, p_full, z_full, &
                0362                       flux, dflux_datmos, factor, dt_xi, myThid, kbot)
                0363
                0364 ! one-step diffusion of a single field
                0365
                0366 real,    intent(in)                          :: delt
                0367 real,    intent(in)   , dimension(:,:,:)     :: xi, t, q, diff, p_half, p_full, z_full
                0368 real,    intent(inout), dimension(:,:)       :: flux
                0369 real,    intent(in)   , dimension(:,:)       :: dflux_datmos
                0370 real,    intent(in)                          :: factor
                0371 real,    intent(inout), dimension(:,:,:)     :: dt_xi
                0372 integer, intent(in)                          :: myThid
                0373
                0374 integer, intent(in)   , dimension(:,:), optional :: kbot
                0375
                0376 real, dimension(size(xi,1),size(xi,2),size(xi,3)  ) :: mu, nu
                0377 real, dimension(size(xi,1),size(xi,2),size(xi,3)-1) :: e, f
                0378
                0379 real, dimension(size(xi,1),size(xi,2))  :: mu_delt_n, nu_n, e_n1,  &
                0380                                            f_delt_n1, delta_xi_n
                0381
                0382 !-----------------------------------------------------------------------
                0383
                0384  call compute_mu    (p_half, mu, myThid)
                0385
                0386  call compute_nu    (diff, p_half, p_full, z_full, t, q, nu, myThid)
                0387
                0388  call vert_diff_down &
                0389      (delt, mu, nu, xi, dt_xi, e, f, mu_delt_n, nu_n, e_n1,  &
                0390       f_delt_n1, delta_xi_n, myThid, kbot)
                0391
                0392  call diff_surface (mu_delt_n, nu_n, e_n1, f_delt_n1,     &
                0393                     dflux_datmos, flux, factor, delta_xi_n, myThid)
                0394
                0395  call vert_diff_up (delt, e, f, delta_xi_n, dt_xi, myThid, kbot)
                0396
                0397 end subroutine vert_diff
                0398
                0399
                0400 !#######################################################################
                0401
                0402 subroutine uv_vert_diff (delt, mu, nu, u, v,  &
                0403                          dtau_datmos, tau_u, tau_v, dt_u, dt_v, dt_t, &
                0404                          dissipative_heat, myThid, kbot )
                0405
                0406 real,    intent(in)                        :: delt
                0407 real,    intent(in)   , dimension(:,:,:)   :: u, v, mu, nu
                0408 real,    intent(in)   , dimension(:,:)     :: dtau_datmos
                0409 real,    intent(inout), dimension(:,:)     :: tau_u, tau_v
                0410 real,    intent(inout), dimension(:,:,:)   :: dt_u, dt_v, dt_t
                0411 real,    intent(out)  , dimension(:,:,:)   :: dissipative_heat
                0412 integer, intent(in)                        :: myThid
                0413
                0414 integer, intent(in)   , dimension(:,:), optional :: kbot
                0415
                0416 real, dimension(size(u,1),size(u,2)) :: mu_delt_n, nu_n, e_n1,    &
                0417                                         f_u_delt_n1, f_v_delt_n1, &
                0418                                         delta_u_n, delta_v_n
                0419
                0420 real, dimension(size(u,1),size(u,2),size(u,3)) :: dt_u_temp, dt_v_temp
                0421
                0422 real, dimension(size(u,1),size(u,2),size(u,3)-1) :: e, f_u, f_v
                0423 integer :: i, j, kb
                0424
                0425 real    :: half_delt, cp_inv
                0426
                0427
                0428 !-----------------------------------------------------------------------
                0429
                0430  half_delt = 0.5*delt
                0431  cp_inv    = 1.0/CP_air
                0432
                0433  if (do_conserve_energy) then
                0434    dt_u_temp = dt_u
                0435    dt_v_temp = dt_v
                0436  endif
                0437
                0438  call vert_diff_down_2 &
                0439      (delt, mu, nu, u, v, dt_u, dt_v, e, f_u, f_v, &
                0440       mu_delt_n, nu_n, e_n1, f_u_delt_n1, f_v_delt_n1,  &
                0441       delta_u_n, delta_v_n, myThid, kbot)
                0442
                0443  call diff_surface (mu_delt_n, nu_n, e_n1, f_u_delt_n1, &
                0444                     dtau_datmos, tau_u, 1.0, delta_u_n, myThid)
                0445  call diff_surface (mu_delt_n, nu_n, e_n1, f_v_delt_n1, &
                0446                     dtau_datmos, tau_v, 1.0, delta_v_n, myThid)
                0447
                0448  call vert_diff_up (delt, e, f_u, delta_u_n, dt_u, myThid, kbot)
                0449  call vert_diff_up (delt, e, f_v, delta_v_n, dt_v, myThid, kbot)
                0450
                0451  if (do_conserve_energy) then
                0452     dt_u_temp = dt_u - dt_u_temp
                0453     dt_v_temp = dt_v - dt_v_temp
                0454     dissipative_heat = - cp_inv*( (u + half_delt*dt_u_temp)*dt_u_temp &
                0455                                  +(v + half_delt*dt_v_temp)*dt_v_temp )
                0456     dt_t = dt_t + dissipative_heat
                0457  else
                0458     dissipative_heat = 0.0
                0459  endif
                0460
                0461 !-----------------------------------------------------------------------
                0462
                0463 end subroutine uv_vert_diff
                0464
                0465 !#######################################################################
                0466
                0467 subroutine tq_vert_diff (delt, mu, nu, t, q,  z_full, &
                0468                          dsens_datmos, devap_datmos, sens, evap, &
                0469                          dt_t, dt_q, myThid,  kbot)
                0470
                0471
                0472 real,    intent(in)                        :: delt
                0473 real,    intent(in)   , dimension(:,:,:)   :: t, q, z_full, mu, nu
                0474 real,    intent(in)   , dimension(:,:)     :: dsens_datmos, devap_datmos
                0475 real,    intent(inout), dimension(:,:)     :: sens, evap
                0476 real,    intent(inout), dimension(:,:,:)   :: dt_t, dt_q
                0477 integer, intent(in)                        :: myThid
                0478
                0479 integer, intent(in)   , dimension(:,:), optional :: kbot
                0480
                0481 real, dimension(size(t,1),size(t,2)) :: mu_delt_n, nu_n,          &
                0482                                         e_n1, f_t_delt_n1, f_q_delt_n1, &
                0483                                         delta_t_n, delta_q_n,      &
                0484                                         flux1, dflux1_dsurf
                0485
                0486 real, dimension(size(t,1),size(t,2),size(t,3)-1) :: e, f_t, f_q
                0487 real, dimension(size(t,1),size(t,2),size(t,3)  ) :: tt
                0488
                0489 integer :: i, j, kb
                0490 real    :: gcp
                0491 !-----------------------------------------------------------------------
                0492
                0493  gcp = GRAV/CP_air
                0494  tt  = t + z_full*gcp
                0495
                0496  call vert_diff_down_2 &
                0497      (delt, mu, nu, tt, q, dt_t, dt_q, e, f_t, f_q,    &
                0498       mu_delt_n, nu_n, e_n1, f_t_delt_n1, f_q_delt_n1, &
                0499       delta_t_n, delta_q_n, myThid, kbot)
                0500
                0501
                0502  call diff_surface (mu_delt_n, nu_n, e_n1, f_t_delt_n1,  &
                0503                     dsens_datmos, sens, CP_air, delta_t_n, myThid)
                0504
                0505  call diff_surface (mu_delt_n, nu_n, e_n1, f_q_delt_n1,  &
                0506                     devap_datmos, evap, 1.0, delta_q_n, myThid)
                0507
                0508  call vert_diff_up (delt, e, f_t, delta_t_n, dt_t, myThid, kbot)
                0509  call vert_diff_up (delt, e, f_q, delta_q_n, dt_q, myThid, kbot)
                0510
                0511
                0512 !-----------------------------------------------------------------------
                0513
                0514 end subroutine tq_vert_diff
                0515
                0516 !#######################################################################
                0517
                0518 subroutine tr_vert_diff (delt, mu, nu, tr, flux, dt_tr, myThid, kbot )
                0519
                0520 real,    intent(in)                        :: delt
                0521 real,    intent(in)   , dimension(:,:,:)   :: mu, nu
                0522 real,    intent(in)   , dimension(:,:,:,:) :: tr
                0523 real,    intent(inout), dimension(:,:,:)   :: flux
                0524 real,    intent(inout), dimension(:,:,:,:) :: dt_tr
                0525 integer, intent(in)                        :: myThid
                0526
                0527 integer, intent(in)   , dimension(:,:), optional :: kbot
                0528
                0529 real, dimension(size(tr,1),size(tr,2)) :: mu_delt_n, nu_n, e_n1
                0530
                0531 real, dimension(size(tr,1),size(tr,2),size(tr,4)) :: f_delt_n1, delta_tr_n
                0532 real, dimension(size(tr,1),size(tr,2)) :: dflux_dtr
                0533
                0534 real, dimension(size(tr,1),size(tr,2),size(tr,3)-1,size(tr,4)) :: f
                0535
                0536 real, dimension(size(tr,1),size(tr,2),size(tr,3)-1) :: e
                0537 integer :: i, j, n, kb, ntr
                0538 character(len=128) :: scheme
                0539 logical, dimension(size(dt_tr,4)) :: skip_tracer_diff
                0540 !-----------------------------------------------------------------------
                0541
                0542  ntr  = size(tr,4) ! number of prognostic tracers
                0543
                0544
                0545  ! setup flags for tracer diffusion
                0546  ! this could be moved to the initialization
                0547  skip_tracer_diff(1:ntr) = .true.
                0548  do n=1,ntr
                0549    ! skip specific humidity (done separately)
                0550      if ( n == sphum ) cycle
                0551    ! skip tracers if diffusion scheme truned off
                0552 ! jmc: --- need to do something ---
                0553 !    if (query_method('diff_vert',MODEL_ATMOS,n,scheme)) then
                0554 !        if(uppercase(trim(scheme)) == 'NONE') cycle
                0555 !    endif
                0556      skip_tracer_diff(n) = .false.
                0557  enddo
                0558
                0559
                0560  dflux_dtr = 0.0
                0561
                0562  call vert_diff_down_n &
                0563      (delt, mu, nu, tr, dt_tr, e, f, mu_delt_n, nu_n,  &
                0564       e_n1, f_delt_n1, delta_tr_n, myThid, skip_tracer_diff, kbot)
                0565
                0566  do n = 1, ntr
                0567    if (skip_tracer_diff(n)) cycle
                0568    call diff_surface (mu_delt_n, nu_n, e_n1, f_delt_n1(:,:,n),  &
                0569                       dflux_dtr, flux(:,:,n), 1.0, delta_tr_n(:,:,n), myThid)
                0570
                0571    call vert_diff_up (delt, e, f(:,:,:,n), delta_tr_n(:,:,n),  &
                0572                       dt_tr(:,:,:,n), myThid, kbot)
                0573  end do
                0574
                0575 !-----------------------------------------------------------------------
                0576
                0577 end subroutine tr_vert_diff
                0578
                0579 !#######################################################################
                0580
                0581 subroutine vert_diff_down &
                0582       (delt, mu, nu, tr, dt_tr, e, f, mu_delt_n, nu_n,  &
                0583        e_n1, f_delt_n1, delta_tr_n, myThid, kbot)
                0584
                0585 !-----------------------------------------------------------------------
                0586
                0587 real,    intent(in)                         :: delt
                0588 real,    intent(in)    , dimension(:,:,:)   :: mu, nu
                0589 real,    intent(in)    , dimension(:,:,:)   :: tr
                0590 real,    intent(inout) , dimension(:,:,:)   :: dt_tr
                0591 real,    intent(out)   , dimension(:,:,:)   :: e
                0592 real,    intent(out)   , dimension(:,:,:)   :: f
                0593 real,    intent(out)   , dimension(:,:)     :: mu_delt_n, nu_n, e_n1
                0594 real,    intent(out)   , dimension(:,:)     :: f_delt_n1, delta_tr_n
                0595 integer, intent(in)                         :: myThid
                0596
                0597 integer, intent(in),    dimension(:,:), optional :: kbot
                0598
                0599 real, dimension(size(tr,1),size(tr,2),size(tr,3)) :: a, b, c, g
                0600
                0601 integer :: i, j, k, kb, nlev
                0602
                0603 !-----------------------------------------------------------------------
                0604
                0605  call explicit_tend (mu, nu, tr, dt_tr, myThid)
                0606
                0607  call compute_e  (delt, mu, nu, e, a, b, c, g, myThid)
                0608
                0609  call compute_f (dt_tr, b, c, g, f, myThid)
                0610
                0611
                0612  if (present(kbot)) then
                0613     do j=1,size(tr,2)
                0614     do i=1,size(tr,1)
                0615         kb = kbot(i,j)
                0616         mu_delt_n(i,j) =  mu(i,j,kb  )*delt
                0617              nu_n(i,j) =  nu(i,j,kb  )
                0618              e_n1(i,j) =   e(i,j,kb-1)
                0619     enddo
                0620     enddo
                0621     do j=1,size(tr,2)
                0622     do i=1,size(tr,1)
                0623         kb = kbot(i,j)
                0624          f_delt_n1(i,j) =     f(i,j,kb-1)*delt
                0625         delta_tr_n(i,j) = dt_tr(i,j,kb  )*delt
                0626     enddo
                0627     enddo
                0628  else
                0629         nlev = size(mu,3)
                0630         mu_delt_n(:,:) =       mu(:,:,nlev  )*delt
                0631              nu_n(:,:) =       nu(:,:,nlev  )
                0632              e_n1(:,:) =        e(:,:,nlev-1)
                0633         f_delt_n1(:,:) =        f(:,:,nlev-1)*delt
                0634        delta_tr_n(:,:) =    dt_tr(:,:,nlev  )*delt
                0635  endif
                0636
                0637
                0638
                0639 !-----------------------------------------------------------------------
                0640
                0641 end subroutine vert_diff_down
                0642
                0643 !#######################################################################
                0644
                0645 subroutine vert_diff_down_2 &
                0646       (delt, mu, nu, xi_1, xi_2, dt_xi_1, dt_xi_2, e, f_1, f_2, &
                0647        mu_delt_n, nu_n, e_n1, f_1_delt_n1, f_2_delt_n1,         &
                0648        delta_1_n, delta_2_n, myThid, kbot)
                0649
                0650 !-----------------------------------------------------------------------
                0651
                0652 real,    intent(in)                       :: delt
                0653 real,    intent(in)    , dimension(:,:,:) :: mu, nu, xi_1, xi_2
                0654 real,    intent(in)    , dimension(:,:,:) :: dt_xi_1, dt_xi_2
                0655 real,    intent(out)   , dimension(:,:,:) :: e, f_1, f_2
                0656 real,    intent(out)   , dimension(:,:)   :: mu_delt_n, nu_n, e_n1,    &
                0657                                              f_1_delt_n1, f_2_delt_n1, &
                0658                                              delta_1_n, delta_2_n
                0659 integer, intent(in)                       :: myThid
                0660
                0661 integer, intent(in),    dimension(:,:), optional :: kbot
                0662
                0663 real, dimension(size(xi_1,1),size(xi_1,2),size(xi_1,3)) :: a, b, c, g, &
                0664                                                       dt_xi_11, dt_xi_22
                0665
                0666 integer :: i, j, k, kb, nlev
                0667
                0668 !-----------------------------------------------------------------------
                0669
                0670 ! local copy of input
                0671   dt_xi_11 = dt_xi_1
                0672   dt_xi_22 = dt_xi_2
                0673
                0674  call explicit_tend (mu, nu, xi_1, dt_xi_11, myThid)
                0675  call explicit_tend (mu, nu, xi_2, dt_xi_22, myThid)
                0676
                0677  call compute_e (delt, mu, nu, e, a, b, c, g, myThid)
                0678
                0679  call compute_f (dt_xi_11, b, c, g, f_1, myThid)
                0680  call compute_f (dt_xi_22, b, c, g, f_2, myThid)
                0681
                0682  if (present(kbot)) then
                0683     do j=1,size(xi_1,2)
                0684     do i=1,size(xi_1,1)
                0685         kb = kbot(i,j)
                0686         mu_delt_n(i,j)  =      mu(i,j,kb  )*delt
                0687              nu_n(i,j)  =      nu(i,j,kb  )
                0688             e_n1(i,j)  =       e(i,j,kb-1)
                0689      f_1_delt_n1(i,j)  =     f_1(i,j,kb-1)*delt
                0690      f_2_delt_n1(i,j)  =     f_2(i,j,kb-1)*delt
                0691         delta_1_n(i,j)  = dt_xi_11(i,j,kb  )*delt
                0692         delta_2_n(i,j)  = dt_xi_22(i,j,kb  )*delt
                0693     enddo
                0694     enddo
                0695  else
                0696         nlev = size(mu,3)
                0697         mu_delt_n(:,:)  =      mu(:,:,nlev  )*delt
                0698              nu_n(:,:)  =      nu(:,:,nlev  )
                0699             e_n1(:,:)  =       e(:,:,nlev-1)
                0700      f_1_delt_n1(:,:)  =     f_1(:,:,nlev-1)*delt
                0701      f_2_delt_n1(:,:)  =     f_2(:,:,nlev-1)*delt
                0702         delta_1_n(:,:)  = dt_xi_11(:,:,nlev  )*delt
                0703         delta_2_n(:,:)  = dt_xi_22(:,:,nlev  )*delt
                0704  endif
                0705
                0706
                0707
                0708 !-----------------------------------------------------------------------
                0709
                0710 end subroutine vert_diff_down_2
                0711
                0712 !#######################################################################
                0713
                0714 subroutine vert_diff_down_n &
                0715       (delt, mu, nu, tr, dt_tr, e, f, mu_delt_n, nu_n,  &
                0716        e_n1, f_delt_n1, delta_tr_n, myThid, skip, kbot)
                0717
                0718 !-----------------------------------------------------------------------
                0719
                0720 real,    intent(in)                         :: delt
                0721 real,    intent(in)    , dimension(:,:,:)   :: mu, nu
                0722 real,    intent(in)    , dimension(:,:,:,:) :: tr
                0723 real,    intent(inout) , dimension(:,:,:,:) :: dt_tr
                0724 real,    intent(out)   , dimension(:,:,:)   :: e
                0725 real,    intent(out)   , dimension(:,:,:,:) :: f
                0726 real,    intent(out)   , dimension(:,:)     :: mu_delt_n, nu_n, e_n1
                0727 real,    intent(out)   , dimension(:,:,:)   :: f_delt_n1, delta_tr_n
                0728 integer, intent(in)                         :: myThid
                0729
                0730 logical, intent(in),    dimension(:),   optional :: skip
                0731 integer, intent(in),    dimension(:,:), optional :: kbot
                0732
                0733 real, dimension(size(tr,1),size(tr,2),size(tr,3)) :: a, b, c, g
                0734
                0735 integer :: i, j, k, n, kb, nlev, ntr
                0736
                0737 !-----------------------------------------------------------------------
                0738
                0739   ntr = size(tr,4)
                0740
                0741  call compute_e  (delt, mu, nu, e, a, b, c, g, myThid)
                0742
                0743  do n = 1, ntr
                0744    if (present(skip)) then
                0745        if(skip(n)) cycle
                0746    endif
                0747    call explicit_tend (mu, nu, tr(:,:,:,n), dt_tr(:,:,:,n), myThid)
                0748    call compute_f (dt_tr(:,:,:,n), b, c, g, f(:,:,:,n), myThid)
                0749  end do
                0750
                0751
                0752  if (present(kbot)) then
                0753     do j=1,size(tr,2)
                0754     do i=1,size(tr,1)
                0755         kb = kbot(i,j)
                0756         mu_delt_n(i,j) =  mu(i,j,kb  )*delt
                0757              nu_n(i,j) =  nu(i,j,kb  )
                0758             e_n1(i,j) =   e(i,j,kb-1)
                0759     enddo
                0760     enddo
                0761     do n=1,size(tr,4)
                0762     if (present(skip)) then
                0763        if(skip(n)) cycle
                0764     endif
                0765     do j=1,size(tr,2)
                0766     do i=1,size(tr,1)
                0767         kb = kbot(i,j)
                0768        f_delt_n1(i,j,n)  =     f(i,j,kb-1,n)*delt
                0769        delta_tr_n(i,j,n) = dt_tr(i,j,kb  ,n)*delt
                0770     enddo
                0771     enddo
                0772     enddo
                0773  else
                0774     nlev = size(mu,3)
                0775     mu_delt_n(:,:)  =  mu(:,:,nlev  )*delt
                0776          nu_n(:,:)  =  nu(:,:,nlev  )
                0777          e_n1(:,:)  =   e(:,:,nlev-1)
                0778     do n=1,size(tr,4)
                0779       if (present(skip)) then
                0780           if(skip(n)) cycle
                0781       endif
                0782        f_delt_n1(:,:,n) =     f(:,:,nlev-1,n)*delt
                0783       delta_tr_n(:,:,n) = dt_tr(:,:,nlev  ,n)*delt
                0784     enddo
                0785  endif
                0786
                0787
                0788
                0789 !-----------------------------------------------------------------------
                0790
                0791 end subroutine vert_diff_down_n
                0792
                0793 !#######################################################################
                0794
                0795 subroutine diff_surface (mu_delt, nu, e_n1, f_delt_n1,  &
                0796                          dflux_datmos, flux, factor, delta_xi, myThid)
                0797
                0798 !-----------------------------------------------------------------------
                0799
                0800 real, intent(in)   , dimension(:,:) :: mu_delt, nu, e_n1, f_delt_n1,  &
                0801                                        dflux_datmos
                0802 real, intent(inout), dimension(:,:) :: flux, delta_xi
                0803 real, intent(in) :: factor
                0804 integer, intent(in)                 :: myThid
                0805
                0806 !-----------------------------------------------------------------------
                0807
                0808  real, dimension(size(flux,1),size(flux,2)) :: dflux
                0809  real :: fff
                0810
                0811  fff = 1.0/factor
                0812
                0813  dflux    = - nu*(1.0 - e_n1)
                0814  delta_xi = delta_xi + mu_delt*nu*f_delt_n1
                0815
                0816  delta_xi = (delta_xi + mu_delt*flux*fff)/&
                0817                       (1.0 - mu_delt*(dflux + dflux_datmos*fff))
                0818
                0819  flux     = flux + dflux_datmos*delta_xi
                0820
                0821
                0822 !-----------------------------------------------------------------------
                0823
                0824 end subroutine diff_surface
                0825
                0826 !#######################################################################
                0827
                0828 subroutine vert_diff_up (delt, e, f, delta_xi_n, dt_xi, myThid, kbot)
                0829
                0830 !-----------------------------------------------------------------------
                0831
                0832 real,    intent(in)                      :: delt
                0833 real,    intent(in),    dimension(:,:,:) :: e, f
                0834 real,    intent(in) ,   dimension(:,:)   :: delta_xi_n
                0835 real,    intent(out),   dimension(:,:,:) :: dt_xi
                0836 integer, intent(in)                      :: myThid
                0837 integer, intent(in),    dimension(:,:), optional :: kbot
                0838
                0839 integer :: i, j, k, kb, nlev
                0840 !-----------------------------------------------------------------------
                0841
                0842  if (present(kbot)) then
                0843      do j = 1, size(dt_xi,2)
                0844      do i = 1, size(dt_xi,1)
                0845          kb = kbot(i,j)
                0846          dt_xi(i,j,kb) = delta_xi_n(i,j)/delt
                0847          do k = kb -1, 1, -1
                0848            dt_xi(i,j,k) = e(i,j,k)*dt_xi(i,j,k+1) + f(i,j,k)
                0849          end do
                0850      end do
                0851      end do
                0852  else
                0853     nlev = size(dt_xi,3)
                0854     dt_xi(:,:,nlev) = delta_xi_n/delt
                0855     do k = size(dt_xi,3)-1, 1, -1
                0856       dt_xi(:,:,k) = e(:,:,k)*dt_xi(:,:,k+1) + f(:,:,k)
                0857     end do
                0858  endif
                0859
                0860 !-----------------------------------------------------------------------
                0861
                0862 end subroutine vert_diff_up
                0863
                0864 !#######################################################################
                0865
                0866 subroutine compute_e (delt, mu, nu, e, a, b, c, g, myThid)
                0867
                0868 !-----------------------------------------------------------------------
                0869
                0870 real,    intent(in)                       :: delt
                0871 real,    intent(in)    , dimension(:,:,:) :: mu, nu
                0872 real,    intent(out)   , dimension(:,:,:) :: e, a, b, c, g
                0873 integer, intent(in)                       :: myThid
                0874
                0875 integer :: k, nlev
                0876
                0877 !-----------------------------------------------------------------------
                0878
                0879  nlev = size(mu,3)
                0880
                0881  a(:,:,1:nlev-1) = - mu(:,:,1:nlev-1)*nu(:,:,2:nlev)*delt
                0882  a(:,:,nlev    ) =   0.0
                0883  c(:,:,2:nlev  ) = - mu(:,:,2:nlev  )*nu(:,:,2:nlev)*delt
                0884  c(:,:,1       ) =   0.0
                0885
                0886  b = 1.0 - a - c
                0887
                0888  e(:,:,1)   =   - a(:,:,1)/b(:,:,1)
                0889  do  k= 2, nlev - 1
                0890     g(:,:,k) = 1.0/(b(:,:,k) + c(:,:,k)*e(:,:,k-1))
                0891     e(:,:,k) = - a(:,:,k)*g(:,:,k)
                0892  enddo
                0893
                0894 !-----------------------------------------------------------------------
                0895
                0896 end subroutine compute_e
                0897
                0898 !#######################################################################
                0899
                0900 subroutine compute_f (dt_xi, b, c, g, f, myThid)
                0901
                0902 !-----------------------------------------------------------------------
                0903 real,    intent(in)    , dimension(:,:,:) :: dt_xi, b, c, g
                0904 real,    intent(out)   , dimension(:,:,:) :: f
                0905 integer, intent(in)                       :: myThid
                0906
                0907 integer :: k
                0908 !-----------------------------------------------------------------------
                0909
                0910  f(:,:,1) =   dt_xi(:,:,1)/b(:,:,1)
                0911
                0912  do  k = 2, size(b,3)-1
                0913     f(:,:,k) = (dt_xi(:,:,k) - c(:,:,k)*f(:,:,k-1))*g(:,:,k)
                0914  enddo
                0915
                0916 !-----------------------------------------------------------------------
                0917
                0918 end subroutine compute_f
                0919
                0920 !#######################################################################
                0921
                0922 subroutine explicit_tend (mu, nu, xi, dt_xi, myThid)
                0923
                0924 !-----------------------------------------------------------------------
                0925
                0926 real,    intent(in)    , dimension(:,:,:) :: mu, nu, xi
                0927 real,    intent(inout) , dimension(:,:,:) :: dt_xi
                0928 integer, intent(in)                       :: myThid
                0929
                0930 real, dimension(size(xi,1),size(xi,2),size(xi,3)) :: fluxx
                0931
                0932 integer :: k, nlev
                0933
                0934 !-----------------------------------------------------------------------
                0935
                0936  nlev = size(mu,3)
                0937
                0938  fluxx(:,:,1)      = 0.0
                0939  fluxx(:,:,2:nlev) = nu(:,:,2:nlev)*(xi(:,:,2:nlev) - xi(:,:,1:nlev-1))
                0940
                0941  dt_xi(:,:,1:nlev-1) = dt_xi(:,:,1:nlev-1) +  &
                0942     mu(:,:,1:nlev-1)*(fluxx(:,:,2:nlev) - fluxx(:,:,1:nlev-1))
                0943  dt_xi(:,:,nlev)     = dt_xi(:,:,nlev) - mu(:,:,nlev)*fluxx(:,:,nlev)
                0944
                0945 !-----------------------------------------------------------------------
                0946
                0947 end subroutine explicit_tend
                0948
                0949 !#######################################################################
                0950
                0951 subroutine compute_mu (p_half, mu, myThid)
                0952
                0953 !-----------------------------------------------------------------------
                0954 real,    intent(in)    , dimension(:,:,:) :: p_half
                0955 real,    intent(out)   , dimension(:,:,:) :: mu
                0956 integer, intent(in)                       :: myThid
                0957
                0958 integer :: nlev
                0959 !-----------------------------------------------------------------------
                0960
                0961 nlev = size(mu,3)
                0962
                0963 mu(:,:,1:nlev) = GRAV / (p_half(:,:,2:nlev+1) -p_half(:,:,1:nlev))
                0964
                0965 !-----------------------------------------------------------------------
                0966
                0967 end subroutine compute_mu
                0968
                0969
                0970 !#######################################################################
                0971
                0972 subroutine compute_nu (diff, p_half, p_full, z_full, t, q, nu, myThid)
                0973
                0974 !-----------------------------------------------------------------------
                0975 real,    intent(in)    , dimension(:,:,:) :: diff, p_half, p_full, &
                0976                                              z_full, t, q
                0977 real,    intent(out)   , dimension(:,:,:) :: nu
                0978 integer, intent(in)                       :: myThid
                0979
                0980 real, dimension(size(t,1),size(t,2),size(t,3)) :: rho_half, tt
                0981 integer :: nlev
                0982 !-----------------------------------------------------------------------
                0983
                0984 nlev = size(nu,3)
                0985
                0986 if ( use_virtual_temp_vert_diff ) then
                0987   tt = t * (1.0 + d608*q)           ! virtual temperature
                0988 else
                0989   tt = t ! Take out virtual temperature effect here to mimic supersource
                0990 endif
                0991
                0992 rho_half(:,:,2:nlev) =  &         ! density at half levels
                0993       2.0*p_half(:,:,2:nlev)/(RDGAS*(tt(:,:,2:nlev)+tt(:,:,1:nlev-1)))
                0994
                0995 if(do_mcm_plev) then
                0996   nu(:,:,2:nlev) = GRAV*rho_half(:,:,2:nlev)*rho_half(:,:,2:nlev)*diff(:,:,2:nlev)/ &
                0997                     (p_full(:,:,2:nlev)-p_full(:,:,1:nlev-1))
                0998 else
                0999   nu(:,:,2:nlev) = rho_half(:,:,2:nlev)*diff(:,:,2:nlev) /  &
                1000                     (z_full(:,:,1:nlev-1)-z_full(:,:,2:nlev))
                1001 endif
                1002 !-----------------------------------------------------------------------
                1003
                1004 end subroutine compute_nu
                1005
                1006 !#######################################################################
                1007
                1008 end module vert_diff_mod
                1009