File indexing completed on 2018-03-02 18:37:43 UTC
view on githubraw file Latest commit add29e06 on 2018-01-31 20:35:05 UTC
b2ea1d2979 Jean*0001 module diffusivity_mod
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012 use gcm_params_mod, only: gcm_LEN_MBUF, gcm_SQZ_R, gcm_stdMsgUnit
0013 use constants_mod, only : grav, vonkarm, cp_air, rdgas, rvgas
0014
0015
0016
0017
0018
0019
0020 use monin_obukhov_mod, only : mo_diff
0021
0022 implicit none
0023 private
0024
0025
0026
0027
0028 public diffusivity, pbl_depth, molecular_diff
0029
0030
0031
0032
0033
0034
0035
0036
0037
0038
0039
0040
0041
0042
0043
0044
0045
0046
0047
0048
0049
0050
0051
0052
0053
0054
0055
0056
0057
0058
0059
0060
0061
0062
0063
0064
0065
0066
0067
0068
0069
0070
0071
0072
0073
0074
0075
0076
0077
0078
0079
0080
0081
0082
0083
0084
0085
0086
0087
0088
0089
0090
0091
0092
0093
0094
15388b052e Jean*0095
b2ea1d2979 Jean*0096
0097
0098
0099
0100
0101
0102 logical :: fixed_depth = .false.
0103 logical :: do_virtual_non_mcm = .false.
0104 real :: depth_0 = 5000.0
0105 real :: frac_inner = 0.1
0106 real :: rich_crit_pbl = 1.0
0107 real :: entr_ratio = 0.2
0108 real :: parcel_buoy = 2.0
0109 real :: znom = 1000.0
0110 logical :: free_atm_diff = .false.
0111 logical :: free_atm_skyhi_diff = .false.
0112 logical :: pbl_mcm = .false.
0113 real :: rich_crit_diff = 0.25
0114 real :: mix_len = 30.
0115 real :: rich_prandtl = 1.00
0116 real :: background_m = 0.0
0117 real :: background_t = 0.0
0118 logical :: ampns = .false.
0119
0120 real :: ampns_max = 1.0E20
0121
0122
0123
0124 namelist /diffusivity_nml/ fixed_depth, depth_0, frac_inner,&
0125 rich_crit_pbl, entr_ratio, parcel_buoy,&
0126 znom, free_atm_diff, free_atm_skyhi_diff,&
0127 pbl_mcm, rich_crit_diff, mix_len, rich_prandtl,&
0128 background_m, background_t, ampns, ampns_max, &
0129 do_virtual_non_mcm
0130
0131
0132
0133
0134
0135 real :: small = 1.e-04
0136 real :: gcp = grav/cp_air
0137 logical :: init = .false.
0138 real :: beta = 1.458e-06
0139 real :: rbop1 = 110.4
0140 real :: rbop2 = 1.405
0141
0142 real, parameter :: d608 = (rvgas-rdgas)/rdgas
0143
0144
0145
0146
0147
0148 subroutine diffusivity_init(myThid)
0149
0150 integer, intent(in) :: myThid
0151 integer :: unit, ierr, io
0152
0153 integer :: iUnit
0154
0155
0156
0157
0158
0159
0160
0161
0162 CALL BARRIER(myThid)
0163 IF ( myThid.EQ.1 ) THEN
0164
0165 WRITE(msgBuf,'(A)') 'DIFFUSIVITY_INIT: opening data.atm_gray'
0166 CALL PRINT_MESSAGE( msgBuf, gcm_stdMsgUnit, gcm_SQZ_R, myThid )
0167 CALL OPEN_COPY_DATA_FILE( &
0168 'data.atm_gray', 'DIFFUSIVITY_INIT', &
0169 iUnit, &
0170 myThid )
0171
0172 READ(UNIT=iUnit,NML=diffusivity_nml)
0173 WRITE(msgBuf,'(A)') &
0174 'DIFFUSIVITY_INIT: finished reading data.atm_gray'
0175 CALL PRINT_MESSAGE( msgBuf, gcm_stdMsgUnit, gcm_SQZ_R, myThid )
0176
15388b052e Jean*0177 #ifdef SINGLE_DISK_IO
b2ea1d2979 Jean*0178 CLOSE(iUnit)
15388b052e Jean*0179 #else
0180 CLOSE(iUnit,STATUS='DELETE')
0181 #endif /* SINGLE_DISK_IO */
b2ea1d2979 Jean*0182
0183
0184
0185
0186
0187
0188
0189
0190
0191
0192 if (frac_inner .le. 0. .or. frac_inner .ge. 1.) &
0193 CALL PRINT_ERROR( 'diffusivity_init'// &
0194 'frac_inner must be between 0 and 1', myThid)
0195
0196
0197 if (rich_crit_pbl .lt. 0.) &
0198 CALL PRINT_ERROR( 'diffusivity_init'// &
0199 'rich_crit_pbl must be greater than or equal to zero', myThid )
0200
0201
0202 if (entr_ratio .lt. 0.) &
0203 CALL PRINT_ERROR( 'diffusivity_init'// &
0204 'entr_ratio must be greater than or equal to zero', myThid )
0205
0206
0207 if (znom .le. 0.) &
0208 CALL PRINT_ERROR( 'diffusivity_init'// &
0209 'znom must be greater than zero', myThid )
0210
0211
0212 if (.not.free_atm_diff .and. free_atm_skyhi_diff)&
0213 CALL PRINT_ERROR( 'diffusivity_init'// &
0214 'free_atm_diff must be set to true if '//&
0215 'free_atm_skyhi_diff = .true.', myThid )
0216
0217
0218
0219 if (rich_crit_diff .le. 0.) &
0220 CALL PRINT_ERROR( 'diffusivity_init'// &
0221 'rich_crit_diff must be greater than zero', myThid )
0222
0223
0224 if (mix_len .lt. 0.) &
0225 CALL PRINT_ERROR( 'diffusivity_init'// &
0226 'mix_len must be greater than or equal to zero', myThid )
0227
0228
0229 if (rich_prandtl .lt. 0.) &
0230 CALL PRINT_ERROR( 'diffusivity_init'// &
0231 'rich_prandtl must be greater than or equal to zero', myThid )
0232
0233
0234 if (background_m .lt. 0.) &
0235 CALL PRINT_ERROR( 'diffusivity_init'// &
0236 'background_m must be greater than or equal to zero', myThid )
0237
0238
0239 if (background_t .lt. 0.) &
0240 CALL PRINT_ERROR( 'diffusivity_init'// &
0241 'background_t must be greater than or equal to zero', myThid )
0242
0243
0244 if (ampns_max .lt. 1.) &
0245 CALL PRINT_ERROR( 'diffusivity_init'// &
0246 'ampns_max must be greater than or equal to one', myThid )
0247
0248
0249 if (ampns .and. .not. free_atm_skyhi_diff) &
0250 CALL PRINT_ERROR( 'diffusivity_init'// &
0251 'ampns is only valid when free_atm_skyhi_diff is & also true', myThid )
0252
0253
0254
0255
0256
0257
0258
0259
0260
0261
0262
0263 init = .true.
0264
0265 ENDIF
0266 CALL BARRIER(myThid)
0267
0268 return
0269 end subroutine diffusivity_init
0270
0271
0272
0273 subroutine diffusivity(t, q, u, v, p_full, p_half, z_full, z_half, &
0274 u_star, b_star, h, k_m, k_t, myThid, kbot)
0275
0276 real, intent(in), dimension(:,:,:) :: t, q, u, v
0277 real, intent(in), dimension(:,:,:) :: p_full, p_half
0278 real, intent(in), dimension(:,:,:) :: z_full, z_half
0279 real, intent(in), dimension(:,:) :: u_star, b_star
0280 real, intent(inout), dimension(:,:,:) :: k_m, k_t
0281 real, intent(out), dimension(:,:) :: h
0282 integer, intent (in) :: myThid
0283 integer, intent(in), optional, dimension(:,:) :: kbot
0284
0285 real, dimension(size(t,1),size(t,2),size(t,3)) :: svcp,z_full_ag, &
0286 k_m_save, k_t_save
0287 real, dimension(size(t,1),size(t,2),size(t,3)+1):: z_half_ag
0288 real, dimension(size(t,1),size(t,2)) :: z_surf
0289 integer :: i,j,k,nlev,nlat,nlon
0290
0291 if(.not.init) call diffusivity_init(myThid)
0292
0293 nlev = size(t,3)
0294
0295 k_m_save = k_m
0296 k_t_save = k_t
0297
0298
0299 if (present(kbot)) then
0300 nlat = size(t,2)
0301 nlon = size(t,1)
0302 do j=1,nlat
0303 do i=1,nlon
0304 z_surf(i,j) = z_half(i,j,kbot(i,j)+1)
0305 enddo
0306 enddo
0307 else
0308 z_surf(:,:) = z_half(:,:,nlev+1)
0309 end if
0310
0311
0312 do k = 1, nlev
0313
0314 z_full_ag(:,:,k) = z_full(:,:,k) - z_surf(:,:)
0315 z_half_ag(:,:,k) = z_half(:,:,k) - z_surf(:,:)
0316
0317 if (do_virtual_non_mcm) then
0318 svcp(:,:,k) = t(:,:,k)*(1. + d608*q(:,:,k)) + gcp*(z_full_ag(:,:,k))
0319 else
0320 svcp(:,:,k) = t(:,:,k) + gcp*(z_full_ag(:,:,k))
0321 endif
0322
0323 end do
0324 z_half_ag(:,:,nlev+1) = z_half(:,:,nlev+1) - z_surf(:,:)
0325
0326 if(fixed_depth) then
0327 h = depth_0
0328 else
0329 call pbl_depth(svcp,u,v,z_full_ag,u_star,b_star,h, myThid, kbot=kbot)
0330 end if
0331
0332 if(pbl_mcm) then
0333 call diffusivity_pbl_mcm (u,v, t, p_full, p_half, &
0334 z_full_ag, z_half_ag, h, k_m, k_t, myThid)
0335 else
0336 call diffusivity_pbl (svcp, u, v, z_half_ag, h, u_star, b_star,&
0337 k_m, k_t, myThid, kbot=kbot)
0338 end if
0339 if(free_atm_diff) &
0340 call diffusivity_free (svcp, u, v, z_full_ag, z_half_ag, h, k_m, k_t, myThid)
0341
0342 k_m = k_m + k_m_save
0343 k_t = k_t + k_t_save
0344
0345
0346
0347 if(entr_ratio .gt. 0. .and. .not. fixed_depth) &
0348 call diffusivity_entr(svcp,z_full_ag,h,u_star,b_star,k_m,k_t, myThid)
0349
0350
0351 if(background_m.gt.0.0) k_m = max(k_m,background_m)
0352 if(background_t.gt.0.0) k_t = max(k_t,background_t)
0353
0354 return
0355 end subroutine diffusivity
0356
0357
0358
0359 subroutine pbl_depth(t, u, v, z, u_star, b_star, h, myThid, kbot)
0360
0361 real, intent(in) , dimension(:,:,:) :: t, u, v, z
0362 real, intent(in) , dimension(:,:) :: u_star,b_star
0363 real, intent(out), dimension(:,:) :: h
0364 integer, intent(in) :: myThid
0365 integer,intent(in) , optional, dimension(:,:) :: kbot
0366
0367 real, dimension(size(t,1),size(t,2),size(t,3)) :: rich
0368 real, dimension(size(t,1),size(t,2)) :: ws,k_t_ref,&
0369 h_inner,tbot
0370 real :: rich1, rich2,&
0371 h1,h2,svp,t1,t2
0372 integer, dimension(size(t,1),size(t,2)) :: ibot
0373 integer :: i,j,k,nlon,&
0374 nlat, nlev
0375
0376 nlev = size(t,3)
0377 nlat = size(t,2)
0378 nlon = size(t,1)
0379
0380
0381 if (present(kbot)) then
0382 ibot(:,:) = kbot
0383 do j = 1,nlat
0384 do i = 1,nlon
0385 tbot(i,j) = t(i,j,ibot(i,j))
0386 enddo
0387 enddo
0388 else
0389 ibot(:,:) = nlev
0390 tbot(:,:) = t(:,:,nlev)
0391 end if
0392
0393
0394 do k = 1,nlev
0395 rich(:,:,k) = z(:,:,k)*grav*(t(:,:,k)-tbot(:,:))/tbot(:,:)&
0396 /(u(:,:,k)*u(:,:,k) + v(:,:,k)*v(:,:,k) + small )
0397 end do
0398
0399
0400
0401
0402
0403 h_inner(:,:)=frac_inner*znom
0404
0405 ws = max(small,ws/vonkarm/h_inner)
0406
0407 do j = 1, nlat
0408 do i = 1, nlon
0409
0410
0411 if (b_star(i,j).le.0.) then
0412
0413 h1 = z(i,j,ibot(i,j))
0414 h(i,j) = h1
0415 rich1 = rich(i,j,ibot(i,j))
0416 do k = ibot(i,j)-1, 1, -1
0417 rich2 = rich(i,j,k)
0418 h2 = z(i,j,k)
0419 if(rich2.gt.rich_crit_pbl) then
0420 h(i,j) = h2 + (h1 - h2)*(rich2 - rich_crit_pbl)&
0421 /(rich2 - rich1 )
0422 go to 10
0423 endif
0424 rich1 = rich2
0425 h1 = h2
0426 enddo
0427
0428
0429 else
0430
0431 svp = tbot(i,j)*(1.+ &
0432 (parcel_buoy*u_star(i,j)*b_star(i,j)/grav/ws(i,j)) )
0433 h1 = z(i,j,ibot(i,j))
0434 h(i,j) = h1
0435 t1 = tbot(i,j)
0436 do k = ibot(i,j)-1 , 1, -1
0437 h2 = z(i,j,k)
0438 t2 = t(i,j,k)
0439 if (t2.gt.svp) then
0440 h(i,j) = h2 + (h1 - h2)*(t2 - svp)/(t2 - t1 )
0441 go to 10
0442 end if
0443 h1 = h2
0444 t1 = t2
0445 enddo
0446
0447 end if
0448 10 continue
0449 enddo
0450 enddo
0451
0452 return
0453 end subroutine pbl_depth
0454
0455
0456
0457 subroutine diffusivity_pbl(t, u, v, z_half, h, u_star, b_star, &
0458 k_m, k_t, myThid, kbot)
0459
0460 real, intent(in) , dimension(:,:,:) :: t, u, v, z_half
0461 real, intent(in) , dimension(:,:) :: h, u_star, b_star
0462 real, intent(inout) , dimension(:,:,:) :: k_m, k_t
0463 integer, intent (in) :: myThid
0464 integer, intent(in) , optional, dimension(:,:) :: kbot
0465
0466 real, dimension(size(t,1),size(t,2)) :: h_inner, k_m_ref,&
0467 k_t_ref, factor
0468 real, dimension(size(t,1),size(t,2),size(t,3)+1) :: zm
0469 real :: h_inner_max
0470 integer :: i,j, k, kk, nlev
0471
0472 nlev = size(t,3)
0473
0474
0475
0476
0477 zm = z_half
0478 if (present(kbot)) then
0479 where(zm < 0.)
0480 zm = 0.
0481 end where
0482 end if
0483
0484 h_inner = frac_inner*h
0485 h_inner_max = maxval(h_inner)
0486
0487 kk = nlev
0488 do k = 2, nlev
0489 if( minval(zm(:,:,k)) < h_inner_max) then
0490 kk = k
0491 exit
0492 end if
0493 end do
0494
0495 k_m = 0.0
0496 k_t = 0.0
0497
0498
0499
0500 myThid )
0501
0502 do k = 2, nlev
0503 where(zm(:,:,k) >= h_inner .and. zm(:,:,k) < h)
0504 factor = (zm(:,:,k)/h_inner)* &
0505 (1.0 - (zm(:,:,k) - h_inner)/(h - h_inner))**2
0506 k_m(:,:,k) = k_m_ref*factor
0507 k_t(:,:,k) = k_t_ref*factor
0508 end where
0509 end do
0510
0511 return
0512 end subroutine diffusivity_pbl
0513
0514
0515
0516 subroutine diffusivity_pbl_mcm(u, v, t, p_full, p_half, z_full, z_half, &
0517 h, k_m, k_t, myThid)
0518
0519 real, intent(in) , dimension(:,:,:) :: u, v, t, z_full, z_half
0520 real, intent(in) , dimension(:,:,:) :: p_full, p_half
0521 real, intent(in) , dimension(:,:) :: h
0522 real, intent(inout) , dimension(:,:,:) :: k_m, k_t
0523 integer, intent(in) :: myThid
0524
0525 integer :: k, nlev
0526 real, dimension(size(z_full,1),size(z_full,2)) :: elmix, htcrit
0527 real, dimension(size(z_full,1),size(z_full,2)) :: delta_u, delta_v, delta_z
0528
0529 real :: htcrit_ss
0530 real :: h_ss
0531 real, dimension(size(z_full,1),size(z_full,2)) :: sig_half, z_half_ss, elmix_ss
0532
0533
0534
0535
0536
0537
0538
0539
0540
0541 real :: tsfc = 288.16
0542 real :: salaps = -6.5e-3
0543
0544 nlev = size(z_full,3)
0545
0546 k_m = 0.
0547
0548 h_ss = depth_0
0549 htcrit_ss = frac_inner*h_ss
0550
0551 do k = 2, nlev
0552
0553
0554
0555
0556
0557
0558
0559
0560
0561
0562
0563
0564
0565 sig_half = p_half(:,:,k)/p_half(:,:,nlev+1)
0566 z_half_ss = -rdgas * .5*(tsfc+tsfc*(sig_half**(-rdgas*salaps/grav))) * alog(sig_half)/grav
0567
0568
0569 elmix_ss = 0.
0570
0571 where (z_half_ss < htcrit_ss .and. z_half_ss > 0.)
0572 elmix_ss = vonkarm*z_half_ss
0573 endwhere
0574 where (z_half_ss >= htcrit_ss .and. z_half_ss < h_ss)
0575 elmix_ss = vonkarm*htcrit_ss*(h_ss-z_half_ss)/(h_ss-htcrit_ss)
0576 endwhere
0577
0578 delta_z = rdgas*0.5*(t(:,:,k)+t(:,:,k-1))*(p_full(:,:,k)-p_full(:,:,k-1))/&
0579 (grav*p_half(:,:,k))
0580 delta_u = u(:,:,k-1) - u(:,:,k)
0581 delta_v = v(:,:,k-1) - v(:,:,k)
0582
0583 k_m(:,:,k) = elmix_ss * elmix_ss *&
0584 sqrt(delta_u*delta_u + delta_v*delta_v)/delta_z
0585
0586 end do
0587
0588 k_t = k_m
0589
0590 return
0591 end subroutine diffusivity_pbl_mcm
0592
0593
0594
0595 subroutine diffusivity_free(t, u, v, z, zz, h, k_m, k_t, myThid)
0596
0597 real, intent(in) , dimension(:,:,:) :: t, u, v, z, zz
0598 real, intent(in) , dimension(:,:) :: h
0599 real, intent(inout) , dimension(:,:,:) :: k_m, k_t
0600 integer, intent(in) :: myThid
0601
0602 real, dimension(size(t,1),size(t,2)) :: dz, b, speed2, rich, fri, &
0603 alpz, fri2
0604 integer :: k
0605
0606 do k = 2, size(t,3)
0607
0608
0609
0610
0611
0612 dz = z(:,:,k-1) - z(:,:,k)
0613 b = grav*(t(:,:,k-1)-t(:,:,k))/t(:,:,k)
0614 speed2 = (u(:,:,k-1) - u(:,:,k))**2 + (v(:,:,k-1) - v(:,:,k))**2
0615 rich= b*dz/(speed2+small)
0616 rich = max(rich, 0.0)
0617
0618 if (free_atm_skyhi_diff) then
0619
0620
0621
0622
0623
0624 where (rich(:,:) >= rich_crit_diff)
0625 fri2(:,:) = 0.0
0626 elsewhere
0627 fri2(:,:) = (1.0 - rich/rich_crit_diff)**2
0628 endwhere
0629 endif
0630
0631
0632
0633
0634
0635 if (ampns) then
0636 alpz(:,:) = MIN ( (1. + 1.e-04*(dz(:,:)**1.5)), ampns_max)
0637 rich(:,:) = rich(:,:) / alpz(:,:)
0638 endif
0639
0640
0641
0642
0643
0644
0645 fri(:,:) = (1.0 - rich/rich_crit_diff)**2
0646
0647
0648
0649
0650
0651
0652
0653
0654
0655
0656
0657 if (free_atm_skyhi_diff) then
0658
0659
0660
0661
0662
0663 if (ampns) then
0664 where (rich < rich_crit_diff .and. zz(:,:,k) > h)
0665 k_m(:,:,k) = mix_len*mix_len*sqrt(speed2)*fri(:,:)* &
0666 ( 1. + 1.e-04*(dz(:,:)**1.5))/dz
0667 k_t(:,:,k) = k_m(:,:,k)* (0.1 + 0.9*fri2(:,:))
0668 end where
0669 else
0670 where (rich < rich_crit_diff .and. zz(:,:,k) > h)
0671 k_m(:,:,k) = mix_len*mix_len*sqrt(speed2)*fri(:,:)/dz
0672 k_t(:,:,k) = k_m(:,:,k)* (0.1 + 0.9*fri2(:,:))
0673 end where
0674 endif
0675 else
0676
0677
0678
0679
0680
0681 where (rich < rich_crit_diff .and. zz(:,:,k) > h)
0682 k_t(:,:,k) = mix_len*mix_len*sqrt(speed2)*fri(:,:)/dz
0683 k_m(:,:,k) = k_t(:,:,k)*rich_prandtl
0684 end where
0685 end if
0686 end do
0687
0688 end subroutine diffusivity_free
0689
0690
0691
0692 subroutine molecular_diff ( temp, press, k_m, k_t, myThid)
0693
0694 real, intent(in), dimension (:,:,:) :: temp, press
0695 real, intent(inout), dimension (:,:,:) :: k_m, k_t
0696 integer, intent(in) :: myThid
0697
0698 real, dimension (size(temp,1), size(temp,2)) :: temp_half, &
0699 rho_half, rbop2d
0700 integer :: k
0701
0702
0703
0704 do k=2,size(temp,3)
0705 temp_half(:,:) = 0.5*(temp(:,:,k) + temp(:,:,k-1))
0706 rho_half(:,:) = press(:,:,k)/(rdgas*temp_half(:,:) )
0707 rbop2d(:,:) = beta*temp_half(:,:)*sqrt(temp_half(:,:))/ &
0708 (rho_half(:,:)*(temp_half(:,:)+rbop1))
0709 k_m(:,:,k) = rbop2d(:,:)
0710 k_t(:,:,k) = rbop2d(:,:)*rbop2
0711 end do
0712
0713 k_m(:,:,1) = 0.0
0714 k_t(:,:,1) = 0.0
0715
0716 end subroutine molecular_diff
0717
0718
0719
0720 subroutine diffusivity_entr(t, z, h, u_star, b_star, k_m, k_t, myThid)
0721
0722 real, intent(in) , dimension(:,:,:) :: t, z
0723 real, intent(in) , dimension(:,:) :: h, u_star, b_star
0724 real, intent(inout) , dimension(:,:,:) :: k_m, k_t
0725 integer, intent(in) :: myThid
0726
0727 integer :: k, nlev
0728
0729 nlev=size(t,3)
0730
0731 do k = 2,nlev
0732 where (b_star .gt. 0. .and. z(:,:,k-1) .gt. h .and. &
0733 z(:,:,k) .le. h)
0734 k_t(:,:,k) = (z(:,:,k-1)-z(:,:,k))*entr_ratio*t(:,:,k)* &
0735 u_star*b_star/grav/max(small,t(:,:,k-1)-t(:,:,k))
0736 k_m(:,:,k) = k_t(:,:,k)
0737 end where
0738 enddo
0739 end subroutine diffusivity_entr
0740
0741
0742
0743 end module diffusivity_mod
0744