24 implicit none ;
private
26 #include <MOM_memory.h>
28 public mixedlayer_restrat
29 public mixedlayer_restrat_init
30 public mixedlayer_restrat_register_restarts
39 real :: ml_restrat_coef
43 real :: ml_restrat_coef2
48 logical :: mle_use_pbl_mld
51 real :: mle_mld_decay_time
52 real :: mle_mld_decay_time2
53 real :: mle_density_diff
57 real :: mle_mld_stretch
59 logical :: mle_use_mld_ave_bug
60 logical :: debug = .false.
64 real,
dimension(:,:),
pointer :: &
65 mld_filtered => null(), &
66 mld_filtered_slow => null()
70 integer :: id_urestrat_time = -1
71 integer :: id_vrestrat_time = -1
72 integer :: id_uhml = -1
73 integer :: id_vhml = -1
74 integer :: id_mld = -1
75 integer :: id_rml = -1
76 integer :: id_udml = -1
77 integer :: id_vdml = -1
78 integer :: id_uml = -1
79 integer :: id_vml = -1
84 character(len=40) :: mdl =
"MOM_mixed_layer_restrat"
91 subroutine mixedlayer_restrat(h, uhtr, vhtr, tv, forces, dt, MLD, VarMix, G, GV, US, CS)
95 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(inout) :: h
96 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)),
intent(inout) :: uhtr
98 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)),
intent(inout) :: vhtr
102 real,
intent(in) :: dt
103 real,
dimension(:,:),
pointer :: mld
108 if (.not.
associated(cs))
call mom_error(fatal,
"MOM_mixedlayer_restrat: "// &
109 "Module must be initialized before it is used.")
112 call mixedlayer_restrat_bml(h, uhtr, vhtr, tv, forces, dt, g, gv, us, cs)
114 call mixedlayer_restrat_general(h, uhtr, vhtr, tv, forces, dt, mld, varmix, g, gv, us, cs)
117 end subroutine mixedlayer_restrat
120 subroutine mixedlayer_restrat_general(h, uhtr, vhtr, tv, forces, dt, MLD_in, VarMix, G, GV, US, CS)
125 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(inout) :: h
126 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)),
intent(inout) :: uhtr
128 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)),
intent(inout) :: vhtr
132 real,
intent(in) :: dt
133 real,
dimension(:,:),
pointer :: MLD_in
138 real :: uhml(SZIB_(G),SZJ_(G),SZK_(G))
139 real :: vhml(SZI_(G),SZJB_(G),SZK_(G))
140 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)) :: &
143 real,
dimension(SZI_(G),SZJ_(G)) :: &
151 real :: rho_ml(SZI_(G))
162 real :: Ihtot,Ihtot_slow
168 real :: uDml(SZIB_(G))
169 real :: vDml(SZI_(G))
170 real :: uDml_slow(SZIB_(G))
171 real :: vDml_slow(SZI_(G))
172 real :: utimescale_diag(SZIB_(G),SZJ_(G))
173 real :: vtimescale_diag(SZI_(G),SZJB_(G))
175 real :: uDml_diag(SZIB_(G),SZJ_(G)), vDml_diag(SZI_(G),SZJB_(G))
176 integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
177 real,
dimension(SZI_(G)) :: rhoSurf, deltaRhoAtKm1, deltaRhoAtK
178 real,
dimension(SZI_(G)) :: dK, dKm1
179 real,
dimension(SZI_(G)) :: pRef_MLD
180 real,
dimension(SZI_(G)) :: rhoAtK, rho1, d1, pRef_N2
181 real :: aFac, bFac, ddRho
182 real :: hAtVel, zpa, zpb, dh, res_scaling_fac, I_l_f
183 logical :: proper_averaging, line_is_empty, keep_going, res_upscale
185 real :: PSI, PSI1, z, BOTTOP, XP, DD
188 psi1(z) = max(0., (1. - (2.*z+1.)**2 ) * (1. + (5./21.)*(2.*z+1.)**2) )
189 bottop(z) = 0.5*(1.-sign(1.,z+0.5))
190 xp(z) = max(0., min(1., (-z-0.5)*2./(1.+2.*cs%MLE_tail_dh) ) )
191 dd(z) = (1.-3.*(xp(z)**2)+2.*(xp(z)**3))**(1.+2.*cs%MLE_tail_dh)
192 psi(z) = max( psi1(z), dd(z)*bottop(z) )
194 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
195 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
197 if (.not.
associated(tv%eqn_of_state))
call mom_error(fatal,
"MOM_mixedlayer_restrat: "// &
198 "An equation of state must be used with this module.")
199 if (.not.
associated(varmix) .and. cs%front_length>0.)
call mom_error(fatal,
"MOM_mixedlayer_restrat: "// &
200 "The resolution argument, Rd/dx, was not associated.")
202 if (cs%MLE_density_diff > 0.)
then
206 dk(:) = 0.5 * h(:,j,1)
207 call calculate_density(tv%T(:,j,1), tv%S(:,j,1), pref_mld, rhosurf, is-1, ie-is+3, tv%eqn_of_state)
212 dk(:) = dk(:) + 0.5 * ( h(:,j,k) + h(:,j,k-1) )
214 deltarhoatkm1(:) = deltarhoatk(:)
215 call calculate_density(tv%T(:,j,k), tv%S(:,j,k), pref_mld, deltarhoatk, is-1, ie-is+3, tv%eqn_of_state)
216 deltarhoatk(:) = deltarhoatk(:) - rhosurf(:)
218 ddrho = deltarhoatk(i) - deltarhoatkm1(i)
219 if ((mld_fast(i,j)==0.) .and. (ddrho>0.) .and. &
220 (deltarhoatkm1(i)<cs%MLE_density_diff) .and. (deltarhoatk(i)>=cs%MLE_density_diff))
then
221 afac = ( cs%MLE_density_diff - deltarhoatkm1(i) ) / ddrho
222 mld_fast(i,j) = dk(i) * afac + dkm1(i) * (1. - afac)
227 mld_fast(i,j) = cs%MLE_MLD_stretch * mld_fast(i,j)
228 if ((mld_fast(i,j)==0.) .and. (deltarhoatk(i)<cs%MLE_density_diff)) &
229 mld_fast(i,j) = dk(i)
232 elseif (cs%MLE_use_PBL_MLD)
then
233 if (.not.
associated(mld_in))
call mom_error(fatal,
"MOM_mixedlayer_restrat: "// &
234 "Argument MLD_in was not associated!")
235 do j = js-1, je+1 ;
do i = is-1, ie+1
236 mld_fast(i,j) = (cs%MLE_MLD_stretch * gv%m_to_H) * mld_in(i,j)
239 call mom_error(fatal,
"MOM_mixedlayer_restrat: "// &
240 "No MLD to use for MLE parameterization.")
244 if (cs%MLE_MLD_decay_time>0.)
then
246 call hchksum(cs%MLD_filtered,
'mixed_layer_restrat: MLD_filtered',g%HI,haloshift=1,scale=gv%H_to_m)
247 call hchksum(mld_in,
'mixed_layer_restrat: MLD in',g%HI,haloshift=1)
249 afac = cs%MLE_MLD_decay_time / ( dt + cs%MLE_MLD_decay_time )
250 bfac = dt / ( dt + cs%MLE_MLD_decay_time )
251 do j = js-1, je+1 ;
do i = is-1, ie+1
255 cs%MLD_filtered(i,j) = max( mld_fast(i,j), bfac*mld_fast(i,j) + afac*cs%MLD_filtered(i,j) )
256 mld_fast(i,j) = cs%MLD_filtered(i,j)
261 if (cs%MLE_MLD_decay_time2>0.)
then
263 call hchksum(cs%MLD_filtered_slow,
'mixed_layer_restrat: MLD_filtered_slow',g%HI,haloshift=1,scale=gv%H_to_m)
264 call hchksum(mld_fast,
'mixed_layer_restrat: MLD fast',g%HI,haloshift=1,scale=gv%H_to_m)
266 afac = cs%MLE_MLD_decay_time2 / ( dt + cs%MLE_MLD_decay_time2 )
267 bfac = dt / ( dt + cs%MLE_MLD_decay_time2 )
268 do j = js-1, je+1 ;
do i = is-1, ie+1
272 cs%MLD_filtered_slow(i,j) = max( mld_fast(i,j), bfac*mld_fast(i,j) + afac*cs%MLD_filtered_slow(i,j) )
273 mld_slow(i,j) = cs%MLD_filtered_slow(i,j)
276 do j = js-1, je+1 ;
do i = is-1, ie+1
277 mld_slow(i,j) = mld_fast(i,j)
281 udml(:) = 0.0 ; vdml(:) = 0.0
282 udml_slow(:) = 0.0 ; vdml_slow(:) = 0.0
284 g_rho0 = gv%g_Earth*us%L_to_m**2*us%s_to_T**2 / gv%Rho0
285 h_neglect = gv%H_subroundoff
286 dz_neglect = gv%H_subroundoff*gv%H_to_Z
287 proper_averaging = .not. cs%MLE_use_MLD_ave_bug
288 if (cs%front_length>0.)
then
290 i_l_f = 1./cs%front_length
292 res_upscale = .false.
309 htot_fast(i,j) = 0.0 ; rml_av_fast(i,j) = 0.0
310 htot_slow(i,j) = 0.0 ; rml_av_slow(i,j) = 0.0
315 h_avail(i,j,k) = max(i4dt*g%areaT(i,j)*(h(i,j,k)-gv%Angstrom_H),0.0)
318 call calculate_density(tv%T(:,j,k),tv%S(:,j,k),p0,rho_ml(:),is-1,ie-is+3,tv%eqn_of_state)
319 line_is_empty = .true.
321 if (htot_fast(i,j) < mld_fast(i,j))
then
323 if (proper_averaging) dh = min( h(i,j,k), mld_fast(i,j)-htot_fast(i,j) )
324 rml_av_fast(i,j) = rml_av_fast(i,j) + dh*rho_ml(i)
325 htot_fast(i,j) = htot_fast(i,j) + dh
326 line_is_empty = .false.
328 if (htot_slow(i,j) < mld_slow(i,j))
then
329 dh = min( h(i,j,k), mld_slow(i,j)-htot_slow(i,j) )
330 rml_av_slow(i,j) = rml_av_slow(i,j) + dh*rho_ml(i)
331 htot_slow(i,j) = htot_slow(i,j) + dh
332 line_is_empty = .false.
335 if (line_is_empty) keep_going=.false.
340 rml_av_fast(i,j) = -(g_rho0*rml_av_fast(i,j)) / (htot_fast(i,j) + h_neglect)
341 rml_av_slow(i,j) = -(g_rho0*rml_av_slow(i,j)) / (htot_slow(i,j) + h_neglect)
346 call hchksum(h,
'mixed_layer_restrat: h',g%HI,haloshift=1,scale=gv%H_to_m)
347 call hchksum(forces%ustar,
'mixed_layer_restrat: u*',g%HI,haloshift=1,scale=us%Z_to_m*us%s_to_T)
348 call hchksum(mld_fast,
'mixed_layer_restrat: MLD',g%HI,haloshift=1,scale=gv%H_to_m)
349 call hchksum(rml_av_fast,
'mixed_layer_restrat: rml',g%HI,haloshift=1, scale=us%m_to_Z)
358 do j=js,je ;
do i=is-1,ie
359 u_star = us%s_to_T*0.5*(forces%ustar(i,j) + forces%ustar(i+1,j))
360 absf = 0.5*us%s_to_T*(abs(g%CoriolisBu(i,j-1)) + abs(g%CoriolisBu(i,j)))
362 if (res_upscale) res_scaling_fac = &
363 ( sqrt( 0.5 * ( g%dxCu(i,j)**2 + g%dyCu(i,j)**2 ) ) * i_l_f ) &
364 * min( 1., 0.5*( varmix%Rd_dx_h(i,j) + varmix%Rd_dx_h(i+1,j) ) )
369 h_vel = 0.5*((htot_fast(i,j) + htot_fast(i+1,j)) + h_neglect) * gv%H_to_Z
370 mom_mixrate = (0.41*9.8696)*u_star**2 / &
371 (absf*h_vel**2 + 4.0*(h_vel+dz_neglect)*u_star)
372 timescale = 0.0625 * (absf + 2.0*mom_mixrate) / (absf**2 + mom_mixrate**2)
373 timescale = timescale * cs%ml_restrat_coef
374 if (res_upscale) timescale = timescale * res_scaling_fac
375 udml(i) = timescale * g%mask2dCu(i,j)*g%dyCu(i,j)* &
376 g%IdxCu(i,j)*(rml_av_fast(i+1,j)-rml_av_fast(i,j)) * (h_vel**2 * gv%Z_to_H)
378 h_vel = 0.5*((htot_slow(i,j) + htot_slow(i+1,j)) + h_neglect) * gv%H_to_Z
379 mom_mixrate = (0.41*9.8696)*u_star**2 / &
380 (absf*h_vel**2 + 4.0*(h_vel+dz_neglect)*u_star)
381 timescale = 0.0625 * (absf + 2.0*mom_mixrate) / (absf**2 + mom_mixrate**2)
382 timescale = timescale * cs%ml_restrat_coef2
383 if (res_upscale) timescale = timescale * res_scaling_fac
384 udml_slow(i) = timescale * g%mask2dCu(i,j)*g%dyCu(i,j)* &
385 g%IdxCu(i,j)*(rml_av_slow(i+1,j)-rml_av_slow(i,j)) * (h_vel**2 * gv%Z_to_H)
387 if (udml(i) + udml_slow(i) == 0.)
then
388 do k=1,nz ; uhml(i,j,k) = 0.0 ;
enddo
390 ihtot = 2.0 / ((htot_fast(i,j) + htot_fast(i+1,j)) + h_neglect)
391 ihtot_slow = 2.0 / ((htot_slow(i,j) + htot_slow(i+1,j)) + h_neglect)
392 zpa = 0.0 ; zpb = 0.0
396 hatvel = 0.5*(h(i,j,k) + h(i+1,j,k))
398 zpa = zpa - (hatvel * ihtot)
399 a(k) = a(k) - psi(zpa)
401 if (a(k)*udml(i) > 0.0)
then
402 if (a(k)*udml(i) > h_avail(i,j,k)) udml(i) = h_avail(i,j,k) / a(k)
403 elseif (a(k)*udml(i) < 0.0)
then
404 if (-a(k)*udml(i) > h_avail(i+1,j,k)) udml(i) = -h_avail(i+1,j,k) / a(k)
409 hatvel = 0.5*(h(i,j,k) + h(i+1,j,k))
411 zpb = zpb - (hatvel * ihtot_slow)
412 b(k) = b(k) - psi(zpb)
414 if (b(k)*udml_slow(i) > 0.0)
then
415 if (b(k)*udml_slow(i) > h_avail(i,j,k) - a(k)*udml(i)) &
416 udml_slow(i) = max( 0., h_avail(i,j,k) - a(k)*udml(i) ) / b(k)
417 elseif (b(k)*udml_slow(i) < 0.0)
then
418 if (-b(k)*udml_slow(i) > h_avail(i+1,j,k) + a(k)*udml(i)) &
419 udml_slow(i) = -max( 0., h_avail(i+1,j,k) + a(k)*udml(i) ) / b(k)
423 uhml(i,j,k) = a(k)*udml(i) + b(k)*udml_slow(i)
424 uhtr(i,j,k) = uhtr(i,j,k) + uhml(i,j,k)*dt
428 utimescale_diag(i,j) = timescale
429 udml_diag(i,j) = udml(i)
434 do j=js-1,je ;
do i=is,ie
435 u_star = us%s_to_T*0.5*(forces%ustar(i,j) + forces%ustar(i,j+1))
436 absf = 0.5*us%s_to_T*(abs(g%CoriolisBu(i-1,j)) + abs(g%CoriolisBu(i,j)))
438 if (res_upscale) res_scaling_fac = &
439 ( sqrt( 0.5 * ( g%dxCv(i,j)**2 + g%dyCv(i,j)**2 ) ) * i_l_f ) &
440 * min( 1., 0.5*( varmix%Rd_dx_h(i,j) + varmix%Rd_dx_h(i,j+1) ) )
445 h_vel = 0.5*((htot_fast(i,j) + htot_fast(i,j+1)) + h_neglect) * gv%H_to_Z
446 mom_mixrate = (0.41*9.8696)*u_star**2 / &
447 (absf*h_vel**2 + 4.0*(h_vel+dz_neglect)*u_star)
448 timescale = 0.0625 * (absf + 2.0*mom_mixrate) / (absf**2 + mom_mixrate**2)
449 timescale = timescale * cs%ml_restrat_coef
450 if (res_upscale) timescale = timescale * res_scaling_fac
451 vdml(i) = timescale * g%mask2dCv(i,j)*g%dxCv(i,j)* &
452 g%IdyCv(i,j)*(rml_av_fast(i,j+1)-rml_av_fast(i,j)) * (h_vel**2 * gv%Z_to_H)
454 h_vel = 0.5*((htot_slow(i,j) + htot_slow(i,j+1)) + h_neglect) * gv%H_to_Z
455 mom_mixrate = (0.41*9.8696)*u_star**2 / &
456 (absf*h_vel**2 + 4.0*(h_vel+dz_neglect)*u_star)
457 timescale = 0.0625 * (absf + 2.0*mom_mixrate) / (absf**2 + mom_mixrate**2)
458 timescale = timescale * cs%ml_restrat_coef2
459 if (res_upscale) timescale = timescale * res_scaling_fac
460 vdml_slow(i) = timescale * g%mask2dCv(i,j)*g%dxCv(i,j)* &
461 g%IdyCv(i,j)*(rml_av_slow(i,j+1)-rml_av_slow(i,j)) * (h_vel**2 * gv%Z_to_H)
463 if (vdml(i) + vdml_slow(i) == 0.)
then
464 do k=1,nz ; vhml(i,j,k) = 0.0 ;
enddo
466 ihtot = 2.0 / ((htot_fast(i,j) + htot_fast(i,j+1)) + h_neglect)
467 ihtot_slow = 2.0 / ((htot_slow(i,j) + htot_slow(i,j+1)) + h_neglect)
468 zpa = 0.0 ; zpb = 0.0
472 hatvel = 0.5*(h(i,j,k) + h(i,j+1,k))
474 zpa = zpa - (hatvel * ihtot)
475 a(k) = a(k) - psi( zpa )
477 if (a(k)*vdml(i) > 0.0)
then
478 if (a(k)*vdml(i) > h_avail(i,j,k)) vdml(i) = h_avail(i,j,k) / a(k)
479 elseif (a(k)*vdml(i) < 0.0)
then
480 if (-a(k)*vdml(i) > h_avail(i,j+1,k)) vdml(i) = -h_avail(i,j+1,k) / a(k)
485 hatvel = 0.5*(h(i,j,k) + h(i,j+1,k))
487 zpb = zpb - (hatvel * ihtot_slow)
488 b(k) = b(k) - psi(zpb)
490 if (b(k)*vdml_slow(i) > 0.0)
then
491 if (b(k)*vdml_slow(i) > h_avail(i,j,k) - a(k)*vdml(i)) &
492 vdml_slow(i) = max( 0., h_avail(i,j,k) - a(k)*vdml(i) ) / b(k)
493 elseif (b(k)*vdml_slow(i) < 0.0)
then
494 if (-b(k)*vdml_slow(i) > h_avail(i,j+1,k) + a(k)*vdml(i)) &
495 vdml_slow(i) = -max( 0., h_avail(i,j+1,k) + a(k)*vdml(i) ) / b(k)
499 vhml(i,j,k) = a(k)*vdml(i) + b(k)*vdml_slow(i)
500 vhtr(i,j,k) = vhtr(i,j,k) + vhml(i,j,k)*dt
504 vtimescale_diag(i,j) = timescale
505 vdml_diag(i,j) = vdml(i)
509 do j=js,je ;
do k=1,nz ;
do i=is,ie
510 h(i,j,k) = h(i,j,k) - dt*g%IareaT(i,j) * &
511 ((uhml(i,j,k) - uhml(i-1,j,k)) + (vhml(i,j,k) - vhml(i,j-1,k)))
512 enddo ;
enddo ;
enddo
516 if (query_averaging_enabled(cs%diag))
then
517 if (cs%id_urestrat_time > 0)
call post_data(cs%id_urestrat_time, utimescale_diag, cs%diag)
518 if (cs%id_vrestrat_time > 0)
call post_data(cs%id_vrestrat_time, vtimescale_diag, cs%diag)
519 if (cs%id_uhml > 0)
call post_data(cs%id_uhml, uhml, cs%diag)
520 if (cs%id_vhml > 0)
call post_data(cs%id_vhml, vhml, cs%diag)
521 if (cs%id_MLD > 0)
call post_data(cs%id_MLD, mld_fast, cs%diag)
522 if (cs%id_Rml > 0)
call post_data(cs%id_Rml, rml_av_fast, cs%diag)
523 if (cs%id_uDml > 0)
call post_data(cs%id_uDml, udml_diag, cs%diag)
524 if (cs%id_vDml > 0)
call post_data(cs%id_vDml, vdml_diag, cs%diag)
526 if (cs%id_uml > 0)
then
527 do j=js,je ;
do i=is-1,ie
528 h_vel = 0.5*((htot_fast(i,j) + htot_fast(i+1,j)) + h_neglect)
529 udml_diag(i,j) = udml_diag(i,j) / (0.01*h_vel) * g%IdyCu(i,j) * (psi(0.)-psi(-.01))
531 call post_data(cs%id_uml, udml_diag, cs%diag)
533 if (cs%id_vml > 0)
then
534 do j=js-1,je ;
do i=is,ie
535 h_vel = 0.5*((htot_fast(i,j) + htot_fast(i,j+1)) + h_neglect)
536 vdml_diag(i,j) = vdml_diag(i,j) / (0.01*h_vel) * g%IdxCv(i,j) * (psi(0.)-psi(-.01))
538 call post_data(cs%id_vml, vdml_diag, cs%diag)
544 call diag_update_remap_grids(cs%diag)
546 end subroutine mixedlayer_restrat_general
550 subroutine mixedlayer_restrat_bml(h, uhtr, vhtr, tv, forces, dt, G, GV, US, CS)
554 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(inout) :: h
555 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)),
intent(inout) :: uhtr
557 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)),
intent(inout) :: vhtr
561 real,
intent(in) :: dt
564 real :: uhml(SZIB_(G),SZJ_(G),SZK_(G))
565 real :: vhml(SZI_(G),SZJB_(G),SZK_(G))
566 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)) :: &
569 real,
dimension(SZI_(G),SZJ_(G)) :: &
573 real :: Rho0(SZI_(G))
591 real :: uDml(SZIB_(G))
592 real :: vDml(SZI_(G))
593 real :: utimescale_diag(SZIB_(G),SZJ_(G))
594 real :: vtimescale_diag(SZI_(G),SZJB_(G))
597 real :: uDml_diag(SZIB_(G),SZJ_(G)), vDml_diag(SZI_(G),SZJB_(G))
600 integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, nkml
601 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
602 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB ; nkml = gv%nkml
604 if (.not.
associated(cs))
call mom_error(fatal,
"MOM_mixedlayer_restrat: "// &
605 "Module must be initialized before it is used.")
606 if ((nkml<2) .or. (cs%ml_restrat_coef<=0.0))
return
608 udml(:) = 0.0 ; vdml(:) = 0.0
610 g_rho0 = gv%g_Earth*us%L_to_m**2*us%s_to_T**2 / gv%Rho0
611 use_eos =
associated(tv%eqn_of_state)
612 h_neglect = gv%H_subroundoff
613 dz_neglect = gv%H_subroundoff*gv%H_to_Z
615 if (.not.use_eos)
call mom_error(fatal,
"MOM_mixedlayer_restrat: "// &
616 "An equation of state must be used with this module.")
631 htot(i,j) = 0.0 ; rml_av(i,j) = 0.0
634 call calculate_density(tv%T(:,j,k),tv%S(:,j,k),p0,rho0(:),is-1,ie-is+3,tv%eqn_of_state)
636 rml_av(i,j) = rml_av(i,j) + h(i,j,k)*rho0(i)
637 htot(i,j) = htot(i,j) + h(i,j,k)
638 h_avail(i,j,k) = max(i4dt*g%areaT(i,j)*(h(i,j,k)-gv%Angstrom_H),0.0)
643 rml_av(i,j) = (g_rho0*rml_av(i,j)) / (htot(i,j) + h_neglect)
653 do j=js,je;
do i=is-1,ie
654 h_vel = 0.5*(htot(i,j) + htot(i+1,j)) * gv%H_to_Z
656 u_star = us%s_to_T*0.5*(forces%ustar(i,j) + forces%ustar(i+1,j))
657 absf = 0.5*us%s_to_T*(abs(g%CoriolisBu(i,j-1)) + abs(g%CoriolisBu(i,j)))
661 mom_mixrate = (0.41*9.8696)*u_star**2 / &
662 (absf*h_vel**2 + 4.0*(h_vel+dz_neglect)*u_star)
663 timescale = 0.0625 * (absf + 2.0*mom_mixrate) / (absf**2 + mom_mixrate**2)
665 timescale = timescale * cs%ml_restrat_coef
668 udml(i) = timescale * g%mask2dCu(i,j)*g%dyCu(i,j)* &
669 g%IdxCu(i,j)*(rml_av(i+1,j)-rml_av(i,j)) * (h_vel**2 * gv%Z_to_H)
671 if (udml(i) == 0)
then
672 do k=1,nkml ; uhml(i,j,k) = 0.0 ;
enddo
674 i2htot = 1.0 / (htot(i,j) + htot(i+1,j) + h_neglect)
679 hx2 = (h(i,j,k) + h(i+1,j,k) + h_neglect)
680 a(k) = (hx2 * i2htot) * (2.0 - 4.0*(z_topx2+0.5*hx2)*i2htot)
681 z_topx2 = z_topx2 + hx2
682 if (a(k)*udml(i) > 0.0)
then
683 if (a(k)*udml(i) > h_avail(i,j,k)) udml(i) = h_avail(i,j,k) / a(k)
685 if (-a(k)*udml(i) > h_avail(i+1,j,k)) udml(i) = -h_avail(i+1,j,k)/a(k)
689 uhml(i,j,k) = a(k)*udml(i)
690 uhtr(i,j,k) = uhtr(i,j,k) + uhml(i,j,k)*dt
694 udml_diag(i,j) = udml(i)
695 utimescale_diag(i,j) = timescale
700 do j=js-1,je ;
do i=is,ie
701 h_vel = 0.5*(htot(i,j) + htot(i,j+1)) * gv%H_to_Z
703 u_star = us%s_to_T*0.5*(forces%ustar(i,j) + forces%ustar(i,j+1))
704 absf = 0.5*us%s_to_T*(abs(g%CoriolisBu(i-1,j)) + abs(g%CoriolisBu(i,j)))
708 mom_mixrate = (0.41*9.8696)*u_star**2 / &
709 (absf*h_vel**2 + 4.0*(h_vel+dz_neglect)*u_star)
710 timescale = 0.0625 * (absf + 2.0*mom_mixrate) / (absf**2 + mom_mixrate**2)
712 timescale = timescale * cs%ml_restrat_coef
715 vdml(i) = timescale * g%mask2dCv(i,j)*g%dxCv(i,j)* &
716 g%IdyCv(i,j)*(rml_av(i,j+1)-rml_av(i,j)) * (h_vel**2 * gv%Z_to_H)
717 if (vdml(i) == 0)
then
718 do k=1,nkml ; vhml(i,j,k) = 0.0 ;
enddo
720 i2htot = 1.0 / (htot(i,j) + htot(i,j+1) + h_neglect)
725 hx2 = (h(i,j,k) + h(i,j+1,k) + h_neglect)
726 a(k) = (hx2 * i2htot) * (2.0 - 4.0*(z_topx2+0.5*hx2)*i2htot)
727 z_topx2 = z_topx2 + hx2
728 if (a(k)*vdml(i) > 0.0)
then
729 if (a(k)*vdml(i) > h_avail(i,j,k)) vdml(i) = h_avail(i,j,k) / a(k)
731 if (-a(k)*vdml(i) > h_avail(i,j+1,k)) vdml(i) = -h_avail(i,j+1,k)/a(k)
735 vhml(i,j,k) = a(k)*vdml(i)
736 vhtr(i,j,k) = vhtr(i,j,k) + vhml(i,j,k)*dt
740 vtimescale_diag(i,j) = timescale
741 vdml_diag(i,j) = vdml(i)
745 do j=js,je ;
do k=1,nkml ;
do i=is,ie
746 h(i,j,k) = h(i,j,k) - dt*g%IareaT(i,j) * &
747 ((uhml(i,j,k) - uhml(i-1,j,k)) + (vhml(i,j,k) - vhml(i,j-1,k)))
748 enddo ;
enddo ;
enddo
753 if (cs%id_uhml > 0 .or. cs%id_vhml > 0) &
755 call pass_var(h, g%domain, to_west+to_south+omit_corners, halo=1)
756 call diag_update_remap_grids(cs%diag)
759 if (query_averaging_enabled(cs%diag) .and. &
760 ((cs%id_urestrat_time > 0) .or. (cs%id_vrestrat_time > 0)))
then
761 call post_data(cs%id_urestrat_time, utimescale_diag, cs%diag)
762 call post_data(cs%id_vrestrat_time, vtimescale_diag, cs%diag)
764 if (query_averaging_enabled(cs%diag) .and. &
765 ((cs%id_uhml>0) .or. (cs%id_vhml>0)))
then
767 do j=js,je ;
do i=isq,ieq ; uhml(i,j,k) = 0.0 ;
enddo ;
enddo
768 do j=jsq,jeq ;
do i=is,ie ; vhml(i,j,k) = 0.0 ;
enddo ;
enddo
770 if (cs%id_uhml > 0)
call post_data(cs%id_uhml, uhml, cs%diag)
771 if (cs%id_vhml > 0)
call post_data(cs%id_vhml, vhml, cs%diag)
772 if (cs%id_MLD > 0)
call post_data(cs%id_MLD, htot, cs%diag)
773 if (cs%id_Rml > 0)
call post_data(cs%id_Rml, rml_av, cs%diag)
774 if (cs%id_uDml > 0)
call post_data(cs%id_uDml, udml_diag, cs%diag)
775 if (cs%id_vDml > 0)
call post_data(cs%id_vDml, vdml_diag, cs%diag)
778 end subroutine mixedlayer_restrat_bml
782 logical function mixedlayer_restrat_init(Time, G, GV, US, param_file, diag, CS, restart_CS)
783 type(time_type),
intent(in) :: time
788 type(
diag_ctrl),
target,
intent(inout) :: diag
795 real :: flux_to_kg_per_s
797 # include "version_variable.h"
802 call get_param(param_file, mdl,
"MIXEDLAYER_RESTRAT", mixedlayer_restrat_init, &
803 "If true, a density-gradient dependent re-stratifying "//&
804 "flow is imposed in the mixed layer. Can be used in ALE mode "//&
805 "without restriction but in layer mode can only be used if "//&
806 "BULKMIXEDLAYER is true.", default=.false.)
807 if (.not. mixedlayer_restrat_init)
return
809 if (.not.
associated(cs))
then
810 call mom_error(fatal,
"mixedlayer_restrat_init called without an "// &
811 "associated control structure.")
815 cs%MLE_MLD_decay_time = -9.e9
816 cs%MLE_density_diff = -9.e9
817 cs%MLE_tail_dh = -9.e9
818 cs%MLE_use_PBL_MLD = .false.
819 cs%MLE_MLD_stretch = -9.e9
821 call get_param(param_file, mdl,
"DEBUG", cs%debug, default=.false., do_not_log=.true.)
822 call get_param(param_file, mdl,
"FOX_KEMPER_ML_RESTRAT_COEF", cs%ml_restrat_coef, &
823 "A nondimensional coefficient that is proportional to "//&
824 "the ratio of the deformation radius to the dominant "//&
825 "lengthscale of the submesoscale mixed layer "//&
826 "instabilities, times the minimum of the ratio of the "//&
827 "mesoscale eddy kinetic energy to the large-scale "//&
828 "geostrophic kinetic energy or 1 plus the square of the "//&
829 "grid spacing over the deformation radius, as detailed "//&
830 "by Fox-Kemper et al. (2010)", units=
"nondim", default=0.0)
834 call get_param(param_file, mdl,
"FOX_KEMPER_ML_RESTRAT_COEF2", cs%ml_restrat_coef2, &
835 "As for FOX_KEMPER_ML_RESTRAT_COEF but used in a second application "//&
836 "of the MLE restratification parameterization.", units=
"nondim", default=0.0)
837 call get_param(param_file, mdl,
"MLE_FRONT_LENGTH", cs%front_length, &
838 "If non-zero, is the frontal-length scale used to calculate the "//&
839 "upscaling of buoyancy gradients that is otherwise represented "//&
840 "by the parameter FOX_KEMPER_ML_RESTRAT_COEF. If MLE_FRONT_LENGTH is "//&
841 "non-zero, it is recommended to set FOX_KEMPER_ML_RESTRAT_COEF=1.0.",&
842 units=
"m", default=0.0)
843 call get_param(param_file, mdl,
"MLE_USE_PBL_MLD", cs%MLE_use_PBL_MLD, &
844 "If true, the MLE parameterization will use the mixed-layer "//&
845 "depth provided by the active PBL parameterization. If false, "//&
846 "MLE will estimate a MLD based on a density difference with the "//&
847 "surface using the parameter MLE_DENSITY_DIFF.", default=.false.)
848 call get_param(param_file, mdl,
"MLE_MLD_DECAY_TIME", cs%MLE_MLD_decay_time, &
849 "The time-scale for a running-mean filter applied to the mixed-layer "//&
850 "depth used in the MLE restratification parameterization. When "//&
851 "the MLD deepens below the current running-mean the running-mean "//&
852 "is instantaneously set to the current MLD.", units=
"s", default=0.)
853 call get_param(param_file, mdl,
"MLE_MLD_DECAY_TIME2", cs%MLE_MLD_decay_time2, &
854 "The time-scale for a running-mean filter applied to the filtered "//&
855 "mixed-layer depth used in a second MLE restratification parameterization. "//&
856 "When the MLD deepens below the current running-mean the running-mean "//&
857 "is instantaneously set to the current MLD.", units=
"s", default=0.)
858 if (.not. cs%MLE_use_PBL_MLD)
then
859 call get_param(param_file, mdl,
"MLE_DENSITY_DIFF", cs%MLE_density_diff, &
860 "Density difference used to detect the mixed-layer "//&
861 "depth used for the mixed-layer eddy parameterization "//&
862 "by Fox-Kemper et al. (2010)", units=
"kg/m3", default=0.03)
864 call get_param(param_file, mdl,
"MLE_TAIL_DH", cs%MLE_tail_dh, &
865 "Fraction by which to extend the mixed-layer restratification "//&
866 "depth used for a smoother stream function at the base of "//&
867 "the mixed-layer.", units=
"nondim", default=0.0)
868 call get_param(param_file, mdl,
"MLE_MLD_STRETCH", cs%MLE_MLD_stretch, &
869 "A scaling coefficient for stretching/shrinking the MLD "//&
870 "used in the MLE scheme. This simply multiplies MLD wherever used.",&
871 units=
"nondim", default=1.0)
872 call get_param(param_file, mdl,
"MLE_USE_MLD_AVE_BUG", cs%MLE_use_MLD_ave_bug, &
873 "If true, do not account for MLD mismatch to interface positions.",&
879 if (gv%Boussinesq)
then ; flux_to_kg_per_s = gv%Rho0
880 else ; flux_to_kg_per_s = 1. ;
endif
882 cs%id_uhml = register_diag_field(
'ocean_model',
'uhml', diag%axesCuL, time, &
883 'Zonal Thickness Flux to Restratify Mixed Layer',
'kg s-1', conversion=flux_to_kg_per_s, &
884 y_cell_method=
'sum', v_extensive=.true.)
885 cs%id_vhml = register_diag_field(
'ocean_model',
'vhml', diag%axesCvL, time, &
886 'Meridional Thickness Flux to Restratify Mixed Layer',
'kg s-1', conversion=flux_to_kg_per_s, &
887 x_cell_method=
'sum', v_extensive=.true.)
888 cs%id_urestrat_time = register_diag_field(
'ocean_model',
'MLu_restrat_time', diag%axesCu1, time, &
889 'Mixed Layer Zonal Restratification Timescale',
's')
890 cs%id_vrestrat_time = register_diag_field(
'ocean_model',
'MLv_restrat_time', diag%axesCv1, time, &
891 'Mixed Layer Meridional Restratification Timescale',
's')
892 cs%id_MLD = register_diag_field(
'ocean_model',
'MLD_restrat', diag%axesT1, time, &
893 'Mixed Layer Depth as used in the mixed-layer restratification parameterization',
'm')
894 cs%id_Rml = register_diag_field(
'ocean_model',
'ML_buoy_restrat', diag%axesT1, time, &
895 'Mixed Layer Buoyancy as used in the mixed-layer restratification parameterization', &
896 'm s2', conversion=us%m_to_Z)
897 cs%id_uDml = register_diag_field(
'ocean_model',
'udml_restrat', diag%axesCu1, time, &
898 'Transport stream function amplitude for zonal restratification of mixed layer',
'm3 s-1')
899 cs%id_vDml = register_diag_field(
'ocean_model',
'vdml_restrat', diag%axesCv1, time, &
900 'Transport stream function amplitude for meridional restratification of mixed layer',
'm3 s-1')
901 cs%id_uml = register_diag_field(
'ocean_model',
'uml_restrat', diag%axesCu1, time, &
902 'Surface zonal velocity component of mixed layer restratification',
'm s-1')
903 cs%id_vml = register_diag_field(
'ocean_model',
'vml_restrat', diag%axesCv1, time, &
904 'Surface meridional velocity component of mixed layer restratification',
'm s-1')
907 if (cs%MLE_MLD_decay_time>0. .or. cs%MLE_MLD_decay_time2>0.)
then
909 (gv%m_to_H_restart /= 0.0) .and. (gv%m_to_H_restart /= gv%m_to_H))
then
910 h_rescale = gv%m_to_H / gv%m_to_H_restart
911 do j=g%jsc,g%jec ;
do i=g%isc,g%iec
912 cs%MLD_filtered(i,j) = h_rescale * cs%MLD_filtered(i,j)
916 if (cs%MLE_MLD_decay_time2>0.)
then
917 if (
query_initialized(cs%MLD_filtered_slow,
"MLD_MLE_filtered_slow", restart_cs) .and. &
918 (gv%m_to_H_restart /= 0.0) .and. (gv%m_to_H_restart /= gv%m_to_H))
then
919 h_rescale = gv%m_to_H / gv%m_to_H_restart
920 do j=g%jsc,g%jec ;
do i=g%isc,g%iec
921 cs%MLD_filtered_slow(i,j) = h_rescale * cs%MLD_filtered_slow(i,j)
927 if (
associated(cs%MLD_filtered))
call pass_var(cs%MLD_filtered, g%domain)
929 end function mixedlayer_restrat_init
932 subroutine mixedlayer_restrat_register_restarts(HI, param_file, CS, restart_CS)
940 logical :: mixedlayer_restrat_init
943 call get_param(param_file, mdl,
"MIXEDLAYER_RESTRAT", mixedlayer_restrat_init, &
944 default=.false., do_not_log=.true.)
945 if (.not. mixedlayer_restrat_init)
return
948 if (
associated(cs))
call mom_error(fatal, &
949 "mixedlayer_restrat_register_restarts called with an associated control structure.")
952 call get_param(param_file, mdl,
"MLE_MLD_DECAY_TIME", cs%MLE_MLD_decay_time, &
953 default=0., do_not_log=.true.)
954 call get_param(param_file, mdl,
"MLE_MLD_DECAY_TIME2", cs%MLE_MLD_decay_time2, &
955 default=0., do_not_log=.true.)
956 if (cs%MLE_MLD_decay_time>0. .or. cs%MLE_MLD_decay_time2>0.)
then
958 allocate(cs%MLD_filtered(hi%isd:hi%ied,hi%jsd:hi%jed)) ; cs%MLD_filtered(:,:) = 0.
959 vd = var_desc(
"MLD_MLE_filtered",
"m",
"Time-filtered MLD for use in MLE", &
960 hor_grid=
'h', z_grid=
'1')
963 if (cs%MLE_MLD_decay_time2>0.)
then
965 allocate(cs%MLD_filtered_slow(hi%isd:hi%ied,hi%jsd:hi%jed)) ; cs%MLD_filtered_slow(:,:) = 0.
966 vd = var_desc(
"MLD_MLE_filtered_slow",
"m",
"c Slower time-filtered MLD for use in MLE", &
967 hor_grid=
'h', z_grid=
'1')
971 end subroutine mixedlayer_restrat_register_restarts