File indexing completed on 2018-03-02 18:37:48 UTC
view on githubraw file Latest commit add29e06 on 2018-01-31 20:35:05 UTC
b2ea1d2979 Jean*0001 module vert_diff_mod
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011 use constants_mod, only: GRAV, RDGAS, RVGAS, CP_air
0012
0013
0014
0015
0016
0017
0018
0019
0020 implicit none
0021 private
0022
0023
0024
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
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
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
0057
0058
0059
0060
0061 real, parameter :: d608 = (RVGAS-RDGAS)/RDGAS
0062
0063
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
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
0080
0081 CALL BARRIER(myThid)
0082 IF ( myThid.EQ.1 ) THEN
0083
0084
0085
0086
0087
0088 sphum = 2
0089
0090
0091
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
0107
0108
0109
0110
0111
0112
0113
0114 do_init = .false.
0115
0116 endif
0117
0118
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
0167
0168
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
0179 diff_m, diff_t, p_half, p_full, &
0180 z_full, tau_u, tau_v, dtau_datmos, &
0181
0182 dt_u, dt_v, dt_t, dt_q, &
0183
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
0197 real, intent(in) , dimension(:,:) :: dtau_datmos
0198 real, intent(inout), dimension(:,:) :: tau_u, tau_v
0199
0200 real, intent(inout), dimension(:,:,:) :: dt_u, dt_v, dt_t
0201 real, intent(in), dimension(:,:,:) :: dt_q
0202
0203 real, intent(out) , dimension(:,:,:) :: dissipative_heat
0204
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
0228
0229
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
0236
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
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
0247 call compute_nu (diff_t, p_half, p_full, z_full, t, q, nu, myThid)
0248
0249
0250
0251
0252
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
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
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
0320
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
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
0539 logical, dimension(size(dt_tr,4)) :: skip_tracer_diff
0540
0541
0542 ntr = size(tr,4)
0543
0544
0545
0546
0547 skip_tracer_diff(1:ntr) = .true.
0548 do n=1,ntr
0549
0550 if ( n == sphum ) cycle
0551
0552
0553
0554
0555
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
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)
0988 else
0989 tt = t
0990 endif
0991
0992 rho_half(:,:,2:nlev) = &
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