1 |
jmc |
1.2 |
! $Header: /u/gcmpack/MITgcm/pkg/atm_phys/diffusivity_mod.F90,v 1.1 2013/05/08 22:14:14 jmc Exp $ |
2 |
jmc |
1.1 |
! $Name: $ |
3 |
|
|
|
4 |
|
|
module diffusivity_mod |
5 |
|
|
|
6 |
|
|
!======================================================================= |
7 |
|
|
! |
8 |
|
|
! DIFFUSIVITY MODULE |
9 |
|
|
! |
10 |
|
|
! Routines for computing atmospheric diffusivities in the |
11 |
|
|
! planetary boundary layer and in the free atmosphere |
12 |
|
|
! |
13 |
|
|
!======================================================================= |
14 |
|
|
|
15 |
|
|
use gcm_params_mod, only: gcm_LEN_MBUF, gcm_SQZ_R, gcm_stdMsgUnit |
16 |
|
|
use constants_mod, only : grav, vonkarm, cp_air, rdgas, rvgas |
17 |
|
|
|
18 |
|
|
!use fms_mod, only: error_mesg, FATAL, file_exist, & |
19 |
|
|
! check_nml_error, open_namelist_file, & |
20 |
|
|
! mpp_pe, mpp_root_pe, close_file, stdlog, & |
21 |
|
|
! write_version_number |
22 |
|
|
|
23 |
|
|
use monin_obukhov_mod, only : mo_diff |
24 |
|
|
|
25 |
|
|
implicit none |
26 |
|
|
private |
27 |
|
|
|
28 |
|
|
! public interfaces |
29 |
|
|
!======================================================================= |
30 |
|
|
|
31 |
|
|
public diffusivity, pbl_depth, molecular_diff |
32 |
|
|
|
33 |
|
|
!======================================================================= |
34 |
|
|
|
35 |
|
|
! form of iterfaces |
36 |
|
|
|
37 |
|
|
!======================================================================= |
38 |
|
|
! subroutine diffusivity (t, q, u, v, p_full, p_half, z_full, z_half, |
39 |
|
|
! u_star, b_star, h, k_m, k_t) |
40 |
|
|
|
41 |
|
|
! input: |
42 |
|
|
|
43 |
|
|
! t : real, dimension(:,:,:) -- (:,:,pressure), third index running |
44 |
|
|
! from top of atmosphere to bottom |
45 |
|
|
! temperature (K) |
46 |
|
|
! |
47 |
|
|
! q : real, dimension(:,:,:) |
48 |
|
|
! water vapor specific humidity (nondimensional) |
49 |
|
|
! |
50 |
|
|
! u : real, dimension(:,:) |
51 |
|
|
! zonal wind (m/s) |
52 |
|
|
! |
53 |
|
|
! v : real, dimension(:,:,:) |
54 |
|
|
! meridional wind (m/s) |
55 |
|
|
! |
56 |
|
|
! z_full : real, dimension(:,:,: |
57 |
|
|
! height of full levels (m) |
58 |
|
|
! 1 = top of atmosphere; size(p_half,3) = surface |
59 |
|
|
! size(z_full,3) = size(t,3) |
60 |
|
|
! |
61 |
|
|
! z_half : real, dimension(:,:,:) |
62 |
|
|
! height of half levels (m) |
63 |
|
|
! size(z_half,3) = size(t,3) +1 |
64 |
|
|
! z_half(:,:,size(z_half,3)) must be height of surface! |
65 |
|
|
! (if you are not using eta-model) |
66 |
|
|
! |
67 |
|
|
! u_star: real, dimension(:,:) |
68 |
|
|
! friction velocity (m/s) |
69 |
|
|
! |
70 |
|
|
! b_star: real, dimension(:,:) |
71 |
|
|
! buoyancy scale (m/s**2) |
72 |
|
|
|
73 |
|
|
! (u_star and b_star can be obtained by calling |
74 |
|
|
! mo_drag in monin_obukhov_mod) |
75 |
|
|
|
76 |
|
|
! output: |
77 |
|
|
|
78 |
|
|
! h : real, dimension(:,:,) |
79 |
|
|
! depth of planetary boundary layer (m) |
80 |
|
|
! |
81 |
|
|
! k_m : real, dimension(:,:,:) |
82 |
|
|
! diffusivity for momentum (m**2/s) |
83 |
|
|
! |
84 |
|
|
! defined at half-levels |
85 |
|
|
! size(k_m,3) should be at least as large as size(t,3) |
86 |
|
|
! only the returned values at |
87 |
|
|
! levels 2 to size(t,3) are meaningful |
88 |
|
|
! other values will be returned as zero |
89 |
|
|
! |
90 |
|
|
! k_t : real, dimension(:,:,:) |
91 |
|
|
! diffusivity for temperature and scalars (m**2/s) |
92 |
|
|
! |
93 |
|
|
! |
94 |
|
|
!======================================================================= |
95 |
|
|
|
96 |
|
|
!--------------------- version number ---------------------------------- |
97 |
|
|
|
98 |
jmc |
1.2 |
character(len=128) :: version = '$Id: diffusivity_mod.F90,v 1.1 2013/05/08 22:14:14 jmc Exp $' |
99 |
jmc |
1.1 |
character(len=128) :: tag = '$Name: $' |
100 |
|
|
|
101 |
|
|
!======================================================================= |
102 |
|
|
|
103 |
|
|
! DEFAULT VALUES OF NAMELIST PARAMETERS: |
104 |
|
|
|
105 |
|
|
logical :: fixed_depth = .false. |
106 |
|
|
logical :: do_virtual_non_mcm = .false. ! applies only to non-'mcm' pbl scheme |
107 |
|
|
real :: depth_0 = 5000.0 |
108 |
|
|
real :: frac_inner = 0.1 |
109 |
|
|
real :: rich_crit_pbl = 1.0 |
110 |
|
|
real :: entr_ratio = 0.2 |
111 |
|
|
real :: parcel_buoy = 2.0 |
112 |
|
|
real :: znom = 1000.0 |
113 |
|
|
logical :: free_atm_diff = .false. |
114 |
|
|
logical :: free_atm_skyhi_diff = .false. |
115 |
|
|
logical :: pbl_mcm = .false. |
116 |
|
|
real :: rich_crit_diff = 0.25 |
117 |
|
|
real :: mix_len = 30. |
118 |
|
|
real :: rich_prandtl = 1.00 |
119 |
|
|
real :: background_m = 0.0 |
120 |
|
|
real :: background_t = 0.0 |
121 |
|
|
logical :: ampns = .false. ! include delta z factor in |
122 |
|
|
! defining ri ? |
123 |
|
|
real :: ampns_max = 1.0E20 ! limit to reduction factor |
124 |
|
|
! applied to ri due to delta z |
125 |
|
|
! factor |
126 |
|
|
|
127 |
|
|
namelist /diffusivity_nml/ fixed_depth, depth_0, frac_inner,& |
128 |
|
|
rich_crit_pbl, entr_ratio, parcel_buoy,& |
129 |
|
|
znom, free_atm_diff, free_atm_skyhi_diff,& |
130 |
|
|
pbl_mcm, rich_crit_diff, mix_len, rich_prandtl,& |
131 |
|
|
background_m, background_t, ampns, ampns_max, & |
132 |
|
|
do_virtual_non_mcm |
133 |
|
|
|
134 |
|
|
!======================================================================= |
135 |
|
|
|
136 |
|
|
! OTHER MODULE VARIABLES |
137 |
|
|
|
138 |
|
|
real :: small = 1.e-04 |
139 |
|
|
real :: gcp = grav/cp_air |
140 |
|
|
logical :: init = .false. |
141 |
|
|
real :: beta = 1.458e-06 |
142 |
|
|
real :: rbop1 = 110.4 |
143 |
|
|
real :: rbop2 = 1.405 |
144 |
|
|
|
145 |
|
|
real, parameter :: d608 = (rvgas-rdgas)/rdgas |
146 |
|
|
|
147 |
|
|
contains |
148 |
|
|
|
149 |
|
|
!======================================================================= |
150 |
|
|
|
151 |
|
|
subroutine diffusivity_init(myThid) |
152 |
|
|
|
153 |
|
|
integer, intent(in) :: myThid |
154 |
|
|
integer :: unit, ierr, io |
155 |
|
|
|
156 |
|
|
integer :: iUnit |
157 |
|
|
CHARACTER*(gcm_LEN_MBUF) :: msgBuf |
158 |
|
|
|
159 |
|
|
!------------------- read namelist input ------------------------------- |
160 |
|
|
|
161 |
|
|
! read namelist and copy to logfile |
162 |
|
|
|
163 |
|
|
! _BARRIER |
164 |
|
|
! _BEGIN_MASTER(myThid) |
165 |
|
|
CALL BARRIER(myThid) |
166 |
|
|
IF ( myThid.EQ.1 ) THEN |
167 |
|
|
|
168 |
|
|
WRITE(msgBuf,'(A)') 'DIFFUSIVITY_INIT: opening data.atm_gray' |
169 |
|
|
CALL PRINT_MESSAGE( msgBuf, gcm_stdMsgUnit, gcm_SQZ_R, myThid ) |
170 |
|
|
CALL OPEN_COPY_DATA_FILE( & |
171 |
|
|
'data.atm_gray', 'DIFFUSIVITY_INIT', & |
172 |
|
|
iUnit, & |
173 |
|
|
myThid ) |
174 |
|
|
! Read parameters from open data file |
175 |
|
|
READ(UNIT=iUnit,NML=diffusivity_nml) |
176 |
|
|
WRITE(msgBuf,'(A)') & |
177 |
|
|
'DIFFUSIVITY_INIT: finished reading data.atm_gray' |
178 |
|
|
CALL PRINT_MESSAGE( msgBuf, gcm_stdMsgUnit, gcm_SQZ_R, myThid ) |
179 |
|
|
! Close the open data file |
180 |
jmc |
1.2 |
#ifdef SINGLE_DISK_IO |
181 |
jmc |
1.1 |
CLOSE(iUnit) |
182 |
jmc |
1.2 |
#else |
183 |
|
|
CLOSE(iUnit,STATUS='DELETE') |
184 |
|
|
#endif /* SINGLE_DISK_IO */ |
185 |
jmc |
1.1 |
|
186 |
|
|
! if (file_exist('input.nml')) then |
187 |
|
|
! unit = open_namelist_file () |
188 |
|
|
! ierr=1; do while (ierr /= 0) |
189 |
|
|
! read (unit, nml=diffusivity_nml, iostat=io, end=10) |
190 |
|
|
! ierr = check_nml_error(io,'diffusivity_nml') |
191 |
|
|
! enddo |
192 |
|
|
! 10 call close_file (unit) |
193 |
|
|
|
194 |
|
|
!------------------- dummy checks -------------------------------------- |
195 |
|
|
if (frac_inner .le. 0. .or. frac_inner .ge. 1.) & |
196 |
|
|
CALL PRINT_ERROR( 'diffusivity_init'// & |
197 |
|
|
'frac_inner must be between 0 and 1', myThid) |
198 |
|
|
! call error_mesg ('diffusivity_init', & |
199 |
|
|
! 'frac_inner must be between 0 and 1', FATAL) |
200 |
|
|
if (rich_crit_pbl .lt. 0.) & |
201 |
|
|
CALL PRINT_ERROR( 'diffusivity_init'// & |
202 |
|
|
'rich_crit_pbl must be greater than or equal to zero', myThid ) |
203 |
|
|
! call error_mesg ('diffusivity_init', & |
204 |
|
|
! 'rich_crit_pbl must be greater than or equal to zero', FATAL) |
205 |
|
|
if (entr_ratio .lt. 0.) & |
206 |
|
|
CALL PRINT_ERROR( 'diffusivity_init'// & |
207 |
|
|
'entr_ratio must be greater than or equal to zero', myThid ) |
208 |
|
|
! call error_mesg ('diffusivity_init', & |
209 |
|
|
! 'entr_ratio must be greater than or equal to zero', FATAL) |
210 |
|
|
if (znom .le. 0.) & |
211 |
|
|
CALL PRINT_ERROR( 'diffusivity_init'// & |
212 |
|
|
'znom must be greater than zero', myThid ) |
213 |
|
|
! call error_mesg ('diffusivity_init', & |
214 |
|
|
! 'znom must be greater than zero', FATAL) |
215 |
|
|
if (.not.free_atm_diff .and. free_atm_skyhi_diff)& |
216 |
|
|
CALL PRINT_ERROR( 'diffusivity_init'// & |
217 |
|
|
'free_atm_diff must be set to true if '//& |
218 |
|
|
'free_atm_skyhi_diff = .true.', myThid ) |
219 |
|
|
! call error_mesg ('diffusivity_init', & |
220 |
|
|
! 'free_atm_diff must be set to true if '//& |
221 |
|
|
! 'free_atm_skyhi_diff = .true.', FATAL) |
222 |
|
|
if (rich_crit_diff .le. 0.) & |
223 |
|
|
CALL PRINT_ERROR( 'diffusivity_init'// & |
224 |
|
|
'rich_crit_diff must be greater than zero', myThid ) |
225 |
|
|
! call error_mesg ('diffusivity_init', & |
226 |
|
|
! 'rich_crit_diff must be greater than zero', FATAL) |
227 |
|
|
if (mix_len .lt. 0.) & |
228 |
|
|
CALL PRINT_ERROR( 'diffusivity_init'// & |
229 |
|
|
'mix_len must be greater than or equal to zero', myThid ) |
230 |
|
|
! call error_mesg ('diffusivity_init', & |
231 |
|
|
! 'mix_len must be greater than or equal to zero', FATAL) |
232 |
|
|
if (rich_prandtl .lt. 0.) & |
233 |
|
|
CALL PRINT_ERROR( 'diffusivity_init'// & |
234 |
|
|
'rich_prandtl must be greater than or equal to zero', myThid ) |
235 |
|
|
! call error_mesg ('diffusivity_init', & |
236 |
|
|
! 'rich_prandtl must be greater than or equal to zero', FATAL) |
237 |
|
|
if (background_m .lt. 0.) & |
238 |
|
|
CALL PRINT_ERROR( 'diffusivity_init'// & |
239 |
|
|
'background_m must be greater than or equal to zero', myThid ) |
240 |
|
|
! call error_mesg ('diffusivity_init', & |
241 |
|
|
! 'background_m must be greater than or equal to zero', FATAL) |
242 |
|
|
if (background_t .lt. 0.) & |
243 |
|
|
CALL PRINT_ERROR( 'diffusivity_init'// & |
244 |
|
|
'background_t must be greater than or equal to zero', myThid ) |
245 |
|
|
! call error_mesg ('diffusivity_init', & |
246 |
|
|
! 'background_t must be greater than or equal to zero', FATAL) |
247 |
|
|
if (ampns_max .lt. 1.) & |
248 |
|
|
CALL PRINT_ERROR( 'diffusivity_init'// & |
249 |
|
|
'ampns_max must be greater than or equal to one', myThid ) |
250 |
|
|
! call error_mesg ('diffusivity_init', & |
251 |
|
|
! 'ampns_max must be greater than or equal to one', FATAL) |
252 |
|
|
if (ampns .and. .not. free_atm_skyhi_diff) & |
253 |
|
|
CALL PRINT_ERROR( 'diffusivity_init'// & |
254 |
|
|
'ampns is only valid when free_atm_skyhi_diff is & also true', myThid ) |
255 |
|
|
! call error_mesg ('diffusivity_init', & |
256 |
|
|
! 'ampns is only valid when free_atm_skyhi_diff is & |
257 |
|
|
! & also true', FATAL) |
258 |
|
|
|
259 |
|
|
! endif !end of reading input.nml |
260 |
|
|
|
261 |
|
|
!---------- output namelist to log------------------------------------- |
262 |
|
|
|
263 |
|
|
! call write_version_number(version, tag) |
264 |
|
|
! if ( mpp_pe() == mpp_root_pe() ) write (stdlog(), nml=diffusivity_nml) |
265 |
|
|
|
266 |
|
|
init = .true. |
267 |
|
|
|
268 |
|
|
ENDIF |
269 |
|
|
CALL BARRIER(myThid) |
270 |
|
|
|
271 |
|
|
return |
272 |
|
|
end subroutine diffusivity_init |
273 |
|
|
|
274 |
|
|
!======================================================================= |
275 |
|
|
|
276 |
|
|
subroutine diffusivity(t, q, u, v, p_full, p_half, z_full, z_half, & |
277 |
|
|
u_star, b_star, h, k_m, k_t, myThid, kbot) |
278 |
|
|
|
279 |
|
|
real, intent(in), dimension(:,:,:) :: t, q, u, v |
280 |
|
|
real, intent(in), dimension(:,:,:) :: p_full, p_half |
281 |
|
|
real, intent(in), dimension(:,:,:) :: z_full, z_half |
282 |
|
|
real, intent(in), dimension(:,:) :: u_star, b_star |
283 |
|
|
real, intent(inout), dimension(:,:,:) :: k_m, k_t |
284 |
|
|
real, intent(out), dimension(:,:) :: h |
285 |
|
|
integer, intent (in) :: myThid |
286 |
|
|
integer, intent(in), optional, dimension(:,:) :: kbot |
287 |
|
|
|
288 |
|
|
real, dimension(size(t,1),size(t,2),size(t,3)) :: svcp,z_full_ag, & |
289 |
|
|
k_m_save, k_t_save |
290 |
|
|
real, dimension(size(t,1),size(t,2),size(t,3)+1):: z_half_ag |
291 |
|
|
real, dimension(size(t,1),size(t,2)) :: z_surf |
292 |
|
|
integer :: i,j,k,nlev,nlat,nlon |
293 |
|
|
|
294 |
|
|
if(.not.init) call diffusivity_init(myThid) |
295 |
|
|
|
296 |
|
|
nlev = size(t,3) |
297 |
|
|
|
298 |
|
|
k_m_save = k_m |
299 |
|
|
k_t_save = k_t |
300 |
|
|
|
301 |
|
|
!compute height of surface |
302 |
|
|
if (present(kbot)) then |
303 |
|
|
nlat = size(t,2) |
304 |
|
|
nlon = size(t,1) |
305 |
|
|
do j=1,nlat |
306 |
|
|
do i=1,nlon |
307 |
|
|
z_surf(i,j) = z_half(i,j,kbot(i,j)+1) |
308 |
|
|
enddo |
309 |
|
|
enddo |
310 |
|
|
else |
311 |
|
|
z_surf(:,:) = z_half(:,:,nlev+1) |
312 |
|
|
end if |
313 |
|
|
|
314 |
|
|
!compute density profile, and heights relative to surface |
315 |
|
|
do k = 1, nlev |
316 |
|
|
|
317 |
|
|
z_full_ag(:,:,k) = z_full(:,:,k) - z_surf(:,:) |
318 |
|
|
z_half_ag(:,:,k) = z_half(:,:,k) - z_surf(:,:) |
319 |
|
|
|
320 |
|
|
if (do_virtual_non_mcm) then |
321 |
|
|
svcp(:,:,k) = t(:,:,k)*(1. + d608*q(:,:,k)) + gcp*(z_full_ag(:,:,k)) |
322 |
|
|
else |
323 |
|
|
svcp(:,:,k) = t(:,:,k) + gcp*(z_full_ag(:,:,k)) |
324 |
|
|
endif |
325 |
|
|
|
326 |
|
|
end do |
327 |
|
|
z_half_ag(:,:,nlev+1) = z_half(:,:,nlev+1) - z_surf(:,:) |
328 |
|
|
|
329 |
|
|
if(fixed_depth) then |
330 |
|
|
h = depth_0 |
331 |
|
|
else |
332 |
|
|
call pbl_depth(svcp,u,v,z_full_ag,u_star,b_star,h, myThid, kbot=kbot) |
333 |
|
|
end if |
334 |
|
|
|
335 |
|
|
if(pbl_mcm) then |
336 |
|
|
call diffusivity_pbl_mcm (u,v, t, p_full, p_half, & |
337 |
|
|
z_full_ag, z_half_ag, h, k_m, k_t, myThid) |
338 |
|
|
else |
339 |
|
|
call diffusivity_pbl (svcp, u, v, z_half_ag, h, u_star, b_star,& |
340 |
|
|
k_m, k_t, myThid, kbot=kbot) |
341 |
|
|
end if |
342 |
|
|
if(free_atm_diff) & |
343 |
|
|
call diffusivity_free (svcp, u, v, z_full_ag, z_half_ag, h, k_m, k_t, myThid) |
344 |
|
|
|
345 |
|
|
k_m = k_m + k_m_save |
346 |
|
|
k_t = k_t + k_t_save |
347 |
|
|
|
348 |
|
|
!NOTE THAT THIS LINE MUST FOLLOW DIFFUSIVITY_FREE SO THAT ENTRAINMENT |
349 |
|
|
!K's DO NOT GET OVERWRITTEN IN DIFFUSIVITY_FREE SUBROUTINE |
350 |
|
|
if(entr_ratio .gt. 0. .and. .not. fixed_depth) & |
351 |
|
|
call diffusivity_entr(svcp,z_full_ag,h,u_star,b_star,k_m,k_t, myThid) |
352 |
|
|
|
353 |
|
|
!set background diffusivities |
354 |
|
|
if(background_m.gt.0.0) k_m = max(k_m,background_m) |
355 |
|
|
if(background_t.gt.0.0) k_t = max(k_t,background_t) |
356 |
|
|
|
357 |
|
|
return |
358 |
|
|
end subroutine diffusivity |
359 |
|
|
|
360 |
|
|
!======================================================================= |
361 |
|
|
|
362 |
|
|
subroutine pbl_depth(t, u, v, z, u_star, b_star, h, myThid, kbot) |
363 |
|
|
|
364 |
|
|
real, intent(in) , dimension(:,:,:) :: t, u, v, z |
365 |
|
|
real, intent(in) , dimension(:,:) :: u_star,b_star |
366 |
|
|
real, intent(out), dimension(:,:) :: h |
367 |
|
|
integer, intent(in) :: myThid |
368 |
|
|
integer,intent(in) , optional, dimension(:,:) :: kbot |
369 |
|
|
|
370 |
|
|
real, dimension(size(t,1),size(t,2),size(t,3)) :: rich |
371 |
|
|
real, dimension(size(t,1),size(t,2)) :: ws,k_t_ref,& |
372 |
|
|
h_inner,tbot |
373 |
|
|
real :: rich1, rich2,& |
374 |
|
|
h1,h2,svp,t1,t2 |
375 |
|
|
integer, dimension(size(t,1),size(t,2)) :: ibot |
376 |
|
|
integer :: i,j,k,nlon,& |
377 |
|
|
nlat, nlev |
378 |
|
|
|
379 |
|
|
nlev = size(t,3) |
380 |
|
|
nlat = size(t,2) |
381 |
|
|
nlon = size(t,1) |
382 |
|
|
|
383 |
|
|
!assign ibot, compute tbot (virtual temperature at lowest level) |
384 |
|
|
if (present(kbot)) then |
385 |
|
|
ibot(:,:) = kbot |
386 |
|
|
do j = 1,nlat |
387 |
|
|
do i = 1,nlon |
388 |
|
|
tbot(i,j) = t(i,j,ibot(i,j)) |
389 |
|
|
enddo |
390 |
|
|
enddo |
391 |
|
|
else |
392 |
|
|
ibot(:,:) = nlev |
393 |
|
|
tbot(:,:) = t(:,:,nlev) |
394 |
|
|
end if |
395 |
|
|
|
396 |
|
|
!compute richardson number for use in pbl depth of neutral/stable side |
397 |
|
|
do k = 1,nlev |
398 |
|
|
rich(:,:,k) = z(:,:,k)*grav*(t(:,:,k)-tbot(:,:))/tbot(:,:)& |
399 |
|
|
/(u(:,:,k)*u(:,:,k) + v(:,:,k)*v(:,:,k) + small ) |
400 |
|
|
end do |
401 |
|
|
|
402 |
|
|
!compute ws to be used in evaluating parcel buoyancy |
403 |
|
|
!ws = u_star / phi(h_inner,u_star,b_star) . To find phi |
404 |
|
|
!a call to mo_diff is made. |
405 |
|
|
|
406 |
|
|
h_inner(:,:)=frac_inner*znom |
407 |
|
|
call mo_diff(h_inner, u_star, b_star, ws, k_t_ref, myThid ) |
408 |
|
|
ws = max(small,ws/vonkarm/h_inner) |
409 |
|
|
|
410 |
|
|
do j = 1, nlat |
411 |
|
|
do i = 1, nlon |
412 |
|
|
|
413 |
|
|
!do neutral or stable case |
414 |
|
|
if (b_star(i,j).le.0.) then |
415 |
|
|
|
416 |
|
|
h1 = z(i,j,ibot(i,j)) |
417 |
|
|
h(i,j) = h1 |
418 |
|
|
rich1 = rich(i,j,ibot(i,j)) |
419 |
|
|
do k = ibot(i,j)-1, 1, -1 |
420 |
|
|
rich2 = rich(i,j,k) |
421 |
|
|
h2 = z(i,j,k) |
422 |
|
|
if(rich2.gt.rich_crit_pbl) then |
423 |
|
|
h(i,j) = h2 + (h1 - h2)*(rich2 - rich_crit_pbl)& |
424 |
|
|
/(rich2 - rich1 ) |
425 |
|
|
go to 10 |
426 |
|
|
endif |
427 |
|
|
rich1 = rich2 |
428 |
|
|
h1 = h2 |
429 |
|
|
enddo |
430 |
|
|
|
431 |
|
|
!do unstable case |
432 |
|
|
else |
433 |
|
|
|
434 |
|
|
svp = tbot(i,j)*(1.+ & |
435 |
|
|
(parcel_buoy*u_star(i,j)*b_star(i,j)/grav/ws(i,j)) ) |
436 |
|
|
h1 = z(i,j,ibot(i,j)) |
437 |
|
|
h(i,j) = h1 |
438 |
|
|
t1 = tbot(i,j) |
439 |
|
|
do k = ibot(i,j)-1 , 1, -1 |
440 |
|
|
h2 = z(i,j,k) |
441 |
|
|
t2 = t(i,j,k) |
442 |
|
|
if (t2.gt.svp) then |
443 |
|
|
h(i,j) = h2 + (h1 - h2)*(t2 - svp)/(t2 - t1 ) |
444 |
|
|
go to 10 |
445 |
|
|
end if |
446 |
|
|
h1 = h2 |
447 |
|
|
t1 = t2 |
448 |
|
|
enddo |
449 |
|
|
|
450 |
|
|
end if |
451 |
|
|
10 continue |
452 |
|
|
enddo |
453 |
|
|
enddo |
454 |
|
|
|
455 |
|
|
return |
456 |
|
|
end subroutine pbl_depth |
457 |
|
|
|
458 |
|
|
!======================================================================= |
459 |
|
|
|
460 |
|
|
subroutine diffusivity_pbl(t, u, v, z_half, h, u_star, b_star, & |
461 |
|
|
k_m, k_t, myThid, kbot) |
462 |
|
|
|
463 |
|
|
real, intent(in) , dimension(:,:,:) :: t, u, v, z_half |
464 |
|
|
real, intent(in) , dimension(:,:) :: h, u_star, b_star |
465 |
|
|
real, intent(inout) , dimension(:,:,:) :: k_m, k_t |
466 |
|
|
integer, intent (in) :: myThid |
467 |
|
|
integer, intent(in) , optional, dimension(:,:) :: kbot |
468 |
|
|
|
469 |
|
|
real, dimension(size(t,1),size(t,2)) :: h_inner, k_m_ref,& |
470 |
|
|
k_t_ref, factor |
471 |
|
|
real, dimension(size(t,1),size(t,2),size(t,3)+1) :: zm |
472 |
|
|
real :: h_inner_max |
473 |
|
|
integer :: i,j, k, kk, nlev |
474 |
|
|
|
475 |
|
|
nlev = size(t,3) |
476 |
|
|
|
477 |
|
|
!assign z_half to zm, and set to zero any values of zm < 0. |
478 |
|
|
!the setting to zero is necessary so that when using eta model |
479 |
|
|
!below ground half levels will have zero k_m and k_t |
480 |
|
|
zm = z_half |
481 |
|
|
if (present(kbot)) then |
482 |
|
|
where(zm < 0.) |
483 |
|
|
zm = 0. |
484 |
|
|
end where |
485 |
|
|
end if |
486 |
|
|
|
487 |
|
|
h_inner = frac_inner*h |
488 |
|
|
h_inner_max = maxval(h_inner) |
489 |
|
|
|
490 |
|
|
kk = nlev |
491 |
|
|
do k = 2, nlev |
492 |
|
|
if( minval(zm(:,:,k)) < h_inner_max) then |
493 |
|
|
kk = k |
494 |
|
|
exit |
495 |
|
|
end if |
496 |
|
|
end do |
497 |
|
|
|
498 |
|
|
k_m = 0.0 |
499 |
|
|
k_t = 0.0 |
500 |
|
|
|
501 |
|
|
call mo_diff(h_inner , u_star, b_star, k_m_ref , k_t_ref, myThid ) |
502 |
|
|
call mo_diff(zm(:,:,kk:nlev), u_star, b_star, k_m(:,:,kk:nlev), k_t(:,:,kk:nlev), & |
503 |
|
|
myThid ) |
504 |
|
|
|
505 |
|
|
do k = 2, nlev |
506 |
|
|
where(zm(:,:,k) >= h_inner .and. zm(:,:,k) < h) |
507 |
|
|
factor = (zm(:,:,k)/h_inner)* & |
508 |
|
|
(1.0 - (zm(:,:,k) - h_inner)/(h - h_inner))**2 |
509 |
|
|
k_m(:,:,k) = k_m_ref*factor |
510 |
|
|
k_t(:,:,k) = k_t_ref*factor |
511 |
|
|
end where |
512 |
|
|
end do |
513 |
|
|
|
514 |
|
|
return |
515 |
|
|
end subroutine diffusivity_pbl |
516 |
|
|
|
517 |
|
|
!======================================================================= |
518 |
|
|
|
519 |
|
|
subroutine diffusivity_pbl_mcm(u, v, t, p_full, p_half, z_full, z_half, & |
520 |
|
|
h, k_m, k_t, myThid) |
521 |
|
|
|
522 |
|
|
real, intent(in) , dimension(:,:,:) :: u, v, t, z_full, z_half |
523 |
|
|
real, intent(in) , dimension(:,:,:) :: p_full, p_half |
524 |
|
|
real, intent(in) , dimension(:,:) :: h |
525 |
|
|
real, intent(inout) , dimension(:,:,:) :: k_m, k_t |
526 |
|
|
integer, intent(in) :: myThid |
527 |
|
|
|
528 |
|
|
integer :: k, nlev |
529 |
|
|
real, dimension(size(z_full,1),size(z_full,2)) :: elmix, htcrit |
530 |
|
|
real, dimension(size(z_full,1),size(z_full,2)) :: delta_u, delta_v, delta_z |
531 |
|
|
|
532 |
|
|
real :: htcrit_ss |
533 |
|
|
real :: h_ss |
534 |
|
|
real, dimension(size(z_full,1),size(z_full,2)) :: sig_half, z_half_ss, elmix_ss |
535 |
|
|
|
536 |
|
|
! htcrit_ss = height at which mixing length is a maximum (75m) |
537 |
|
|
|
538 |
|
|
! h_ss = height at which mixing length vanishes (4900m) |
539 |
|
|
! elmix_ss = mixing length |
540 |
|
|
|
541 |
|
|
! Define some constants: |
542 |
|
|
! salaps = standard atmospheric lapse rate (K/m) |
543 |
|
|
! tsfc = idealized global mean surface temperature (15C) |
544 |
|
|
real :: tsfc = 288.16 |
545 |
|
|
real :: salaps = -6.5e-3 |
546 |
|
|
|
547 |
|
|
nlev = size(z_full,3) |
548 |
|
|
|
549 |
|
|
k_m = 0. |
550 |
|
|
|
551 |
|
|
h_ss = depth_0 |
552 |
|
|
htcrit_ss = frac_inner*h_ss |
553 |
|
|
|
554 |
|
|
do k = 2, nlev |
555 |
|
|
|
556 |
|
|
! TK mods 8/13/01: (code derived from SS) |
557 |
|
|
! Compute the height of each half level assuming a constant |
558 |
|
|
! standard lapse rate using the above procedure. |
559 |
|
|
! WARNING: These should be used with caution. They will |
560 |
|
|
! have large errors above the tropopause. |
561 |
|
|
|
562 |
|
|
! In order to determine the height, the layer mean temperature |
563 |
|
|
! from the surface to that level is required. A surface |
564 |
|
|
! temperature of 15 deg Celsius and a standard lapse rate of |
565 |
|
|
! -6.5 deg/km will be used to estimate an average temperature |
566 |
|
|
! profile. |
567 |
|
|
|
568 |
|
|
sig_half = p_half(:,:,k)/p_half(:,:,nlev+1) |
569 |
|
|
z_half_ss = -rdgas * .5*(tsfc+tsfc*(sig_half**(-rdgas*salaps/grav))) * alog(sig_half)/grav |
570 |
|
|
|
571 |
|
|
!compute mixing length as in SS (no geographical variation) |
572 |
|
|
elmix_ss = 0. |
573 |
|
|
|
574 |
|
|
where (z_half_ss < htcrit_ss .and. z_half_ss > 0.) |
575 |
|
|
elmix_ss = vonkarm*z_half_ss |
576 |
|
|
endwhere |
577 |
|
|
where (z_half_ss >= htcrit_ss .and. z_half_ss < h_ss) |
578 |
|
|
elmix_ss = vonkarm*htcrit_ss*(h_ss-z_half_ss)/(h_ss-htcrit_ss) |
579 |
|
|
endwhere |
580 |
|
|
|
581 |
|
|
delta_z = rdgas*0.5*(t(:,:,k)+t(:,:,k-1))*(p_full(:,:,k)-p_full(:,:,k-1))/& |
582 |
|
|
(grav*p_half(:,:,k)) |
583 |
|
|
delta_u = u(:,:,k-1) - u(:,:,k) |
584 |
|
|
delta_v = v(:,:,k-1) - v(:,:,k) |
585 |
|
|
|
586 |
|
|
k_m(:,:,k) = elmix_ss * elmix_ss *& |
587 |
|
|
sqrt(delta_u*delta_u + delta_v*delta_v)/delta_z |
588 |
|
|
|
589 |
|
|
end do |
590 |
|
|
|
591 |
|
|
k_t = k_m |
592 |
|
|
|
593 |
|
|
return |
594 |
|
|
end subroutine diffusivity_pbl_mcm |
595 |
|
|
|
596 |
|
|
!======================================================================= |
597 |
|
|
|
598 |
|
|
subroutine diffusivity_free(t, u, v, z, zz, h, k_m, k_t, myThid) |
599 |
|
|
|
600 |
|
|
real, intent(in) , dimension(:,:,:) :: t, u, v, z, zz |
601 |
|
|
real, intent(in) , dimension(:,:) :: h |
602 |
|
|
real, intent(inout) , dimension(:,:,:) :: k_m, k_t |
603 |
|
|
integer, intent(in) :: myThid |
604 |
|
|
|
605 |
|
|
real, dimension(size(t,1),size(t,2)) :: dz, b, speed2, rich, fri, & |
606 |
|
|
alpz, fri2 |
607 |
|
|
integer :: k |
608 |
|
|
|
609 |
|
|
do k = 2, size(t,3) |
610 |
|
|
|
611 |
|
|
!---------------------------------------------------------------------- |
612 |
|
|
! define the richardson number. set it to zero if it is negative. save |
613 |
|
|
! a copy of it for later use (rich2). |
614 |
|
|
!---------------------------------------------------------------------- |
615 |
|
|
dz = z(:,:,k-1) - z(:,:,k) |
616 |
|
|
b = grav*(t(:,:,k-1)-t(:,:,k))/t(:,:,k) |
617 |
|
|
speed2 = (u(:,:,k-1) - u(:,:,k))**2 + (v(:,:,k-1) - v(:,:,k))**2 |
618 |
|
|
rich= b*dz/(speed2+small) |
619 |
|
|
rich = max(rich, 0.0) |
620 |
|
|
|
621 |
|
|
if (free_atm_skyhi_diff) then |
622 |
|
|
!--------------------------------------------------------------------- |
623 |
|
|
! limit the standard richardson number to between 0 and the critical |
624 |
|
|
! value (rich2). compute the richardson number factor needed in the |
625 |
|
|
! eddy mixing coefficient using this standard richardson number. |
626 |
|
|
!--------------------------------------------------------------------- |
627 |
|
|
where (rich(:,:) >= rich_crit_diff) |
628 |
|
|
fri2(:,:) = 0.0 |
629 |
|
|
elsewhere |
630 |
|
|
fri2(:,:) = (1.0 - rich/rich_crit_diff)**2 |
631 |
|
|
endwhere |
632 |
|
|
endif |
633 |
|
|
|
634 |
|
|
!--------------------------------------------------------------------- |
635 |
|
|
! if ampns is activated, compute the delta z factor. define rich |
636 |
|
|
! including this factor. |
637 |
|
|
!--------------------------------------------------------------------- |
638 |
|
|
if (ampns) then |
639 |
|
|
alpz(:,:) = MIN ( (1. + 1.e-04*(dz(:,:)**1.5)), ampns_max) |
640 |
|
|
rich(:,:) = rich(:,:) / alpz(:,:) |
641 |
|
|
endif |
642 |
|
|
|
643 |
|
|
!--------------------------------------------------------------------- |
644 |
|
|
! compute the richardson number factor to be used in the eddy |
645 |
|
|
! mixing coefficient. if ampns is on, this value includes it; other- |
646 |
|
|
! wise it does not. |
647 |
|
|
!--------------------------------------------------------------------- |
648 |
|
|
fri(:,:) = (1.0 - rich/rich_crit_diff)**2 |
649 |
|
|
|
650 |
|
|
!--------------------------------------------------------------------- |
651 |
|
|
! compute the eddy mixing coefficients in the free atmosphere ( zz |
652 |
|
|
! > h). in the non-ampns case, values are obtained only when the |
653 |
|
|
! standard richardson number is sub-critical; in the ampns case values |
654 |
|
|
! are obtained only when the richardson number computed with the |
655 |
|
|
! ampns factor is sub critical. when the ampns factor is activated, |
656 |
|
|
! it is also included in the mixing coefficient. the value of mixing |
657 |
|
|
! for temperature, etc. is reduced dependent on the ri stability |
658 |
|
|
! factor calculated without the ampns factor. |
659 |
|
|
!--------------------------------------------------------------------- |
660 |
|
|
if (free_atm_skyhi_diff) then |
661 |
|
|
|
662 |
|
|
!--------------------------------------------------------------------- |
663 |
|
|
! this is the skyhi-like formulation -- possible ampns factor, ratio |
664 |
|
|
! of k_m to k_t defined based on computed stability factor. |
665 |
|
|
!--------------------------------------------------------------------- |
666 |
|
|
if (ampns) then |
667 |
|
|
where (rich < rich_crit_diff .and. zz(:,:,k) > h) |
668 |
|
|
k_m(:,:,k) = mix_len*mix_len*sqrt(speed2)*fri(:,:)* & |
669 |
|
|
( 1. + 1.e-04*(dz(:,:)**1.5))/dz |
670 |
|
|
k_t(:,:,k) = k_m(:,:,k)* (0.1 + 0.9*fri2(:,:)) |
671 |
|
|
end where |
672 |
|
|
else |
673 |
|
|
where (rich < rich_crit_diff .and. zz(:,:,k) > h) |
674 |
|
|
k_m(:,:,k) = mix_len*mix_len*sqrt(speed2)*fri(:,:)/dz |
675 |
|
|
k_t(:,:,k) = k_m(:,:,k)* (0.1 + 0.9*fri2(:,:)) |
676 |
|
|
end where |
677 |
|
|
endif |
678 |
|
|
else |
679 |
|
|
|
680 |
|
|
!--------------------------------------------------------------------- |
681 |
|
|
! this is the non-skyhi-like formulation -- no ampns factor, ratio |
682 |
|
|
! of k_m to k_t defined by rich_prandtl. |
683 |
|
|
!--------------------------------------------------------------------- |
684 |
|
|
where (rich < rich_crit_diff .and. zz(:,:,k) > h) |
685 |
|
|
k_t(:,:,k) = mix_len*mix_len*sqrt(speed2)*fri(:,:)/dz |
686 |
|
|
k_m(:,:,k) = k_t(:,:,k)*rich_prandtl |
687 |
|
|
end where |
688 |
|
|
end if |
689 |
|
|
end do |
690 |
|
|
|
691 |
|
|
end subroutine diffusivity_free |
692 |
|
|
|
693 |
|
|
!======================================================================= |
694 |
|
|
|
695 |
|
|
subroutine molecular_diff ( temp, press, k_m, k_t, myThid) |
696 |
|
|
|
697 |
|
|
real, intent(in), dimension (:,:,:) :: temp, press |
698 |
|
|
real, intent(inout), dimension (:,:,:) :: k_m, k_t |
699 |
|
|
integer, intent(in) :: myThid |
700 |
|
|
|
701 |
|
|
real, dimension (size(temp,1), size(temp,2)) :: temp_half, & |
702 |
|
|
rho_half, rbop2d |
703 |
|
|
integer :: k |
704 |
|
|
|
705 |
|
|
!--------------------------------------------------------------------- |
706 |
|
|
|
707 |
|
|
do k=2,size(temp,3) |
708 |
|
|
temp_half(:,:) = 0.5*(temp(:,:,k) + temp(:,:,k-1)) |
709 |
|
|
rho_half(:,:) = press(:,:,k)/(rdgas*temp_half(:,:) ) |
710 |
|
|
rbop2d(:,:) = beta*temp_half(:,:)*sqrt(temp_half(:,:))/ & |
711 |
|
|
(rho_half(:,:)*(temp_half(:,:)+rbop1)) |
712 |
|
|
k_m(:,:,k) = rbop2d(:,:) |
713 |
|
|
k_t(:,:,k) = rbop2d(:,:)*rbop2 |
714 |
|
|
end do |
715 |
|
|
|
716 |
|
|
k_m(:,:,1) = 0.0 |
717 |
|
|
k_t(:,:,1) = 0.0 |
718 |
|
|
|
719 |
|
|
end subroutine molecular_diff |
720 |
|
|
|
721 |
|
|
!======================================================================= |
722 |
|
|
|
723 |
|
|
subroutine diffusivity_entr(t, z, h, u_star, b_star, k_m, k_t, myThid) |
724 |
|
|
|
725 |
|
|
real, intent(in) , dimension(:,:,:) :: t, z |
726 |
|
|
real, intent(in) , dimension(:,:) :: h, u_star, b_star |
727 |
|
|
real, intent(inout) , dimension(:,:,:) :: k_m, k_t |
728 |
|
|
integer, intent(in) :: myThid |
729 |
|
|
|
730 |
|
|
integer :: k, nlev |
731 |
|
|
|
732 |
|
|
nlev=size(t,3) |
733 |
|
|
|
734 |
|
|
do k = 2,nlev |
735 |
|
|
where (b_star .gt. 0. .and. z(:,:,k-1) .gt. h .and. & |
736 |
|
|
z(:,:,k) .le. h) |
737 |
|
|
k_t(:,:,k) = (z(:,:,k-1)-z(:,:,k))*entr_ratio*t(:,:,k)* & |
738 |
|
|
u_star*b_star/grav/max(small,t(:,:,k-1)-t(:,:,k)) |
739 |
|
|
k_m(:,:,k) = k_t(:,:,k) |
740 |
|
|
end where |
741 |
|
|
enddo |
742 |
|
|
end subroutine diffusivity_entr |
743 |
|
|
|
744 |
|
|
!======================================================================= |
745 |
|
|
|
746 |
|
|
end module diffusivity_mod |
747 |
|
|
|